草薙の研究ログ

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

混合ガンマ分布モデルをデータにフィットさせる

Rのmixtoolsパッケージでは,混合ガンマ分布モデルもデータへフィットさせることができる。

#λ= .25,.75,α - 4, 10,β= 1, .2が正解 
set.seed(0)
dat<-c(rgamma(100,4,1),rgamma(300,10,.2))
hist(dat,breaks=20, main="",col="lightblue")

model<-gammamixEM(dat)
model$gamma.pars
model$lambda

まあまあ,て感じ。これ結構混ざっちゃうと,うまくいかないかもね。

混合ガンマ分布の確率密度関数は,

dmixgamma<-function(x,lambda,a1,a2,b1,b2){
 y<-lambda*dgamma(x,a1,b1)+(1-lambda)*dgamma(x,a2,b2)
 y
}

これで絵を描くとこう。

hist(dat,breaks=20, main="",col="lightblue",freq=F,ylim=c(0,.06),xlim=c(0,100))
x<-seq(0,100,.01)
lines(x,dmixgamma(x,0.25,3.05,11.54,1/1.23,1/4.37),lwd=2,col=2)


f:id:kusanagik:20170126152344p:plain



…これ,単変量で,しかも前もって2つの混合だってわかっているんだったら,EMアルゴリズムでなくても,この程度の単純さなら普通の(?)L-BFGSとかでいけそうだよね。
そしたら結局はfitdist,mle,bbmle,optimとかでも同じことできそう。
指数正規合成分布とガンマ分布の混合とか,そういうのも自由にできるはず。おもったより簡単ぽい。

単変量混合正規分布モデルをデータにフィットさせる

外国語教育研究では2つの山とか3つの山とかの分布になっているのを見ることがある。こういうときは,混合分布モデルをデータにフィットさせるといいかもだ。Rではmclustもいいけど,mixtoolsというパッケージがある。

#準備
library(mixtools)

#数値例の作成
#ここでは混合比(λ)が1:2,母平均(μ)が50, 120,母標準偏差(σ)が10, 20が正解
set.seed(92)
huta<-c(rnorm(100,50,10),rnorm(200,120,20))

#可視化して確認
par(mfrow=c(1,2))
hist(huta,col="lightgray",xlab="Value",main="")
plot(ecdf(huta))

f:id:kusanagik:20170124191930p:plain

#mixtoolsのEMアルゴリズム関数を使って母数を推定
fit<-normalmixEM(huta)
summary(fit)

summary of normalmixEM object:
           comp 1    comp 2
lambda   0.669647  0.330353
mu     117.657305 50.067067
sigma   20.090447  9.126634
loglik at estimate:  -1428.916 

#ここではだいたい正解だった
#対数尤度
fit$loglik
#AIC (自由母数は5個でいいよね?)
fitaic<--2*fit$loglik+2*5
fitaic

#誤差をブートストラップで
boot.se(fit)

#ヒストグラムに確率密度関数を合わせる
hist(huta,col="lightgray",xlab="Value",main="",freq=F,ylim=c(0,.015),breaks=20)

#その前に推定した母数をもつ混合正規分布の確率密度関数を定義する
dnormmix<-function(x){
y<-fit$lambda[1]*dnorm(x,fit$mu[1],fit$sigma[1])+fit$lambda[2]*dnorm(x,fit$mu[2],fit$sigma[2])
y
}

x<-seq(0,200,0.1)
lines(x,dnormmix(x),lwd=2,col=2)


f:id:kusanagik:20170124193111p:plain

うむ。間便な混合正規分布確率密度関数も書いておく。これで好きにお絵かきできる。

dnormmix<-function(x,lambda, mu1, mu2, sigma1, sigma2){
y<-lambda*dnorm(x,mu1,sigma1)+(1-lambda)*dnorm(x,mu2,sigma2)
y
}

(なぜか混合分布モデルは自分の青春だった,って感じがする)

任意の累積分布関数を仮定した一標本コルモゴロフ・スミルノフ検定

#正規分布の場合
set.seed(0)
dat<-rnorm(100,0,1)
ks.test(dat,"pnorm")

#特定の平均と標準偏差をもつ正規分布
set.seed(0)
dat<-rnorm(100,0,1)
ks.test(dat,"pnorm",1,1)

#ガンマ分布
set.seed(0)
dat<-rgamma(100,2,3)
ks.test(dat,"pgamma",2,3)

#ワイブル分布
set.seed(0)
dat<-rweibull(100,2,3)
ks.test(dat,"pweibull",2,3)

コルモゴロフ・スミルノフ検定で大事なのは,(よく正規性の検定に使われるけど)帰無仮説は任意の累積分布関数であって,その帰無仮説は自由に立てられる(一般に正規性の検定で使われる方法は,標本から推定した母平均と母標準偏差を項にもつ正規分布の累積分布関数を帰無仮説として設定している,ということ)。帰無仮説を棄却できないとき,「標本がその累積分布関数に適合していないとはいえない」だけであって,積極的に任意の累積分布関数から発生した標本だとか,母集団の分布はまさにこれだとか,ちょっと厳しくいえば正規分布に従っている,とか,そういうように語気強くいうべきではないかも。もちろん,正規性を示すひとつの証拠であることは間違いないけど。

自分で用意した分布をデータへフィットさせる

