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

　こんにちは。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枚分)

2019/12/15追記：この記事の内容をさらに発展させた内容が書籍になりました。→「たのしいベイズモデリング2：事例で拓く研究のフロンティア」の第1章

目次
1. 目的と概要
2. エンゼルが当たる確率の生成メカニズムとしての多項分布モデル
3. 使用するデータおよびStanとRのコード
4. 金と銀のエンゼルが当たる確率の同時推定の結果
5. 銀のエンゼルは金のエンゼルの何倍出やすいか

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のコード

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

• チョコボール統計
  ※「第107回 チョコボール計測」までのデータを使用しました（番外編のデータは未使用）。
  →【通常期間】金：0, 銀：13，総数：289
  →【2倍キャンペーン期間】金：1，銀：(0)，総数：114
• テレビ番組「月曜から夜ふかし」（日本テレビ，2015年9月15日放送）
  →【通常期間】金：1, 銀：42，総数：915
• Youtube動画「【検証】チョコボール1000個開封して金のエンゼル、銀のエンゼルの確率を調べてみた」
  →【通常期間】金：2, 銀：40，総数：1000
• Youtube動画「金のエンゼル確率2倍はウソだった？チョコボール1000個開封して確かめた結果…」
  →【2倍キャンペーン期間】金：0, 銀：(0)，総数：1000

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

表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) も示しています。

図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の事後分布を見れば分かります。

図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枚分の価値があるということになります。金のエンゼルが当たっても，すぐにおもちゃのカンヅメと交換してしまうのはもったいないのかもしれません。

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

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

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

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

2019/12/15追記：この記事の内容をさらに発展させた内容が書籍になりました。→「たのしいベイズモデリング2：事例で拓く研究のフロンティア」の第1章

目次
1. 目的と概要
2. 金のエンゼル2倍キャンペーン
3. 使用させて頂くデータ
4. 事前分布
5. StanとRのコード
6. 分析結果

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倍キャンペーン期間中は銀のエンゼルが出ないようになっていたようです。

• チョコボール統計
  ※「第107回 チョコボール計測」までのデータを使用しました（番外編のデータは未使用）。
  →【通常期間】金：0, 銀：13，総数：289
  →【2倍キャンペーン期間】金：1，銀：(0)，総数：114
• テレビ番組「月曜から夜ふかし」（日本テレビ，2015年9月15日放送）
  →【通常期間】金：1, 銀：42，総数：915
• Youtube動画「【検証】チョコボール1000個開封して金のエンゼル、銀のエンゼルの確率を調べてみた」
  →【通常期間】金：2, 銀：40，総数：1000
• Youtube動画「金のエンゼル確率2倍はウソだった？チョコボール1000個開封して確かめた結果…」
  →【2倍キャンペーン期間】金：0, 銀：(0)，総数：1000

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

表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つのベータ分布を図示すると以下のようになります。

図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の事後分布です (破線で事前分布も示しています)。

図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を乗算すると，幾何分布と呼ばれる確率分布になります)。

図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追記：最後の分析の結果に誤りがあったのでグラフと数値を訂正し，少し加筆しました。

認知心理学への実践：データ生成メカニズムのベイズモデリング[スライド紹介とWS感想]

　こんにちは。mutopsyです。昨日（2018年6月30日），専修大学にて，広島ベイズ塾第三回ワークショップ「心理学者のためのベイズ統計学：モデリングの実際と，モデル選択・評価」が開催されました。本ワークショップですが，告知後ほんの数日で200名以上の予約参加者が集まりまして，当日のYouTubeのライヴ配信もアクティブ人数が常に40〜60人程度となるほどの盛り上がりを見せました。会場にお越しの皆様，配信を視聴してくださった皆様，そして一緒にワークショップの運営をしてくださった皆様，本当にありがとうございました。ベイズ塾生として発表した僕にとってもただただ有益で幸せな時間でした。（他の皆様の発表もめちゃくちゃ勉強になりました！）

　午後からは，心理学の分野ごとにベイズモデリングの実践例を紹介するセッションが行われましたが，そのセッションで僕は「認知心理学への実践：データ生成メカニズムのベイズモデリング」というタイトルの発表を行いました。発表で使用したスライドの一部をSlideShareに公開しています（↓）。

認知心理学への実践：データ生成メカニズムのベイズモデリング by @mutopsy

