MCMC再訪(1)：Taglibro de H：SSブログ

	ブログをはじめるログイン

MCMC再訪(1)　[統計] [編集]

秋の某統計高座のテキスト作成のため、もう一度MCMCを最初からやってみる。

例題1:
正規分布することがわかっているある母集団から、X = (5.89, 5.69, 4.71, 6.73, 4.90, 3.34, 4.60, 4.00) という標本が得られたとき、その母集団の平均と標準偏差をベイズ推定する。

コード

##
## Sample 1
##

library(R2WinBUGS)

## data
X <- c(5.89, 5.69, 4.71, 6.73,
       4.90, 3.34, 4.60, 4.00)

## model
model <- function() {
  for (i in 1:N) {
    X[i] ~ dnorm(mu, tau);
  }
  
  ## priors
  mu ~ dnorm(0.0, 1.0E-6);
  tau ~ dgamma(1.0E-3, 1.0E-3);
  sigma <- 1/sqrt(tau);
}
write.model(model, "sample1.bug")

## Wine
WINE <- "/opt/local/bin/wine"
WINEPATH <- "/opt/local/bin/winepath"
bugs.dir <- paste(Sys.getenv("HOME"),
                  ".wine/drive_c/Program Files/WinBUGS14",
                  sep = "/")

## number of chains
n.chains <- 3

## initial values
init1 <- list(mu = -10, tau = 0.1)
init2 <- list(mu =   0, tau = 1)
init3 <- list(mu =  10, tau = 10)

## MCMC
post.samp <- bugs(data = list(X = X, N = length(X)),
                  inits = list(init1, init2, init3),
                  parameters.to.save = c("mu", "sigma"),
                  model = "sample1.bug",
                  n.chains = n.chains,
                  n.iter = 21000,
                  n.burnin = 1000,
                  n.thin = 10,
                  debug = FALSE,
                  bugs.directory = bugs.dir,
                  working.directory = getwd(),
                  clearWD = TRUE,
                  useWINE = TRUE,
                  newWINE = TRUE,
                  WINE = WINE,
                  WINEPATH = WINEPATH)

## show results
print(post.samp)
plot(post.samp)

## mcmc object
post.samp.mcmc <- mcmc.list(lapply(1:(post.samp$n.chain),
                                   function(chain)
                                   as.mcmc(post.samp$sims.array[,chain,])))
summary(post.samp.mcmc)
plot(post.samp.mcmc)

## Rhat (Gelman-Rubin's convergence diagnostic)
gelman.diag(post.samp.mcmc)

## mode
sapply(1:2, function(v) {
  br <- seq(min(post.samp$sims.array[,,v]),
            max(post.samp$sims.array[,,v]),
            length = length(post.samp$sims.array[,,v]) / 10)
  loc <- which.max(hist(post.samp$sims.array[,,v],
                        breaks = br,
                        plot = FALSE)$counts)
  hist(post.samp$sims.array[,,v], breaks = br,
       plot = FALSE)$mids[loc]
  }
)

解説

library(R2WinBUGS)

R2WingBUGSパッケージを使用する。

X <- c(5.89, 5.69, 4.71, 6.73,
       4.90, 3.34, 4.60, 4.00)

データを用意する。

model <- function() {
  for (i in 1:N) {
    X[i] ~ dnorm(mu, tau);
  }
  
  ## priors
  mu ~ dnorm(0.0, 1.0E-6);
  tau ~ dgamma(1.0E-3, 1.0E-3);
  sigma <- 1/sqrt(tau);
}
write.model(model, "sample1.bug")

モデルを用意して"sample1.bug"というファイル名でファイルに落としておく。BUGS言語のdnorm(mu, tau)関数では、Rのdnorm()とは異なり、引数はmu=平均、tau=分散の逆数であることに注意。モデルは、X[i] (i = 1 ..N)がNormal(mu, 1/sqrt(tau))に従うとし、muおよびtauの事前分布は無情報事前分布としている。

WINE <- "/opt/local/bin/wine"
WINEPATH <- "/opt/local/bin/winepath"

Wineを使うための設定。wineとwinepathのパスを指定する。Windowsから使うときは不要。

bugs.dir "- paste(Sys.getenv("HOME"),
                  ".wine/drive_c/Program Files/WinBUGS14",
                  sep = "/")

bugs.dirはWinBUGSをインストールしてあるディレクトリを指定する。

n.chains <- 3

マルコフ鎖の数を3とする。

init1 <- list(mu = -10, tau = 0.1)
init2 <- list(mu =   0, tau = 1)
init3 <- list(mu =  10, tau = 10)

初期値の設定。マルコフ鎖の数だけ用意する。後段のbugs()関数でinits = NULLとすると自動的に初期値を設定してくれるが、複雑なモデルの場合には、不適当な初期値が選ばれていきなり計算が失敗することもある。

post.samp <- bugs(data = list(X = X, N = length(X)),
                  inits = list(init1, init2, init3),
                  parameters.to.save = c("mu", "sigma"),
                  model = "sample1.bug",
                  n.chains = n.chains,
                  n.iter = 21000,
                  n.burnin = 1000,
                  n.thin = 10,
                  debug = FALSE,
                  bugs.directory = bugs.dir,
                  working.directory = getwd(),
                  clearWD = TRUE,
                  useWINE = TRUE,
                  newWINE = TRUE,
                  WINE = WINE,
                  WINEPATH = WINEPATH)

