サイコロ問題

MikuHatsune2015-05-24

推定の章
 
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のコード部分は別ファイルに記述して、テキストとして読み込むのがかっこいいスタイルらしいが、コピペで一発作動してほしいからベタに書いた。
推定するべきはサイコロの各目が出る確率pが、長さ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が出やすそうな確率がちょっとだけ上がる。
これを上記の条件っぽく拡張して書いたのだが、pの確率で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の数値とそれっぽく似る。
 
車両の保有台数がべき乗分布に比例 PMF(x)\propto (\frac{1}{x})^\alpha する、とすると、上限値の設定による誤差がかなり解消される。

# 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

となり、なんとなく幅をもってもやっと推定ができる。