　認知心理学の考え方とベイズモデリングの相性の良さについて冒頭で説明した後で，認知心理学におけるモデリングの実践例を3つ紹介し，そこから分かるベイズモデリングの利点を列挙していくという流れで発表をしました。実践例として，(1) 自由再生実験で得られる系列位置曲線の生成メカニズムをモデル化した記憶のSIMPLEモデルの紹介，(2) 心的回転に関する既存のモデルをベイズ流に拡張した事例の紹介，それから(3)心理学の概念モデルを確率モデルに翻訳し，直接的なモデル検証を試みた事例の紹介を行いました。それぞれのモデルの詳細には立ち入りませんでしたが，「ベイズモデリングを使えばこんなこともできるんだ！」という実感を少しでもお伝えできたのであれば幸いです。（なお，実践例の(2)と(3)は僕自身の未発表のアイディアとデータを含む内容でしたので，上記のスライドでは非公開にしております。）

　この発表で最も伝えたかったのは，実践例の中身ではなく，冒頭で説明した認知心理学とベイズモデリングの付き合い方に関する哲学じみた部分です。僕が本格的にベイズの勉強を始めたのは1年と数か月前ですが，ベイズを学び始めた当初はしばらくの間，僕の研究分野（認知心理学）においてベイズを用いる利点がいまいち見い出せずにおりました。ここ数年，ベイズは心理学でもずいぶん流行っているのですが，情報がかなり錯綜している印象で，ある種の「神話」じみた言説を未だに多く見かけます（「ベイズ万能論」のような風潮？）。僕自身，そうした「神話」がきっかけでベイズの扉を開くことができた訳ですが，ベイズをきちんと学ぶにつれてそういった神話に疑問を抱くようになり，どんどん幻想が崩れていくうちに，結局ベイズの何が良いのかが分からなくなってしまいました。けれど，ベイズ塾や学会・研究会・SNSを通じて出会った多くの方々とともに勉強を続け，Stanで遊び続けてきた結果，ようやく1つの答えにたどり着くことができました。それが冒頭の内容であり，この発表の芯の部分です。つまり今回の僕の発表は，ベイジアンとしての僕の（現時点での）集大成でもある訳です。これまでの学びを自分なりに整理できたという意味でも，このワークショップは僕にとって本当に意味のあるものでした。

Stanでポン・デ・リングを作ってみた

　こんにちは。mutopsyです。この記事は，Stan Advent Calendar 2017の25日目のエントリー記事です（実際にこの記事を書いたのは12月30日ですが）。

　みなさんはドーナツは何が好きですか？僕はフレンチクルーラーとエンゼルクリームがお気に入りですが，たまーにポン・デ・リングなんかも食べたくなります。このブログを訪れてくださった皆様の中には「Stanでポン・デ・リングが作れたら嬉しいのになあ……」とお考えになられた経験のある方もおられるでしょう。本日，なんとなくTogetterまとめを見ていたら，こんなまとめを見つけました。

→「おっぱい関数選手権」の名古屋代表の名大生さん、1つの数式だけでポン・デ・リングを表示することに成功「天才の所業」

　このまとめは，ポン・デ・リングの形状（三次元と二次元の両方）がたった1つの数式で表現可能であることを発見した@CHARTMANqさんのツイートをまとめたものです。このまとめを読んで甚く感動した僕は，さっそくStanで実装してみることにしました。Stanを乱数生成器として用いて，二次元のポン・デ・リング分布からのMCMCサンプルを発生させるという試みです。三次元でもできますが，とりあえず二次元ポン・デ・リング不等式を用いてみました。

　二次元ポン・デ・リング不等式は，極座標を用いて(r-2|cos4θ|-10)(r+|cos4θ|-6)≤0と表せます（上記まとめの@CHARTMANqさんのツイートをご参照ください）。極座標というのは，原点からの距離 (r) と角度 (θ) を使って空間上（ここでは平面上）の点の位置を表す座標系です。この不等式を満たすrとθの領域が 二次元のポン・デ・リングを表します。

　二次元ポン・デ・リング分布からの乱数（MCMCサンプル）を発生させるStanコードの例は以下の通りです。

Stanコード01: 二次元のポン・デ・リング分布から乱数を発生させるコード

parameters{
  real<lower=0> r;  //極座標のr
  real<lower=0,upper=2*pi()> theta;  //極座標のθ(0°〜360°の範囲)
}