WinBUGSを呼び出す。

print(post.samp)

結果の表示

Inference for Bugs model at "sample1.bug", fit using WinBUGS,
 3 chains, each with 21000 iterations (first 1000 discarded), n.thin = 10
 n.sims = 6000 iterations saved
         mean  sd 2.5%  25%  50%  75% 97.5% Rhat n.eff
mu        5.0 0.5  4.1  4.7  5.0  5.3   5.9    1  4400
sigma     1.2 0.4  0.7  1.0  1.1  1.4   2.1    1  6000
deviance 25.2 2.3 23.0 23.6 24.5 26.1  31.5    1  6000

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = Dbar-Dhat)
pD = 2.2 and DIC = 27.4
DIC is an estimate of expected predictive error (lower deviance is better).

plot(post.samp)

結果を図示する。

post.samp.mcmc <- mcmc.list(lapply(1:(post.samp$n.chain),
                                   function(chain)
                                   as.mcmc(post.samp$sims.array[,chain,])))

結果をmcmcオブジェクトに変換する。

summary(post.samp.mcmc)

サマリーを表示する。


Iterations = 1:2000
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 2000 

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

           Mean     SD Naive SE Time-series SE
mu        4.988 0.4626 0.005972       0.006057
sigma     1.217 0.3820 0.004932       0.004982
deviance 25.237 2.2658 0.029251       0.027407

2. Quantiles for each variable:

           2.5%     25%    50%    75%  97.5%
mu        4.061  4.7150  4.992  5.269  5.891
sigma     0.721  0.9552  1.139  1.394  2.145
deviance 23.020 23.6200 24.550 26.110 31.490

plot(post.samp.mcmc)

mcmcオブジェクトのグラフ表示。軌跡と密度が表示される。

gelman.diag(post.samp.mcmc)

Gelman-Rubinの収束診断(Rhat)を求める。

Potential scale reduction factors:

         Point est. 97.5% quantile
mu                1           1.00
sigma             1           1.00
deviance          1           1.00

Multivariate psrf

1

sapply(1:2, function(v) {
  br <- seq(min(post.samp$sims.array[,,v]),
            max(post.samp$sims.array[,,v]),
            length = length(post.samp$sims.array[,,v]) / 10)
  loc <- which.max(hist(post.samp$sims.array[,,v],
                        breaks = br,
                        plot = FALSE)$counts)
  hist(post.samp$sims.array[,,v], breaks = br,
       plot = FALSE)$mids[loc]
  }
)

1–2番目のパラメータについて、だいたいの最頻値を求める。

[1] 4.9334574 0.9833632

タグ：MCMC winbugs

2009-08-27 05:53 nice!(0) コメント(8) トラックバック(1)
共通テーマ：学問

nice! 0

大変わかりやすい説明に感謝いたします。
ずいぶん前の記事への質問で失礼します。
5.89, 5.69, 4.71, 6.73, 4.90, 3.34, 4.60, 4.00
の平均とSDを普通に（エクセルで）求めると4.9825, 1.0868となります。これは最尤推定値と同じものだと思うのですが、こちらの計算では一番最後の部分の、1-2番目のパラメータについてのだいたいの最頻値がこれに相当するはず、という解釈であっていますでしょうか？　だいたいでなくもっと厳密に計算すると4.9825, 1.0868に近づいていくのでしょうか？
by Sunahara (2012-10-01 23:56)

完全な定常分布に収束して、事前分布がimproperな一様分布だとそうなるとおもいます。

ただ、ベイズ推定では、推定する対象は事後分布そのものなので、最頻値はあまり問題にしないような。
by hiroki (2012-10-02 05:27)

お返事ありがとうございます。
　実はある問題に対して単純なGLM、変量効果を加えたGLMM、変量効果が空間的に自己相関を示すモデル（CAR）でcoefficientsの値を比べて、予測モデルとしてはGLMで十分ですよと持って行きたいと思っているのです。私の知る限り空間自己相関モデルはベイズでないと解けませんが、最尤推定値と事後分布の平均値を比べるのは厳密にはおかしい。GLMからベイズでやる手ありますが、最終的に持って行きたい結論は「単純なモデルで簡単に予測できます」なのでなるべく簡単に見せたい。ではWinBUGSで最頻値を求める方法はないかと情報を探してこちらにたどり着きました。
　こちらのコードを参照して、また例題としてこちらの問題でMCMCのchainを長くしたりして最頻値が最尤推定値になるか試してみます。
by Sunahara (2012-10-02 11:11)

まあたしかに、GLMで片がつく問題をわざわざMCMC推定するようなことは通常は必要ないですね。
by hiroki (2012-10-02 20:01)

でも雑誌のレフェリーを納得させるためには通常必要ないこともやらざるをえない状況なのです。
by Sunahara (2012-10-02 22:11)

なんとそんなことが……。
by hiroki (2012-10-03 05:29)

あ、具体的にこういうことを要求されているわけではないのですが、こちらがGLMで十分と思っていることに対して空間自己相関を入れろと言われ、それは必要ないとわかってもらうための根拠が必要なのです。
by Sunahara (2012-10-03 11:40)

ああ、空間構造のあるデータでしたら、空間自己相関を入れたモデルを求められる可能性がありますね。
by hiroki (2012-10-03 20:36)