「おまえは今まで買ったチョコボールの数をおぼえているのか?」
……こんなことを聞かれたら,あなたならどう答えますか?「覚えとらんわ」と一蹴するのは簡単ですが,答えられないのはちょっと悔しい。という訳で,この記事では僕が今までに買ったチョコボールの箱数の推定を試みます。
※この記事は,Stan Advent Calendar 2019の3日目のエントリー記事です。
1. チョコボールとエンゼルについて
改めまして,こんにちは。mutopsyです。冒頭でも説明した通り,この記事では僕が今までに買ったチョコボールの箱数を推定することを目指します。
チョコボールは森永製菓から発売されている国民的チョコレート菓子の1つです。チョコボールの箱の「くちばし」と呼ばれる開口部には,ごくまれにエンゼルのイラストが印刷されており,銀のエンゼルなら5枚,金のエンゼルなら1枚で「おもちゃのカンヅメ」と交換することができます。金のエンゼルと銀のエンゼルが当たる確率は企業秘密とのことですが,以前12,283箱分のチョコボールデータを用いて推定してみたところ,銀のエンゼルが当たる確率は約3.60% (3.13〜4.10%程度),金のエンゼルが当たる確率は約0.05% (0.02〜0.12%程度) であると推定されました(「たのしいベイズモデリング2:事例で拓く研究のフロンティア」の第1章;過去の記事も参照:[その1] [その2])。今回の分析では,銀のエンゼルが当たる確率についての情報を利用します。
1. 利用可能な情報とモデル
当然ながら,今までに買ったチョコボールの箱数を推定するためには何らかの情報が必要です。今回僕が注目したのは,今までに当てた銀のエンゼルの数です。これまでに買ったチョコボールの総数そのものはさすがに覚えていませんが,銀のエンゼルが当たるというイベントはそれなりに記憶に残る事象であり,かつ生起頻度も少ないので,それぐらいなら覚えているのも難しくないでしょう。さらにいえば,僕は以前に銀のエンゼルを5枚集めておもちゃのカンヅメと交換してもらったこともあり,そしておそらくそれ以来銀のエンゼルを見かけていないので,これまでに当てた銀のエンゼル数が5枚であることがはっきりしています。そういう訳で,今回の分析では「今までに当てた銀のエンゼルの数は5枚」という情報と「銀のエンゼルが当たる確率は約3.60%」という情報を使って今までに当てたチョコボールの箱数の推定を試みます。ちなみに金のエンゼルを当てたことはありませんが,金のエンゼルが当たる確率は非常に低いため,「今までに当てた金のエンゼルの数は0枚」という情報はあまり役に立ちません。簡単のためにも今回の分析では金のエンゼルの情報については用いないことにします(反対に,「金のエンゼルをn (≧1) 枚当てたことがある」という情報は有用なので積極的に利用すべきです)。
銀のエンゼルが当たる確率は約3.60%と述べましたが,先に紹介した書籍ではベイズ推定を用いて銀のエンゼルが当たる確率θsilverの確率分布 (事後分布)を計算しました (3.60%という値は事後分布のMAP推定値)。せっかくなので今回は点推定値ではなく分布の持つ情報をフル活用しましょう。この確率分布は以下のベータ分布とほぼ一致します (※書籍の分析では弱情報事前分布を用いたので厳密には一致しませんがほとんど無視できる程度の差です)。
この分布は,9755 (= 351 + 9404) 箱のチョコボールから351枚の銀のエンゼルを見つけたときに二項分布モデルを用いて推定したθsilverの分布に相当します (一様事前分布を想定した場合)。チョコボールをν箱開封したときに銀のエンゼルがNsilver枚得られる確率は,二項分布を用いて
と表現することが可能です。このνが今回推定したいパラメータですね。ただし,νは自然数をとる離散パラメータなので,離散パラメータを直接推定することができないStanを使う場合は一工夫必要になります。そこでこの記事では,まず始めに離散パラメータを推定することができるJAGSを用いた推定法を紹介し,続いてStanで同様の分析を行う方法を説明します。そして最後に,数値積分を用いて正確な事後分布を得る方法も簡単に紹介します。ちなみに今回はNsilver = 5の場合の推定結果を示しますが,他の値を使うことももちろんできますし,Nsilverを確率分布で表現できるようにモデルを拡張することも可能です (例えば,「正確には覚えていないが3〜5枚のいずれかのはず」といった記憶の曖昧さも表現できます)。
3-1. JAGSによる推定
JAGSはStanと違って離散パラメータを直接推定することができるので,上記のモデルをそのままコードに起こすだけで分析を実行することができます。分析に使用するJAGSコードは以下の通り。
angels_postmodel_jags.txt
model{
#モデル
N_silver ~ dbinom(theta_silver, nu)
#銀のエンゼルが当たる確率の事前分布
theta_silver ~ dbeta(alpha, beta)
#今までに買ったチョコボールの数の事前分布(離散一様分布)
nu ~ dcat(p_each[])
for (m in 1:Max_nu){
p_each[m] <- 1/Max_nu
}
}
ベータ分布と二項分布を用いて上記のモデルをそのまま書いています。Stanとは異なり変数の宣言は不要です。νの事前分布として十分幅の広い離散一様分布を指定しています。今回は上限 (Max_nu) を1000とします (データとしてRから渡す)。以下のRコードを実行すれば走ります。
JAGSで分析するためのRコード
#以下の環境で動作確認済み
#JAGS4.3.0, R2jags0.5-7, R3.5.2, Windows 10 Pro
#JAGSを下記URLからダウンロードしインストールする必要があります。
# http://mcmc-jags.sourceforge.net/
library(rjags)
library(R2jags)
#銀のエンゼルが当たる確率の事後分布 (ベータ分布) のパラメータ
alpha <- 351 + 1
beta <- 9404 + 1
#今回の分析に必要な情報
N_silver = 5 #今までに当てた銀のエンゼルの数
Max_nu = 1000 #今までに買ったチョコボールの数の上限 (十分大きな値にしておく)
#データの作成
data.post <- list(N_silver = N_silver, Max_nu = Max_nu,
alpha = alpha, beta = beta)
#モデルの読み込みとサンプリング
fit.jags <- jags(data=data.post, model.file="angels_postmodel_jags.txt", parameters.to.save=c("nu"),
n.chains=4, n.iter=10500, n.burnin=500, n.thin = 1, DIC=TRUE)
#MAPやHDIといった要約統計量の推定精度がどのくらい必要かに応じてn.iterを変える。
#n.iter=250500ぐらいまで増やせばそれなりの精度になります (時間はかかりますが)。
maphdi(fit.jags) #MAPとHDIを計算するための自作関数;過去記事参照
#結果は以下の通り:
# MAP Lower95 Upper95 Lower99 Upper99
#deviance 3.531314 3.438523 7.376154 3.438523 10.26979
#nu 141.986992 46.000000 296.000000 36.000000 369.00000
MAPとHDIの計算部分は以前の記事で紹介した自作関数を使用しています。νの事後分布を可視化すると以下のようになります。