transformed parameters{
  real<upper=0> D;  //上限を指定することで不等式を表現
  D = (r - 2*fabs(cos(4*theta)) - 10)*(r + fabs(cos(4*theta))-6); //二次元ポン・デ・リング（不）等式
}

generated quantities{
  real x;
  real y;
  x = cos(theta) * r;  //直交座標のxを計算
  y = sin(theta) * r;  //直交座標のyを計算
}

　今回はStanを乱数生成のためだけに用いているのでdataブロックは必要ありません。また，事前分布として一様分布を用いるため，modelブロックも不要です。Stanに慣れていると変な感じがするかもしれません。rとθの事前分布（一様分布）のうち，ポン・デ・リング不等式を満たすパラメータのみをMCMCサンプルとして生成することにより，二次元ポン・デ・リング一様分布からの乱数を発生させます。このモデルを実行するためのRコードは以下の通りです。

Rコード01: 上記モデルを走らせるためのRコード

library(rstan)
data <- list(Dummy = 1) #ダミーデータを作る（実際には使わない）
stanmodel <- stan_model(file="stancode/Ponde_2d.stan") #上記Stanコードを読み込む
fit <- sampling(
  stanmodel,data=data,seed=123,
  chains=4,
  iter=10000
)

　上記のコードではダミーデータを渡していますが，実際にはこのデータに意味はありません。Stanコードを走らせてみたところ，Rhatは全て1.01未満となり問題なく収束しました。各パラメータ (r, θ, x, y) の事後分布は以下のようになります。

　周辺分布だけを見ても，なんとなくポン・デ・リングっぽくありませんか。θの分布の凹凸やxとyの二峰分布具合はまさにポン・デ・リングですよね。ではいよいよxとyの同時分布を見てみましょう。

　どう見てもポン・デ・リングです。本当にありがとうございました。ぜひ皆さんもStanでポン・デ・リングを作ってみてください。もちろん数式が分かればポン・デ・リング以外のものも作れます。お試しあれ。それでは皆様，よいお年を！

　（「別にStanじゃなくてもよくね？」とつっこんではいけません。）

実験心理学者のためのベイジアンモデリング──正答率データをどうやって分析するか──

　こんにちは。mutopsyです。この記事は，Stan Advent Calendar 2017の14日目のエントリー記事です。この記事では，実験心理学でよく扱われる正答率・誤答率（あるいは正答数・誤答数）データの分析を例に，よくある素朴な分析（平均値の点推定・区間推定およびt検定）の背後にあるモデルを可視化した上で，そのモデルを改善した別の3つのモデルを紹介します。また，情報量規準によるモデル比較も行います。そして最後に，実験心理学におけるベイジアンモデリングの有用性について述べたいと思います。

目次
1. はじめに：正答率データを素朴に分析してみる
2. モデル1：制約のない正規分布モデル
3. モデル2：μの範囲を制限した正規分布モデル
4. モデル3：θを固定した二項分布モデル
5. モデル4：θが個人ごとに異なる二項分布モデル
6. 4つのモデルを比較する
7. まとめ

1. はじめに：正答率データを素朴に分析してみる

　改めまして，こんにちは。mutopsyこと武藤拓之です。大阪大学で認知心理学の研究をしています。この記事では，実験心理学でもそれ以外の分野でもお馴染みの，正答率（正答数）データを分析する方法について考えてみたいと思います。昨日はフォン・ミーゼス分布に関する若干マニアックな記事を投稿しましたが，この記事ではもう少しだけ需要のありそうな話をするつもりです。

　さて，いきなりですが，下のグラフ（図1）をご覧ください。このグラフは，筆者が過去に行った記憶に関する実験のデータを要約したものです¹。この実験では，実験参加者さんに2枚の画像を順番に見てもらい，2枚目の画像に写っている対象が1枚目の画像と同じであったか異なっていたかを回答してもらいました。この試行を繰り返すと，実験参加者さん一人ひとりの正答数のデータを得ることができます。正答数を合計試行数で割れば正答率になるので，正答率のデータが実験参加者の人数分だけ手に入ると考えることもできます。図1のグラフの青い点は，実験参加者一人ひとりの正答率をプロットしたものです。この実験では難易度の異なる2つの条件が設定されており，条件1と条件2にはそれぞれ異なる実験参加者が割り当てられたと考えてください。なお，この記事の内容は，正答数と合計試行数が分かっているデータであればどんなデータにも当てはまる話なので，読者の皆様にとって分かりやすいデータをイメージして頂ければ良いかと思います。

