ベイズマックスへようこそ。
このブログでは,心理学と統計が大好きなmutopsyが,
統計に関するあれやこれやを自由気ままに投稿します。

2019年01月09日

二項分布モデルと幾何分布モデルで推定した確率パラメータの事後分布の比較

 こんにちは。mutopsyです。この記事では,同一のカウントデータを二項分布でモデリングした場合と幾何分布でモデリングしたときに,推定される確率パラメータの事後分布が一致することを説明します。間違いがあったらご指摘ください。

 次の2つの例を考えてみましょう。

  1. N箱のチョコボールを購入したら銀のエンゼルが1枚当たった。銀のエンゼルが当たる確率θ1を推定せよ。
  2. 銀のエンゼルが当たるまでチョコボールを購入したらN箱目で当たった。銀のエンゼルが当たる確率θ2を推定せよ。

θ1θ2はどちらも試行にかかわらず一定の値をとり,かつ各試行は独立であると仮定します。ここではチョコボールを例に挙げていますが,2つの事象 (e.g., 当たるか外れるか) のカウントデータであれば何でも構いません。どちらの文も,「銀のエンゼルがN箱中1箱当たった」という同じ結果を述べていますが,データの収集方法に違いがあります。1番目のケースではNがあらかじめ決まっていた状況で銀のエンゼルが1回当たったという結果が得られていますが,2番目のケースでは銀のエンゼルが1回当たることがあらかじめ決まっていた状況でN箱を要したという結果が得られているという違いです。この違いは確率パラメータの推定結果に影響を与えるのでしょうか。

 この例のようにNの決め方 (決まり方) が異なるケースでは,異なる統計モデルを想定するのが自然です。1番目の例のように,Nが所与のときに銀のエンゼルが当たった枚数 (K = 1) のデータが得られた場合には,二項分布を用いて

1 ~ Binomial(N, θ1)

と表現することができます。一方,2番目の例のように,銀のエンゼルが当たった枚数 (K = 1) が所与のときにNというデータが得られた場合には,幾何分布を用いて

N ~ Geometric(θ2)

と表現できます。幾何分布は,同一の確率パラメータを持つベルヌーイ分布に従う試行を独立に繰り返し行ったときに初めて成功する (今回の例では銀のエンゼルが当たる) のに要する回数の分布です。二項分布では銀のエンゼルが当たった回数 (K = 1) が説明されるのに対し,幾何分布では銀のエンゼルが当たるまでに要した試行回数 (N) が説明されるのがポイントです。

 では,2つのモデルでθ1θ2の推定結果がどうなるのかを確認してみましょう。今回は,N = {10, 20, 50, 100} の場合のθ1θ2をStanでベイズ推定してみます。4通りのN と2つのモデル (Is_Censored = 0なら二項分布モデル,Is_Censored = 1なら幾何分布モデル) の組み合わせを用いて8個のθをまとめて推定します。Stanコードは次の通り。rstan2.18.2,Rtools3.5,R3.5.2,Windows10で実行しました。

Stanコード01: 二項分布モデルと幾何分布モデルの比較
data {
  int N_ID; //標本の数 = 8
  int N[N_ID]; //試行回数 = {10,20,50,100,10,20,50,100}
  int K[N_ID]; //成功回数 (今回は1で固定;2以上の場合は負の二項分布になる)
  int Is_Censored[N_ID]; //打ち切られたデータか否か = {0,0,0,0,1,1,1,1}
}

parameters{
  real<lower=0,upper=1> theta[N_ID]; //当たる確率
}

model{
  for (i in 1:N_ID){
    if (Is_Censored[i] == 0){
      K[i] ~ binomial(N[i],theta[i]); //二項分布
    } else if (Is_Censored[i] == 1){
      N[i]-K[i] ~ neg_binomial(K[i],theta[i]/(1-theta[i])); //負の二項分布 (K = 1のとき幾何分布と一致)
    }
  }
}

幾何分布は負の二項分布の特殊なケースなのでモデル上は負の二項分布を使って表現しています (現時点ではStanには幾何分布用の関数はありません)。neg_binomialの第2引数にはθそのものではなくθのオッズを与える必要があります。また,左辺にはNではなくN-1を与えます (Stanで使われている負の二項分布の定義に合わせるため)。このモデルを走らせるためのRコードは以下の通り。

Rコード01: MCMCの実行
library(rstan)
N_ID <- 8
N <- c(10,20,50,100,10,20,50,100)
K <- rep(1,N_ID)
Is_Censored <- c(0,0,0,0,1,1,1,1)

#コンパイル
stanmodel <- stan_model(file="binomgeom.stan")

#データの準備
data <- list(N_ID = N_ID, N = N,
             K = K, Is_Censored = Is_Censored)

#サンプリング
fit <- sampling(stanmodel, data=data, seed=123,
                 chains=4, iter=250500, warmup=500, thin=1)

 各々の場合のθの事後分布とMAP推定値は以下の通り。

posterior_binomgeom.png
図1. 幾何分布モデルと二項分布モデルで推定したθの事後分布とそのMAP推定値。

二項分布モデルで推定したθ1を半透明の赤色,幾何分布モデルで推定したθ2を半透明の青色で描画していますが,ほとんど重なっていて紫色にしか見えません。MAPも,多少の誤差はありますがほぼ一致しています。どうやら,どちらのモデルを使っても推定結果は一致ししそうです。

 両者が一致することは,事後分布を解析的に求めれば分かります。まず,二項分布の確率質量関数は,Nを所与として

posterior_binomgeom.png

と表されます。K = 1のとき,

posterior_binomgeom.png

となります。データが1つしかないときのθ1の尤度関数は

posterior_binomgeom.png

と表せます。最初のNθ1とは無関係な定数なので消去しています。最後に,θ1について無情報事前分布を仮定して尤度関数を正規化することで,次の事後分布が得られます。

posterior_binomgeom.png

 次に幾何分布ですが,確率質量関数は

posterior_binomgeom.png

と表されます。二項分布モデルの場合と同じように,これをθ2の関数 (尤度関数) とみなして正規化すると,

posterior_binomgeom.png

という事後分布が得られます。先ほどの二項分布モデルの事後分布と一致していることが確認できると思います。実はこの事後分布はどちらも (α, β) = (2, N) のベータ分布と一致します。二項分布と幾何分布の確率質量関数の違いは組み合わせ (C) を掛けるか否かの違いのみで,この部分は確率パラメータに依存しないため,事後分布を計算する過程で消去されるという訳です。ちなみに,データが複数ある場合やKが2以上の場合 (つまり幾何分布ではなく負の二項分布を使う場合) にも事後分布は同じパラメータを持つベータ分布 (α = 1 + 成功回数,β = 1 + 失敗回数) となり一致します (計算は省略)。

 という訳で,二項分布モデルを使っても幾何分布モデルを使っても確率パラメータの推定結果は一致することが分かりました。とはいえ,データの持つ情報 (データの収集方法の違い) を正確に表現にするためにもモデル上は両者を区別したほうが良さそうです。

