分布の推定
1 はじめに
データから母集団の確率分布を推定することは重要である。
正規近似ができることを理論的に知っていても、あてはまりがどの程度かを知り、必要に応じて標本数を増やすなどの措置を講ずることは調査研究には必要と思う。
2 fitdist関数
草彅の研究ログ様の説明をまとめてみます。
library(fitdistrplus)## Loading required package: MASS## Loading required package: survival2.1 数値例の作成
ガンマ分布
rgamma (n, shapeは回数, scale = 1/rateは周期,
rateは振動数)
単位時間中にある事象が発生する平均回数をλとする。（周期：\(\mu=\dfrac{1}{\lambda}\) としてもよい）
その事象のk回目が起こるまでの時間をXとするとき、Xはガンマ分布に従う。
確率密度関数：\(f\left(x| k,\lambda \right) =\dfrac{\lambda ^{k}}{\Gamma \left( k\right)}x^{k-1}e^{-\lambda x}=\dfrac{1}{\mu ^{k}\Gamma \left( k\right)}x^{k-1}e^{-\frac{x}{\mu}}\).
パラメータ： 回数 k （形状尺度） 周期 \(\mu\)または振動数 \(\lambda\)（母数尺度） .
期待値： \(k\mu=\dfrac{k}{\lambda}=\dfrac{2}{3}=0.67\)
#周期1/3のイベントが2回起きるまでの時間Xの分布
dat<-rgamma(300,2,3)
#分布の確認
hist(dat)
#mean(dat)
plot(ecdf(dat))#データから経験分布関数を計算する関数 
- Empirical
は経験分布を用いて推定（データから推定）したもの。.
- QQプロットは、（片方の\(q\)分位数 , もう片方の\(q\)分位数）をプロットしたもの。
- PPプロットは、それぞれの分布関数 \(F_1,F_2\)に対し、\((F_1(x),~F_2(x))\)をプロットしたもの。
注意
\(0<q<1\)
は媒介変数としての役割を果たすので、曲線の媒介変数表示をしていることになります。
- ２つの分布が一致しているかを確認したいときは、Q-QプロットまたはP-Pプロットが
y=x に近いかで判断
- ２つの分布が同じ種類の分布なのかを確認したいときは、Q-Qプロットが直線かどうかで判断
具体例で学ぶ数学様「Q-QプロットとP-Pプロットの意味と違い」参照
#これをガンマ分布に当てはめてみる
fit<-fitdist(dat,"gamma")
plot(fit)
summary(fit)## Fitting of the distribution ' gamma ' by maximum likelihood
## Parameters :
## estimate Std. Error
## shape 2.237434 0.1708232
## rate 3.543868 0.3031775
## Loglikelihood: -117.0049 AIC: 238.0098 BIC: 245.4174
## Correlation matrix:
## shape rate
## shape 1.0000000 0.8924357
## rate 0.8924357 1.0000000- Loglikelihood: 対数尤度
- AIC: \(AIC=-2\log L+2k\)
を最小にするモデルを選ぶ
- BIC: \(AIC=-2\log L+k\log n\) を最小にするモデルを選ぶ
参考サイト.
3 適合度比較.
データを，ガンマ分布、ワイブル分布，対数正規分布に当てはめて，適合度を比較する。
当然、ガンマ分布への当てはまりが良いはずであるが。
#フィットさせる
gam<-fitdist(dat,"gamma")
weib<-fitdist(dat,"weibull")
logn<-fitdist(dat,"lnorm")
#比較（リストで入れる）
gofstat(list(gam, weib, logn), fitnames=c("gamma", "Weibull", "lognormal"))## Goodness-of-fit statistics
## gamma Weibull lognormal
## Kolmogorov-Smirnov statistic 0.04053029 0.02331442 0.0827493
## Cramer-von Mises statistic 0.07882477 0.03370755 0.5155087
## Anderson-Darling statistic 0.49834184 0.32320174 2.9287695
##
## Goodness-of-fit criteria
## gamma Weibull lognormal
## Akaike's Information Criterion 238.0098 238.1026 271.1472
## Bayesian Information Criterion 245.4174 245.5102 278.5548赤池情報量基準 AIC をみると、ガンマ分布とワイブル分布では当てはまりに差がない。
3.1 Goodness-of-fit.
- Kolmogorov-Smirnov statistic.
九州大学 工学部地球環境工学科 様によれば、有意水準5%の境界値は 1.36. - Cramer-von Mises statistic.
統計科学研究所様によれば、有意水準 5%の境界値は 0.46136 - Anderson-Darling statistic.
econoshift.comアンダーソン・ダーリング検定の5%の境界値は 2.492
3.2 二項分布
二項分布 \(B(n,~p)\)
は、nが大きいとき正規分布 \(N(np,np(1-p))\) で近似される。
Rの二項乱数を用いてシミュレーションする。
rbinom(n, size, prob) は、
確率probで起こる試行をsize回行ったときの二項乱数、nは乱数の個数である。
library(fitdistrplus)3.2.1 サイズ30.
正規分布へ当てはめ
theta=1/6;n=30
x <- rbinom(500, n, theta)
fit <- fitdist(x,"norm","mle")
summary(fit)## Fitting of the distribution ' norm ' by maximum likelihood
## Parameters :
## estimate Std. Error
## mean 5.01400 0.09288492
## sd 2.07697 0.06567949
## Loglikelihood: -1074.924 AIC: 2153.849 BIC: 2162.278
## Correlation matrix:
## mean sd
## mean 1.000000e+00 -4.161371e-09
## sd -4.161371e-09 1.000000e+00#gofstat(fit)
plot(fit)
ポアソン分布へ当てはめ.
theta=1/6;n=30
x <- rbinom(500, n, theta)
fit <- fitdist(x,"pois","mle")
summary(fit)## Fitting of the distribution ' pois ' by maximum likelihood
## Parameters :
## estimate Std. Error
## lambda 5.12 0.1011929
## Loglikelihood: -1076.088 AIC: 2154.176 BIC: 2158.391#gofstat(fit)
plot(fit)
library(fitdistrplus)
theta=1/6;n=30
x <- rbinom(500, n, theta)
fit <- fitdist(x,"binom",fix.arg=list(size = 30), start=list(prob=0.1))
summary(fit)## Fitting of the distribution ' binom ' by maximum likelihood
## Parameters :
## estimate Std. Error
## prob 0.1678668 0.00305155
## Fixed parameters:
## value
## size 30
## Loglikelihood: -1067.048 AIC: 2136.097 BIC: 2140.311gofstat(fit)## Chi-squared statistic: 4.569746
## Degree of freedom of the Chi-squared distribution: 6
## Chi-squared p-value: 0.6000542
## Chi-squared table:
## obscounts theocounts
## <= 2 53.00000 49.93520
## <= 3 62.00000 67.23527
## <= 4 96.00000 91.55309
## <= 5 95.00000 96.03918
## <= 6 79.00000 80.72523
## <= 7 47.00000 55.83347
## <= 8 41.00000 32.38207
## > 8 27.00000 26.29650
##
## Goodness-of-fit criteria
## 1-mle-binom
## Akaike's Information Criterion 2136.097
## Bayesian Information Criterion 2140.311plot(fit)
適合度比較.
データを，二項分布、正規分布，ポアソン分布に当てはめて，適合度を比較する。
当然、二項分布への当てはまりが良いはずであるが。
#フィットさせる
binom <- fitdist(x,"binom",fix.arg=list(size = 30),
start=list(prob = 1/6))
norm<-fitdist(x,"norm")
pois<-fitdist(x,"pois")
#比較（リストで入れる）
gofstat(list(binom,norm, pois), fitnames=c("binom","norm", "pois"))## Chi-squared statistic: 4.569654 39.03417 8.012722
## Degree of freedom of the Chi-squared distribution: 6 5 6
## Chi-squared p-value: 0.6000664 2.337526e-07 0.2371727
## Chi-squared table:
## obscounts theo binom theo norm theo pois
## <= 2 53 49.93438 34.89021 60.82625
## <= 3 62 67.23460 45.55659 69.17814
## <= 4 96 91.55263 73.09124 87.09528
## <= 5 95 96.03916 92.96762 87.72237
## <= 6 79 80.72561 93.74794 73.62831
## <= 7 47 55.83400 74.94727 52.97031
## <= 8 41 32.38253 47.50107 33.34481
## > 8 27 26.29708 37.29805 35.23452
##
## Goodness-of-fit criteria
## binom norm pois
## Akaike's Information Criterion 2136.097 2143.176 2143.763
## Bayesian Information Criterion 2140.311 2151.606 2147.9783.2.2 サイズ100.
theta=1/6;n=100
x <- rbinom(500, n, theta)
fit <- fitdist(x,"norm","mle")
summary(fit)## Fitting of the distribution ' norm ' by maximum likelihood
## Parameters :
## estimate Std. Error
## mean 16.654000 0.1649138
## sd 3.687585 0.1166116
## Loglikelihood: -1361.955 AIC: 2727.91 BIC: 2736.34
## Correlation matrix:
## mean sd
## mean 1.000000e+00 3.935334e-08
## sd 3.935334e-08 1.000000e+00plot(fit)
gofstat(fit)## Goodness-of-fit statistics
## 1-mle-norm
## Kolmogorov-Smirnov statistic 0.08311536
## Cramer-von Mises statistic 0.39244489
## Anderson-Darling statistic 2.20997737
##
## Goodness-of-fit criteria
## 1-mle-norm
## Akaike's Information Criterion 2727.91
## Bayesian Information Criterion 2736.34適合度検定 gofstat() には、
- Kolmogorov-Smirnov statistic.
- Cramer-von Mises statistic.
- Anderson-Darling statistic.
が出力されているが、検出力には差がありそうである。
Shapiro-Wilk 検定は出力されていない。
3.2.3 適合度比較.
適合度を比較する。
#フィットさせる
binom<-fitdist(x,"binom",fix.arg=list(size = 100),
start=list(prob = 1/6))
norm<-fitdist(x,"norm")
pois<-fitdist(x,"pois")
#比較（リストで入れる）
gofstat(list(binom,norm, pois), fitnames=c("binom","norm", "pois"))## Chi-squared statistic: 11.86516 26.10621 20.58888
## Degree of freedom of the Chi-squared distribution: 10 9 10
## Chi-squared p-value: 0.2941821 0.001962808 0.02414994
## Chi-squared table:
## obscounts theo binom theo norm theo pois
## <= 11 32 39.09411 31.30340 48.86985
## <= 12 25 26.08507 20.42727 27.80452
## <= 13 48 35.28337 28.70386 35.61973
## <= 14 46 43.81270 37.49102 42.37222
## <= 15 54 50.19333 45.51677 47.04446
## <= 16 43 53.28240 51.36569 48.96740
## <= 17 49 52.60814 53.88061 47.97077
## <= 18 45 48.47271 52.53512 44.38362
## <= 19 50 41.80193 47.61290 38.90342
## <= 20 39 33.82908 40.11042 32.39487
## <= 22 39 44.22881 54.26943 45.13852
## > 22 30 31.30835 36.78351 40.53060
##
## Goodness-of-fit criteria
## binom norm pois
## Akaike's Information Criterion 2719.906 2727.91 2729.309
## Bayesian Information Criterion 2724.120 2736.34 2733.524サイズが大きくなると、正規分布でも良い近似になる。
4 研究.
日新火災海上保険株式会社様「R と Excel を用いた分布推定の実践例」によれば、
分布推定に伴うリスク.
データから推定した分布に従って将来予測を行うと、
推定分布と真の分布との差に起因する誤差が生じる.
- パラメータリスク.
パラメータの推定誤差に起因するリスク.
- モデルリスク.
モデル(確率分布の種類)の選択誤りに起因するリスク.
これらのリスクを計量する方法を考える.
4.1 参考サイト.
pythonを使った分布推定の解説サイトを含め、今後、勉強したいサイトを列挙しておく。