図1.　この記事で分析に用いる正答率データの平均値とその95%信頼区間。青い点は個人ごとの正答率。

　このデータのcsvファイルはこちらからダウンロードできます。→ accuracy_data.csv

　このデータは，各行が実験参加者1人分のデータに対応しており，以下の5つの変数を含んでいます。

• ID: 実験参加者のID（1から24までの整数）
• N.Corr: 正解した試行数（1から6までの整数）
• N.Total: 合計試行数（6で固定）
• Corr.Rate: 正答率（「正解した試行数」を「合計試行数」で割った値）
• Group: 条件群（1または2）

　実験心理学の文脈では，このようなデータが得られたら，まずは平均値の点推定値と誤差範囲を計算するのが一般的だと思います。ここでは各条件における平均正答率に関心があるので，正答率の算術平均（＝母平均の点推定値）と母平均の95%信頼区間を計算してみました。これらの値は先ほどの図1に折れ線グラフとエラーバーで示されています。赤色の四角形が各条件の平均正答率，エラーバーがその95%信頼区間です。信頼区間ではなく標準誤差をプロットする場合も多いですが，このようにして算術平均とエラーバーをグラフに描くのが一般的かと思います。データの傾向を確認したら，例えば2つの条件間で平均正答率に差があるか否かを検定したり（Welchのt検定），あるいはある条件の平均正答率がチャンスレベル（ここでは50%）よりも有意に高いと言えるかどうかを検定したり（一標本のt検定），研究仮説に合った分析を行います。実際にこのデータに対してWelchのt検定を行うと，p = .024 という結果が得られ，条件2の平均正答率は条件1の平均正答率よりも平均して23.6%高い，という結論が得られます。

　ただ，改めて図1のグラフを見ると，気になる点がいくつかあります。まず，条件2の平均正答率の95%信頼区間をよく見ると，区間の上側が1.0 (=100%) の値よりも上にあります。正答率は理論上0%から100%の範囲内に収まるはずなのに，平均正答率の母数が100%を超過しているかもしれないという奇妙な推定結果が得られてしまいました。また，実験参加者一人ひとりの正答率の散らばり方も気になります。例えば，条件1でも条件2でも正答率が100%の実験参加者がそれなりにいるように見えます。少なくとも正規分布に従っているようには見えません。たとえデータそのものが正規分布に従っていなくても，サンプルサイズが十分に大きければ平均値の分布は中心極限定理により正規分布に近付きますが，今回のように高々 n = 12 の場合には微妙な気もします。それに，やはり平均値の信頼区間が100%をまたいでしまうのは何となく気持ち悪いです。

　機械的に信頼区間を求めたりt検定を行ったり，あるいは分散分析したり回帰分析したり，ということを行ってしまうと，これらの分析が結局何をやっているのかが見えにくくなってしまいます。そこでまず，上述した分析の背後にあるモデルを可視化してみようと思います。そうすれば，先の分析の違和感の正体が掴めるはずです。

2. モデル1：制約のない正規分布モデル

　この節では，1節で説明したような素朴な分析が何をやっているのかを可視化するために，先の分析の背後にあるモデルをベイズ的に表現してみようと思います。具体的には平均値の区間推定の背後にあるモデルをStanで記述してみます（平均値差の推定については6節で触れます）。「頻度主義 vs. ベイジアン」という枠組みで考えてしまうと違和感を持たれるかもしれませんが，以下のモデルはいわゆる頻度主義の枠組みで想定されているモデルと（ほぼ）等価なモデルであり，後で確認するように，推定結果も（ほぼ）一致します。なお，この記事ではパラメータの解釈に関する頻度主義とベイズの違いについては触れません。前節で説明したようなよく使われる分析方法の背後にあるモデルは，以下に示す「モデル1：制約のない正規分布モデル」として記述することができます²。特に重要な行を赤色でハイライトしています。

Stanコード01: 【モデル1】制約のない正規分布モデル

data {
  int N;  //データの総数
  real CorrRate[N];  //正答率
  int Group[N];  //実験条件（群）
}

parameters {
  real mu[2];  //平均正答率 (μ)
  real<lower=0> sigma[2];  //標準偏差 (σ)
}

