草薙の研究ログ

英語の先生をやってます。

単変量の分布母数の推定からあれこれ:fitdistrplusパッケージが便利

ある確率分布を手持ちの単変量データに当てはめ,その母数を最尤推定するという方法があって,Rでは通常MASSパッケージのfitdistr関数を使うのがお手軽なんだけど,fitdistrplusという便利なパッケージが出ていてこれがとてもいいかも。ま,結局はoptimにかぶせているだけだが。

例えば,

  • さまざまな可視化
  • ブートストラップ
  • 任意の母数を固定した状態での別の母数の推定
  • 最尤推定以外の方法
  • 適合度比較

みたいなのをかなり手軽にやってくれる。いちいち自分で書かなくてもいい。

可視化

たとえば,可視化に関してならば,

#数値例の作成
library(fitdistrplus)
dat<-rgamma(300,2,3)

#分布の確認
hist(dat)
plot(ecdf(dat))

#これをガンマ分布に当てはめてみる
fit<-fitdist(dat,"gamma")
plot(fit)
summary(fit)


f:id:kusanagik:20170116191733p:plain


ほら,便利。これだけ。
左上は,ヒストグラムに推定された母数の値を項にもつ確率密度関数
右上はQ-Qプロット。この例だと,何故か大きい値があんま当てはまっていないことがわかる。
左下は経験累積密度関数と推定された母数の値の累積密度関数。右下はP-Pプロット。

もちろん,summaryも見やすい。こんな感じで普通欲しいのは一気に出してくれる。

Fitting of the distribution ' gamma ' by maximum likelihood 
Parameters : 
      estimate Std. Error
shape 1.965111  0.1488415
rate  2.841928  0.2450182
Loglikelihood:  -156.0496   AIC:  316.0992   BIC:  323.5068 
Correlation matrix:
          shape      rate
shape 1.0000000 0.8785199
rate  0.8785199 1.0000000


対数尤度プロットなる機能もある。これは,2つの母数の組み合わせの尤度をヒートマップみたいにしたのね。
尤度関数の2次元版,みたいなもんだろうか。ガンマは二母数だからね。

llplot(fit)


f:id:kusanagik:20170116195113p:plain

ブートストラップ

ブートストラップで誤差や信頼区間の検討もお手軽にできる。

fit<-fitdist(dat,"gamma")

#B=1000でパラメトリック
bfit1<-bootdist(fit,bootmethod="param",niter=1000)

#B=500でノンパラメトリック
bfit2<-bootdist(fit,bootmethod="nonparam",niter=500)

#結果の表示
summary(bfit1)
summary(bfit2)


ここではブートストラップした値の中央値と2.5%と97.5%点を出してくれる。←これはパーセンタイル法でのブートストラップ信頼区間でもある。

結果を色々見るには,

#パーセンタイル法での信頼区間
bfit1$CI

#ブートストラップ推定値
boot.shape<-bfit1$estim[,1]
boot.shape
boot.rate<-bfit1$estim[,2]
boot.rate

#これをヒストグラムで描いて,中央値,パーセンタイル信頼区間を描き入れる。

par(mfrow=c(1,2))
hist(boot.shape,col="lightgray")
abline(v=quantile(boot.shape,c(0.025,.5,.975)),col=2)
hist(boot.rate,col="lightgray")
abline(v=quantile(boot.rate,c(0.025,.5,.975)),col=2)


f:id:kusanagik:20170116195204p:plain


片方の母数を固定してもう片方の母数を最尤推定

やりたかったのこれ。ある理論的な見地から,片方の母数,例えば形状母数を任意の値に固定した状態で尺度母数を推定したいとか,これって結構外国語教育研究的によくあること。こういうときは,

#形状母数を無理やり3に固定(ホントは2)
fitfix<-fitdist(dat,"gamma",method="mle",fix.arg=list("shape"=3))
#結果を見る
summary(fitfix)

Fitting of the distribution ' gamma ' by maximum likelihood 
Parameters : 
     estimate Std. Error
rate 4.525801    0.15086
Fixed parameters:
      value
shape     3
Loglikelihood:  -149.1229   AIC:  300.2459   BIC:  303.9496 

#1つだから尺度母数のみの対数尤度曲線が出る
llplot(fitfix)


f:id:kusanagik:20170116200537p:plain

最尤推定以外の推定方法

最尤推定以外の方法も実装されている。methodに入れてやればいい。

  • mle 最尤推定。もちろん初期値も与えられるし,optimの方法も指定できる
  • mme モーメント法
  • qme 分位点法
  • mge 最大適合度法 よく知らない
fitfist(dat,"gamma",method="mme")
fitfist(dat,"gamma",method="qme")
fitfist(dat,"gamma",method="mge")

適合度比較

件のデータを,ワイブル,対数正規分布にも当てはめてみて,適合度を比較する。

#フィットさせる
gam<-fitdist(dat,"gamma")
weib<-fitdist(dat,"gamma")
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.04626154 0.04626154 0.08208341
Cramer-von Mises statistic   0.05906861 0.05906861 0.43467475
Anderson-Darling statistic   0.40584196 0.40584196 2.93641451

Goodness-of-fit criteria
                                  gamma  Weibull lognormal
Aikake's Information Criterion 278.3222 278.3222  321.8015
Bayesian Information Criterion 285.7298 285.7298  329.2090


こんな感じ。便利ね。コルモゴロフ・スミルノフやアンダーソン・ダーリング検定ってのはわかる。残りのクラメール・フォン=ミーゼスってのは,こんなこと?なるほど。


Cramér–von Mises criterion - Wikipedia


…ま,最初はこうやってモデリングについて慣れていくことが大事で,こういうのがまさに出発点って感じだなあと最近しみじみ。