So-net無料ブログ作成

MCMCの勉強(1) [統計]

今さら感はあるが、MCMC (Markov Chain Monte Carlo; マルコフ連鎖モンテカルロ)を使えるようになろうと、まずは簡単な例から試してみた。


手始めに、正規乱数から生成した標本の平均と標準偏差を推定してみる。

やはりRを使用。MCMCpackパッケージを あらかじめインストールしておいて、呼び出す。MCMCpack中のMCMCmetrop1R()関数を利用して、メトロポリス法によるMCMC推定をおこなう。

library(MCMCpack)

乱数系列を初期化。

set.seed(1)

平均10、標準偏差3の乱数を1000個生成して、xに入れる。

m <- 10
s <- 3
x <- rnorm(1000, m, s)
MCMC推定に使用する関数を用意する。betaは要素数2のベクトル。beta[1]が平均、beta[2]が標準偏差で、betaを推定する。関数の返り値は対数尤度。標準偏差は非負なので、負になったときは対数尤度をマイナス無限大とする。
llnormfun <- function(beta, x) {
	if (beta[2] < 0.0)
		l <- -Inf
	else
		l <- sum(log(dnorm(x, mean=beta[1], sd=beta[2])))
	return(l)
}
推定値の初期値を平均0、標準偏差6とし、繰り返し10000回、初期値依存として捨てる数を500として、MCMC推定を実行。
post.samp <- MCMCmetrop1R(llnormfun, theta.init=c(0,6), 
	x=x, 
	thin=1, mcmc=10000, burnin=500, 
	tune=1.0, 
	verbose=500, logfun=TRUE, optim.maxit=100)
結果をグラフで表示する。
plot(post.samp)

結果を数値で表示。
summary(post.samp)
結果は以下の通り。
Iterations = 501:10500
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

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

      Mean      SD  Naive SE Time-series SE
[1,] 9.967 0.09639 0.0009639       0.003111
[2,] 3.109 0.06957 0.0006957       0.002424

2. Quantiles for each variable:

      2.5%   25%   50%    75%  97.5%
var1 9.776 9.902 9.970 10.033 10.157
var2 2.974 3.061 3.108  3.157  3.248
実際の値を表示する。
> c(mean(x), sd(x))
結果
[1] 9.965056 3.104748
標準偏差の初期値を小さくすると、エラーになる。
post.samp <- MCMCmetrop1R(llnormfun, theta.init=c(0,1), 
 	x=x, 
 	thin=1, mcmc=10000, burnin=500, 
 	tune=1.0, 
 	verbose=500, logfun=TRUE, optim.maxit=100)
結果
Warning message:
Mode and Hessian were not found with call to optim().
Sampling proceeded anyway. 
 in: MCMCmetrop1R(llnormfun, theta.init = c(0, 1), x = x, thin = 1,  
Hessian from call to optim() not negative definite.
Sampling (as specified) cannot proceed.
エラー:Check data and fun() and call MCMCmetrop1R() again. 
tuneの値を小さくすると、初期値の影響が長く残る。
post.samp <- MCMCmetrop1R(llnormfun, theta.init=c(0,6), 
	x=x, 
	thin=1, mcmc=10000, burnin=500, 
	tune=0.1, 
	verbose=500, logfun=TRUE, optim.maxit=100)

plot(post.samp)

このあたりは 計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 に解説があったところ。
タグ:MCMC
nice!(0)  コメント(1)  トラックバック(2) 

nice! 0

コメント 1

ととろ

MCMCPackを動かしている日本語サイトを
初めてみました。
同じ事をやらせてみたのですが、案外と
早いのに、びっくり。
自分にも使えるかも?と言う気分になり
ました。
by ととろ (2006-08-19 23:39) 

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

トラックバック 2