posted by mutopsy at 15:22 | 統計

2019年01月04日

StanやJAGSの計算結果からパラメータのMAP推定値とHDIをまとめて得るためのRの自作関数

 こんにちは。mutopsyです。この記事では,StanやJAGSのfitオブジェクトを渡すとパラメータのMAP推定値とHDI (最高密度区間) をまとめて返してくれる自作関数を紹介します。

 StanやJAGSでパラメータの推定結果の要約を出力すると,通常は事後分布のEAP (i.e., 平均) や分位点の値が出力されますが,事後分布のMAP推定値やHDIが知りたい場合もありますよね。そんなわけで,そういう関数を作ってみました (もう既にあるかもしれないけど知ったこっちゃない)。自作関数といっても,HDIの計算にはHDIntervalというパッケージを利用しますが。関数は以下の通り。

Rコード01: MCMCサンプルからMAP推定値とHDIを計算する関数
library(HDInterval)
maphdi <- function(fit,credMass=0.99,include=TRUE,pars=NA){
  if(credMass<0|credMass>1) print("警告:credMassの値が正しくありません。")
  Class <- class(fit)[1]
  if(Class=="stanfit"){
    e <- as.data.frame(rstan::extract(fit))
  }else if(Class=="rjags"){
    e <- as.data.frame(fit$BUGSoutput$sims.list)
  }else if(Class=="data.frame"){
    e <- fit
  }else{
    print("エラー:未対応なオブジェクトです。")
    stop()
  }

  if(is.na(pars[1])) pars <- colnames(e)
  if (include==TRUE){
    e <- e[,pars]
  }else{
    e <- e[,!colnames(e) %in% pars]
  }
  if(length(pars)==1){
    MAP <- density(e)$x[which.max(density(e)$y)]
    HDI95 <- hdi(e)
    HDIAny <- hdi(e,credMass = credMass)
    Result <- data.frame(MAP = MAP, Lower95 = HDI95[1], Upper95 = HDI95[2],
                         LowerAny = HDIAny[1],UpperAny = HDIAny[2])
    rownames(Result) <- pars
  }else{
    MAP <- apply(e,2,function(z) density(z)$x[which.max(density(z)$y)])
    HDI95 <- apply(e,2,HDInterval::hdi)
    HDIAny <- apply(e,2,function(x) HDInterval::hdi(x,credMass = credMass))
    Result <- data.frame(MAP = MAP, Lower95 = HDI95[1,], Upper95 = HDI95[2,],
                         LowerAny = HDIAny[1,],UpperAny = HDIAny[2,])
  }

  colnames(Result)[4:5] <- c(
    paste("Lower",round(credMass*100,2),sep=""),
    paste("Upper",round(credMass*100,2),sep="")
  )
  return(Result)
  }

 試しにRからStanのfitオブジェクトを渡してみます。データとモデルは何でもよいので適当なものを渡します。

Rコード02: 関数を使ってみる
>fit.eg <- sampling(stanmodel, data=data, seed=123,
                 chains=4, iter=10500, warmup=500, thin=1)
> maphdi(fit.eg)
              MAP    Lower95    Upper95    Lower99    Upper99
beta.1  1.8502036  1.3540767  2.3442492  1.1987086  2.5187496
beta.2  0.6478270  0.5188389  0.7805164  0.4767612  0.8295815
beta.3  0.7103451  0.5996639  0.8201469  0.5627496  0.8524817
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.8891769 -0.2392025
sigma   0.3135227  0.2814639  0.3548575  0.2728315  0.3686322
lp__   97.7211912 93.6796922 99.1662181 91.7998798 99.3115815

 こんな風に,Stanのfitオブジェクトをmaphdi()に引数として渡すだけで各パラメータのMAP推定値と95% HDI・99% HDIを計算して出力することができます。Stan2.18, rstan2.18.2, Rtools3.5, R3.5.2, Windows10で動作確認済み。JAGSにも対応しています (JAGS4.3.0, R2jags0.5-7で確認済み)。

 また,fitオブジェクトの代わりに,MCMCサンプルが格納されたデータフレームを渡してもMAP推定値とHDIを返してくれます。

Rコード03: MCMCサンプルを直接渡す
> e <- as.data.frame(rstan::extract(fit.eg))
> head(e)
    beta.1    beta.2    beta.3     beta.4     sigma     lp__
1 2.194799 0.5785987 0.6255656 -0.4086454 0.3094092 97.46736
2 1.726557 0.6609587 0.7787883 -0.7231942 0.3058624 97.31513
3 2.053161 0.6021407 0.6534262 -0.3992853 0.3425830 96.40984
4 1.941408 0.6503287 0.6399921 -0.4178822 0.2976153 97.93040
5 2.518750 0.4448252 0.7676119 -0.7428853 0.3209268 91.90027
6 1.758346 0.6890341 0.6522970 -0.4046868 0.2985565 97.78386
> maphdi(e)
              MAP    Lower95    Upper95    Lower99    Upper99
beta.1  1.8502036  1.3540767  2.3442492  1.1987086  2.5187496
beta.2  0.6478270  0.5188389  0.7805164  0.4767612  0.8295815
beta.3  0.7103451  0.5996639  0.8201469  0.5627496  0.8524817
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.8891769 -0.2392025
sigma   0.3135227  0.2814639  0.3548575  0.2728315  0.3686322
lp__   97.7211912 93.6796922 99.1662181 91.7998798 99.3115815

 引数のcredMassに0〜1の範囲の値を与える事で,99% HDIの代わりに任意の確度のHDIを出力することもできます。95% HDIは必ず出力されるようになっています。

Rコード04: 任意のHDIを出力
> maphdi(e,credMass = 0.50)
              MAP    Lower95    Upper95    Lower50    Upper50
beta.1  1.8502036  1.3540767  2.3442492  1.7011924  2.0385575
beta.2  0.6478270  0.5188389  0.7805164  0.6060299  0.6964437
beta.3  0.7103451  0.5996639  0.8201469  0.6729364  0.7490473
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.6429855 -0.4713128
sigma   0.3135227  0.2814639  0.3548575  0.3030061  0.3280369
lp__   97.7211912 93.6796922 99.1662181 96.7857008 98.5912001

 引数のparsにパラメータの名前を与えれば,特定のパラメータの結果のみを出力することもできます。また,include = FALSE とすると,指定したパラメータ以外の結果を出力することもできます。

Rコード05: 一部のパラメータのみ計算

> maphdi(fit.eg,pars=c("beta.1","beta.2"))
            MAP   Lower95   Upper95   Lower99   Upper99
beta.1 1.850204 1.3540767 2.3442492 1.1987086 2.5187496
beta.2 0.647827 0.5188389 0.7805164 0.4767612 0.8295815
> maphdi(fit.eg,pars=c("beta.3"))
             MAP   Lower95   Upper95   Lower99   Upper99
