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)