model {
  //モデル
  for (n in 1:N){
    CorrRate[n] ~ normal(mu[Group[n]],sigma[Group[n]]);
    //↑正答率データが正規分布に従う
  }
  
  //事前分布
  for (k in 1:2){
    mu[k] ~ uniform(-9999999,9999999); 
    sigma[k] ~ uniform(0,9999999);
  }
}

generated quantities{
  real diff_mu;
  diff_mu = mu[2] - mu[1];  //2群の平均正答率の差
}

　このモデルでは，条件1と条件2の正答率データがそれぞれ独立な正規分布から生成されると仮定しています。正規分布は平均値μと標準偏差σという2つのパラメータを持ちますが，今回の分析ではこの平均値μの分布が求めたい平均正答率の分布に相当します。以下，コードの解説です。

• 1-5行目: dataブロックでは分析に用いるデータを宣言します。ここではデータの総数（今回はデータが24行なのでN=24），正答率，実験条件（群）を使用します。
• 7-10行目: parametersプロックでは，このモデルで推定するパラメータを宣言します。ここでは正規分布のパラメータμとσを条件ごとに推定します。σは非負の値なので下限を0に指定しておきます（これを指定しないとうまく推定できません）。
• 15行目: 正規分布から一人ひとりの正答率データが生成されることを表現しています。この部分がこのモデルの核となります。
• 21-22行目: 事前分布を指定しています。この部分を省略した場合，Stanは十分に幅の広い一様分布を事前分布として採用します。なので本当はこの2行は不要なのですが，モデルの仮定を示すためにあえて明示的に書いています。伝統的な分析による推定結果は，-∞から+∞の範囲の一様分布³を事前分布としてベイズ推定した時の結果と一致します。ここで指定しているuniform(-9999999,9999999)は，-9999999から+9999999までの範囲の一様分布を表します。つまり，平均値がいくつになるか全く分からないという事前の信念を表現しています。本当はuniform(-Inf,Inf)とでも書きたいところですが，Stanではこのような書き方は（おそらく）できないので，あえて極端に範囲の広い事前分布を指定してみました。
• 26-29行目: 条件間の平均値差を計算しています。単にmu[0]とmu[1]の引き算をしているだけです。

　このモデルに先ほどのデータを当てはめて，推定したμの事後モード（MAP推定値）と95%最高密度区間を求めると，以下の図2のような結果が得られます。個人のデータのプロットはローデータそのものなので，値は図1と厳密に同じです（描画の都合上，横位置は若干ずれています）。折れ線グラフで表された平均値（μ）の点推定値（MAP推定値）とその95%最高密度区間は上記のコードで推定したものです。推定上の誤差が多少あるとはいえ，図1とほとんど同じ結果になっていることが分かるかと思います。

図2.　制約のない正規分布モデル（モデル1）で推定した平均正答率。エラーバーは95%最高密度区間を表す。

　改めてこのモデルを見てみると，まず事前分布の箇所が気になるかと思います。正答率の平均値は理論上0から1の範囲に収まるはずなのに，事前知識としてほぼ-∞から+∞の範囲の一様分布を用いるのは明らかに不自然です。素朴な信頼区間の計算やt検定の背後にはこのようなモデルが想定されているのです。そこで次のモデルでは，「平均正答率（μ）は0から1までの値しか取らない」という事前知識をパラメータの事前分布として反映させてみようと思います。

3. モデル2：μの範囲を制限した正規分布モデル

　モデル2は，正規分布の平均値パラメータμについて「0から1までの値しか取らない」という事前知識を与えたモデルとなります。先ほどのモデル1を僅かに修正するだけでこのモデルを作ることができます。以下のコードの赤色でハイライトされた部分のみがモデル1と異なる部分です。

Stanコード02: 【モデル2】μの範囲を制限した正規分布モデル

data {
  int N;  //データの総数
  real CorrRate[N];  //正答率
  int Group[N];  //実験条件（群）
}

parameters {
  real<lower=0, upper=1> mu[2];  //平均正答率 (μ)
  real<lower=0> sigma[2];  //標準偏差 (σ)
}

model {
  //モデル
  for (n in 1:N){
    CorrRate[n] ~ normal(mu[Group[n]],sigma[Group[n]]);
    //↑正答率データが正規分布に従う
  }
  
  //事前分布
  for (k in 1:2){
    mu[k] ~ uniform(0,1); 
    sigma[k] ~ uniform(0,9999999);
  }
}