beta.3 0.7103451 0.5996639 0.8201469 0.5627496 0.8524817
> maphdi(fit.eg,include=FALSE,pars=c("beta.1","beta.2"))
              MAP    Lower95    Upper95    Lower99    Upper99
beta.3  0.7103451  0.5996639  0.8201469  0.5627496  0.8524817
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.8891769 -0.2392025
sigma   0.3135227  0.2814639  0.3548575  0.2728315  0.3686322
lp__   97.7211912 93.6796922 99.1662181 91.7998798 99.3115815

 以上です。Enjoy!

posted by mutopsy at 19:10 | Rに関するTips

2019年01月01日

既知の確率質量関数からMAP推定値と最高密度区間 (HDI) を計算するためのRの自作関数

 あけましておめでとうございます。mutopsyです。年末にノロウイルスらしきものに感染してしまいましたがぶじ元気に2019年を迎えることができました。今年も何卒よろしくお願いいたします。新年早々ではありますが,必要に駆られてRの自作関数を作ったので共有します。

 既に分かっている確率質量関数を使ってMAP (最大事後確率) 推定値とHDI (highest density interval; 最高密度区間) を求める関数です。実現値として離散値kを返す既知の確率質量関数 (e.g., 二項分布) で事後分布が表現できる場合に,そのような事後分布のMAPとHDIを計算するのに使えます。ちなみに連続値を返す確率密度関数の場合でも,離散化すればこの関数を使って近似的にMAP推定値とHDIを求めることができるはずです。

 2019年1月2日追記:MAPとなる値やHDIの範囲が2つ以上ある場合 (例えば2峰分布など) にも対応できるように関数を修正しました。結果の出力の仕方も微修正。

Rコード01: 既知の確率質量関数からMAP推定値とHDIを計算する関数
hdimass <- function(prob.k,start.k=0,credMass=0.95){
  #警告メッセージ
  if(sum(prob.k)>1) print("警告:確率の和が1を越えています。")
  if(credMass<0|credMass>1) print("警告:credMassの値が正しくありません。")
  
  #HDIの計算の下準備
  credMass1m <- 1 - credMass
  accum<-1:length(prob.k)
  for(k in 1:length(prob.k)){
    if(k ==1){
      accum[k] <- sort(prob.k)[k]
    }else{
      accum[k] <- accum[k-1] + sort(prob.k)[k]
    }
  }
  Interval <- grep(TRUE,prob.k>max(sort(prob.k)[accum<=(min(abs(accum-credMass1m))+credMass1m)]))
  
  #HDIが複数ある場合やHDIが存在しない場合 (全範囲がHDI) を考慮
  if (length(Interval) <= 1){
    Check <- 0
    Interval <- 1:length(prob.k)
  } else{
    Check <- Interval[2:length(Interval)]-Interval[1:(length(Interval)-1)]-1
  }
  if(sum(Check)==0){
    Lower <- min(Interval) + (start.k-1)
    Upper <- max(Interval) + (start.k-1)
  } else{
    Lower <- c(min(Interval),Interval[grep(TRUE,Check!=0)+1]) + (start.k-1)
    Upper <- c(Interval[grep(TRUE,Check!=0)],max(Interval)) + (start.k-1)
  }
  
  #実際にその範囲の値となる確率
  ActualProb <- sum(prob.k[Interval])
  
  #MAPの計算
  MAP <- grep(TRUE,prob.k == max(prob.k)) + (start.k-1)
  
  #結果をリストにまとめて返す
  Result <- list(MAP=MAP,Lower=Lower,Upper=Upper,ActualProb=ActualProb)
  return(Result)
}

 第1引数 (prob.k) としてそれぞれのkが実現する確率のベクトルを,第2引数 (start.k) としてkの最小値 (デフォルトはstart.k=0) を与えることで,MAPと95% HDIの上限・下限を格納したリストが返り値として得られます。また,実際にHDIの範囲内の値をとる確率も返してくれます (離散分布なので95% HDIであっても厳密には95%にはなりませんが,95%に最も近い値になるはずです)。第1引数 (prob.k) には現実的な長さのベクトルを与えれば十分です。例えば N = 20,θ = 0.5 の二項分布であれば,長さが20もあれば総和がほぼ1になるので計算上は問題ありません (以下の例では一貫して長さ500のベクトルを与えていますがもっと少なくても結果は変わりません)。ちなみに第3引数であるcreadMassに(0,1)の範囲の値を与えれば (例えばcredMass = 0.99),異なる確度のHDIを出力することもできます (99% HDIなど)。

 計算の中身を簡単に説明します。まず最初にfor文で確率質量 (prob.k)を昇順にソートし,順番に累積確率を計算しています。次に,累積確率が閾値 (95% HDIなら5%)に最も近い点を探し,各々のkがHDIに含まれるか否かを論理値ベクトルで表現します。この論理値ベクトルを使って,(credMass)% HDIに含まれるkの値を抽出し,Intervalという内部の変数に格納しています。HDIが1つしかない場合は,Intervalの最小値と最大値ををHDIの下限・上限とみなしています。HDIが2つ以上ある場合は,Interval内で数値が飛び飛びになっている部分を見つけてその変化点をHDIの境界とみなしています。MAP推定値は最大事後確率を実現するkの値を求めることで計算しています。

 実際にこの関数を使って,さまざまなパラメータを持つ二項分布のMAPとHDIを計算してみました。

Rコード02: パラメータと引数を変えながら二項分布で試してみた結果
> #(1) 成功確率0.5,試行回数20の二項分布 (95% HDI)
> hdimass(dbinom(0:500,20,0.5),start.k = 0)
$MAP
[1] 10

$Lower
[1] 6

$Upper
[1] 14

$ActualProb
[1] 0.9586105

> #(2) 成功確率0.2,試行回数20の二項分布 (95% HDI)
> hdimass(dbinom(0:500,20,0.2),start.k = 0)
$MAP
[1] 4

$Lower
[1] 1

$Upper
[1] 7

$ActualProb
[1] 0.9563281

> #(3) 成功確率0.6,試行回数100の二項分布 (95% HDI)
> hdimass(dbinom(0:500,100,0.6),start.k = 0)
$MAP
[1] 60

$Lower
[1] 51

$Upper
[1] 69

$ActualProb
[1] 0.948118

> #(4) 成功確率0.6,試行回数100の二項分布 (99% HDI)
> hdimass(dbinom(0:500,100,0.6),start.k = 0,credMass = 0.99)
$MAP
[1] 60

$Lower
[1] 48

$Upper
[1] 72

$ActualProb
[1] 0.9896389

> #(5) 成功確率0.5,試行回数20,k>=7の切断二項分布 (95% HDI)
> hdimass(dbinom(7:500,20,0.5)/sum(dbinom(7:500,20,0.5)),start.k = 7,credMass = 0.95)
$MAP
[1] 10

