推定の章
4, 6, 8, 12, 20面サイコロがひとつずつあって、どれを振ったかわからないが1回振って6が出たらしい。
このとき、N面サイコロが選ばれる確率を考える。
side <- c(4, 6, 8, 12, 20) n <- 6 p <- ifelse(side < n, 0, 1/side) ps <- post_p <- p/sum(p)
[1] 0.0000000 0.3921569 0.2941176 0.1960784 0.1176471
となって結果は一致する。
その後、何回か繰り返して6,8,7,7,5,4が出たとすると
trial <- c(6, 8, 7, 7, 5, 4) for(roll in trial){ p <- ifelse(side < roll, 0, 1/side) p <- post_p*p post_p <- p/sum(p) ps <- rbind(ps, post_p) } colnames(ps) <- side cols <- rev(grey(seq(0.1, 1, length=nrow(ps)))) matplot(t(ps), type="l", lwd=3, lty=1, col=cols, xlab="Side of die", ylab="Probability", ylim=c(0, 1), xaxt="n") axis(1, at=seq(side), labels=side)
post_p
[1] 0.000000000 0.000000000 0.943248454 0.055206128 0.001545418
となって8面サイコロが振られた確率が高くなる。
これをrstanでやってみよう。とりあえず、6面サイコロを1回振って6が出たという簡単なシミュレーションを行う。
本来ならstanのコード部分は別ファイルに記述して、テキストとして読み込むのがかっこいいスタイルらしいが、コピペで一発作動してほしいからベタに書いた。
推定するべきはサイコロの各目が出る確率が、長さNの単体 simplex に収まっていて、これにしたがって多項分布 categorical からサンプリングされる。
library(rstan) library(coda) code1 <- ' data { int<lower = 1> N; int<lower = 1> die; } parameters { simplex[N] p; } model { die ~ categorical(p); } ' # 入力データ p1 <- list(N=6, die=6) # 推定 fit <- stan(model_code = code1, data = p1, iter = 1000, chains = 4)
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat p[1] 0.14 0.00 0.12 0.00 0.05 0.11 0.21 0.45 1477 1 p[2] 0.14 0.00 0.12 0.01 0.04 0.11 0.21 0.45 1310 1 p[3] 0.15 0.00 0.12 0.00 0.05 0.11 0.21 0.47 1538 1 p[4] 0.14 0.00 0.12 0.00 0.05 0.11 0.21 0.44 1626 1 p[5] 0.14 0.00 0.12 0.00 0.05 0.11 0.20 0.46 1475 1 p[6] 0.28 0.00 0.16 0.04 0.16 0.26 0.39 0.63 1545 1 lp__ -15.23 0.08 1.96 -20.34 -16.11 -14.81 -13.80 -12.78 539 1
6面サイコロを高々1回振って6が出たので、6が出やすそうな確率がちょっとだけ上がる。
これを上記の条件っぽく拡張して書いたのだが、の確率でN面サイコロを選んで、という条件分岐がうまく書けない。
# これはエラー code1 <- ' data { int<lower = 1> N; int<lower = 1> die; int side[N]; } parameters { simplex[N] p; } transformed parameters { simplex[max(side)] q[N]; for(i in 1:N) for(j in 1:max(side)) q[i, j] <- if_else(j <= side[i], 1.0/side[i], 0); } model { die ~ categorical(q[2]); } ' # 入力データ p1 <- list(N=length(dice), side=side, die=6) # 推定 fit <- stan(model_code = code1, data = p1, iter = 1000, chains = 4)
というわけで諦めている。increment_log_prob を使えと神からのお告げがあったが保留。
機関車問題に移る。
1000台の車両を持っている会社のうち、60番というナンバリングの車両を見かけたとしたら、この会社が保有している車両の数は
hypos <- seq(1000) N <- 60 p <- ifelse(hypos < N, 0, 1/hypos) post_p <- p1 <- p/sum(p) plot(hypos, post_p, type="l", lwd=3, xlab="No. trains", ylab="Probability") sum(hypos*post_p) # 事後確率分布の平均
[1] 333.4199
となって計算上は合う。
その後、30番目、90番目の車両を見たら、確率は更新されて
# three trains N <- 60 p <- ifelse(hypos < N, 0, 1/hypos) post_p <- p/sum(p) Ns <- c(30, 90) for(N in Ns){ p <- ifelse(hypos < N, 0, 1/hypos) post_p <- post_p*p post_p <- post_p/sum(post_p) } sum(hypos*post_p) # 事後確率分布の平均
[1] 164.3056
となる。上限を500もしくは2000にしてもThinkBayesの数値とそれっぽく似る。
車両の保有台数がべき乗分布に比例 する、とすると、上限値の設定による誤差がかなり解消される。
# power law alpha <- 1 prior_dis <- hypos^(-alpha)/sum(hypos^(-alpha)) N <- 60 p <- ifelse(hypos < N, 0, 1/hypos)*prior_dis post_p <- p2 <- p/sum(p) Ns <- c(30, 90) for(N in Ns){ p <- ifelse(hypos < N, 0, 1/hypos) p <- p/sum(p) post_p <- post_p*p post_p <- post_p/sum(post_p) } sum(hypos*post_p) # 事後確率分布の平均
[1] 133.2752
matplot(cbind(p1, p2), type="l", lwd=3, lty=1, xlab="No. trains", ylab="Probability") legend("topright", legend=c("uniform", "power law"), lwd=3, lty=1, col=1:2)
信用区間を求めると
mapply(function(x) head(which(cumsum(post_p)>x), 1), c(0.05, 0.95))
[1] 91 242
となり、なんとなく幅をもってもやっと推定ができる。