草薙の研究ログ

英語教育関係。でも最近は統計(特にR)ネタが中心。

MCMCを使って指数正規合成分布(ex-Gaussian)の母数を推定

RのMCMCpackにはMCMCmetrop1Rっていう関数があって,これは任意(自作)の対数尤度の関数をいれてMCMCでサンプリングすることができる。なので,結構手軽にMCMCを使ってデータに好きな分布を当てはめることが可能。

ここでは(まったくそんなことはしなくてもいいんだけど),指数正規合成分布の母数(μ,σ,τ)をMCMCを使って推定してみる。

#準備
library(retimes)
library(MCMCpack)

#データの例
dat<-rexgauss(1000,3000,100,800)

#関数の準備
llf<-function(beta,x){
sum(log(dexgauss(x,beta[1],beta[2],beta[3])))
}

#MCMCしてみる,初期値は全部適当,burninとかmcmcの数とかデフォルトのまま
m<-MCMCmetrop1R(llf,theta.init=c(2000,100,100),x=dat)

#結果を見てみる
summary(m)

Iterations = 501:20500
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 20000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

       Mean    SD Naive SE Time-series SE
[1,] 3005.7 13.00  0.09189         0.3102
[2,]  110.8 10.92  0.07724         0.2642
[3,]  779.5 27.79  0.19647         0.6677

2. Quantiles for each variable:

        2.5%    25%    50%    75%  97.5%
var1 2979.96 2997.0 3005.7 3014.3 3031.0
var2   90.49  103.5  110.5  117.8  133.7
var3  726.38  760.6  778.9  798.1  836.6

plot(m)


f:id:kusanagik:20170128185732p:plain