generated quantities{
  real diff_mu;
  diff_mu = mu[2] - mu[1];  //2群の平均正答率の差
  }
}

• 8行目: parametersブロックで平均値パラメータを宣言する時に，あらかじめ下限0と上限1と指定しておきます。
• 21行目: 平均値パラメータの事前分布として，0から1の範囲の一様分布を指定しています。parametersブロックで平均値パラメータの下限と上限を指定しているので，この行を省略しても自動的に範囲[0,1]の一様分布が事前分布として指定されます。しかし，この上限と下限の制約がパラメータの事前分布であるということを強調するために，あえてここでは明示的に書いています。

　このモデルによる推定結果は以下の図3のようになります。目論見通り，条件2の区間推定幅の上端が100%を越えなくなりました。なお，この記事の下部にある図6にこの記事で用いた全てのモデルの結果を示しているので，比較の際にはこちらをご覧いただくと分かりやすいかと思います。

図3.　μの範囲を限定した正規分布モデル（モデル2）で推定した平均正答率。エラーバーは95%最高密度区間を表す。

　さて，平均値の範囲の問題は一応解決しましたが，そもそも正規分布を用いたモデルで本当に良いのだろうかという疑問は残ります。次の節では，正規分布ではなく二項分布を用いたモデルについて説明します。

4. モデル3：θを固定した二項分布モデル

　表が出る確率がθであるコインをn回投げた時に表が出る回数kは，pとnをパラメータに持つ二項分布で表現できます。もちろんコイン投げ以外のケースにも二項分布は使えます。例えば，記憶実験のある試行で正答する確率をθ，合計試行数をnと置けば，正答する回数kはやはり二項分布に従います（二項分布を使ったモベイズ推定については過去の記事でも紹介しています）。二項分布はnが大きくなるにつれて正規分布に近づきますが，今回分析に使用しているデータの合計試行数はたったの6なので，正規分布ではなく二項分布を使ったモデルのほうが適切であるように思われます（正規分布を用いるということは，二項分布モデルにおいてnが∞であると仮定することと一致します）。また，kとnは正の整数で，k > nとなることはなく，確率θは0から1の範囲に収まります。つまり，二項分布を用いる以上はモデル2のようにわざわざ範囲を制限する必要がありません（モデル上では制限をかけないと収束しませんが，この制限は事前の信念に基づく制限ではなく，モデルから自然に要求される制限であるという点で異なります）。

　それでは二項分布を使ったモデルについて考えてみましょう。上述したように，ある試行で正答する確率θ，合計試行数nの二項分布から正答回数が生じるモデルを想定しますが，まずは，θが条件のみによって決定するというモデルを考えます。つまり，条件1における正答確率をθ₁と置くと，条件1の実験参加者は，どの参加者であっても任意の試行においてθ₁の確率で正答し，1-θ₁の確率で誤答すると考えます。このパラメータθは平均正答率ではなく各試行で共通の正答確率を表すという点でモデル1とモデル2のμとは厳密には異なります。ただし，このモデルでは同じ条件の参加者でθが共通しているという仮定を置いているため，各条件における正答確率θの参加者間平均は正答確率θと一致します（θ * 12 / 12 = θ）。そのため，ここではθをμと同じように解釈することにします（研究仮説において両者を区別する必要が特になければ問題ありません）。このモデルのStanコードは以下のようになります。推定するパラメータがmuからthetaに変わっていることと，データが正規分布ではなく二項分布に従うと仮定していることが先のモデルとの違いです。

Stanコード03: 【モデル3】θを固定した二項分布モデル

data {
  int N;  //データの総数
  int CorrResp[N];  //正答数
  int TotalResp[N]; //合計試行数
  int Group[N]; //実験条件（群）
}

parameters {
  real<lower=0, upper=1> theta[2]; //正答確率 (θ)
}

model {
  //モデル
  for (n in 1:N){
    CorrResp[n] ~ binomial(TotalResp[n],theta[Group[n]]);
    //↑正答数データが二項分布に従う
  }
  
  //事前分布
  for (k in 1:2){
    theta[k] ~ uniform(0,1);
  }
}

generated quantities{
  real diff_theta;
  diff_theta = theta[2] - theta[1];  //2群の正答確率の差
}