$Lower
[1] 7

$Upper
[1] 13

$ActualProb
[1] 0.9388129

> #(6) MAPとHDIが2つずつある混合二項分布 (95% HDI)
> hdimass(dbinom(0:500,100,0.3) * 0.5 + dbinom(0:500,100,0.7) * 0.5)
$MAP
[1] 30 70

$Lower
[1] 22 61

$Upper
[1] 39 78

$ActualProb
[1] 0.9501802
#↑95% HDIは[22,39]または[61,78]の範囲

 いい感じに計算できていると思います。

 最後に,この関数を作った経緯を少しだけ説明します。僕が今行っている分析で離散パラメータの推定を行っているのですが,MCMCの結果からMAP推定値とHDIを十分高い精度で計算するためには発生させるMCMCサンプルの数 (JAGSのn.iterの値) をかなり増やさなければならないことが分かりました。これは現実的ではなかったので,数値積分を利用して事後分布を近似的に計算し,その事後分布からMAP推定値とHDIを計算することにしました。その副産物として今回の関数ができたという訳です。Rを使って事後分布を近似的に求める方法についてもいつか記事で解説するかもしれません。

 という訳で,需要があるかどうかも怪しい誰得関数の紹介から2019年が始まりました。Enjoy!

posted by mutopsy at 16:52 | Rに関するTips

2018年12月22日

銀のエンゼルが当たる確率のベイズ推定:金のエンゼルの何倍当たりやすいか?

 こんにちは。mutopsyです。この記事は,Stan Advent Calendar 2018の22日目のエントリー記事です。この記事では,森永製菓が発売しているおなじみのチョコレート菓子「チョコボール」で金のエンゼルが当たる確率と銀のエンゼルが当たる確率を,多項分布モデルを用いて同時推定し,銀のエンゼルが金のエンゼルの何倍出やすいかを明らかにすることを試みます。ちなみにこの記事は,以前書いた記事「金のエンゼルが当たる確率のベイズ推定」の続編という位置づけですが,前回の記事を読んでいなくても理解できるように書いたつもりです。前回と同様,チョコボール50周年を記念して行われた「金のエンゼル2倍キャンペーン」のデータを有効活用してベイズ推定を行います。

 ちなみに結論から言えば以下のような推定結果が得られました:

  • 金のエンゼルが当たる確率は0.02%─0.26% (MAPは0.09%;前回の結果のreproduce)
  • 銀のエンゼルが当たる確率は3.28%─5.51% (MAPは4.29%)
  • 珍しさという点で,金のエンゼル1枚は銀のエンゼル約11─159枚分の価値がある (MAPは33枚分)

【お詫び】当初の予定では,帰無仮説検定の代わりとしてのベイズファクターについて初学者向けの解説記事を書く予定でおりましたが,不測の事態が重なり準備が間に合わなかったため,急きょ内容を変更させて頂きました。楽しみにしていた方には申し訳ございません。いつかちゃんと書きます!

1. 目的と概要

 改めまして,こんにちは。チョコボールと早めのMCMCが効いたのか,先週から患っていた風邪が治りつつあるmutopsyです。前回の記事「金のエンゼルが当たる確率のベイズ推定」では,先人たちが収集したチョコボールデータを用いて二項分布モデルによるなんちゃってメタ分析を行い,金のエンゼルが当たる確率がMAPで約0.09% (99% HDI [0.02%, 0.26%]) であることを明らかにしました。続編となる今回の記事では,金のエンゼルが当たる確率と銀のエンゼルが当たる確率を同時に推定し,さらに銀のエンゼルが金のエンゼルの何倍出やすいのかを区間推定した結果を報告します

 1つの推定方法としては,前回と同じように二項分布モデルを用いて銀のエンゼルが出る確率 (金のエンゼルが出ず,はずれでもない確率) を推定するという手が考えられます。その推定結果を,前回得られた金のエンゼルが当たる確率の推定結果の点推定値 (0.09%) で除することで,銀のエンゼルが金のエンゼルの何倍出やすいかを計算することができます。ところが,この方法では「99%の確率で,銀のエンゼルは金のエンゼルよりも〇〜△倍出やすい」といった幅をもった推定 (区間推定) を行うことは困難です。そこで本記事では,金のエンゼルが当たる確率p_goldと銀のエンゼルが当たる確率p_silverを1つのモデルで同時に推定し,さらにこれらの比 ratio = p_gold / p_silver の生成量を計算します。このようなモデルを用いることで,金のエンゼルが当たる確率と銀のエンゼルが当たる確率の不確実性に関する情報 を保持したまま,その比の区間推定を行うことが可能です。このように,パラメータの事後分布が持つ情報を損失することなく差や比などの生成量の事後分布を簡単に推定できるのは,ベイズ推定の強みの1つです。

2. エンゼルが当たる確率の生成メカニズムとしての多項分布モデル

 では,金のエンゼルが当たる確率と銀のエンゼルが当たる確率を同時に推定するにはどのようなモデルを用いれば良いのでしょうか。前回の記事では,あるチョコボールの箱のクチバシを確認したときに「金のエンゼルが当たる事象」と「金のエンゼルが当たらない事象」という2つの事象が起こり得るというベルヌーイ分布を仮定し,この試行を独立に複数回繰り返したときに金のエンゼルが当たる回数が二項分布に従うことを利用しました。しかし,実際には「金のエンゼルが当たる事象」「銀のエンゼルが当たる事象」「いずれのエンゼルも当たらない事象」という3つの事象が存在します。このように起こり得る事象が3つ以上存在する場合は,ベルヌーイ分布ではなく次のようにカテゴリカル分布を仮定するのが自然です。

Y ~ Categorical(p)
p = {p_gold, p_silver, (1 - p_gold - p_silver)}

ここで,Yはいずれの事象が観測されたかを表す名義変数 (e.g., 0 = はずれ,1 = 銀のエンゼルが当たった,2 = 金のエンゼルが当たった),p はそれぞれの事象が生起する確率を成分とするベクトルを表します (e.g., p = {1%, 5%, 94%}; 和は必ず100%になる)。いずれのエンゼルも当たらない確率は,金のエンゼルが当たる確率p_goldと銀のエンゼルが当たる確率p_silverを用いて1 - p_gold - p_silverと表現できます。要するに,ベルヌーイ分布を3値以上に拡張した分布がカテゴリカル分布です。もちろん事象の数が4つ以上の場合もカテゴリカル分布が利用できます (e.g., ダイスを振っていずれの目が出るか)。カテゴリカル分布に従う試行を独立にN回繰り返したときにそれぞれの事象が生起する回数は多項分布に従います。

Count ~ Multinomial(N, p)
Count = {Count_Gold, Count_Silver, Count_Miss}
p = {p_gold, p_silver, (1 - p_gold - p_silver)}

