==== 確率分布のパラメータの最尤推定 ====

MASSパッケージに収められている関数fitdistrは、「互いに独立に同一の確率分布に従う確率変数の実現値」であるデータから確率分布のパラメータの最尤推定値を算出してくれる。

|関数|fitdistr|
|目的|確率分布のパラメータの最尤推定値を算出する|
|書式|fitdistr(データ, 確率分布, 初期値もしくは初期値リスト)|
|返り値(表示)|パラメータ名と最尤推定値と(漸近)標準誤差|
|返り値(詳細)|estimate:最尤推定値のリスト, sd:標準誤差のリスト, vcov:漸近共分散行列, n:サンプルサイズ, loglik:対数尤度の最大値|

データはc()で囲んだ数値列(リストという)で与える。

<code>
x = c(1,3,4,2,3,4,3,2,1)
</code>

確率分布は次の一覧から選んで指定する。

|確率分布|与え方|パラメータ数|初期値を与えるべきか否か|
|ベータ分布|densfun="beta"|2|start=c(shape1=2,shape2=1)初期値必須|
|コーシー分布(自由度1のt分布)|densfun="cauchy"|2|start=list(location=,scale=)|
|カイ２乗分布|densfun="chi-squared"|2|不明|
|指数分布|densfun="exponential"|1|初期値不要|
|ガンマ分布|densfun="gamma"|2|start=list(shape=,rate=)|
|幾何分布|densfun="geometric"|1|初期値不要|
|対数正規分布|densfun="log-normal"|2|初期値不要|
|対数正規分布|densfun="lognormal"|2|初期値不要|
|ロジスティック分布|densfun="logistic"|2|start=list(location=,scale=)|
|負の二項分布|densfun="negative binomial"|2|start=list(size=,mu=)|
|正規分布|densfun="normal"|2|初期値不要|
|ポアソン分布|densfun="Poisson"|1|初期値不要|
|t分布|densfun="t"|2|不明|
|ワイブル分布|densfun="weibull"|2|start=list(shape=,scale=)|

初期値はパラメータの数分の長さの数字列を名前付きリストで与える。

正規分布に従うデータから、正規分布のパラメータを推定するには、次のようにする。
<code>
fitdistr(c(1,2,1,3,4,1,2,3,1),densfun="normal")
</code>

ワイブル分布のパラメータを推定するには、次の１行を実行する。
<code>
fitdistr(c(1,2,1,3,4,1,2,3,1),densfun="weibull")
</code>
データを別のオブジェクトに保存してから、関数fitdistrに渡しても良い。
<code>
x = c(1,2,1,3,4,1,2,3,1)
fitdistr(x,densfun="weibull")
</code>
初期値を指定する場合には、確率分布ごとのパラメータ名つきのリストで与える。
<code>
fitdistr(c(1,2,1,3,4,1,2,3,1),densfun="weibull",start=list(shape=1,scale=1))
</code>

<code>
library(MASS)
</code>

[[r:survival:data.frame]]のようにデータを与えてあるとする。関数fitdistrを用いて、ワイブル分布を仮定し、その２つの母数を推定してみる。
<code>
fitdistr(X$time, densfun="weibull")
</code>
その結果、次の数値を得る。
<code>
     shape       scale  
  1.6822280   6.8219220 
 (0.4279334) (1.3515357)
</code>
ワイブル分布には形状母数shapeと尺度母数scaleがある。それぞれの点推定値は1.6822280と6.8219220であり、標準誤差は括弧の中に表示される。

これらは、次のようにしても表示される。
<code>
X.fit.weibull <- fitdistr(X$time, densfun="weibull")
names(X.fit.weibull)
X.fit.weibull$estimate
X.fit.weibull$sd
X.fit.weibull$loglik
</code>
最後の数値は、最尤推定法で最大化した対数尤度関数の最大値である。この値を参考に、寿命分布を選ぶこともある。

対数正規分布も推定しておく。
<code>
X.fit.lognormal <- fitdistr(X$time, densfun="lognormal")
names(X.fit.lognormal)
X.fit.lognormal$estimate
X.fit.lognormal$sd
X.fit.lognormal$loglik
</code>

同様に、ガンマ分布。
<code>
X.fit.gamma <- fitdistr(X$time, densfun="gamma")
names(X.fit.gamma)
X.fit.gamma$estimate
X.fit.gamma$sd
X.fit.gamma$loglik
</code>