• 3-4行目: 今度は正答率ではなく正答数と合計試行数をデータとして与えます。今回のデータでは合計試行数は6で固定ですが，この値が実験参加者によって違っていてもこのモデルは走ります。
• 9行目: 正規分布の時とは推定するパラメータが違うので，名前もmuではなくthetaにしています。thetaの範囲は0から1に指定しておかないとうまく収束しません。
• 15行目: 二項分布から正答数データが生成されるモデルを表現しています。

　このモデルを当てはめると以下の図4のような結果が得られます（他のモデルとの比較は図6参照）。縦軸が平均正答率（μ）ではなく正答確率（θ）になっていることに注意してください。パラメータのMAP推定値はモデル1・モデル2とあまり変わっていませんが，これまでのモデルよりも条件1における区間推定幅がかなり狭くなっていることが読み取れます。また，当然ながら区間推定幅は100%をまたぎません。

図4.　θを固定した二項分布モデル（モデル3）で推定した正答確率。エラーバーは95%最高密度区間を表す。

　この節では二項分布を使ったモデルを紹介しましたが，「θが実験参加者にかかわらず一定である」という仮定はなんとなく直感に反するように思えます。記憶成績が良い人も悪い人もいる，というように，θにも個人差があると考える方が自然でしょう。そこで，次の節ではθが人によって異なるモデルの一例を示します。

5. モデル4：θが個人ごとに異なる二項分布モデル

　この節では，二項分布のパラメータθが人によって異なることを仮定したモデルの一例を紹介します。以下にそのStanコードを示します。赤色でハイライトされた行はモデル3と異なる箇所です。

Stanコード04: 【モデル4】θが個人ごとに異なる二項分布モデル

data {
  int N;  //データの総数
  int CorrResp[N];  //正答数
  int TotalResp[N]; //合計試行数
  int Group[N]; //実験条件（群）
}

parameters {
  real<lower=0, upper=1> theta[N]; //正答確率 (θ)
}

model {
  //モデル
  for (n in 1:N){
    CorrResp[n] ~ binomial(TotalResp[n],theta[n]);
    //↑正答数データが二項分布に従う
  }
  
  //事前分布
  for (n in 1:N){
    theta[n] ~ uniform(0,1);
  }
}

generated quantities{
  real diff_mu;
  real mu[2];
  mu[1] = 0;
  mu[2] = 0;
  for (n in 1:N) mu[Group[n]] = mu[Group[n]] + theta[n]/12;
  diff_mu = mu[2] - mu[1];
}

• 9行目: 前節のモデル3ではパラメータthetaは条件によって異なると仮定していため，2つのthetaを推定すれば良かったのですが，今回のモデルでは実験参加者の人数分だけthetaを推定します。
• 15行目: 条件ごとのthetaではなく実験参加者ごとのthetaを二項分布のパラメータに指定します。
• 20-21行目: 事前分布も全てのthetaに対して指定します。
• 27-30行目: 推定したthetaの平均値muを条件ごとに計算しています。このmuは正規分布のパラメータμではなく，ただの算術平均です（thetaが分布するのでmuも分布します）。

　このモデルでパラメータを推定した結果を以下の図5に示します（モデルの比較は図6を参照）。縦軸は正答確率 (θ) の参加者間平均です。モデル3に比べると，条件2のパラメータの推定値が低くなっています。

図5.　θが個人ごとに異なる二項分布モデル（モデル4）で推定した正答確率の平均。エラーバーは95%最高密度区間を表す。

6. 4つのモデルを比較する

　さて，これまで4つのモデルを使って正答率データのモデリングを行ってきましたが，結局どのモデルを使うのが良いのでしょうか。図6は4つのモデルによる推定結果を並べたものですが，モデルによって結果が異なっている様子が見て取れます。

図6.　4つのモデルによる正答率（正答確率）の推定結果の比較。

　ここではモデル比較のために，情報量規準のWAICを用いることにします。WAICは予測の良さを表す指標で，値が小さいモデルほど真のモデルに近いことを意味します⁴。それぞれのモデルのWAICは図6に示した通りです⁵。この結果からは，モデルの良さは「モデル1 ＜ モデル2 ＜ モデル3 ＜ モデル4」であると言えそうです。やはり，素朴な分析の背後にあるモデルと等価な正規分布モデルはあまり適切ではないようです。この中ではθが人によって異なる二項分布モデルが最も良いモデルであるといえそうです。