fitdistrplusパッケージのfitdist関数は,dnorm, pnormのように,dとpの関数が定義されていればどんなものでも指定することができる。自作でもいいってこと。ま,でも殆どの場合,自分程度が思いつくような分布はすでに用意されているってのがR。

たとえば,レイリー分布の場合。レイリー分布はVGAMパッケージに含まれている。レイリー分布のパラミタはひとつ(形状母数)。

#パッケージの読み込み
library(VGAM)

#例の作成
dat<-rrayleigh(1000,300)

#fitdist関数によるフィット
#初期値を指定してやらなければならない。とりあえずここではデータの中央値を入れてやる。
fit<-fitdist(dat,"rayleigh",start=list(scale=median(dat)))
summary(fit)

Fitting of the distribution ' rayleigh ' by maximum likelihood 
Parameters : 
      estimate Std. Error
scale 296.7956   4.692907
Loglikelihood:  -6652.58   AIC:  13307.16   BIC:  13312.07 

plot(fit)

f:id:kusanagik:20170117133218p:plain


よし。

じゃ,次は指数正規合成分布。これはretimesパッケージにある。ただし,dexgaussとrexgaussしかないので,pexgaussを自分で作ってしまう。

#パッケージ読み込み
library(retimes)

#pexgaussを作成
pexgauss<-function(x,mu,sigma,tau){
tau<-1/tau
u<-tau*(x-mu)
v<-tau*sigma
phi1<-pnorm(u,0,v)
phi2<-pnorm(u,v^2,v)
e<-exp(1)
p<-phi1-e^(-u+(v^2)/2+log(phi2))
p
}

#数値例の作成
dat<-rexgauss(1000,1000,500,1000)

#fitdist関数によるフィット
#初期値を指定してやらなければならない。とりあえずここではデータの正解を入れてやる。
fit<-fitdist(dat,"exgauss",start=list(mu=1000,sigma=500,tau=1000))
summary(fit)


ちなみに,retimesパッケージにはtimefitっていう専用の関数がある。ま,ほぼおなじ推定値が帰ってくる。

timefit(dat)


なんかもうoptimやりゃいいって気がしてくるな。

単変量の分布母数の推定からあれこれ: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


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

最尤推定した母数のもとでの確率密度曲線をヒストグラムに描き足す

えっと,リクエストがあったのでここに書く。(SPSSならそれっぽい曲線もつけてくれるのにRはそんなこともできないのか?といわれた)

まずはこんなデータがあるとしよう。形状母数が3,尺度母数1のガンマ分布にしたがう300個。

set.seed(0)
dat<-rgamma(300,3,1)

これをヒストグラムにする。

hist(dat,col="lightblue",main="")

f:id:kusanagik:20170111175015p:plain

まずはガンマ分布を仮定して,母数を最尤推定する。

library(MASS)
fitdistr(dat,densfun="gamma")

ま,3.01と1.03とかが推定されたと。
この母数のもとでの関数をヒストグラムに描き足したい。

まずはhist関数に,freq=Fを入れる。上に少し余裕を持たせるためにylimを指定しよう。

hist(dat,col="lightblue",main="",freq=F,ylim=c(0,.4))

f:id:kusanagik:20170111174956p:plain

すると縦軸が密度になる。

この上に関数を描けばいい。

たとえば,

x<-seq(0,9,.001)
lines(x,dgamma(x,3.01,1.03),col="blue")

f:id:kusanagik:20170111180126p:plain

はい。上手く当てはまっている感ある。

これで,別の関数も試してみたらいいかもだ。ワイブルと正規分布いこう。

#ここまでの作図
hist(dat,col="lightblue",main="",freq=F,ylim=c(0,.4))
x<-seq(0,9,.001)
lines(x,dgamma(x,3.01,1.03),col="blue")

#ワイブルからパラメーターを取る
w.shape<-fitdistr(dat,densfun="weibull")$estimate[[1]]
w.scale<-fitdistr(dat,densfun="weibull")$estimate[[2]]

#ワイブルを描き足す
lines(x,dweibull(x,w.shape,w.scale),col="red")

#正規分布からパラメーターを取る
mu<-fitdistr(dat,densfun="normal")$estimate[[1]]
sigma<-fitdistr(dat,densfun="normal")$estimate[[2]]

#正規分布を描き足す
lines(x,dnorm(x,mu,sigma),col="green")

f:id:kusanagik:20170111180744p:plain


緑がダメ,みたいな感じがわかる。

ちなみにAICを比較してみたり…

aics<-numeric(3)
aics[1]<-AIC(fitdistr(dat,densfun="gamma"))
aics[2]<-AIC(fitdistr(dat,densfun="weibull"))
aics[3]<-AIC(fitdistr(dat,densfun="normal"))
aics


ほい。

負の二項分布に関する2つのパラメーター化の方法

負の二項分布(negative binomial distribution)は,(a)その母数を成功回数(ないしサイズ母数)rと,成功確率pとする場合(こっちに親しみ)と,(b)サイズ母数rと平均μとするときの二種類があるんだそうだ。パラメーター化の方法が違うってのは,ガンマ分布にもそういうことあるよね(片方はなぜかベイジアンが好むというやつ)。

ただし,pは,r/(r+μ)なのですぐもとまる。

Rのnbinom関数系は両方イケる。

手を抜いてMASSパッケージのfitdistr関数に最尤推定をお願いすると,俄然(b)の方でもとめてくる。

忘れないようにメモ。