多項分布は二項分布を3値以上に拡張したものに相当します。なお,Stanでモデルを記述するときにはNは省略します (Countの成分の和から自動的に計算されるため)。

 続いて,前回の記事と同様,「金のエンゼル2倍キャンペーン」について考えます。チョコボール誕生50周年を記念して,2017年11月〜2018年2月の期間,金のエンゼルが当たる確率が通常の2倍になるというキャンペーンが行われました。また,この期間は銀のエンゼルは出現しないことが明言されていました。つまり,金のエンゼル2倍キャンペーン中は,金のエンゼルが当たる確率が2倍になり,銀のエンゼルが当たる確率が0倍になる,と表現することが可能です。そこで先ほどのモデルを少し修正し,

Count ~ Multinomial(N, p)
Count = {Count_Gold, Count_Silver, Count_Miss}
p = {p_gold × W_Gold, p_silver × W_Silver, (残り)}

と表現します。通常時には (W_Gold, W_Silver) = (1,1),金のエンゼル2倍キャンペーン時には (W_Gold, W_Silver) = (2,0) という重みを付与することで,キャンペーン時の確率を表現することができます。さらに,銀のエンゼルが金のエンゼルの何倍当たりやすいかを表す生成量ratioは,

ratio = p_gold / p_silver

と表現することができます。

 最後に,金のエンゼルが当たる確率と銀のエンゼルが当たる確率の事前分布について考えましょう。ベルヌーイ分布や二項分布のパラメータpの事前分布としては前回の記事でも使用したベータ分布がよく用いられますが,カテゴリカル分布や多項分布のパラメータpの事前分布としてはディリクレ分布が用いられることが多いです。ディリクレ分布はベータ分布を多値に拡張したものと考えることができます。ただし,今回の記事ではディリクレ分布は使用せず,金のエンゼルが当たる確率と銀のエンゼルが当たる確率のそれぞれについて事前分布を考えることにします (これらの事前分布が決まれば,いずれのエンゼルも当たらない確率の事前分布も自動的に決まります)。

 まず,金のエンゼルが当たる確率の事前分布に関しては,前回の記事と同様,「今まで100箱以上のチョコボールを確認したが,金のエンゼルは1度も当たったことがない」という筆者の経験を反映させて,(α, β) = (1, 101) のベータ分布を考えようと思いますが,今回はそれに加えて,「金のエンゼルが当たる確率は50%未満である」という常識的な情報も付与した切断ベータ分布を金のエンゼルが当たる確率の事前分布として仮定します。

p_gold ~ Beta(1, 101)[0.0, 0.5]

この切断ベータ分布は,値が0.5以上となる確率が0になるようにベータ分布を修正したものです。もっとも,もともとBeta(1, 101) から0.5以上の値が得られる確率はもともとほぼ0%なので,Beta(1, 101)を事前分布として用いるのとほとんど変わりません。

 続いて銀のエンゼルが当たる確率の事前分布ですが,こちらに関してはこれといった事前知識がありません。銀のエンゼルは筆者自身何度か当てたことがあり,おもちゃのカンヅメと交換してもらったこともあるのですが,どのくらい当たりやすいのかについては実際のところよく分かりません。ただ,常識的に考えて50%よりも高いということは無さそうですので,

p_silver ~ Uniform(0.0, 0.5)

という一様分布を事前分布として仮定することにします。ところで,金のエンゼル・銀のエンゼルともに確率の上限が50%となるように事前分布を設定しましたが,これにより,金のエンゼルが当たる確率と銀のエンゼルが当たる確率の和が100%を超えてしまうのを未然に防ぐことができます。確率の和が100%を超えたり0%を下回ったりすることを許容するモデルでは推定がうまくいかないことがあるので,このように何らかの工夫を施すことが必要になる場合があります。他の方法としては,各事象の生起確率の事前分布をディリクレ分布などを用いて同時に設定したり,softmax関数やmin関数・max関数などを使って確率の和が1になるように補正したりする方法などがあります。Stanであれば,simplex型で変数を宣言することで,成分の和が1になるベクトルとしてパラメータを宣言することもできます (今回のモデルではsimplex型ではうまくいかないのでvecter型を指定します)。

3. 使用するデータおよびStanとRのコード

 前回の記事と同様,分析で使用するデータは以下の出典から拝借いたしました。

 これらの結果を次のようなデータとしてまとめました。

表1. 今回の分析対象となるチョコボールデータ
DataID Count_Gold Count_Silver Count_Miss W_Gold W_Silver
1 0 13 276 1 1
2 1 0 113 2 0
3 1 42 872 1 1
4 2 40 958 1 1
5 0 0 1000 2 0

2018/12/23追記:「月曜から夜ふかし」のデータは,あらかじめ開封する箱の数を決めていたのではなく,金のエンゼルが当たるまで箱を開封し続けた結果得られたデータです。このような場合には二項分布や多項分布を用いるのは適切ではなく,本当は幾何分布を用いるべきなのですが,今回は目をつむることにします。この辺りの事情を考慮したモデルはまた別の機会にお伝えします。)

 Stanの仕様に合わせて,今回は試行の総数は使用せず,「金のエンゼルが当たった回数 (Count_Gold)」,「銀のエンゼルが当たった回数 (Count_Silver)」,「いずれのエンゼルも当たらなかった回数 (Count_Miss)」を基本データとして渡します。また,金のエンゼル2倍キャンペーン期間の確率の重みをW_GoldとW_Slverの2つの変数を使って与えておきます (キャンペーン期間中は金のエンゼルが当たる確率が2倍,銀のエンゼルが当たる確率が0倍になる)。

 前節で説明したモデルを実装するためのStanコードは以下の通りです。

Stanコード01: 金のエンゼルが当たる確率と銀のエンゼルが当たる確率の同時推定
data {
  int N;  //データセットの数
  int Count_Gold[N];  //金のエンゼルが当たった回数
  int Count_Silver[N];  //銀のエンゼルが当たった回数
  int Count_Miss[N]; //いずれのエンゼルも当たらなかった回数
  int W_Gold[N]; //金のエンゼルが当たる確率の重み (1か2)
  int W_Silver[N]; //銀のエンゼルが当たる確率の重み (1か0)
}

transformed data{
  int Count[N,3];
  Count[,1] = Count_Gold;
  Count[,2] = Count_Silver;
  Count[,3] = Count_Miss;
}

parameters{
  real<lower=0,upper=0.5> p_gold; //金のエンゼルが当たる確率
  real<lower=0,upper=0.5> p_silver; //銀のエンゼルが当たる確率
}

transformed parameters{
  vector[3] p[N];
  for (n in 1:N){
    p[n,1] = p_gold * W_Gold[n];
    p[n,2] = p_silver * W_Silver[n];
    p[n,3] = 1 - (p[n,1] + p[n,2]);
  }
}

model{
  //モデル
  for (n in 1:N){
    Count[n] ~ multinomial(p[n]);
  }
  
  //事前分布
  p_gold ~ beta(1,101);
}