ということで,MAP推定値で見ると,僕はどうやら今まで約139箱のチョコボールを購入したようです。体感的にもそんなもんかなといったところです。ただし事後分布の幅は広めなので注意が必要です。また,iterationの数が少ないと点推定値はかなりぶれます。ちなみに今までに買ったチョコボールの数が100個を超える確率を計算したところ約84%となりました。
3-2. Stanによる推定
前節ではJAGSによる推定方法を紹介しましたが,今回のモデルの場合は実は負の二項分布からの乱数を発生させるだけでνの事後分布を得ることが可能です。Stanを使うまでもないのですが,例えば銀のエンゼルが当たる確率をデータから推定した上でその値を使ってνを推定したい,という場合は全部Stanで書いた方が楽なので,そういった場合を念頭に置いてStanによる推定法を紹介します。今回のモデルでは
という二項分布を用いていますが,Nsilverが定数でνが確率変数のとき,νの事後分布は負の二項分布を用いて
と表現することができます (負の二項分布は,成功率がpのときにr回成功するまでに要する試行数が従う分布です)。したがって,既知のNsilverとθsilverを用いれば,MCMCするまでもなくνの事後分布を得ることができます。ただし,Stanのneg_binomial関数では2番目の引数が成功確率ではなくオッズ (成功確率と失敗確率の比) になっているので書き方を工夫する必要があります。Stanコードは以下の通りです。
angels_postmodel_stan.stan
data{
int alpha;
int beta;
int N_silver;
}
parameters{
real theta_silver;
}
model{
theta_silver ~ beta(alpha,beta);
}
generated quantities{
int nu;
nu = neg_binomial_rng(N_silver+1,theta_silver/(1-theta_silver)) + N_silver;
}
モデル記述部ではθsilverをベータ分布から発生させているだけです。生成量として負の二項分布から乱数を発生させればνの事後分布を得ることができます。上記コードを実行するためのRコードは以下の通り。ただし前処理の部分はJAGSコードと同じなので省略しています。
Stanで分析するためのRコード
#以下の環境で動作確認済み
#rstan2.18.2, Rtools3.5, R3.5.2, Windows 10 Pro
library(rstan)
... #JAGSのときと同じ前処理
stanmodel_post <- stan_model(file="angels_postmodel_stan.stan")
fit.stan <- sampling(stanmodel_post, data=data.post, seed=123,
chains=4, iter=250500, warmup=500, thin=1
) #JAGSよりも速いのでiterを多めにしている。
round(maphdi(fit.stan),2)
#結果は以下の通り:
# MAP Lower95 Upper95 Lower99 Upper99
#theta_silver 0.04 0.03 0.04 0.03 0.04
#nu 135.62 46.00 296.00 32.00 367.00
#lp__ -1514.99 -1516.88 -1514.96 -1518.29 -1514.96
そしてνの事後分布は以下の通り。