　これまでは触れてきませんでしたが，各モデルにおける条件間の（平均）正答率の差（diff_muまたはdiff_theta）の推定結果も見てみましょう。以下の図7に結果を示しました。これを見ると，どのモデルにおいても（平均）正答率の差の区間推定幅には0が含まれていませんが，幅の広さにはだいぶ違いが見られます。モデル1の95%最高密度区間はWelchのt検定の計算過程で用いられる95%信頼区間と一致します。どうやら二項分布モデルのほうが推定精度が良さそうです。

図7.　4つのモデルによる正答率の条件間差の推定結果。

　ただし，可能なモデルは他にもあるので（例えばベータ二項分布モデルのような，θに分布を仮定するモデルや，感度と反応バイアスを区別する信号検出理論モデル，項目困難度と能力を区別する項目反応理論モデルなど），今回のモデル4があらゆるモデルの中で最良であるとは限りません。研究仮説に応じて適切なモデルを選択することも必要でしょう。

7. まとめ

　この記事では，素朴な分析（平均値の点推定・区間推定およびt検定）の背後にあるモデルを可視化した上で，そのモデルを修正した別の3つのモデルを紹介しました。また，モデルごとの推定結果の違いを概観し，WAICを用いたモデル比較も行いました。

　この記事で最も主張したかったことは，分析の背後にあるモデルを意識することの重要性についてです。機械的に分析を実行することに慣れてしまうと，明らかにデータに合わないモデルを無意識に用いてしまう危険性があります。Stanを使ってベイズ的にモデルを記述すれば，そのモデルが何を仮定しているのかは簡単に分かりますし，この記事で実際にお見せしたように，モデルの修正も容易に行えます。この記事で紹介したモデル1-4のパラメータは，ベイズを用いないで最尤法で解くことも可能だとは思いますが，ベイズ推定で用いられるMCMCアルゴリズムは強力かつ汎用的で，またStanを使えばモデルの記述も修正も拡張も比較的簡単に行えるという点で，非常に一般性が高いです。共通の文法を用いてあらゆるモデルを記述しパラメータを効率的に推定できるベイジアンモデリングはとても魅力的であるように思います。ベイズというと，パラメータの解釈における頻度主義との違いが強調されることが多いですが，個人的にはモデリングの柔軟性こそがベイズの良さであるように感じています。あなたもベイジアンモデラーになってみませんか。

脚注

¹本当はもっと多くの条件を設定していたのですが，ここでは2つの条件のデータのみを抽出しています。また，実際にはこの2条件は対応のあるデータなのですが，分かりやすさのために，ここでは対応のないデータとして扱うことにします。とはいえこの記事で書かれている内容は，もっと条件数の多い実験で得られたデータや対応のあるデータに対してももちろん適用可能です（Stanコードは適宜修正しなければなりませんが）。

²分かりやすさを優先し，各実験参加者の正答率が平均μ（＝平均正答確率）の正規分布に従うモデル，すなわち平均正答確率を直接推定するモデルを書きました。しかし，正答率ではなく正答数が正規分布に従うモデルを書き，generated quantitiesブロックで平均正答数を合計試行数で割って平均正答率を計算するほうが，モデルの適用範囲は広がります。また，渡すデータが後の二項分布モデル（モデル3・モデル4）と同じになるのでそのままWAICを計算できるようになります。

³厳密に言えば，この関数は全区間で積分しても1にならないので確率密度関数ではない。

⁴情報量規準は万能ではなく，あくまでも1つの判断基準でしかないのでご注意ください。

⁵モデル比較のためにWAICを計算する際には，どのモデルにも同じデータを与える必要があります。上述のモデル1とモデル2では正答率データを渡していたのに対し，モデル3とモデル4では正答数と合計試行数のデータを渡していたので，そのままWAICを計算すると，モデル1とモデル2のWAICが極端に小さくなってしまいます。そのため，モデル1とモデル2のWAICを計算する際は，正答確率ではなく正答数が正規分布に従うモデルに修正しなければなりません。この時，修正したモデルには合計試行数のデータを渡す必要はありません（正規分布モデルは合計試行数が∞の二項分布モデルと等しいため）。このことは関西学院大学の清水先生に教えて頂きました！

2018年1月25日：csvデータへのリンク切れを修正しました。