generated quantities{
  real ratio;
  ratio = p_silver / p_gold; //銀のエンゼルは金のエンゼルの何倍当たりやすいか
}

金のエンゼルが当たる確率と銀のエンゼルが当たる確率はどちらも50%以下である,という事前分布の情報は,パラメータを宣言するときに「upper = 0.5」と上限を指定することで表現できます。

 続いて,データを上記Stanコードに渡すためのRコードを示します。

Rコード01: チョコボールデータをStanに渡す
library(rstan)
d <- read.csv("チョコボールデータ.csv")
stanmodel <- stan_model(file="Stancode/chocoball_multinomial.stan")
data <- list(N = nrow(d), Count_Gold = d$Count_Gold, Count_Silver = d$Count_Silver,
             Count_Miss = d$Count_Miss, W_Gold = d$W_Gold, W_Silver = d$W_Silver)

fit2 <- sampling(stanmodel, data=data, seed=123,
                 chains=4, iter=50500, warmup=500, thin=1)

 分析にはrstan 2.16.2およびR 3.4.2を使用しました。長さ50,500のチェインを4本発生させ,ウォームアップ期間を500とし,得られた200,000個のMCMCサンプルを用いて事後分布を近似しました。Rhatは1.01未満で,目視からも十分に収束していることが確認できました。

4. 金と銀のエンゼルが当たる確率の同時推定の結果

 それでは,金のエンゼルが当たる確率と銀のエンゼルが当たる確率の推定結果をお示しします。ついでに,いずれのエンゼルも当たらない確率の事後分布も示します。次の図は,各々の確率の事後分布です。MAPと99%最高密度区間 (HDI) も示しています。

multinomial_result.png
図1. 金のエンゼルが当たる確率p_gold,銀のエンゼルが当たる確率p_silver,およびいずれのエンゼルも当たらない確率の事後分布。MAPおよび99%最高密度区間も記してある。

 まず,金のエンゼルが当たる確率のMAPは0.086% (99% HDI[0.017%, 0.259%],95% HDI [0.027%, 0.208%]) でした。前回の記事の推定結果と比べると僅かに値が小さいですが,これは推定上の誤差であると思われます (もしかしたら切断ベータ分布を事前分布に用いた影響もあるかもしれません)。とはいえ小数第2位までの値は一致しているので,前回の結果が再現できたと言っても良いでしょう。続いて,今回新たに得られた結果として,銀のエンゼルが当たる確率のMAPは4.290% (99% HDI[3.281%, 5.511%], 95% HDI[3.506%, 5.201%]),いずれのエンゼルも当たらない確率のMAPは95.60% (99% HDI[94.38%, 96.62%], 95% HDI[94.66%, 96.37%]) であると推定されました。つまり結論としては,99%の確率で,金のエンゼルが当たる確率は0.02%─0.26%,銀のエンゼルが当たる確率は3.28%─5.51%,いずれのエンゼルも当たらない確率は94.38%─96.62%の範囲にあると主張することができそうです。

5. 銀のエンゼルは金のエンゼルの何倍出やすいか

 最後に,「銀のエンゼルは金のエンゼルの何倍当たりやすいか」という疑問に答えることを試みます。この値は生成量のratioの事後分布を見れば分かります。

goldangel_posterior.png
図2. 銀のエンゼルが金のエンゼルの何倍当たりやすいかを表す倍率ratioの事後分布。MAPと99%最高密度区間も示してある。

この結果から,MAPで見ると,銀のエンゼルは金のエンゼルの32.6倍当たりやすいということが分かります。ただし,99%最高密度区間は11.0─158.6倍 (95% HDIは13.6─102.4倍) とかなり広い点には注意が必要です。ところで,おもちゃのカンヅメは金のエンゼル1枚もしくは銀のエンゼル5枚と交換できますが,金のエンゼル1枚の価値は本当に銀のエンゼル5枚分と同等なのでしょうか?これを検証するために,MCMCサンプルを使って「銀のエンゼルが当たる確率は金のエンゼルが当たる確率の5倍よりも大きい」という事象の生起確率を計算したところ100%となりました (2万個のMCMCサンプルは全て5よりも大きな値でした)。したがって,少なくとも珍しさという点において,金のエンゼルは銀のエンゼル5枚分以上の価値があると考えて間違いはなさそうです。より具体的に言えば,金のエンゼル1枚は銀のエンゼル約11─159枚分の価値があるということになります。金のエンゼルが当たっても,すぐにおもちゃのカンヅメと交換してしまうのはもったいないのかもしれません。

 なお,前回の記事でも同じことが言えますが,確信区間の広さからも分かるように,エンゼルが当たる確率についてより正確な推定を行うにはまだまだデータが足りません。さらに多くのチョコボールデータを収集することができれば,同じモデルを使ってより精度の高い推定を行うことが可能となります。観測とデータ収集は科学の基本。データをどんどん蓄積していけば,エンゼルたちの微笑みの秘密にもっともっと近づけるはずです。筆者も貢献できるように頑張ります。ちなみにチョコボールデータを用いた分析のアイディアはもう少し残っているので (離散パラメータの推定に関するネタと,もう少し実用的なネタ),また機会があれば報告します。

posted by mutopsy at 00:00 | 統計

2018年12月18日

金のエンゼルが当たる確率のベイズ推定

 こんにちは。mutopsyです。この記事は,Stan Advent Calendar 2018の18日目のエントリー記事です。この記事では,森永製菓が発売しているおなじみのチョコレート菓子「チョコボール」で金のエンゼルが当たる確率の推定方法を例に,Stanによるベイズ推定のいろはをお伝えすることを目指します。チョコボール50周年を記念して行われた「金のエンゼル2倍キャンペーン」のデータを有効活用しつつ,事前分布としてベータ分布を用いた二項分布モデルで金のエンゼルが当たる確率を推定します。以前,ポン・デ・リングのネタ記事を書いたことがありますが,今度はチョコボールのガチ記事です。

2018/12/22追記続編にあたる記事を書きました。→「銀のエンゼルが当たる確率のベイズ推定:金のエンゼルの何倍当たりやすいか?

