宜しくお願いします。

針が上になっている状態を表、下になっている状態を裏とする。

これはRを用いて計算が可能である。

X.raw <- c(1,0,1,0,0,1,1,1,0,1,0,1,1,0,1,0,0,0,1,1,0,0,1,1,1,0,0,1,0,0,1,0)

n <- length(X.raw)
p.hat <- sum(X.raw)/n
var.p.hat <- p.hat*(1-p.hat)/n
cat("\n", "Point Estimate of p:", p.hat, "\n", 
    "Sample Variance Estimate of p:", var.p.hat, "\n", 
    '95% Confidence Interval:', p.hat-pnorm(1-0.25)*sqrt(var.p.hat), 
    "<=p<=", p.hat-pnorm(0.25)*sqrt(var.p.hat), "\n", "\n")

ベイズ法に関する計算、特に事後分布の算出は、WinBUGSというソフトウェアを用いるのが簡単である。

DATA list(n=32, x=c(1,0,1,0,0,1,1,1,0,1,0,1,1,0,1,0,0,0,1,1,0,0,1,1,1,0,0,1,0,0,1,0))

モデルと事前分布と初期値の指定。

MODEL model{
  for ( i in 1:n ) {
    x[i] ~ dbin(p,1)
  }
  p ~ dbeta(2,3) # objective prior 
  # p ~ dbeta(1,1)  # flat prior
}
INIT list {p=0.5}

ここでベータ分布はパラメータを変えると、次のように動く。

ここから先の手順は、多少複雑でこの手引きに頼るとよい。手順だけ記すと、次の通り。

1. 事後分布を得るためのWinBUGSのコードを読み込ませる (〔Model〕→〔Specification〕)
2. モデルを読み込ませる (上で書いたコードのうちmodelをマウスで指定してから、〔Specification Tool〕→〔check model〕、成功すると ``model is syntactically correctと表示される) - データを読み込ませる (上で書いたコードのうち、データの記述に相当するlistをマウスで指定してから、〔Specification Tool〕→〔load data〕、成功すると ``data loaded)
3. 走らせるマルコフ連鎖の本数を指定する (〔Specification Tool〕→〔num of chains〕)
4. コンパイルする (〔Specification Tool〕→〔compile〕、成功すると ``model compiled) - 初期値を指定する場合、読み込む (〔Specification Tool〕→〔load inits〕、マルコフ連鎖の本数だけ実行する必要あり、不足分があると ``this chain contains uninitialized variables.)
5. 初期値が不足している分は生成させる (〔Specification Tool〕→〔gen inits〕)
6. 事後分布を得るためにマルコフ連鎖モンテカルロ(MCMC)を走らせる準備をする (〔Model〕→〔Update〕、〔Inference〕→〔Sample〕
7. モニターしたい変数をSample Monitor Toolに登録する (〔Sample Monitor Tool〕→〔node〕にすべての変数を、ひとつひとつ手で入力し、〔set〕)
8. まずは1000歩、走らせてみる (〔Update Tool〕でupdatesに1000とタイプしてから、〔update〕) これはバーンインと呼ぶ準備である。多次元空間を一度に一次元方向ずつだけ移動する点の列なので、初期値の影響がなくなるぐらいまでまず走らせる (変数の数にもよるが、1000程度以上か)
9. 次にまた9000歩、走らせてみてから、historyを眺めて、定常状態になったことを確認する。定常分布に収束できるとは限らないので、複数のマルコフ連鎖を走らせるほうが良い (chains>=2)、 各変数の履歴(history)を見て、定常分布に至っていることを確認する (ホワイトノイズに見えることが望ましい) などの注意を憶えておく。
10. 本格運転 (〔Update Tool〕でマルコフ連鎖の歩を進め、〔Sample Monitor Tool〕で様子を眺める)
11. 十分な長さだけ、マルコフ連鎖モンテカルロを走らせる ($50000$, $100000$, などともかく十分な長さ)
12. 履歴(history)、直近の履歴(trace)、密度関数(density)、分位点の推移(quantiles)、などを眺めて、定常分布となっていることを確認 (〔Sample Monitor Tool〕の〔node〕欄に ``*'' を入力してから、各ボタンをクリックするのが便利) し、自己相関(auto cor)も確認し(妙なパターンがなく、滑らかに減少していれば良い、速く減少していればもっと良いがMCMCなので仕方ない)てから、事後分布を拾う
13. 定常と判断した範囲(beg〜end)までについてstatsを表示させる (〔Sample Monitor Tool〕でbegの値を1ではなく、定常と判断した歩数に変更)

点推定は、モーメント法・最尤法と同じ。データの与え方も同じ。下記のコードで、ブートストラップによる分散の推定と、信頼区間の構成まで行える。

n <- length(X.raw)
B <- 100000
p.boot <- NULL
for( b in c(1:B) ) {
  X.boot <- sample(X.raw,size=n, replace=TRUE)
  p.boot <- append(p.boot, sum(X.boot)/n)
}
hist(p.boot)
quantile(p.boot, probs=c(0.025, 0.975))
cat("\n", "Point Estimate of p:", p.hat, "\n", 
    "Bootstrap Variance Estimate of p:", var(p.boot), "\n", 
    '95% Bootstrap Confidence Interval:', quantile(p.boot, probs=0.025), 
    "<=p<=", quantile(p.boot, probs=0.975), "\n", "\n")

これは点推定の方法は正しくても、分布の仮定がそれほど正しくないかもしれない時に、一番最初の信頼区間とは異なった範囲を示す。

ただし、ブートスラップ反復回数 B の設定には要注意。

# Beta(1,1)
plot(c(1:99)/100,dbeta(c(1:99)/100,1,1), xlab="p", ylab="prior probability", type="l", lty=1, xlim=c(0,1), ylim=c(0,4.5))
# Beta(3,2)
lines(c(1:99)/100,dbeta(c(1:99)/100,3,2), xlab="p", ylab="prior probability", type="l", lty=2)
# Beta(2,3)
lines(c(1:99)/100,dbeta(c(1:99)/100,2,3), xlab="p", ylab="prior probability", type="l", lty=3)
# Beta(10,5)
lines(c(1:99)/100,dbeta(c(1:99)/100,10,5), xlab="p", ylab="prior probability", type="l", lty=4)
# Beta(5,10)
lines(c(1:99)/100,dbeta(c(1:99)/100,5,10), xlab="p", ylab="prior probability", type="l", lty=5)
# Beta(15,5)
lines(c(1:99)/100,dbeta(c(1:99)/100,15,5), xlab="p", ylab="prior probability", type="l", lty=6)
# Beta(5,15)
lines(c(1:99)/100,dbeta(c(1:99)/100,5,15), xlab="p", ylab="prior probability", type="l", lty=7)