多少誤差はありますが,JAGSの結果とほぼ一致していることが分かります。
3-3. 数値積分による推定
これまでは,二項分布のパラメータνを推定する方法として(1)JAGSで離散パラメータとして推定する方法と(2)Stanを使って (あるいは使わなくても) 負の二項分布から乱数を発生させる方法を紹介しました。ただし,これらの方法は乱数を用いた方法のためiterationが小さいと(少ないと)推定結果が安定せず,たとえばMAP推定値を1の位まで正確に求めるためにはiterationをかなり大きく設定する必要があります。そこでこの節では,乱数に頼らずに事後分布を数値積分で求める方法を簡単に紹介します。いろいろなやり方があると思いますが,ここでは高校数学で習う区分求積法を用いたコードを紹介します。Rコードは以下の通り。
数値積分で推定するためのRコード
Max_nu = 1000
#θがベータ分布に従うときの二項分布モデル(r|ν,θ)の,nの尤度を返す関数を定義 (区分求積法でθを消去)
cal.lik.nu <- function(nu,r,alpha,beta,bin.n=2000){
output.lik.nu<-0
bin.w <- 1/bin.n
for (x in 0:(bin.n-1)){
p.t <- x/bin.n+bin.w/2
output.lik.nu <- output.lik.nu + bin.w * dbinom(r,nu,p.t) * dbeta(p.t,alpha,beta)
}
return(output.lik.nu)
}
#全てのnについて尤度を計算しベクトルとして格納
lik.nu <- 1:Max_nu
for(n in 1:Max_nu){
lik.nu[n] <- cal.lik.nu(n,N_Silver,prior_pars$alpha,prior_pars$beta)
}
#尤度を正規化して確率に直す (nの事前分布に離散一様分布を仮定した場合)
prob.nu <- lik.nu / sum(lik.nu)
#MAPとHDIの計算と格納 ※参照:http://bayesmax.sblo.jp/article/185297395.html
accum<-1:length(prob.nu)
for(i in 1:length(prob.nu)){
if (i ==1){
accum[i] <- sort(prob.nu)[i]
} else{
accum[i] <- accum[i-1] + sort(prob.nu)[i]
}
}
Interval95 <- grep(TRUE,prob.nu>max(sort(prob.nu)[accum<=(min(abs(accum-0.05))+0.05)]))
Interval99 <- grep(TRUE,prob.nu>max(sort(prob.nu)[accum<=(min(abs(accum-0.01))+0.01)]))
HDI95 <- c(min(Interval95),max(Interval95))
HDI99 <- c(min(Interval99),max(Interval99))
Result.intg <- data.frame(MAP=which.max(prob.nu),Lower99=HDI99[1],
Lower95=HDI95[1],Upper95=HDI95[2],Upper99=HDI99[2],
row.names = "n_total")
#MAPとHDIの表示
print(Result.intg) #MCMCサンプルから計算した値よりも正確
# MAP Lower99 Lower95 Upper95 Upper99
#n_total 138 34 50 300 369
#今までに買ったチョコボールの数が100を超える確率
print(sum(prob.nu[101:length(prob.nu)])) #MCMCサンプルから計算した値よりも正確
#[1] 0.8419658
事後分布は次の通り。

この方法を使うと高速かつ正確に事後分布を求めることができます。iterationをものすごく増やせばJAGSやStanでも同じ結果が得られます。モデルによってはMCMCより数値積分したり解析的に解いたりしたほうが楽な場合もあるということですね。
4. まとめ
という訳で,僕が今までに買ったチョコボールの数を推定することにぶじ成功しました。離散パラメータを推定するためにJAGSを使うもよし,Stanで乱数を発生させるのもよし,数値積分に頼るのもよし。手段はどうあれ目的は達成できましたね。これでもう大丈夫。「おまえは今まで買ったチョコボールの数をおぼえているのか?」と聞かれたら,僕ならこう答えます。
「おぼえていないが どうやら50〜300箱のようだ──ッ」
※この記事の内容はもともと「たのしいベイズモデリング2」の第1章に含めるつもりでいたのですが,紙面の都合で泣く泣くカットせざるを得なかったためアドカレのネタとして使い回すことにしました。という裏話。