MATLAB Note/統計/確率分布
Last-modified: 2013-01-19 (土) 03:07:52
確率分布(確率密度関数) †
- 確率分布の特性と分類に関して、さまざまな確率分布(滋賀大中川先生) に総合的な解説があります。
- 以下、Statistics Toolbox の関数を使用します。
- 参考:Statistics Toolbox でサポートされている確率分布
- Rでも同じことができるはずです。R Note/統計/確率分布 を参照して下さい。
連続変数(Continuous variable) †
- 連続的に変化する値の出現頻度を扱う分布。
正規分布(Normal distribution)、ガウス分布(Gaussian distribution) †
- μ = 0, σ = 1 の正規分布(標準正規分布)を作ってみます。
- まずは以下を実行します。
%パラメータの宣言 mu = 0; sigma = 1;
- この分布に従う乱数を10000個(10000×1の配列で)作るには、関数 normrnd を使います。
data = normrnd(mu, sigma, [10000 1]); % 乱数生成 bar(data); xlim([0,10000]); title('normrnd'); % プロット
- ヒストグラムで書いてみます。
[histdata, position] = hist(data); % ヒストグラム bar(position, histdata); xlim([-5,5]); title('normrnd'); % プロット
- ヒストグラムの頂点を結んだ線が確率密度曲線(の近似)になります。より厳密な正規分布の確率密度曲線は、関数 normpdf でプロットできます。ここでは x座標 -5 ~ 5 までの区間(0.01刻み)で関数生成&プロットしてみます。
linescale = [-5:0.01:5]; % x座標のスケール pdfdata = normpdf(linescale, mu, sigma); % 確率密度関数生成 plot(linescale, pdfdata); xlim([-5,5]); title('normpdf'); % プロット
- ヒストグラムの各バーを、x座標の小さいものから順に足しあわせてプロットしていけば累積密度曲線(の近似)になります。より厳密な正規分布の累積密度曲線は、関数 normcdf でプロットできます。
histsum = cumsum(histdata); % ヒストグラムの累積和 cdfdata = normcdf(linescale, mu, sigma); % 累積密度関数生成 subplot(2,1,1); bar(position, histsum); xlim([-5,5]); title('normcdf'); subplot(2,1,2); plot(linescale, cdfdata); xlim([-5,5]); % プロット
- まずは以下を実行します。
- 得られたデータセットが正規分布に従うと仮定して、データからμとσの値を推定することもできます。
- データのフィッティングには関数 normfit を使います。
- 上で作ったμ = 0, σ = 1 の乱数正規分布のデータを使ってパラメータの推定を行います(参考:Statistics Toolbox/正規分布による近似)。
[muhat,sigmahat,muci,sigmaci] = normfit(data) % パラメータ推定 pdfdata = normpdf(linescale, mu, sigma); % 推定したパラメータから確率密度関数生成 [histdata, position] = hist(data); % 乱数正規分布のヒストグラム histdata = histdata / (max(histdata) / max(pdfdata)); % ヒストグラムのスケール調整 figure, bar(position, histdata); xlim([-5,5]); % プロット hold on plot(linescale, pdfdata,'r','LineWidth',2); xlim([-5,5]); title('normfit'); hold off
- ここでは暫定的に、ヒストグラムの頂点と確率密度関数の頂点が一致するようにスケール調整をしていますが、本当ならヒストグラムの面積が1になるように調整するべきです。Statistics Toolbox/正規分布による近似を参考にして下さい。
- 以下が出力されます。muhatが推定した平均(μ)、sigmahatが推定した標準偏差(σ)です。muci、sigmaciは、それぞれ95%信頼区間です。これらは最尤推定の関数 mle によって計算されます。*3
muhat = 0.0017 sigmahat = 0.9915 muci = -0.0178 0.0211 sigmaci = 0.9779 1.0054
- なお、データが正規分布であることを仮定できるなら、ガウス混合分布モデルの混合数1でも同様のフィッティングができます。
options = statset('Display','final'); obj = gmdistribution.fit(data,1,'Options',options); % 混合数1のGMM pdfdata = normpdf(linescale, obj.mu, obj.Sigma); % GMMパラメータから確率密度関数生成 figure, bar(position, histdata); xlim([-5,5]); % プロット hold on plot(linescale, pdfdata,'r','LineWidth',2); xlim([-5,5]); title('normfit'); hold off
- ちなみに、正規分布の確率密度関数の二次導関数をメキシカン・ハット関数といい*4、ウェーブレット変換 や自己組織化マップ の近傍強化関数*5などに使われています。
plot( -1 * diff(diff(pdfdata)) ); title('Mexican hat function');
- ガウス混合分布モデル
- 複数の正規分布の混合とみなせるデータセットから、元の各正規分布のパラメータを推定する。
ガンマ分布(Gamma distribution) †
- ガンマ分布 - Wikipedia
- 形状母数 k、尺度母数 θ から求まる
- 左に歪んだL字分布。
- 【引用】 電子部品の寿命分布や通信工学におけるトラフィックの待ち時間分布に応用される。*6*7
- 蝸牛基底膜の聴覚細胞の応答関数に用いられる(ガンマチャープフィルタ、和歌山大入野先生)
- k = 3, θ = 2.0 のガンマ分布を作ってみます。
- まずは以下を実行します。
%パラメータの宣言 k = 3; theta = 2;
- この分布に従う乱数を10000個(10000×1の配列で)作るには、関数 gamrnd を使います。
data = gamrnd(k, theta, [10000 1]); % 乱数生成 bar(data); xlim([0,10000]); title('gamrnd'); % プロット
- ヒストグラムで書いてみます。
[histdata, position] = hist(data); % ヒストグラム bar(position, histdata); xlim([0,30]); title('gamrnd'); % プロット
- より厳密なガンマ分布の確率密度曲線は、関数 gampdf でプロットできます。
linescale = [0:0.01:30]; % x座標のスケール pdfdata = gampdf(linescale, k, theta); % 確率密度関数生成 plot(linescale, pdfdata); xlim([0,30]); title('gampdf'); % プロット
- 累積密度曲線は、関数 gamcdf でプロットできます(図は省略)。
histsum = cumsum(histdata); % ヒストグラムの累積和 cdfdata = gamcdf(linescale, k, theta); % 累積密度関数生成 subplot(2,1,1); bar(position, histsum); xlim([0,30]); title('gamcdf'); subplot(2,1,2); plot(linescale, cdfdata); xlim([0,30]); % プロット
- まずは以下を実行します。
ワイブル分布(Weibull distribution) †
ロジスティック分布(Logistic distribution) †
- ロジスティック分布 - Wikipedia
- 第一変数μ 、 第二変数 s から求まる
- 左右対称の分布。
- 【引用】 正規分布と同様に対称なS字(シグモイド)型の分布関数、釣鐘型の確率密度関数を持ち一見して両者は類似しているが、ロジスティック分布の方が裾が長く密度関数は平均から離れても下がりにくい*11
- 【引用】 ロジスティック分布の累積分布関数に見る曲線はロジスティック曲線と呼ばれ、分布そのものよりこちらの曲線の方が様々な場面で現れています*12
- データが直線で回帰できないときは、ロジスティック関数(曲線)で回帰モデルを作る(ロジスティック回帰)
- ロジスティック関数の逆関数を ロジット関数 といい、こちらも良く用いられる。
- ロジスティック関数のさまざまな応用例は ロジスティック分布(Logistic distribution) を参照
ベータ分布(Beta distribution) †
- Beta distribution - Wikipedia
- 第一変数 p、第二変数 q から求まる (英語版Wikipedia ではαとβ、以下は pとqで統一します)
- パラメータ次第で左右対称、左に歪んだL字分布、右に歪んだJ字分布を表現できる。
- 【引用】 独立に一様分布 U(0,1) に従う p+q-1 個の確率変数を大きさの順に並べ替えたとき,小さい方から p 番め(大きい方からは q 番目)の確率変数 X の分布がベータ分布 B(p,q) となる*13
- 【引用】 工学的な実用例はあまりない*14とされますが、様々な形状の分布を表現できるため統計モデル(一般化線形モデル)の確率分布として注目されています。*15
- p = 2.5, q = 1.5 のベータ分布を作ってみます。
- まずは以下を実行します。
%パラメータの宣言 p = 2.5; q = 1.5;
- この分布に従う乱数を10000個(10000×1の配列で)作るには、関数 betarnd を使います(図は省略)。
data = betarnd(p, q, [10000 1]); % 乱数生成 bar(data); xlim([0,10000]); title('betarnd'); % プロット
- ヒストグラムで書いてみます。
[histdata, position] = hist(data); % ヒストグラム bar(position, histdata); xlim([0,1.2]); title('betarnd'); % プロット
- x座標の値の範囲は0~1
- より厳密なベータ分布の確率密度曲線は、関数 betapdf でプロットできます。
linescale = [0:0.01:1.2]; % x座標のスケール pdfdata = betapdf(linescale, p, q); % 確率密度関数生成 plot(linescale, pdfdata); xlim([0,1.2]); title('betapdf'); % プロット
- まずは以下を実行します。
- ベータ分布の確率密度曲線の平均と分散は、関数 betastat で求めることができます。
[M, V] = betastat(p, q)
- p = 2.5, q = 1.5 のときの平均は0.6250、分散は0.0469です。これは、以下の式でも求めることができます(参考:ベータ分布の性質(関西学院高丹羽先生))。
M = p / (p + q) V = (p * q) / ((p + q) ^ 2 * (p + q + 1))
- 累積密度曲線は、関数 betacdf でプロットできます(図は省略)。
histsum = cumsum(histdata); % ヒストグラムの累積和 cdfdata = betacdf(linescale, p, q); % 累積密度関数生成 subplot(2,1,1); bar(position, histsum); xlim([0,1.2]); title('gamcdf'); subplot(2,1,2); plot(linescale, cdfdata); xlim([0,1.2]); % プロット
- 得られたデータセットがベータ分布に従うと仮定して、データからpとqの値を推定することもできます。
- データのフィッティングには関数 betafit を使います。
- 上で作ったp = 2.5, q = 1.5 の乱数ベータ分布をテストデータとして与えて、このパラメータを推定できるかどうか試してみます。
[phat, pci] = betafit(data) % パラメータ推定 pdfdata = betapdf(linescale, phat(1), phat(2)); % 確率密度関数生成 [histdata, position] = hist(data); % 乱数ベータ分布のヒストグラム histdata = histdata / (max(histdata) / max(pdfdata)); % ヒストグラムのスケール調整 figure, bar(position, histdata); xlim([0,1.2]); % プロット hold on plot(linescale, pdfdata, 'r','LineWidth',2); xlim([0,1.2]); title('betafit'); hold off
- ここでは暫定的に、ヒストグラムの頂点と確率密度関数の頂点が一致するようにスケール調整をしていますが、本当ならヒストグラムの面積が1になるように調整するべきです。Statistics Toolbox/正規分布による近似を参考にして下さい。
- 以下が出力されます。phat(1)が推定したpパラメータ、phat(2)が推定したqパラメータです。pci(:,1)、pci(:,2)は、それぞれpとqの95%信頼区間です。これらは最尤推定の関数 mle によって計算されます。
phat = 2.5035 1.4856 pci = 2.4353 1.4479 2.5717 1.5233
- p = 2.5035、q = 1.4856 になりました。おおむね正しく推定できています。
- ベータ混合分布モデル*16
- 複数のベータ分布の混合とみなせるデータセットから、元の各ベータ分布のパラメータを推定する。
- ガウス混合分布モデルでも、ある程度の分布の歪みには頑健に対応できるとのこと。
- 複数のベータ分布の混合とみなせるデータセットから、元の各ベータ分布のパラメータを推定する。
ディリクレ分布(Dirichlet distribution) †
- ディリクレ分布 - Wikipedia
- 【引用】ベータ分布を多変量に拡張して一般化した形をしており、そのため多変量ベータ分布とも呼ばれる*17
- ディリクレ分布とベータ分布の関係に関する議論
- ベータ分布とディリクレ分布(ryamadaの遺伝学・遺伝統計学メモ)
- 「統計学関連なんでもあり」BBS
- Infinite Mixture Model(ディリクレ分布(など)に基づく混合分布モデル)
離散変数(Categorical variable) †
- 離散的に変化する値(二値のみ、整数のみ等)の出現頻度を扱う分布。
- 【引用】 サイコロを投げるような場合における確率関数(略)。確率変数は離散的な値だけを取ることができます。このような場合を離散型分布といいます*18
二項分布(Binomial distribution) †
- n = 20, p = 0.4 の二項分布を作ってみます。
- ヒストグラムで書いてみます。
[histdata, position] = hist(data); % ヒストグラム bar(position, histdata); xlim([0,20]); title('binornd'); % プロット
- ヒストグラムの頂点を結んだ線が確率密度曲線(の近似)になります。より厳密な二項分布の確率密度曲線は、関数 binopdf でプロットできます。ここでは x座標 0 ~ 20 までの区間について「0.1刻み」と「1刻み」の二通りで関数生成&プロットしてみます。
linescale1 = [0:0.1:20]; % x座標のスケール(0.1刻み) pdfdata1 = binopdf(linescale1, n, p); % 確率密度関数生成 linescale2 = [0:1:20]; % x座標のスケール(1刻み) pdfdata2 = binopdf(linescale2, n, p); % 確率密度関数生成 subplot(2,1,1); plot(linescale1, pdfdata1); xlim([0,20]); title('binopdf'); subplot(2,1,2); plot(linescale2, pdfdata2); xlim([0,20]); % プロット
- このように、離散分布の関数は、x座標に連続値(小数点を含む値)を指定しても、x座標の値が整数のとき以外は0を出力します。
- ヒストグラムの各バーを、x座標の小さいものから順に足しあわせてプロットしていけば累積密度曲線(の近似)になります。より厳密な二項分布の累積密度曲線は、関数 binocdf でプロットできます。
histsum = cumsum(histdata); % ヒストグラムの累積和 cdfdata1 = binocdf(linescale1, n, p); % 累積密度関数生成(0.1刻み) cdfdata2 = binocdf(linescale2, n, p); % 累積密度関数生成(1刻み) subplot(3,1,1); bar(position, histsum); xlim([0,20]); title('binocdf'); subplot(3,1,2); plot(linescale1, cdfdata1); xlim([0,20]); subplot(3,1,3); plot(linescale2, cdfdata2); xlim([0,20]); % プロット
ポアソン分布 †
- 【引用】 1 回ごとの事象の生起確率 p が一定であるとき,実験を繰り返し行うことを ベルヌーイ試行 というが,独立に実験を n 回繰り返したとき,x 回事象が生じる確率は二項分布である。ポアソン分布は、二項分布の中でも特に出現確率の低い事象を扱うとき、よくあてはまる。*21
- 【引用】 製品中の不良品の個数,一定時間内に電話がかかってくる回数などがあげられる。*22
- 値に上限と下限をもつような、「閉じた尺度」のデータはポアソン分布に近づく。
- 0%~100%、0点~100点など。こういったデータは分布の中心が端に近づくにつれて分布の裾野がつまり(天井効果、フロア効果)、L字型やJ字型のポアソン分布を作りやすい。
- 片方にのみ開いた尺度(負値を持たないデータ、年収など)も、裾野が引っ張られてポアソン分布を作る。
- 値に上限と下限をもつような、「閉じた尺度」のデータはポアソン分布に近づく。
ベルヌーイ分布 †
分布形状の測定 †
- 以下、Statistics Toolbox / 形状の測度 に基づきます。
尖度 kurtosis と 歪度 skewness †
- 尖度:分布の頂点がどれだけ尖っているか。
- 歪度:分布の左右がどれだけ非対称か。
- 【引用】 確率分布の分布特性を示すためには期待値および分散が通常用いられるが、分布型の違いを示す指標の1つに3次モーメント(3乗の期待値)と4次モーメント(4乗の期待値)がある。(略)平均と分散の影響を除くように標準化された3次モーメントと標準化された4次モーメントを考える。(略)標準正規確率変数の分布に歪みはなく、0を中心に左右対称であるから歪度は0である。歪度の符号により正の歪み、負の歪みを持つ分布といわれる。(略)Z の4次モーメント E (Z4) - 3 は尖度 β2 と呼ばれる*24
- 正規分布(左右対称の分布)では歪度0、L字型分布は歪度+、J字型分布は歪度-
- 上で作った 正規分布 と ガンマ分布 の尖度、歪度を比べてみます。
mu = 0; sigma = 1; dataNorm = normrnd(mu, sigma, [10000 1]); % 正規分布に従う10000個の乱数 k = 3; theta = 2; dataGam = gamrnd(k, theta, [10000 1]); % ガンマ分布に従う10000個の乱数 k = kurtosis([dataNorm dataGam]) % 尖度(1列目:正規分布、2列目:ガンマ分布) y = skewness([dataNorm dataGam]) % 歪度(1列目:正規分布、2列目:ガンマ分布)
- 引数のベクトルは確率密度曲線ではなく、それぞれの分布に従う乱数(ヒストグラムを書く前の生データ)であることに注意して下さい。
- 結果
(正規) (ガンマ) k = 2.9537 5.4342 y = -0.0096 1.1623
- 正規分布の尖度は3、歪度は0に近くなっています。
- ガンマ分布の尖度は正規分布より大きく、歪度は正の値になっている(L字に近い)ことがわかります(これは k=3, θ=2.0 のガンマ分布の場合なので、他のパラメータを指定すれば尖度、歪度は変わります)。
- なお、歪度の大きい分布の形状を正規分布に近づけるための変換手法が提案されています。量的分析法 勉強会/確率分布/標本分布が正規分布ではないが、ANOVAを使いたい場合 を参照して下さい。
中心モーメント moment †
百分位数値(パーセンタイル)prctile と 分位数(クォンタイル) quantile †
標準化z得点 zscore †
分布形状の検定 †
- あるデータ群が正規分布に従っているかどうかを検定したいときは、1標本コルモゴルフ-スミルノフ検定 を使います。
- 2つのデータ群の分布が一致しているかどうかを検定したいときは、2標本コルモゴルフ-スミルノフ検定 を使います。
- その他
- Matlabによる初歩の検定(理研 最上先生) でも、kstwo.m を公開なさっています。
- Rによるコルモゴロフ・スミルノフ検定の実例
- QQプロット(正規性検定)
*1 確率と統計 3.3.2 連続型分布(静岡理工科大 菅沼先生) より引用
*2 確率と統計 3.3.2 連続型分布(静岡理工科大 菅沼先生) を参考にしました。
*3 最尤推定によるベータ分布のパラメータ推定に関して、以下のウェブサイトがたいへん詳しいです:ベータ分布の最尤推定(Rで実装) BIC (ベイズ情報量基準) その3 また脱線でベータ分布の最尤推定
*4 参考:短時間フーリエ変換と連続ウェーブレット変換(信州大井澤先生)
*5 内積型自己組織化マップでは、正規化によって抑制型学習を行うので、強化関数は普通の正規分布です。
*6 ガンマ分布 - Wikipedia より引用
*7 ポアソン分布、ワイブル分布、ベータ分布、ガンマ分布の特性と工学的な応用については、ポワソン配置とワイブル分布(雑然か整然か) や ポアソン過程・指数分布・ガンマ分布・ワイブル分布・ベータ分布 に詳細な説明あり。
*8 ワイブル分布 - Wikipedia より引用
*9 ワイブル分布 - Wikipedia による
*10 ポアソン分布、ワイブル分布、ベータ分布、ガンマ分布の特性と工学的な応用については、ポワソン配置とワイブル分布(雑然か整然か) や ポアソン過程・指数分布・ガンマ分布・ワイブル分布・ベータ分布 に詳細な説明あり。
*11 ロジスティック分布 - Wikipedia による
*12 ロジスティック分布(Logistic distribution) より引用
*13 ベータ分布の性質 より引用
*14 ベータ分布の性質 より引用
*15 例えば、Smithson M, Verkuilen J, "A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables," Psychol Methods, 11(1), 54-71, 2006. などで議論されています。(理研のIさんが教えて下さいました。ありがとうございました。)
*16 参考:Ji Y, Wu C, Liu P, Wang J, Coombes KR, "Applications of beta-mixture models in bioinformatics," Bioinformatics, Vol.21, No.9, pp.2118-2122, 2005. 著者ウェブサイト からRスクリプトのダウンロードが出来る。
*17 ディリクレ分布 - Wikipedia より引用
*18 確率と統計 3.3.1 離散型分布(静岡理工科大 菅沼先生) より引用
*19 二項分布 - Wikipediaより引用
*20 二項分布(群馬大青木先生)より引用
*21 離散分布と連続分布 の ポアソン分布(群馬大青木先生)より引用
*22 ポアソン分布 より引用
*23 尖度 - Wikipedia より引用
*24 歪度 - Wikipedia より引用
添付ファイル: plot_beta_4.png 3304件 [詳細] plot_norm_GMM.png 3185件 [詳細] plot_norm_5.png 3264件 [詳細] plot_norm_mexican.png 3381件 [詳細] plot_beta_3.png 3433件 [詳細] plot_beta_2.png 3195件 [詳細] plot_bino_4.png 3126件 [詳細] plot_bino_3.png 3193件 [詳細] plot_bino_2.png 3105件 [詳細] plot_bino_1.png 3122件 [詳細] plot_gam_3.png 3309件 [詳細] plot_norm_4.png 3373件 [詳細] plot_gam_2.png 3271件 [詳細] plot_gam_1.png 3136件 [詳細] plot_norm_3.png 3528件 [詳細] plot_norm_2.png 3255件 [詳細] plot_norm_1.png 3189件 [詳細]