混合ガンマ分布モデルをデータにフィットさせる
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)
…これ,単変量で,しかも前もって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))
#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)
うむ。間便な混合正規分布の確率密度関数も書いておく。これで好きにお絵かきできる。
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)
よし。
じゃ,次は指数正規合成分布。これは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)
ほら,便利。これだけ。
左上は,ヒストグラムに推定された母数の値を項にもつ確率密度関数。
右上は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)
ブートストラップ
ブートストラップで誤差や信頼区間の検討もお手軽にできる。
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)
片方の母数を固定してもう片方の母数を最尤推定
やりたかったのこれ。ある理論的な見地から,片方の母数,例えば形状母数を任意の値に固定した状態で尺度母数を推定したいとか,これって結構外国語教育研究的によくあること。こういうときは,
#形状母数を無理やり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)
最尤推定以外の推定方法
最尤推定以外の方法も実装されている。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="")
まずはガンマ分布を仮定して,母数を最尤推定する。
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))
すると縦軸が密度になる。
この上に関数を描けばいい。
たとえば,
x<-seq(0,9,.001) lines(x,dgamma(x,3.01,1.03),col="blue")
はい。上手く当てはまっている感ある。
これで,別の関数も試してみたらいいかもだ。ワイブルと正規分布いこう。
#ここまでの作図 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")
緑がダメ,みたいな感じがわかる。
ちなみに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
ほい。