1. 目的と概要

 改めまして,こんにちは。先日ぶじに博士論文を提出することができたmutopsyです。面白いことに(?),博士論文を提出したその日から風邪を引いてしまい,まだ熱が引いておりませんが,昔から「風邪にはMCMC」「Stanは百薬の長」とも言いますので,パラメータを収束させつつ体の熱を発散させようという意気込みでこの記事を執筆することにいたします。

 タイトルの通り,この記事では「金のエンゼルが当たる確率」を主題として取り上げます。ご存知ない方はおられないかと思いますが,金のエンゼルとは,森永製菓から発売されているチョコレート菓子「チョコボール」の箱の上部にある「クチバシ」にごくまれに印刷されているもので,いわゆる「当たり」です。これ1枚と,夢がつまった「おもちゃのカンヅメ」を交換することができます。ちなみに銀のエンゼルも存在し,こちらは5枚でおもちゃのカンヅメと交換できますが,今回は銀のエンゼルについては触れないことにします(2018/12/22追記:続編で銀のエンゼルも扱いました。→「銀のエンゼルが当たる確率のベイズ推定:金のエンゼルの何倍当たりやすいか?」)。

 筆者はチョコボールの大ファンで,子どものころから今に至るまで何度となく食べ続けてきました。その甲斐あって,銀のエンゼルには何度かお目にかかり,おもちゃのカンヅメを入手することもできましたが,未だかつて金のエンゼルをこの目で見たことはありません。では一体,金のエンゼルが当たる確率は一体どのくらいなのでしょうか。ベイズ推定という武器を手に入れたいま,この謎に取り掛かろうと思います。

 とはいえもちろん,先行研究は多数存在します。例えばチョコボール統計さんは,日々データを蓄積し,さまざまな分析を通じてチョコボールの神秘を解き明かそうと試みておられます。また,後述するように,テレビやYoutubeの企画などでチョコボールを大量買いした結果を報告している先行研究も存在します。これらの先人の知恵をお借りし,巨人の肩に乗り,金のエンゼルが当たる確率をなるべく正確に推定するということが本記事の主たる目的です。

 基本的なモデルは非常にシンプルです。まず,あるチョコボールの箱のクチバシに金のエンゼルが存在する確率をpと置きます。このpが今回推定したいパラメータです。ここでは,確率pはチョコボールの味や開封者,購入時期の違いなどに依らず一定であると仮定します。さて,いま,目の前にチョコボールの箱があるとします。この箱のクチバシを確認して,金のエンゼルがあれば1,なければ0と記録することにします。なお,ここでは銀のエンゼルがあった場合も0として記録することにします。このとき,金のエンゼルがあったか否かを表す2値変数Yは確率pのベルヌーイ分布に従います。

Y ~ Bernoulli(p)

ベルヌーイ分布は非常に単純な分布で,確率pでY = 1,確率(1-p)でY = 0となることを表現しているだけに過ぎません。このような試行をCountTrial回繰り返したとき,金のエンゼルが当たる回数CountGoldは次の二項分布に従います。

CountGold ~ Binomial(CountTrial, p)

これは,表が出る確率がpであるコインをCountTrial回投げたときに表が出る回数の分布と同じです。以上より,「チョコボールの箱を開けた回数」と「その中で金のエンゼルが見つかった回数」のデータがあれば,金のエンゼルが当たる確率pを推定することができるということです。ところが,次節で紹介するような事情を考慮すると,もう少しモデルを柔軟にしたほうが良さそうであることが分かります。

2. 金のエンゼル2倍キャンペーン

 チョコボール誕生50周年を記念して,2017年11月〜2018年2月の期間,金のエンゼルが当たる確率が通常の2倍になるというキャンペーンが行われました。そこで,このキャンペーン期間中のデータを有効活用できるように,先ほどの二項分布モデルを次のように拡張します。

CountGold ~ Binomial(CountTrial, p × W)

変更箇所はWという重みを表す観測変数です。通常の商品であればW = 1,金のエンゼル2倍キャンペーン対象商品であればW = 2 として確率に重みを付けることにより,このキャンペーンの内容をモデルに反映させることができます。それでは,このモデルを使って金のエンゼルが当たる確率を推定してみましょう。

3. 使用させて頂くデータ

 分析で使用するデータは以下の出典から拝借いたしました。今回は使用しませんが,銀のエンゼルの集計結果も一緒に掲載しています。なお,金のエンゼル2倍キャンペーン期間中は銀のエンゼルが出ないようになっていたようです。

 これらの結果を次のようなデータとしてまとめました。

表1. 今回の分析対象となるチョコボールデータ
DataID CountGold CountSilver CountTrial W
1 0 13 289 1
2 1 0 114 2
3 1 42 915 1
4 2 40 1000 1
5 0 0 1000 2

2018/12/23追記:「月曜から夜ふかし」のデータは,あらかじめ開封する箱の数を決めていたのではなく,金のエンゼルが当たるまで箱を開封し続けた結果得られたデータです。このような場合には二項分布を用いるのは適切ではなく,本当は幾何分布を用いるべきなのですが,今回は目をつむることにします。この辺りの事情を考慮したモデルはまた別の機会にお伝えします。)

 今回の分析ではこれらのデータを使用します。複数の研究で得られたデータを統合して分析するという意味ではいわゆるメタ分析に近いかもしれませんが,今回はデータ収集者などに由来するバイアスは存在しないと仮定しているので,モデル自体はとてもシンプルです。もしかしたら出版バイアスぐらいは考慮したほうが良いのかもしれません(e.g., 金のエンゼルが1つも出なかった場合,どこにも報告されずにその結果がお蔵入りになる,など)。

4. 事前分布

 ベイズ推定を行う際は,あらかじめパラメータの事前分布を設定しておく必要があります。今回の場合であれば,金のエンゼルが当たる確率pの事前分布を決めなければなりません。一番単純な決め方は,次のように一様分布を仮定することでしょうか。

p ~ Uniform(0, 1)
p ~ Beta(1, 1)

つまり,pは確率なので0─1 (i.e., 0%─100%) の範囲に収まるが,実際に確率がどの程度なのかは見当も付かない,という信念を反映させた事前分布となります。もちろんこれでも良いのですが,筆者としてはあまり腑に落ちません。上述した通り,筆者はこれまで幾度となくチョコボールを食べてきましたが,一度も金のエンゼルを見たことがありません。きちんと計測した訳ではありませんが,少なく見積もっても100箱以上のチョコボールは確認したはずです。したがって,金のエンゼルが当たる確率は相当低いに違いありません。少なくとも50%─100%という範囲に収まるという事はあり得なさそうですが,一様分布を仮定した場合,pが0%─50%の範囲に存在するのと同じ程度の確率で50%─100%の範囲に存在するということになってしまいます。筆者のこの信念を事前分布として反映させるために,

p ~ Beta(1, 101)

というベータ分布を事前分布として用いることにします。ベータ分布は2つのパラメータαとβによって形状が決まり,実現値は0─1の範囲をとります (0と1は含まない)。(α, β) = (1, 1) のとき,一様分布と一致します。2つのベータ分布を図示すると以下のようになります。

beta_distribution.png
図1. ベータ分布の説明。

