目次

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

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

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

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

x = c(1,3,4,2,3,4,3,2,1)

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

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

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

正規分布に従うデータから、正規分布のパラメータを推定するには、次のようにする。 

fitdistr(c(1,2,1,3,4,1,2,3,1),densfun="normal")

ワイブル分布のパラメータを推定するには、次の１行を実行する。 

fitdistr(c(1,2,1,3,4,1,2,3,1),densfun="weibull")

データを別のオブジェクトに保存してから、関数fitdistrに渡しても良い。 

x = c(1,2,1,3,4,1,2,3,1)
fitdistr(x,densfun="weibull")

初期値を指定する場合には、確率分布ごとのパラメータ名つきのリストで与える。 

fitdistr(c(1,2,1,3,4,1,2,3,1),densfun="weibull",start=list(shape=1,scale=1))
library(MASS)

data.frameのようにデータを与えてあるとする。関数fitdistrを用いて、ワイブル分布を仮定し、その２つの母数を推定してみる。 

fitdistr(X$time, densfun="weibull")

その結果、次の数値を得る。 

     shape       scale  
  1.6822280   6.8219220 
 (0.4279334) (1.3515357)

ワイブル分布には形状母数shapeと尺度母数scaleがある。それぞれの点推定値は1.6822280と6.8219220であり、標準誤差は括弧の中に表示される。

これらは、次のようにしても表示される。 

X.fit.weibull <- fitdistr(X$time, densfun="weibull")
names(X.fit.weibull)
X.fit.weibull$estimate
X.fit.weibull$sd
X.fit.weibull$loglik

最後の数値は、最尤推定法で最大化した対数尤度関数の最大値である。この値を参考に、寿命分布を選ぶこともある。

対数正規分布も推定しておく。 

X.fit.lognormal <- fitdistr(X$time, densfun="lognormal")
names(X.fit.lognormal)
X.fit.lognormal$estimate
X.fit.lognormal$sd
X.fit.lognormal$loglik

同様に、ガンマ分布。 

X.fit.gamma <- fitdistr(X$time, densfun="gamma")
names(X.fit.gamma)
X.fit.gamma$estimate
X.fit.gamma$sd
X.fit.gamma$loglik