ちなみに今回使用するBeta(1, 101)は,pの事前分布を一様分布とした場合に,100箱のチョコボールを確認したがどのクチバシにも金のエンゼルが印刷されていなかったときのpの事後分布と一致します。これは筆者の実際の経験をうまく反映した分布であるといえます (実際のβはもう少し大きいかもしれない)。この情報は正確にカウントした結果ではなく,あくまでも筆者の信念を反映したものであるという意味で,データではなく事前分布として扱います。

 なお,データを見てから「確率は低そうだ」と判断して事前分布を設定し,その上で同じデータを使って分析してはいけません(情報が重複してしまうため)。もし「100箱開封したが金のエンゼルが1つも当たらなかった」という経験をデータとして扱うのであれば,事前分布には一様分布を用いるのが妥当でしょう(この場合,経験を事前分布として扱った場合と同じ結果が得られます)。もしくは,「当たり」なのだから確率は十分低いはず,という常識に即した信念を,緩い事前分布の形で表現するのも良いでしょう。どの事前分布を選ぶかは分析の目的と分析者の信念次第ですね (例えば,科学的な研究であればより客観的かつ恣意的でない事前分布を採用するべきでしょう)。

5. StanとRのコード

 これまでに説明したモデルを実装するためのStanコードは以下の通りです。

Stanコード01: 金のエンゼルが当たる確率の推定
data {
  int N;  //データセットの数
  int CountGold[N];  //金のエンゼルが当たった数
  int CountTrial[N]; //総数
  int W[N]; //重み (1か2)
}

parameters{
  real<lower=0,upper=1> p; //金のエンゼルが当たる確率
}

model {
  //モデル
  for (n in 1:N){
    CountGold[n] ~ binomial(CountTrial[n],p * W[n]);
  }
  
  //事前分布
  p ~ beta(1,101);
}

generated quantities{
  real p_zero[1000];
  for (i in 1:1000){
    p_zero[i] = (1-p)^i; //i箱見て1箱も金のエンゼルが入っていない確率
  }
}

 最後のgenerated quantitiesブロックでは,i箱のチョコボールを確認して1つも金のエンゼルが入っていない確率p_zeroを計算しています。今回の場合はStan内で生成量を計算するよりも,要約統計量を用いて事後的に生成量を計算したほうが効率的ですが,念のため書き方を示しておきました。

 続いて,データを上記Stanコードに渡すためのRコードを示します。

Rコード01: チョコボールデータをStanに渡す
library(rstan)
d <- read.csv("金のエンゼルデータ.csv")
stanmodel01 <- stan_model(file="Stancode/chocoball.stan")
data <- list(N = nrow(d), CountGold = d$CountGold,
             CountTrial = d$CountTrial, W = d$W)

fit1 <- sampling(stanmodel01,data=data,seed=123,
                 chains=4,iter=10500,warmup=500,thin=1)

 分析にはrstan 2.16.2およびR 3.4.2を使用しました。長さ10,500のチェインを4本発生させ,ウォームアップ期間を500とし,得られた40,000個のMCMCサンプルを用いて事後分布を近似しました。Rhatは1.01未満で,目視からも十分に収束していることが確認できました。

6. 分析結果

 それでは分析結果をお示しします。次の図は,金のエンゼルが当たる確率pの事後分布です (破線で事前分布も示しています)。

goldangel_posterior.png
図2. 金のエンゼルが当たる確率pの事後分布 (破線は事前分布)。MAPおよび99%最高密度区間も記してある。

 この結果から,金のエンゼルが当たる確率の点推定値はおよそ0.091%であると結論づけることができそうです。10,000箱に9箱の割合です。なかなか厳しいですね。また,99%最高密度区間 (HDI) の下限は0.018%,上限は0.262%という結果でした (95%HDIは[0.027%, 0.209%])。つまり,どんなに当たる確率が低かったとしても0.018%未満であるということはなさそうですし,どんなに確率を高く見積もったとしても0.262%を上回ることはなさそうです。

 最後に,チョコボールをn箱空けた時に一度も金のエンゼルが出ない確率の事後予測分布を示します。これは先ほどのStanコードのgenerated quantitiesブロックで計算した生成量に基づく分析です。開封する箱の数を横軸に取ると,グラフは次のようになります (ちなみに,縦軸の値に金のエンゼルが当たる確率pを乗算すると,幾何分布と呼ばれる確率分布になります)。

goldangel_posterior.png
図3. 金のエンゼルが一度も当たらない確率のMAPと95%最高密度区間を,開封する箱の数ごとに示したグラフ。

例えば,100回チョコボールを開封して金のエンゼルが1つも出ない確率はMAPで91.3% (95%HDI[81.1%, 97.3%]) となります。200回でも83.1% (95%HDI[65.4%, 94.3%]) です。たまにチョコボールを食べている人がたまたま金のエンゼルを当てられる確率はそこまで高くなさそうですね。筆者がこれまで一度も金のエンゼルを当てたことがないのも不思議ではなさそうです(その情報を使って分析しているので若干トートロジックではありますが)。一方,今回の事後予測分布から,762回チョコボールを開封すれば,金のエンゼルが1つも出ない確率のMAPが50%を下回るようです。つまり,762箱のチョコボールを確認すれば,少なくとも1箱から金のエンゼルが見つかる確率が50%を超えるであろうことが見込まれます。また,3,291箱開封すれば,少なくとも1つは金のエンゼルが見つかる確率が95%を超え,5,059箱で99%を超えることも示されました。本気で金のエンゼルを狙うのであれば,3,000〜5,000箱程度のチョコボールを用意すればほぼ確実でしょう (ただし,確信区間の広さを考慮するとそれでも足りないかもしれません)。ちなみに,金のエンゼル2倍キャンペーン時に1,000箱のチョコボールを開封して金のエンゼルが1つも手に入らない確率はMAPで16.2%,95%最高密度区間はおよそ1.5%─58.6%の範囲となります。先行研究のYoutube動画ではキャンペーン期間中に1,000箱買っても1枚も金のエンゼルが見つからなかったことが報告されていますが,これは十分ありうることと言えるでしょう。

 今回の分析において留意すべき点はいくつかありますが,そもそも金のエンゼルが観測される頻度が極めて少ないことは気に留めるべきところかもしれません。実際,合計3,318箱ものチョコボールを開封しても金のエンゼルはたったの4枚しか出ていません。これだけ起こりにくい現象の生起確率を高い精度で推定するためにはもっともっと多くのデータが必要でしょう。銀のエンゼルが当たる確率についても同様の分析を行うことが可能ですので,結果を比較してみるのも面白いかもしれません (事前分布は変更する必要があります)。本当は「おまけ」として,多項分布を使って金のエンゼルが当たる確率と銀のエンゼルが当たる確率を同時に推定した分析結果も掲載する予定でいたのですが,残念ながら間に合いませんでした。申し訳ございません。いつかリベンジします(2018/12/22追記:リベンジしました→「銀のエンゼルが当たる確率のベイズ推定:金のエンゼルの何倍当たりやすいか?」)。おもちゃのカンヅメに対する人類の挑戦はこれからも続く……。今回は自分のデータではなく人様のデータを拝借させて頂きましたが,筆者自身もチョコボールデータの収集に貢献できればと考えています。けどその前にまず風邪を治します!おやすみなさい!

※2018/12/18 16:55追記:最後の分析の結果に誤りがあったのでグラフと数値を訂正し,少し加筆しました。

posted by mutopsy at 05:00 | 統計