[Statistics]Bayesian Data Analysis Third edition

Link

Home page for the book, "Bayesian Data Analysis"

www.stat.columbia.edu

Part1 Fundamentals of Bayesian Inference

Chapter1 Probability and inference

1.1 The three steps of Baysian data analysis

1. 完全なモデル確率の準備
   ある問題における観測可能な量と観測不可能な量すべてに対する結合分布を設定する。このモデルは、基礎となる科学的問題やデータ収集プロセスに関する知識と整合している必要がある。
2. 観測データに対する条件付
   適切な事後分布（観測されたデータが与えられたときに、最も関心のある未観測量の条件付き確率分布）を計算し、解釈する。
3. モデルの適合度と、結果として得られる事後分布の意味を評価する。
   モデルがデータにどの程度適合しているのか、実質的な結論は妥当であるか、あるいはステップ1のモデリングの仮定に対して結果はどの程度敏感か、こういった疑問に対してモデルの変更あるいは拡張を行い、3つのステップを繰り返すことができる。

ベイジアン思考の主な動機の一つは、統計的結論の常識的な解釈を容易にすることである。例えば、未知の量に対するベイジアン（確率）区間は、その未知の量を高い確率で含んでいると直接的に解釈できる。これは頻度主義者の（信頼）区間とは対照的で、頻度主義者の信頼区間は、反復的な実践で行われる可能性のある類似の推論の系列に対してのみ厳密に解釈できる。

どういうことかというと、

1. ベイジアン統計学のアプローチでは、確率は未知の量（たとえば未来の事象や固定されたパラメータ）についての信念の度合いを表す。したがって、ベイジアンの確率区間（または信頼区間）は、その区間が真の未知の量を含む確率を直接的に示している。
2. 一方、頻度主義の統計学の視点では、確率は長期的な頻度または反復的な実験から得られる結果を表す。したがって、頻度主義者の信頼区間は、同じ実験を何度も繰り返した場合に、その区間が真の未知の値を含む頻度を示している。

この違いにより、ベイジアン統計学は結果を直感的に理解しやすいという特徴がある。なぜなら、ベイジアンの解釈は直接的で、未知の値についての現在の知識と信念を反映しているからである。

1.2 General notation for statistical inference

統計的推論とは、観測されていない量について、数値データから結論を導き出すことである。例えば、新しい抗がん剤の臨床試験は、その新薬を投与された集団の5年生存確率を、標準的な治療を受けている集団のそれと比較するようにデザインされるかもしれない。これらの生存確率は、患者の大規模な集団に言及しており、全集団で実験することは実行不可能であり、倫理的にも容認できない。したがって、真の確率、特にその差についての推論は、患者のサンプルに基づかなければならない。この例では、全集団をどちらか一方の治療にさらすことが可能であったとしても、誰も両方の治療にさらすことは不可能であるため、因果推論-各患者で観察された結果と、もう一方の治療にさらされた場合のその患者の観察されない結果との比較-を評価するためには、統計的推論が依然として必要である。

我々は2種類の推定量（統計的推論が行われる未観測量）を区別する。第一に、潜在的に観測可能な量、例えばプロセスの将来の観測値や、臨床試験の例で治療を受けなかった場合の結果などであり、第二に、直接観測できない量、つまり観測されたデータにつながる仮説的プロセスを支配するパラメータ（例えば回帰係数）である。これら2種類の推定量を区別することは必ずしも正確ではないが、特定の問題に対する統計モデルが現実世界にどのように適合するかを理解する方法として一般的に有用である。

Parameters, data, and predictions

約束として、\(\theta\)を測定不能なベクトル量あるいは興味のあるpopulationのパラメータとする。これは、例えばランダムに振り当てられた治療の生存確率などである。また、\(y\)を観測されたデータとする。これは、上の例だと割付群の生存者の数や死亡者の数である。\(\tilde{y}\)は未知だが観測可能である数量とする。例えば、既存のトライアルにおいて、新たにリクルートされた同様の患者に対する各々の治療効果である。本書全体の約束として一般的に、パラメータにギリシャ文字を用い、観測された・あるいは観測可能なスカラーやベクトル（ときには行列）に対して小文字のローマ字を用い、観察された・あるいは観測可能な行列に対して大文字のローマ字を用いる。行列表記を用いるときは、ベクトルは列ベクトルと考える。例えば、\(u\)が\(n\)次元のベクトルであるとき、\(u^Tu\)はスカラーで、\(uu^T\)は\(n\times n\)の行列になる。

Observational units and variables

多くの統計研究では、\(n\)個の対象物または単位についてデータを収集し、そのデータを\(y = (y_1, \cdots, y_n)\)というベクトルとして書くことができる。臨床試験の例では、患者\(i\)が\(5\)年後に生存していれば\(y_i\)を\(1\)、死亡していれば\(0\)とラベル付けすることができる。複数の変数が各ユニットで測定される場合、各\(y_i\)は実際にはベクトルであり、データセット\(y\)全体は行列（通常\(n\)行を持つ）となる。\(y\)変数は，「結果」と呼ばれ，推論を行うときに，サンプリング・プロセスや母集団の自然変動によって，変数の観測値が別の結果になった可能性を許容したいという意味で，「無作為」とみなされる．

Exchangeability

統計解析の出発（多くは暗黙のうちに）となる仮定が、\(y_i\)の交換可能性である。つまり、$$p(y_1, y_2, \cdots, y_n)$$はインデックスの順序に変わらず不変であるという仮定である。結果に関連する情報が説明変数ではなく、単位指数で伝達される場合は、非交換性モデルが適切であろう。交換可能性という考え方は統計学の基本であり、本書を通して何度も触れることになる。我々は一般的に、交換可能な分布からのデータを、分布\(p(\theta)\)を持つ未知のパラメータ・ベクトル\(\theta\)が与えられたとき、独立かつ同値に分布する（iid: independently and identically distributed）としてモデル化する。臨床試験の例では，生存の未知の確率である\(\theta\)を与えて，結果\(y_i\)をiidとしてモデルするかもしれない。

Explanatory variables

各ユニットには、わざわざランダムとしてモデル化する必要のないオブザベーションがあるのが通常である。臨床試験の例では、そのような変数には、試験の各患者の年齢と以前の健康状態が含まれるかもしれない。我々は，この\(2\)番目のクラスの変数を説明変数，または共変量と呼び，それらを\(x\)とラベルする．\(X\)をランダムとして扱うと、交換可能性の概念は、\((x, y)_i\)の\(n\)個の値の分布が、インデックスの任意の並べ替えによって変化しないことを要求するように拡張できる。インデックスがランダムに割り当てられていると考えることができる十分な関連情報を\(X\)に組み込んだ後、交換可能モデルを仮定することは常に適切である。2つのユニットが同じ\(x\)の値を持つなら、それらの\(y\)の分布も同じであるという意味で、\(x\)を与えられた\(y\)の分布が、調査中のすべてのユニットで同じであるということは、交換可能性の仮定から導かれる。説明変数\(x\)のどれでも、それらをモデルしたい場合は\(y\)カテゴリに移動できる。説明変数（予測変数ともいう）の役割について、調査、実験、観察研究の分析の文脈では第8章で、回帰モデルの文脈では本書の後半で詳細に議論する。

Hierarchical modeling

第\(5\)章とそれ以降の章では、階層モデル（マルチレベルモデルとも呼ばれる）に焦点を当てる。階層モデルでは、各レベルの単位で交換可能性を語ることができる。たとえば、\(2\)つの医療処置が、別々のランダム化実験で、複数の異なる都市の患者に適用されたとする。そのとき，他に情報がなければ，各都市内の患者を交換可能として扱うのが妥当であろ。そして、異なる都市からの結果自体も交換可能なものとして扱う。実際には、都市レベルの説明変数として、前述した個人レベルの説明変数と同様に、各都市について我々が持っているあらゆる関連情報を含めることが理にかなっており、そうすれば、これらの説明変数が与えられた条件付き分布は交換可能であろう。

1.3 Bayesian inference

ベイズ統計学におけるパラメータ\(\theta\)や未観測のデータ\(\tilde{y}\)に対する結論は確率的に記載される。これらの確率的な言及は、観測されたデータ\(y\)で条件付られており、\(p(\theta \mid y)\)あるいは\(p(\tilde{y}\mid y)\)と書かれる。多くの教科書に記載される統計的推測はレトロスペクティブな手順を用いた、真に未知の\(\theta\)に条件付けられたあり得る\(y\)の値に対しての\(\theta\)や\(\tilde{y}\)の推定に基づいており、ベイズ統計学の観測データに対する条件付は、この点で従前の統計解析から逸脱していると言える。このような違いがあるにもかかわらず、多くの単純な分析では、2つの統計的推論のアプローチから表面的には同じような結論が得られることがわかるだろう。しかし、ベイズ法を用いた分析は 、より複雑な問題にも容易に拡張することができる。

Probability notation

\(p(\cdot \mid \cdot)\)は条件付き確率、\(p(\cdot)\)は周辺確率である。しばしば”distribution”と”density”を同一の意味で用いる。文脈に応じて、混乱をさけるためイベント発生確率は\(\text{Pr}(\cdot)\)と示す。例えば、\(\displaystyle \text{Pr}(\theta > 2) = \int_{\theta > 2}{p(\theta)d\theta}\)などのようにである。正規分布のようによく知られた分布をつかつ時には、\(\theta \sim N(\mu, {\sigma}^2)\)のように表記する。ここで\(\mu\)は平均値で、\({\sigma}^2\)は分散である。一貫して、\(N(\mu, {\sigma}^2)\)と書いたときはランダム変数に対してであり、\(N(\theta | \mu, {\sigma}^2)\)と書いた時は密度関数について述べている。

Bayes’s rule

$$p(\theta \mid y) = \frac{p(\theta, y)}{p(y)} = \frac{p(\theta)p(y\mid \theta)}{p(y)} \tag{1.1}\label{eq1.1}$$

比がわかれば十分なので、\(p(y)\)を省略して以下の形で記載されることも多い。

$$p(\theta \mid y ) \propto p(\theta) p(y \mid \theta) \tag{1.2}\label{eq1.2}$$

\((1.2)\)式で重要なのは、\(p(y\mid \theta)\)が\(y\)ではなく\(\theta\)の関数となっていることである。このベイズの定理によって、観測データ\(y\)に基づきパラメータ\(\theta\)の推測を行うことが可能になる。つまり、\(p(\theta, y)\)についてのモデルの構築が\(p(\theta \mid y)\)を要約する計算を行うことによって可能になる。

Prediction

データ\(y\)を考慮する際、分布は未知であるが観測可能な場合、$$p(y) = \int{p(y, \theta)d\theta} = \int{p(\theta)p(y \mid \theta)d\theta} \tag{1.3}\label{eq1.3}$$となる。これは\(y\)の周辺分布marginal distributionと呼ばれるが、より情報量の高い名称として、prior predictive distributionという。つまり、以前の観測によって条件付けられてはいないが、観測可能な量の分布によって予測できるという点で、pridectiveというのである。データ\(y\)を観測した後は、未観測の\(\tilde{y}\)の予測が同様の手順によって可能となる。すなわち、$$\begin{eqnarray}p(\tilde{y}\mid y) & = & \int{p(\tilde{y}, \theta\mid y)d\theta}\\ & = & \int{p(\tilde{y}\mid \theta, y)p(\theta \mid y)d\theta}\\ & = & \int{p(\tilde{y}\mid \theta)p(\theta \mid y)d\theta} \tag{1.4}\label{eq1.4}\end{eqnarray}$$となる。こちらはposterior predictive distributionと呼ばれる。\(1\)行目から\(2\)行目の変形は$$\begin{eqnarray}p(\tilde{y}, \theta \mid y) & = & \frac{p(\tilde{y}, \theta, y)}{p(y)}\\ & = & \frac{p(\tilde{y}\mid \theta ,y)p(\theta, y)}{p(y)}\\ & = & p(\tilde{y}\mid \theta, y)p(\theta\mid y)\end{eqnarray}$$なる変形に基づく。\(2\)行目から\(3\)行目の変形は、与えられた\(\theta\)に対して\(\tilde{y}\)と\(y\)が互いに独立であることに基づく。

Likelihood

選択された確率モデルでベイズのルールを用いるということは、データ\(y\)が\(p(y|\theta)\)を通してのみ、事後推論\eqref{eq1.2}に影響を与えるということであり、この\(p(y|\theta)\)を定数\(y\)に対する\(\theta\)の関数とみなすと、尤度関数と呼ばれる。このように、ベイズ推論は尤度原理と呼ばれるものに従うことになる。尤度原理は合理的であるが、それは特定の分析に採用するモデルまたはモデル群の枠組みの中でのみ成り立つ。実際には、選択したモデルが正しいと確信できることはほとんどない。第6章では、サンプリング分布（データの繰り返される実測を想像すること）がモデルの仮定をチェックする上で重要な役割を果たしうることを説明する。実際、我々の考える応用ベイズ統計学者とは、様々な可能性のあるモデルの下でベイズの法則を適用することを厭わない人である。

Likelihood and odds ratios

あるモデルのもとで、点\(\theta_1\)と\(\theta_2\)で評価される事後密度\(p(\theta|y)\)の比は、\(\theta_2\)に対する\(\theta_1\)の事後オッズと呼ばれる。この概念の最も馴染み部会応用は、\(\theta_2\)が\(\theta_1\)の補数をみなされる離散パラメータの場合である。オッズは確率の代替表現を提供し、ベイズの法則がオッズで表現されるときに特に以下のような単純な形を取るという魅力的な特性を有する。$$\frac{p(\theta_1|y)}{p(\theta_2|y)} = \frac{p(\theta_1)p(y|\theta_1)/p(y)}{p(\theta_1)p(y|\theta_1)/p(y)} = \frac{p(\theta_1)p(y|\theta_1)}{p(\theta_2)p(y|\theta_2)} \tag{1.5}\label{eq1.5}$$つまり、事後オッズは事前オッズに尤度比（\(\frac{p(y|\theta_1)}{p(y|\theta_2)}\)）を掛けたものに等しい。

1.4 Discrete examples: genetics and spell checking

次にベイズの定理を、母集団全体を記述するパラメータの推定ではなく、特定の離散量についての推論が当面の目標である\(2\)つの例で示す。これらの離散的な例では、事前確率、尤度、事後確率を直接見ることができる。

Inference about a genetic status

ヒトの男性にはX染色体とY染色体が1本ずつ、女性にはX染色体が2本あり、それぞれの染色体は片親から受け継がれる。血友病はX染色体連鎖潜性遺伝を示す疾患であり、X染色体上に疾患の原因となる遺伝子を受け継いだ男性が罹患するのに対し、2本のX染色体のうち1本だけに遺伝子を持つ女性は罹患しない。一般に、このような遺伝子を2つ受け継いだ女性は致死的であり、このようなことはまれである。

Prior distribution

罹患した兄弟がいるある女性を考えてみる。つまり、彼女の母親は血友病遺伝子のキャリアであり、「良い」血友病遺伝子と「悪い」血友病遺伝子を1つずつ持っているはずである。彼女の父親は血友病に罹患しておらず、したがって女性自身が血友病遺伝子を持つ確率は半々である。興味のある未知の量、つまり女性の状態には2つの値しかない：女性は遺伝子のキャリアであるか（\(\theta = 1\)）、そうでないか（\(\theta = 0\)）。これまでに提供された情報に基づいて、未知量\(\theta\)の事前分布は単純に\(\displaystyle Pr(\theta = 1) = Pr(\theta = 0) = \frac{1}{2}\)と表現できる

Data model and likelihood

事前情報の更新に用いられるデータは、その女性の息子の罹患状態である。彼女に二人の息子があり、どちらも罹患していないとする。\(y_i = 1\)あるいは\(y_i = 0\)は、それぞれ罹患した息子または罹患していない息子を表すとする。二人の息子の結果は交換可能であり、未知の\(\theta\)を条件として独立である。\(2\)つの独立したデータは、以下の尤度関数を生成する。$$\begin{eqnarray}Pr(y_1 = 0, y_2 = 0|\theta = 1) & = & (0.5)(0.5) & = & 0.25 \\ Pr(y_1 = 0, y_2 = 0|\theta = 0) & = & (1)(1) &= & 1\end{eqnarray}$$これらの式は、もしその女性が保因者であれば、彼女の息子はそれぞれ50％の確率で遺伝子を受け継いで罹患するが、保因者でなければ、彼女の息子が罹患しない確率は1に近いという事実から導かれる（実際には、母親が保因者でなくても罹患する確率は0ではないが、このリスク（突然変異率）は小さいので、この例では無視できる）。

Poserior distribution

ベイズの法則は、データの情報と事前確率を組み合わせるために使用することができる。\(y\)を結合データ\(y_1, y_2\)に用いると、以下のように単純になる。$$\begin{eqnarray}Pr(\theta = 1|y) & = & \frac{p(y|\theta = 1)Pr(\theta = 1)}{p(y|\theta = 1)Pr(\theta = 1)+p(y|\theta = 0)Pr(\theta = 0)}\\ & = & \frac{(0.25)(0.5)}{(0.25)(0.5)+(1.0)(0.5)}\\ & = & \frac{0.125}{0.625}\\ & = & 0.20\end{eqnarray}$$直感的には、ある女性が罹患していない子供を持つ場合、その女性が保因者である可能性が低くなることは明らかであり、ベイズの法則はその補正の程度を決定するための正式なメカニズムを提供する。この結果は、事前オッズと事後オッズで表すこともできる。その女性が保因者である事前オッズは\(\displaystyle \frac{0.5}{0.5} = 1\)である。罹患していない2人の息子に関する情報に基づく尤度比は\(\displaystyle \frac{0.25}{1} = 0.25\)なので、事後オッズは\(1\cdot 0.25 = 0.25\)となる。確率に直すと、先ほどと同じように\(\displaystyle \frac{0.25}{(1+0.25)} = 0.2\)となる。

Adding more data

ベイズ分析の重要な点は、逐次分析が簡単にできることである。例えば、この女性に罹患していない三男が生まれたとしよう。計算全体をやり直す必要はない。むしろ、前の事後分布を新しい事前分布として使い、次のようになる。$$Pr(\theta = 1|y_1, y_2, y_3) = \frac{(0.5)(0.2)}{(0.5)(0.2) + (1)(1-0.2)} = 0.111$$あるいは、三男が罹患していると仮定すれば、女性が保因者である事後確率が\(1\)になることを確認するのは簡単である（ここでも突然変異の可能性は無視する）。

Spelling correction

単語の分類は不確実性を管理する問題である。例えば、誰かが「radom」とタイプしたとする。どう読むべきだろうか？「random」や「radon」などのスペルミスやミスタイプかもしれないし、（この段落で最初に使われたように）意図的に「radom」とタイプしたのかもしれない。「radom」が実際にランダムを意味する確率はどのくらいだろうか？\(y\)をデータ、\(\theta\)をその人が意図的にタイプした単語とすると、次のようになる。$$Pr(\theta|y = ‘radom’) \propto p(\theta)Pr(y = ‘radom’|\theta) \tag{1.6}\label{eq1.6}$$である（\eqref{eq1.2}）。この積が非正規化事後密度である。この場合、簡単のために、意図する単語\(\theta\)について3つの可能性（random、radon、radom）だけを考えるとすると、まず\(\theta\)の3つの値すべてについて非正規化密度を計算し、次に正規化することで、関心のある事後確率を計算することができる。$$p(random|’radom’) = \frac{p(\theta_1)p(‘radom’|\theta_1)}{\sum_{j=1}^{3}{p(\theta_j)p(‘radom’|\theta_j)}}$$ここで、\(\theta_1 = random, \theta_2 = radon, \theta_3 = radom\)とする。事前確率\(p(\theta_j)\)は、最も単純に、ある大規模なデータベースにおけるこれらの単語の頻度から得ることができ、理想的には、手元の問題に適合したものである（例えば、問題の単語がそのような文書に出現している場合、最近の学生の電子メールのデータベースなどを利用できるだろう）。尤度\(p(y|\theta_j)\)はスペルミスやタイピングエラーのモデリングから得ることができる。

Prior distribution

他の文脈がなければ、いくつかのデータベースにおけるこれら3つの単語の相対的な頻度に基づいて事前確率\(p(\theta_j)\)を割り当てることは理にかなっている。以下はGoogleの研究者から提供された確率である。

これらの可能性だけを考えているので、\(3\)つの数字の合計が\(1\)になるように再正規化することもできる（\(\displaystyle p(random) = \frac{760}{760+60.5+3.12}\)など）が、補正しても\eqref{eq1.6}の比例定数に吸収されるだけなので，その必要はない。上の表に戻ると、言語資料に含まれる「radom」の確率がこれほど高いことに驚く。ウィキペディアでこの単語を調べてみると、「ポーランドで最大かつ最高の観客動員数を誇る航空ショー … … ポーランドの設計による半自動9mmパラピストルの非公式名称 … … 」が開催される中規模の都市であることがわかった。私たちが目にした文書では、「ラドム」の相対的な確率は高すぎるように思われる。上記の確率が私たちの用途に適切でないと思われる場合、それは私たちがまだモデルに含めていない事前情報や信念を持っていることを意味する。この点については、まずこの例に対するモデルの意味を理解した後で、話を戻そう。

Likelihood

Googleのスペルミスやタイプミスのモデルから、いくつかの条件付き確率を紹介しよう。

この尤度関数は確率分布ではないことを強調しておく。そうではなく、未知のパラメータ\(\theta\)の\(3\)つの異なる可能性に対応する、3つの異なる確率分布からの特定の結果（’radom’）の条件付き確率の集合である。

つまり、この\(5\)文字の単語が正しく入力される確率は97％、「random」から誤って1文字抜けてこの文字列が得られる確率は0.2％、「radon」の最後の文字をミスタイプして得られる確率はもっと低い、ということである。これらの確率について強い直感はないので、ここではグーグルのエンジニアを信頼することにする。

Posterior distribution

事前確率と尤度を掛け合わせて結合確率を求め、再正規化して事後確率を求める。

したがって、このモデルを条件として、タイプされた単語「radom」が正しい可能性は、「random」のタイプミスである可能性の約2倍であり、「radon」の間違いである可能性は非常に低い。もっと詳しく分析すれば、この3語以外の可能性も含まれるが、基本的な考え方は同じである。

Decision making, model checking, and model improvement

ここからは2つの方向性が考えられる。ひとつは、この単語が正しくタイプされた確率の3分の2を受け入れるか、あるいは単純に「radom」を一発で正解とする方法である。第二の選択肢は、たとえば「radom」はタイプミスのように見えるし、正しい確率の推定値は高すぎると思う、と言ってこの確率に疑問を呈することである。

事後分布の主張に異議を唱える場合、そのモデルがデータに適合していない、あるいはこれまでのモデルに含まれていない追加的な事前情報があると言っている。したがって、事後分布をめぐる論争は、事前分布か尤度のどちらかに追加情報があるという主張に対応しなければならない。

この問題については、尤度を批判する根拠は特にない。一方、事前確率は文脈に大きく依存する。’random’という単語はもちろん我々の統計学に関する著書の中でも頻出であり、’radon’は時々登場するが（第9.4節参照）、’radom’は我々にとって全く新しい単語である。「radom」の確率の高さに対する驚きは、我々の特定の問題に関連する新たな知識を意味する。

このモデルは、事前確率に文脈情報を含めることによって、最も即座に精緻化することができる。例えば、調査対象の文書が統計の本であれば、「ランダム」と入力した可能性が高くなる。xをモデルで使われる文脈情報とすると、ベイズの計算は次のようになる。$$p(\theta|x, y) \propto p(\theta|x)p(y| \theta, x)$$

第一近似として、最後の項を\(p(y|\theta)\)と単純化し、特定のエラーの確率（つまり、意図した単語\(\theta\)が与えられたときに特定の文字列\(y\)をタイプする確率）が文脈に依存しないようにすることができる。これは完全な仮定ではないが、モデル化と計算の負担を減らすことができる。

ベイズ推論における現実的な課題は、データからこれらの確率をすべて推定するモデルを設定することである。その時点で、上に示したように、ベイズの法則を簡単に適用して、手元の問題に対するモデルの意味を決定することができる。

1.5 Probability as a measure of uncertainty

我々はすでに確率密度などの概念を用いており、実際、読者が基本的な確率論にかなりの程度精通していることを前提としている（ただし、セクション1.8では、ベイズ分析でしばしば生じるいくつかの確率計算について簡単な技術的レビューを行う）。しかし、ベイズの枠組みにおける確率の用途は、非ベイズ統計学よりもはるかに広いので、より詳細な統計的例を検討する前に、少なくとも確率の概念の基礎を簡単に検討することが重要である。確率とは、「結果」の集合で定義される数値量であり、非負であり、互いに排他的な結果に対して加法的であり、可能なすべての互いに排他的な結果に対して和が\(1\)になるものである。

ベイズ統計学では、確率は不確実性の基本的な尺度や基準として用いられる。このパラダイムの中では、「明日雨が降る」確率を論じることも、サッカーワールドカップでブラジルが優勝する確率を論じることも、コイン投げで表が出る確率を論じることと同様に正当である。したがって、未知の推定値が特定の値の範囲にある確率を考えることは、100の大きさを持つ既知の固定母集団から10個の項目を無作為に抽出したサンプルの平均値が特定の範囲にある確率を考えるのと同じくらい自然なことなのである。これら2つの確率のうち、最初の確率はデータが取得された後の方がより興味深いのに対して、2番目の確率は事前の方がより適切である。ベイズ法は、ある状況や「自然の状態」（観察できない、あるいはまだ観察されていない）に関する（データに基づく）部分的な知識について、確率を尺度として体系的に記述することを可能にする。指針となる原則は、未知のものに関する知識の状態は確率分布によって記述されるということである。

不確実性の数値的尺度とはどういう意味だろうか？例えばコイン投げで「表」が出る確率は\(\displaystyle \frac{1}{2}\)であると広く合意されている。これはなぜだろうか？一般的に2つの正当性が与えられているようである。

1. 対称性または交換可能性の議論
   $$\text{probability} = \frac{\text{number of favorable cases}}{\text{number of possibilities}}$$同じようにありそうな可能性を仮定する。コイントスの場合、これは本当に物理的な議論であり、コインの落ち方を決定する際に働く力と、トスの初期物理的条件についての仮定に基づいている。
2. 頻度論
   確率=物理的に互いに独立した同一の方法で行われると仮定された、長い連続したトスで得られる相対的な頻度。

上記の議論はどちらもある意味で主観的であり、コインの性質やトスの手順についての判断を必要とし、また、どちらも等しく起こりそうな出来事、同一の測定値、独立性の意味についての意味論的議論を含んでいる。頻度の議論は、同じトスの長い連続という仮定の概念を伴うという点で、ある種の特別な難しさがあるように思われるかもしれない。厳密に考えれば、この観点は、少なくとも概念的には、同一の事象の長い連続の中にたまたま組み込まれていない一回のコイントスに対する確率の記述を許さない。

以下の例は、確率判断がいかに主観的になりやすいかを示している。まず、次のようなコイン実験を考えてみよう。あるコインが両方とも表か両方とも裏のどちらかであり、それ以上の情報は何もないとする。それでもなお、表の確率を語ることができるだろうか？一般的な言い回しでは間違いなく可能である。この新しい確率をどのように評価するかは、あまり明確ではないかもしれないが、おそらく「表」と「裏」というラベルの交換可能性に基づいて、多くの人が同じ値\(\displaystyle \frac{1}{2}\)に同意するだろう。

では、さらにいくつかの例を考えてみよう。サッカーでコロンビアがブラジルと対戦するとしよう。 コロンビアが勝つ確率は？明日雨が降る確率は？明日雨が降った場合、コロンビアが勝つ確率は？指定されたロケットの打ち上げが失敗する確率は？これらの質問はいずれも常識的には妥当であるように思えるが、参照されている確率について強力な度数解釈を考えるのは難しい。しかし、頻度による解釈は通常構築可能であり、これは統計学において非常に有用なツールである。例えば、将来のロケットの打ち上げを、同じ種類の潜在的な打ち上げの母集団からのサンプルとみなし、失敗した過去の打ち上げの頻度を調べることができる（この例の詳細については、本章末の参考文献を参照）。このようなことを科学的に行うことは、確率モデル（あるいは、少なくとも、比較可能な事象の「参照集合」）を作成することを意味し、このことは、問題の結果を交換可能であり、したがって同じように可能性が高いと考えなければならない、単純なコイン投げに類似した状況に私たちを引き戻す。

なぜ確率は不確実性を定量化する合理的な方法なのか？よく次のような理由が挙げられる。

1. 類推によれば、物理的なランダム性は不確実性を誘発するので、ランダムな事象の言語で不確実性を記述することは合理的であると思われる。一般的な会話では、「おそらく」や「ありそうもない」といった用語が多く使われており、科学的推論の問題に、より正式な確率計算を拡張することは、このような用法と矛盾しないように思われる。
2. 公理的または規範的アプローチ：決定理論に関連するこのアプローチは、すべての統計的推論を、利益と損失を伴う意思決定の文脈に置く。そして、合理的な公理（順序性、推移性など）により、不確実性は確率で表現されなければならない。この規範的根拠は示唆に富んでいるが、説得力はない。
3. ベットの一貫性。ある事象に（あなたが）付随する確率\(p\)を、以下のように定義する。イベント\(E\)が発生した場合に\($1\)のリターンと\($p\)を交換（つまりベット）する割合\(p \in [0, 1]\)とする。つまり、\(E\)が発生した場合、あなたは\($(1-p)\)を得るが、\(E\)が発生しなかった場合、あなたは\($p\)を失う。例えば、
   ・コイントス：コイントスを公平な賭けと考えると、次のような偶数オッズが示唆される。\(\displaystyle p=\frac{1}{2}\)
   ・試合のオッズ：ある試合で、AチームがBチームに対して10対1のオッズで勝つことに賭ける場合（つまり、10勝つために1を賭ける）、Aチームが勝つ「確率」は少なくとも1/11である。

首尾一貫性の原則とは、起こりうるすべての出来事に対する確率の割り当てが、賭けによって確実な利益を得ることができないようなものでなければならない、というものである。この原則のもとで構成される確率は、確率論の基本公理を満たさなければならないことが証明できる。ベッティングの根拠には根本的な問題がある。
・正確なオッズは、すべてのイベントにおいて、どちらかの方向に賭けることができるものでなければならない。正確なオッズがわからないのに、どうやってオッズを決めるのか？
・ある人があなたと賭けることを望んでいて、あなたが知らない情報を持っている場合、あなたがその賭けに乗るのは賢明ではないかもしれない。実際には、確率はベッティングの不完全な（必要ではあるが十分ではない）指針である。

これらの考察はすべて、確率が応用統計における不確実性を要約するための合理的なアプローチである可能性を示唆しているが、最終的な証明は応用の成功にある。本書の残りの章では、確率が統計応用における不確実性を扱うための豊かで柔軟な枠組みを提供することを示す。

Subjectivity and objectivity

確率を用いる統計的手法はすべて、世界の数学的理想化に依存するという意味で主観的である。ベイズ法は、事前分布に依存しているため、特に主観的であると言われることがあるが、ほとんどの問題では、モデルの「尤度」と「事前」の両方の部分を特定するために科学的判断が必要である。たとえば，線形回帰モデルは，一般に，回帰パラメータについて仮定されるかもしれない事前分布と少なくとも同じくらい疑わしい．多くの交換可能なユニットが観察されるという意味で、再現性があるときはいつでも、データから確率分布の特徴を推定し、分析をより「客観的」にする余地がある。実験全体が数回繰り返されるのであれば、第5章で議論したように、事前分布のパラメータはデータから推定することができる。しかし、いずれにせよ、科学的判断を必要とする要素、特に分析に含まれるデータの選択、分布に仮定されるパラメトリックな形式、モデルのチェック方法は残る。

1.6 Example: probabilities from football point spreads

経験的データともっともらしい実質的仮定を用いて確率を割り当てる方法の例として、プロ（アメリカン）フットボールの試合における特定の結果の確率を推定する方法を考える。これは確率の割り当ての例であり、ベイズ推論の例ではない。主観的な評価を行う、観測データに基づく経験的確率を用いる、パラメトリック確率モデルを構築する、などである。

Football point spreads and game outcomes

サッカーの専門家は、両チームの実力差を測る尺度として、サッカーの試合ごとにポイントスプレッドを提示する。このポイントスプレッドの意味するところは、お気に入りのチームAが、格下のチームBに4点差以上で勝つという命題は、フェアなベットとみなされるということである。言い換えれば、Aが3.5点差以上で勝つ確率は\(\displaystyle \frac{1}{2}\)である。ポイントスプレッドが整数である場合、チームAはポイントスプレッドよりも多くのポイントで勝つ可能性が、ポイントスプレッドよりも少ないポイントで勝つ（または負ける）可能性と同じくらい高いということになる。ポイントスプレッドの割り当て自体は、確率論的推論における興味深い練習である。1つの解釈は、ポイントスプレッドは、ギャンブル集団のゲームの可能な結果についての信念の分布の中央値であるということである。この例の残りの部分では、ポイントスプレッドを与えられたものとして扱い、それがどのように導かれたかは気にしない。

1981年、1983年、1984年のシーズンに行われた672のプロフットボールゲームのポイントスプレッドと実際のゲーム結果を図1.1にグラフ化した(1982年シーズンの多くは労働争議のため中止された)。散布図の各ポイントは、ポイントスプレッド\(x\)と実際の結果(favoriteのスコアからunderdogのスコアを引いたもの)\(y\)を表示する(ポイントスプレッドが0のゲームでは、「favorite」と「underdog」のラベルはランダムに割り当てられた)。グラフ上の各点の x 座標と y 座標には、複数の点がぴったり重ならないように、小さなランダムなジッターが加えられている。

Assigning probabilities based on observed frequencies

特定の事象に確率を割り当てることは興味深い： Pr(好きなチームが勝つ)、Pr(好きなチームが勝つ｜ポイントスプレッドは3.5ポイント)、Pr(好きなチームがポイントスプレッド以上に勝つ)、Pr(好きなチームがポイントスプレッド以上に勝つ｜ポイントスプレッドは3.5ポイント)などである。新聞を読んだりフットボールの試合を見たりして得た非公式な経験に基づく主観的な確率を報告することもある。有利なチームが試合に勝つ確率は確かに0.5より大きいはずで、おそらく0.6から0.75の間だろう。より複雑な事象になればなるほど、直感や知識が必要になる。より体系的なアプローチは、図1.1のデータに基づいて確率を割り当てることである。同点の試合を2分の1の勝ちと2分の1の負けとして数え、ポイントスプレッドがゼロの試合を無視すると（したがって、お気に入りは存在しない）、次のような経験的推定値が得られる。

・Pr（好きなチームが勝つ）=\(\displaystyle \frac{410.5}{655} = 0.63\)
・Pr（好きなチームが勝つ|\(x = 3.5\)） \(\displaystyle =\frac{36}{59} = 0.61\)
・Pr（好きなチームがスプレッドポイント差以上で勝つ） \(\displaystyle =\frac{308}{655} = 0.47\)
・Pr（好きなチームがスプレッドポイント差以上で勝つ|\(\displaystyle x = 3.5 = \frac{32}{59} = 0.54\)）

これらの経験的な確率の割り当ては、知識豊富なサッカーファンの直感と一致するという点で、どれも理にかなっているように見える。しかし、このような確率の割り当ては、直接関連するデータポイントが少ない事象については問題がある。例えば、8.5ポイント人気はこの3年間に5回中5回勝ったのに対し、9ポイント人気は20回中13回勝った。しかし、現実的には8.5ポイント人気よりも9ポイント人気の方が勝つ確率は高いと予想される。ポイントスプレッド8.5ではサンプルサイズが小さいため、確率の割り当てが不正確になる。パラメトリックモデルを使った別の方法を考える。

A parametric model for the difference between outcome and point spread

図1.2aは、フットボール・データセットのゲームについて、観察されたゲーム結果とポイント・スプレッドの差\(y-x\)を、ポイント・スプレッドに対してプロットしたものである（ここでも両方の座標にランダム・ジッターを加えた）。このプロットは、\(y-x\)の分布を\(x\)から独立したものとしてモデル化することがおおよそ妥当であることを示唆している（練習問題6.10参照）。図1.2bは、すべてのフットボール・ゲームの差\(y-x\)のヒストグラムであり、フィットされた正規密度が重ねられている。このプロットは，確率変数 \(d = y-x\) の周辺分布を正規分布で近似することが妥当かもしれないことを示唆している。\(d\)の\(672\)値の標本平均は\(0.07\)であり、標本標準偏差は\(13.86\)であり、フットボールの試合結果は、平均がポイント・スプレッドに等しく、標準偏差がほぼ\(14\)ポイント（2回のタッチダウン）のほぼ正規分布であることを示唆している。残りの議論では、\(d\)の分布は\(x\)に依存せず、各\(x\)について平均\(0\)、標準偏差\(14\)の正規分布であるとする。つまり、$$d|x \sim \text{N}(0, 14^2)$$であり、図1.2bに示された通りである。実際のデータでもよくあることだが、フットボールのスコアもポイントスプレッドも連続値ではない。

Assigning probabilities using the parametric model

とはいえ、このモデルは事象に確率を割り当てるのに便利な近似値を提供する。\(d\)が平均\(0\)の正規分布を持ち、ポイント・スプレッドから独立している場合、好きなチームがポイント・スプレッドより多く勝つ確率は\(\displaystyle \frac{1}{2}\)であり、ポイント・スプレッドのどの値に対しても条件付きであり、したがって無条件でも同様である。正規モデルによって得られる確率を\(Pr_{norm}\)とすると、\(x\)ポイントのお気に入りがゲームに勝つ確率は、正規モデルを仮定して次のように計算できる。$$Pr_{norm}(y>0|x) = Pr_{norm}(d > -x|x) = 1-\Phi\left(-\frac{x}{14}\right)$$ここで、\(\Phi\)は標準正規累積関数である。例えば、
・\(\displaystyle Pr_{norm}\)（好きなチームが勝つ| \(x = 3.5\)）\(=0.60\)
・\(\displaystyle Pr_{norm}\)（好きなチームが勝つ| \(x = 8.5\)）\(=0.73\)
・\(\displaystyle Pr_{norm}\)（好きなチームが勝つ| \(x = 9.5\)）\(=0.74\)

となる。3.5ポイント人気の確率は、先に示した経験値と一致するが、8.5ポイントおよび9ポイント人気の確率は、小さなサンプルに基づく経験値よりも直感的に理解できる。

1.7 Example: calibration for record linkage

確率をデータから推定する別の例を用いて、確率の本質的な経験的性質（「主観的」でも「個人的」でもない）を強調する。レコード・リンケージとは、異なるデータベースから同じ個人に対応するレコードを識別するためのアルゴリズム技術の使用を指す。レコード・リンケージ技術は様々な場面で使用されている。ここで説明する作業は、米国国勢調査と大規模な国勢調査後の調査との間のレコード・リンケージの文脈で定式化され、最初に適用されたものである。この最初のステップの目標は、過度のエラー率を発生させることなく、コンピュータによって可能な限り多くのレコードを「一致」と宣言することであり、それによって「一致」と宣言されなかったすべてのレコードについて、結果として発生する手作業による処理のコストを回避することである。

Estimating match probabilities empirically

レコード・リンケージの文献では、多変量レコードの個々の情報フィールドに「重み」を割り当て、2つのレコード間の一致の近さを要約する、我々が\(y\)と呼ぶ「スコア」を得る問題に多くの注意が払われてきた。ここでは、これらのルールが選択されたという意味で、このステップが完了したと仮定する。次のステップは一致候補ペアの割り当てであり、レコードの各ペアはそれぞれのデータベースから互いに最もよく一致する可能性のあるものから構成される。そして、指定された重み付けルールにより、一致候補ペアの順序付けが行われる。国勢調査局における動機付け問題では、「一致を宣言する」対「追跡調査へ送る」の二者択一が行われ、レコードが一致を宣言されるカットオフ・スコアが必要である。偽一致率は、偽に一致したペアの数を宣言された一致ペアの数で割ったものとして定義される。このような決定問題で特に重要なのは、マッチしたペアの候補が正しいマッチである確率を、そのスコアの関数として評価する正確な方法である。スコアを確率に変換する簡単な方法は存在するが、これらは極めて不正確で、典型的には極めて楽観的な偽マッチ率の推定につながる。例えば、公称偽マッチ確率が\(10^{-3}\)から\(10^{-7}\)の範囲にあるレコード（つまり、ほぼ確実にマッチするとみなされるペア）を手作業でチェックしたところ、実際の偽マッチ率は1％の範囲に近かった。公称偽マッチ確率が1％の記録では、実際の偽マッチ率は5％であった。

サッカーの例で、ポイント・スプレッドを条件として異なる試合結果の確率を推定するために過去のデータを使ったのと同じようにである。我々のアプローチは、得点\(y\)を用い、\(y\)の関数として試合の確率を経験的に推定することである。

Existing methods for assigning scores to potential matches

第22章で詳述する混合モデリングを用いて、正確なマッチ確率を求める。以前に得られたマッチ候補のスコア分布は、真のマッチのスコア分布と非マッチのスコア分布の「混合」であると考えられる。混合モデルのパラメータはデータから推定される。推定されたパラメータにより、スコアに関する任意の決定閾値に対する偽の一致（真の一致でない一致と宣言されたペア）の確率の推定値を計算することができる。実際に使用された手順では、混合モデルのいくつかの要素（例えば、正規分布の混合を適用できるようにするために必要な最適変換）は、（我々のキャリブレーション手順を適用するデータとは別の）既知の一致状態を持つ「訓練」データを使用して適合されたが、ここではそれらの詳細については説明しない。その代わりに、マッチ状態が未知のデータセットに対して、この方法をどのように使用するかに焦点を当てる。

我々のアプリケーション・データセットでは、マッチの状態はわからない。したがって、我々は、2つの成分分布と、各成分に属する得点の母集団の割合を推定する1つの結合ヒストグラムに直面する。混合モデルの下では、得点の分布は次のように書くことができる。$$p(y) = Pr(match)p(y|match) + Pr(non-match)p(y|non-match) \tag{1.7}\label{eq1.7}$$

混合確率(Pr(match))と、マッチ(p(y|match))と非マッチ(p(y|non-match))の分布のパラメータは、(第22章で説明する)混合モデルのアプローチを、マッチの状態が不明なデータから得られた結合ヒストグラムに適用して推定する。この方法を用いてレコード・リンケージの決定を行うために、決定しきい値（ペアがマッチすると「宣言」されるスコア）の関数として偽マッチ率を与える曲線を構築する。与えられた決定閾値に対して、\eqref{eq1.7}の確率分布を用いて、分布\(p(y|non-match)\)に由来する閾値以上のスコア\(y\)が偽一致となる確率を推定することができる。閾値が低いほど、より多くのペアをマッチとして宣言することになる。より多くのマッチを宣言すると、エラーの割合が増加する。ここで説明するアプローチは、各閾値について客観的なエラー推定を提供するはずである。(次の段落の検証を参照)そうすれば、意思決定者は、より多くのマッチを自動的に宣言する(したがって、事務的な労力を減らす)という目標と、より少ないミスをするという目標の間で、許容できるバランスを提供する閾値を決定することができる。

External validation of the probabilities using test data

上記の方法は、照合状況が判明しているデータを用いて外部で検証された。この方法は1988年の国勢調査の3つの異なる場所からのデータに適用されたので、3つのテストが可能であった。他の2つの結果も同様であった。混合モデルは、テスト地点におけるすべての候補ペアの得点に当てはめられた。図1.4は、閾値が高い（マッチしない）ケースから低い（ほとんどすべての候補ペアをマッチすると宣言する）ケースまで、マッチすると宣言されたケースの割合で、予想される偽マッチ率（および不確実性の境界）を示している。偽マッチの割合は、宣言されたマッチの数の増加関数であり、これは理にかなっている。グラフ上で右方向に進むにつれて、我々はより弱いケースをマッチであると宣言している。

図1.4上の線は、モデルから推定される偽マッチ率の期待割合と95%事後境界を示している（これらの境界は、偽マッチ率が95%事後確率の範囲内にあると推定される範囲を示している）。(これらの境界は、偽マッチ率が95%事後確率で存在する推定範囲を示す。事後区間の概念については次章で詳しく説明する)。グラフの点は実際の偽マッチ率を示しており、これはモデルとよく一致している。特に、モデルは、ほとんどの誤照合を避けるために、90％以下のケースを照合とし、残りの10％程度をあきらめることを推奨しており、点も同様のパターンを示している。

ファイルの大部分をほとんどエラーなくマッチさせることは明らかに可能である。また、偽マッチ率が加速するある時点で、マッチ候補の質は劇的に悪くなる。図1.5は、偽マッチ率が加速する関心領域における較正手順の動作を強調するために、先ほどの表示を拡大鏡で見たものである。予測された偽マッチ率曲線は、観測された偽マッチ率曲線が急峻に上昇するポイントに近いところで上向きに曲がっており、これは校正法の特に心強い特徴である。真の確率に近い予測確率と、真の値を含む有益な区間推定値を提供するという点で、この校正法は優れている。これと比較すると、経験的な校正を行わずに重みを乗じることによって構築されたマッチ確率の元の推定値は、非常に不正確であった。

1.8 Some useful results from probability theory

読者は確率や確率分布に関する初歩的な操作に慣れているものとする。特に、本書の重要な部分を理解するために必要な基本的な確率の背景には、共同密度の操作、単純モーメントの定義、変数の変換、シミュレーションの方法などが含まれる。このセクションでは、これらの前提条件を簡単に復習し、本書の残りの部分で使用される表記上の慣例を明確にする。付録Aは、一般的に使用される確率分布に関する情報である。

Means and variances of conditional distributions

ある変数\(u\)の平均と分散が、\(v\)の条件付き平均と分散によって与えられることがしばしばある。$$E(u) = E(E(u\mid v)) \tag{1.8}\label{eq1.8}$$この式は$$\begin{eqnarray} E(u) & = & \int\!\!\!\int{up(u, v)dudv}\\ & = & \int\!\!\!\int{up(u\mid v)du p(v)dv}\\ & = & \int{E(u\mid v)p(v)dv}\end{eqnarray}$$と証明される。分散については、以下の公式が成り立つ。$$var(u) = E(var(u\mid v)) + var(E(u\mid v)) \tag{1.9}\label{eq1.9}$$これを見ると、条件付き分散の平均の和と条件付き平均の分散の和が、\(u\)の分散になっている。証明は以下のようになる。$$\begin{eqnarray}E(var(u\mid v)) + var(E(u\mid v)) & = & E(E(u^2\mid v)-(E(u\mid v))^2) + E((E(u\mid v))^2) – (E(E(u\mid v)))^2 \\ & = & E(u^2)-E((E(u\mid v))^2) + E((E(u\mid v))^2)-(E(u))^2 \\ & = & E(u^2)-(E(u))^2\\ & = & var(u)\end{eqnarray}$$である。

Transformation of variables

計算のための変数変換はしばしば行われる。確率変数\(p_u(u)\)を\(v = f(u)\)によって変換する。ただし、\(u, v\)のベクトルとしての次元は同一とする。\(p_u(u)\)が離散分布で、\(f\)が全単射の場合、\(p_v(v)\)は以下の式で与えられる。$$p_v(v) = p_u(f^{-1}(v))$$これは自明である。\(f\)が全単射でない場合は、区間を区切り全単射な関数として、\(p_v(v)\)を表せば良い。

\(p_u(u)\)が連続分布で、\(f\)が全単射の場合、$$p_v(v) = \mid J\mid p_u(f^{-1}(v))$$となる。\(J\)は変換\(u = f^{-1}(v)\)に関するヤコビアンである。

よく用いられる変換は、\(u \in (0, 1)\)に対して$$\text{logit}(u) = \log{\left(\frac{u}{1-u}\right)} \tag{1.10}\label{eq1.10}$$である。これは、$$\begin{eqnarray}v & = & \log{\left(\frac{u}{1-u}\right)}\\ e^{v} & = & \frac{u}{1-u}\\ u & = & \frac{e^v}{1+e^{v}}\end{eqnarray}$$であるから、\(\displaystyle \text{logit}^{-1}(v) = \frac{e^v}{1+e^v}\)となる。その他probit関数（正規分布の累計分布関数の逆関数）もよく用いられる。

1.9 Computation and Software

Sampling using the inverse cumulative distribution function

ここでは例として逆累計分布関数を用いた乱数サンプリングについて述べる。一次元の確率密度関数\(p(v)\)に対して、累計分布関数(Cumulative distribution function: CDF)\(F(v*)\)を以下のように定義する。$$\begin{eqnarray}F(v*) & = & Pr(v\leq v*) \\ & = & \begin{cases}\displaystyle \sum_{v\leq v*}{p(v)}\ \ if\ p \ is \ discrete \\ \int_{-\infty}^{v*}{p(v)dv}\ \ if\ p\ is\ continuous.\end{cases}\end{eqnarray}$$

まず、\([0, 1]\)の一様分布から乱数表やコンピュータの乱数機能を用いて、任意の数\(U\)を一つ取り出す。そして、\(v = F^{-1}(U)\)とする。\(F\)は必ずしも全単射でなくても良い。実際、離散分布の場合は累積分布関数は全単射にならない。しかし、\(F^{-1}(U)\)は累積分布関数の全区間での積分が\(1\)になるという条件のもと、一意に定まる。\(v\)は\(p\)からランダムに引き出されており、\(F^{-1}(U)\)が単純な形であれば簡単に計算できる。離散分布では、\(F^{-1}(U)\)はテーブルで表現できる。

例えば指数分布\(\displaystyle f(x: \lambda) = \lambda e^{-\lambda x}\)のような連続分布の場合、累積分布関数は\(F(x: \lambda) = 1-e^{-\lambda x}\)と与えられる（積分すればすぐに分かる）。なので、\(U = F(v)\)となる\(v\)は\(\displaystyle v = -\frac{\log{(1-U)}}{\lambda}\)である。\(1-U\)が一様分布\([0, 1]\)に従うことを思い出すと、指数分布からのランダムサンプリングは\(\displaystyle -\frac{\log{U}}{\lambda}\)で得ることができる。

逆関数法 - Wikipedia

ja.wikipedia.org

何をしているのかというと、累積分布関数が\(F(v*)\)で表されるような乱数を生成したいという状況で、一様分布のサンプリングが可能であれば、そこに\(F^{-1}\)を咬ませることで目的が達成されるということである。指数分布に従うランダムサンプリングの例として、Python3でのcodeを以下に示す。

import numpy as np
from scipy import stats
from matplotlib import pylab as plt
import seaborn as sns
l = 0.5
U = stats.uniform.rvs(size = 1000)
r = -np.log(U)/l
sns.set_palette("Blues")
plt.title("Inveree transform method")
sns.histplot(r)

Chapter 2 Single-parameter models

このChapterでは\(\theta\)は一次元で、よく知られた\(4\)つの分布、すなわち二項分布、正規分布、ポアソン分布そして指数分布に従うものとして扱う。

2.1 Estimating a probability from binomial data

ベルヌーイ試行とは、取りうる結果が成功か失敗かのいずれかで、各試行において成功の確率が同じであるというランダム施行である。二項分布は\(n\)回の交換可能な試行の列から生じるデータ、あるいは各試行が成功と失敗とラベル付けされた\(2\)つの可能な結果のいずれかをもたらす大きな母集団から得られるデータの、自然なモデルを提供する。交換可能なので、モデルは\(n\)回の施行および\(y\)回の成功によって要約される。\(\theta\)を集団全体における成功の確率、あるいは各試行における成功の確率とすると、二項分布モデルは$$p(y\mid \theta) = \text{Bin}(y\mid n, \theta) = \begin{pmatrix}n\\ y\end{pmatrix}{\theta}^{y}{(1-\theta)}^{n-y} \tag{2.1}\label{eq2.1}$$と表現することができる。ここで左辺では\(n\)が式に含まれていないが、試行回数は実験計画の一部として固定されているとみなすからである。ここで論ぜられる確率については、\(n\)に条件付けられると考えられる。

Example. Estimating the probability of a female birth

\eqref{eq2.1}では\(n\)回の施行は、\(\theta\)の値と独立しているものと考えている。つまり、同一家系での多胎産などの場合は、前回の出産の性別が次回の出産の性別に影響を与える可能性が排除できず、試行が交換可能という仮定は成り立たない。Figure2.1は\(\theta = 0.6\)と一定に保ち、試行回数\(n\)と\(y\)を変化させたグラフである。

これをみると、\(\theta\)の事後分布の信頼区間は\(n\)が増えると狭くなっていることが分かる。これもPythonで簡単に再現できる。ただし、下のcodeのpermは二項係数である。

n = 100
y = 60
t = np.linspace(0, 1, 1000)
p = perm(n, y)*pow(t, y)*pow(1-t, n-y)
sns.relplot(t, p)

ベイズ解析を行うために、\(\theta\)の事前分布を指定しないといけない。ここでは単純に一様分布\([0, 1]\)を使う。\eqref{eq1.2}を利用し、\eqref{eq2.1}に当てはめると、\(\theta\)の事後分布は$$p(\theta\mid y) \propto {\theta}^y {(1-\theta)}^{n-y} \tag{2.2}\label{eq2.2}$$と表される。\(n, y\)を固定すると、二項係数\(\displaystyle \begin{pmatrix}n \\ y\end{pmatrix}\)は\(\theta\)に依存せず、\(\theta\)の事後分布を計算する際定数のように扱うことができる。\eqref{eq2.2}は正規化されていない形のベータ関数を用いて、以下のように表すことができる。$$\theta \mid y \sim \text{Beta}(y+1, n-y+1) \tag{2.3}\label{eq2.3}$$なぜベータ分布を用いるのか、一つの理由としては二項分布のパラメータ\(\theta\)と同様に、ベータ分布は\(0\)と\(1\)の区間に制限されていることがあげられる。二つ目に、ベータ関数はパラメータの組み合わせによって一様分布、正規分布、またはU字型分布など様々な形状をとることができることをあげる。三つ目に、ベータ分布が二項分布の共役事前分布(conjugate prior)であることが大きな理由としてあげられる。尤度の共役事前分布とは、与えられた尤度とともに用いたとき、事前分布と同じ関数型を有する事後分布を返す事前分布のことである。簡単に言うと、事前分布としてベータ分布を用い、尤度として二項分布を用いたとき、事後分布として常にベータ分布を得る。共役事前分布についてはその他の組み合わせも存在する。例えば、正規分布はそれ自体の共役事前分布である。

長年に渡って計算上の理由から、ベイジアンは共役事前分布を用いることに制限されてきた。つまり、ベイズ統計学における一般的な問題の解決が、事後分布を解析的に導くことに集約されてきた。事後分布を解析的に導くことができない場合、分析がそこで終了してしまう。MCMCなど適切な計算方法が発達したことで、任意の事後分布に対して解析的なアプローチが可能になったのである。

Historical note: Bayes and Laplace

ここでの計算で$$\begin{eqnarray}p(y) & = & \int_{0}^{1}{\begin{pmatrix}n \\ y\end{pmatrix} {\theta}^y {(1-\theta)}^{n-y}d\theta} \tag{2.5}\label{eq2.5}\\ & = & \frac{1}{n+1}\ \ \ for\ \ \ y = 0, \cdots, n\end{eqnarray}$$とある。用いているのは、ベータ関数の公式$$\int_{0}^{1}{x^m(1-x)^ndx} = \frac{m!n!}{(m+n+1)!}$$である。

Prediction

二項分布の例では、一様分布に従う\(\theta\)に対して、事後分布は\eqref{eq2.5}のように明示的な形であることができた。このモデル下では、可能な\(y\)の取り得る値はすべて等しく、事前に決まっている。事後分布の予測について、より興味があるのは別の\(n\)回の試行よりは、\(1\)回の新しい試行であろう。\(\tilde{y}\)を新しい\(1\)回の結果とすると、これは最初の\(n\)回の試行と交換可能であるから、$$\begin{eqnarray}\text{Pr}(\hat{y} = 1\mid y) & = & \int_{0}^{1}{\text{Pr}(\tilde{y} = 1\mid \theta, y)p(\theta\mid y)d\theta}\\ & = & \int_{0}^{1}{\theta p(\theta\mid y)d\theta}\\ & = & E(\theta\mid y)\\ & = & \frac{y+1}{n+2}\end{eqnarray}$$となる。以下ではこの証明を記す。まず、\(y\)は観測されたデータ、\(\tilde{y}\)はこれから予測されるデータであることに留意する。また、\(\text{Pr}\)と書いた場合、未観測のデータに対する予測確率を表しているのだと考えられる。$$\text{Pr}(\tilde{y} = 1\mid y) = \int_{0}^{1}{\text{Pr}(\tilde{y} = 1 \mid \theta, y)p(\theta \mid y)d\theta}$$であるが、これは\eqref{eq1.4}を用いている。次の\(\text{Pr}(\tilde{y} = 1\mid \theta, y) = \theta\)であるが、\(\theta\)は\(y\)（観測された成功回数）に依らず、新たな\(\tilde{y}\)が成功するかどうかはまさに\(\theta\)なのであるから、成り立っている。最後の\(E(\theta\mid y)\)であるが、\eqref{eq2.3}から\(\theta \mid y\)はベータ分布で表され、ベータ分布について既知の知識\(X\sim \text{Beta}(\alpha, \beta)\)のとき\(\displaystyle E(X) = \frac{\alpha}{\alpha+\beta}\)から、$$\begin{eqnarray}\theta\mid y & \sim & \text{Beta}(y+1, n-y+1)\\ E(\theta\mid y) & = & \frac{y+1}{y+1+n-y+1}\\ & = & \frac{y+1}{n+2} \end{eqnarray}$$となる。この結果はLapalceのrule of succession（連続の法則）として知られている。

Rule of succession - Wikipedia

en.wikipedia.org

例えば、観測された成功が\(0\)回だった場合、次の試行で成功する確率は\(\displaystyle \frac{1}{n+2}\)になるし、観測された成功が\(n\)回だった場合、次の試行で成功する確率は\(\displaystyle \frac{n+1}{n+2}\)になる。

2.2 Posterior as compromise between data and prior information

\eqref{eq1.8}, \eqref{eq1.9}の重要な応用が以下になる。$$E(\theta) = E(E(\theta \mid y)) \tag{2.7}\label{eq2.7}$$および$$\text{var}(\theta) = E(\text{var}(\theta \mid y)) + \text{var}(E(\theta \mid y)) \tag{2.8}\label{eq2.8}$$

\eqref{eq2.7}は\(\theta\)の事前平均は、可能なデータの分布全体にわたるすべての可能な事後平均の平均である。\eqref{eq2.8}はより興味深く、事後分散は可能なデータの分布に対する事後平均の変動に依存する量だけ、平均して事前分散よりも小さくなる。

事前分布が一様分布に従う例では事前平均は\(\displaystyle \frac{1}{2}\)であり、事前分散は\(\displaystyle \frac{1}{12}\)である。これは、$$\begin{eqnarray}E(X) & = & \int_{0}^{1}{xdx} \\ & = & \frac{1}{2}\\ var(X) & = & E(X^2) -E(X)^2 \\ & = & \int_{0}^{1}{x^2dx}-\left(\frac{1}{2}\right)^2\\ & = & \frac{1}{3}-\frac{1}{4}\\ & = & \frac{1}{12}\end{eqnarray}$$から分かる。事後平均\(\displaystyle \frac{y+1}{n+2}\)は事前平均\(\displaystyle \frac{1}{2}\)とサンプル比\(\displaystyle \frac{y}{n}\)の間の妥協点となる。すなわち、\(\displaystyle \frac{1}{2}\)と\(\displaystyle \frac{y}{n}\)の間の点になる（差を作ると、\(\displaystyle \frac{1}{2}-\frac{y+1}{n+2} = \frac{n-2y}{2(n+2)}\)であり、\(\displaystyle \frac{y}{n}-\frac{y+1}{n+2} = \frac{2y-n}{n(n+2)}\)となる）。またデータのサンプルサイズが大きくなるにつれ事前平均の果たす役割が小さくなることも分かる。これはベイズ推定の一般的な特徴で、事後分布は事前情報とデータの間の妥協点に集中し、サンプルサイズが大きくなるにつれ妥協点はデータによって制御される。

2.3 Summarizing posterior inference

事後分布はパラメータ\(\theta\)に関して現時点で必要な情報のすべてを持っている。事後確率分布\(p(\theta\mid y)\)をグラフィカルな形で報告するのが望ましいが、その他平均mean、中央値median、最頻値mode、あるいは四分位など必要に応じて必要な値で結果を報告するのが良い。例えば、\eqref{eq2.3}のベータ分布では、最頻値modeは\(\displaystyle \frac{y}{n}\)は\(\theta\)の不偏推定量を与えることが知られている。

2.4 Informative prior distributions

二項分布の例では、\(\theta\)の事前分布として一様分布を扱った。以降は一般的にどのように事前分布を指定するかを扱う。事前分布に与えることができる2つの基本的な解釈を検討する。母集団の解釈(Population interpretation)は事前分布はパラメーターが取り得る値の母集団を表し、そこから現在関心のある\(\theta\)が引き出されている。より主観的な知識解釈(knowledge interpretation)は\(\theta\)に関する知識および不確実性をその値が事前分布からのランダムな実現と見なすことができるかのように表現する必要がある。新しい産業プロセスでの失敗の確率推定など、多くの問題について、仮想的な場合を除いて現在の\(\theta\)が導き出された\(\theta\)に完全に一致するという母集団はない。通常事前分布にはすべてのもっともらしい\(\theta\)の値を含める必要があるが、データに含まれる\(\theta\)に関する情報が妥当な事前確率の仕様を上回ることが多く、事前分布を実際に真の値の周りに集中させる必要はない。

Binomial example with different prior distributions

\(\theta\)の関数として、\eqref{eq2.1}の尤度が以下の様に表されるとする。$$p(y\mid \theta) \propto {\theta}^a(1-\theta)^b $$この場合、事前分布が同様の形状であるならば、事後分布もまた同様の形状になる。そのような事前分布を以下のように表しておく。$$p(\theta) \propto {\theta}^{\alpha-1}(1-\theta)^{\beta-1}$$これはベータ分布\(\text{Beta}(\alpha, \beta)\)であり、\(\theta \propto \text{Beta}(\alpha, \beta)\)と表現できる。\(p(\theta)\)と\(p(y\mid \theta)\)を見比べると、事前分布は\(\alpha-1\)の成功、\(\beta-1\)の失敗を表していると考えられる。このようなパラメータはハイパーパラメータhyperparameterと言われる。この\(\alpha, \beta\)についても平均、分散などをもった事前分布を仮定することが可能である。ここでは\(\alpha, \beta\)が定まったとして計算を進めると、$$\begin{eqnarray}p(\theta \mid y) & \propto & {\theta}^y(1-\theta)^{n-y}{\theta}^{\alpha-1}{(1-\theta)}^{\beta-1}\\ & = & {\theta}^{y+\alpha-1}(1-\theta)^{n-y+\beta-1}\\ & = & \text{Beta}(\theta \mid \alpha+n, \beta + n-y)\end{eqnarray}$$となる。このように事後分布が事前分布と同様のパラメータ形式を用いて記述されることを、共役conjugateという。ベータ分布は事前分布が二項分布であるときの、共役分布である。一般に共役分布は既知のパラメータ形式を取るという点で数学的に簡便である。もちろん、現実あるいは臨床的な事前の知識が共役分布と矛盾する場合は、無理やり共役分布を用いてはならない。

\(\theta\)の事前分布にベータ分布を用いた場合、将来の成功の事後確率は$$E(\theta \mid y ) = \frac{\alpha+y}{\alpha+\beta+n}$$で表される。この期待値は常にサンプル比である\(\displaystyle \frac{y}{n}\)と事前分布の平均である\(\displaystyle \frac{\alpha}{\alpha+\beta}\)の間にある。このことは、$$\begin{eqnarray}\frac{y}{n}-\frac{\alpha+y}{\alpha+\beta+n} & = & \frac{y(\alpha+\beta+n)-n(\alpha+y)}{n(\alpha+\beta+n)}\\ & = & \frac{y(\alpha+\beta)-n\alpha}{n(\alpha+\beta+n)}\\ \frac{\alpha+y}{\alpha+\beta+n}- \frac{\alpha}{\alpha+\beta} & = & \frac{(\alpha+\beta)(\alpha+y)-\alpha(\alpha+\beta+n)}{(\alpha+\beta)(\alpha+\beta+n)}\\ & = & \frac{y(\alpha + \beta)-n\alpha}{(\alpha+\beta)(\alpha+\beta+n)}\end{eqnarray}$$であることから分かる。事後分布の分散は$$\begin{eqnarray}\text{var}(\theta \mid y) & = & \frac{(\alpha + y)(\beta + n-y)}{(\alpha+\beta+n)^2 (\alpha+\beta+n+1)}\\ & = & \frac{E(\theta \mid y)(1-E(\theta\mid y))}{\alpha+\beta+n+1}\end{eqnarray}$$である。固定された\(\alpha, \beta\)に対して\(n, n-y\)が大きくなると、\(\displaystyle E(\theta\mid y) \approx \frac{y}{n}\)となり、上の分散の式は\(\displaystyle \text{var}(\theta\mid y) \approx \frac{1}{n}\frac{y}{n}\left(1-\frac{y}{n}\right)\)となる。\(n\)が大きいときこれは\(0\)に近づく。つまり、サンプル数が大きくなると事前分布のパラメータは事後分布に対してほとんど影響を与えないことが分かる。後に出てくるが、中心極限定理はベイズ統計学の文脈では以下のように表現される。$$\left(\frac{\theta-E(\theta\mid y)}{\sqrt{\text{var}(\theta\mid y)}}\middle| y\right) \to \text{N}(1, 0)$$と表すことができる。この結果は、しばしば\(n\)（サンプル数）が大きいとき、事後分布を正規分布で代替させることの正当化として用いられる。二項分布のパラメータ\(\theta\)を考える場合、\(\theta\)をlogit変換した方が二項分布を事後分布に用いるときはより正確である。つまり、\(\theta\)の推定を行うより\(\displaystyle \log{\left(\frac{\theta}{1-\theta}\right)}\)の推定を行った方がよい。この変換で\(\theta\)の区間は\([0, 1]\)から\((-\infty, \infty)\)になり、正規分布での推定によりフィットする。

Conjugate prior distributions

共役は正式には以下のように定義される。\(\mathcal{F}\)が尤度（サンプリング分布）\(p(y\mid \theta)\)の類であり、\(\mathcal{P}\)が事前分布\(\theta\)の類であるとする。以下が成り立つとき\(\mathcal{P}\)は\(\mathcal{F}\)と共役である。$$p(\theta \mid y) \in \mathcal{P}\ \ \text{for all}\ \ p(\cdot\mid \theta) \in \mathcal{F}\ \ \text{and} \ \ p(\cdot)\in \mathcal{P}$$すべての分布類として\(\mathcal{P}\)を選択した場合、どのような類のサンプリング分布を使用しても\(\mathcal{P}\)は常に共役となるので、この定義は形式的にあいまいであるといえる。我々が興味あるのは、自然共役分布であり、これは\(\mathcal{P}\)を尤度と同じ関数形式を持つすべての密度の集合と見なすことによって得られる。自然共役の利用は計算上の利点に加え、実際上の利点も大きいが、そのことは後の章でも述べられる。

Conjugate prior distributions, exponential families, and sufficient statistics

指数型分布族と十分統計量に触れてこのセクションを終える。難解なので、飛ばしても良い。類\(\mathcal{F}\)が指数型分布族に含まれるとき、以下の形式を有する。$$p(y_i \mid \theta) = f(y_i) g(\theta)e^{{\phi(\theta)}^T u(y_i)}$$要素\(\phi(\theta), u(y_i)\)は一般的に、\(\theta\)と同様の次元を持つベクトルである。ベクトル\(\phi(\theta)\)は\(\mathcal{F}\)の自然パラメータnatural parameterと呼ばれる。各々が独立で同一な分布に従う観測値\(y = (y_1, y_2, \cdots, y_n)\)に対応する尤度は以下で与えられる。$$p(y \mid \theta) = \left(\prod_{i=1}^{n}{f(y_i)}\right)g(\theta)^{n}\exp{\left(\phi(\theta)^{T}\sum_{i=1}^{n}{u(y_i)}\right)}$$これは上の式をかけ合わせたものである。より簡便に\(\theta\)の関数として書くと、下のような形になる。$$p(y \mid \theta) \propto g(\theta)^{n}e^{\phi(\theta)^T}t(y),\ \ \text{where} \ \ t(y) = \sum_{i=1}^{n}{u(y_i)}$$\(\theta\)の尤度は\(t(y)\)の値を通じてのみデータ\(y\)に依存するので、\(t(y)\)は\(\theta\)の十分統計量sufficient statisticと呼ばれる。事前分布が$$p(\theta)\propto g(\theta)^{\eta}e^{\phi(\theta)^{T}\nu}$$の形であるとき、事後分布は$$p(\theta\mid y)\propto g(\theta)^{\eta + n}e^{\phi(\theta)^T(\nu + t(y))}$$と表すことができる。

一般に、指数方分布族は自然共役事前分布を持つ唯一の族である。それは、特定の不規則な場合を除き、指数型分布のみがすべての\(n\)に対して十分統計量が固定されている唯一の分布であるからである。

Example. Probability of a girl birth given placenta previa

男女比に影響を及ぼす可能性のある要因の具体例として、前置胎盤を考えてみる（前置胎盤とは、胎盤が子宮の低い位置に着床し、胎児が正常な経膣分娩の妨げとなる妊娠の異常な状態である）。ドイツで行われた前置胎盤の性別に関する初期の研究によると、合計980件の出産のうち、437件が女性であった（\(\displaystyle \frac{437}{980} = 0.4459\cdots\)）。このことは、前置胎盤の集団における女性出産の割合が、一般集団における女性出産の割合である0.485未満であるという主張の根拠となるであろうか？

Analysis using a uniform prior distribution.

女児出生確率が一様であるという仮定の下では、事後分布は\(\text{Beta}(438, 544)\)である。\(\theta\)の事後平均は\(0.446\)で、事後標準偏差は\(0.016\)である。正確な事後分位数は、ベータ密度の数値積分を用いて求めることができ、実際にはコンピュータよって計算する。中央値は\(0.446\)で、中央の95%事後区間は\([0.415, 0.477]\)になる。この\(95\%\)事後区間は，計算された事後平均と標準偏差を用いて正規近似を用いると得られる区間と，小数点以下3桁まで一致する。事後分布の近似正規性についてのさらなる議論は、第4章でされる。

多くの状況では、事後密度関数を直接計算することは不可能である。そのような場合、推論を得るために事後分布からのシミュレーションを使用することが有用である。図2.3の最初のヒストグラムは、\(\text{Beta}(438, 544)\)事後分布からの1000回のDrawの分布を示す。1000個の順序付きドローの25番目と976番目のドローを取ることによって得られる\(95\%\)事後区間の推定値は\([0.415, 0.476]\)であり、事後分布からの1000個のドローの中央値は\(0.446\)である。1000回の抽選の標本平均と標準偏差は\(0.445\)と\(0.016\)で、正確な結果とほぼ同じになる。\(95\%\)事後区間の正規近似は、\([0.445 \pm 1.96 – 0.016] = [0.414, 0.476]\)となる。標本が大きいことと、\(\theta\)の分布が\(0\)と\(1\)から離れて集中しているという事実のために，正規近似はこの例ではよく機能している。

すでに述べたように、比率を推定する場合、正規近似は一般にロジット変換、つまり\(\displaystyle \log{\left(\frac{\theta}{1-\theta}\right)}\)を適用することによって改善する。これによって、パラメータ空間を単位区間から実線に変換する。図2.3 の2番目のヒストグラムは、変換されたDrawの分布を示している。1000回のDrawに基づくロジット尺度の推定事後平均と標準偏差はそれぞれ\(-0.220\)と\(0.065\)である。\(\theta\)の\(95\%\)事後区間に対する正規近似は、ロジット尺度上の\(95\%\)区間\([-0.220 \pm 1.96 – 0.065]\)を反転することによって得られ、これは元の尺度上では\([0.414, 0.477]\)となる。ロジット尺度の使用による改善は、標本サイズが小さいか、あるいは\(\theta\)の分布が\(0\)または\(1\)に近い値を含むときに最も顕著である。どのような実際のデータ分析でも、適用される文脈を念頭に置くことが重要である。この例で関心のあるパラメータは、伝統的に出生数の男女比（性比）、つまり\(\displaystyle \frac{1-\theta}{\theta}\)として表現される、男女比の事後分布を3番目のヒストグラムに示す。性比の事後中央値は\(1.24\)で、\(95\%\)事後区間は\([1.10, 1.41]\)である。この事後分布は、通常のヨーロッパ人種の性比である\(1.06\)をはるかに上回る値に集中しており、前置胎盤がある場合に女児が出生する確率が一般集団よりも低いことを示唆している。

Analysis using different conjugate prior distributions.

提案された事前分布に対する\(\theta\)に関する事後推論の感度分析をTable 2.1に示す。表の最初の行は一様事前分布\(\alpha = 1, \beta = 1\)、それ以降の行は一般人口における女性の出生比率である\(0.485\)付近に次第に集中する事前分布を用いている。最初の列は\(\theta\)の事前平均を示し、2番目の列は\(\alpha+\beta\)によって測定される事前情報の量を指数化する。大きな標本に基づく事後推論は、事前分布に特に敏感ではない。事前分布が100または200の出生に相当する情報を含む表の下部でのみ、事後分布が事前分布に顕著に引っ張られる。

Analysis using a nonconjugate prior distribution.

この問題に対する共役ベータ族に代わるものとして、\(0.485\)を中心とするが、真の値が遠くにある可能性を考慮するため、中心から遠く離れて平坦な事前分布を望むかもしれない。図2.4aの区分線形事前分布密度は、この形式の事前分布の例である。この事前分布は平均\(0.493\)、標準偏差\(0.21\)で、\(\alpha+\beta = 5\)のベータ分布の標準偏差に似ている。非正規化事後分布は、\(\theta\)値の格子、\((0.000,0.001,…,1.000)\)で、各点で事前密度と二項尤度を掛け合わせることによって得られる。事後シミュレーションは、\(\theta\)値の離散格子上の分布を正規化することによって得ることができる。図2.4bは、離散的な事後分布からの1000回のDrawを行った後のヒストグラムである。事後分布の中央値は\(0.448\)で、\(95\%\)中央事後区間は\([0.419, 0.480]\)である。事前分布はデータに制圧されているので、これらの結果はベータ分布に基づく表2.1の結果と一致する。このようなグリッド・アプローチをとる場合、粗すぎて事後質量のかなりの部分を歪めてしまうようなグリッドを避けることが重要である。

2.5 Normal distribution with known variance

正規分布は，ほとんどの統計モデリングの基本である。中心極限定理は，あまり分析的に便利でない実際の尤度の近似として，多くの統計問題で正規尤度を使うことを正当化するのに役立つ。また、後の章で見るように、正規分布自体が良いモデル適合を提供しない場合でも、t分布や有限混合分布を含むより複雑なモデルの構成要素として有用である。今のところは、正規モデルが適切であると仮定して、ベイズの結果を単純に作業する。まず1つのデータ点についての結果を導き、次に多数のデータ点がある一般的なケースについての結果を導く。

Likelihood of one data point

単純な例として、データ\(y\)の尤度を以下のように定義する。$$p(y \mid \theta) = \frac{1}{\sqrt{2\pi\sigma}}e^{-\frac{1}{2{\sigma}^2}(y-\theta)^2}$$つまり、\(\theta\)は平均で、\(\sigma\)は我々の既知である分散とする。

Conjugate prior and posterior distributions

\(\theta\)の関数として考えると、上の尤度は二次形式の指数となるので、共役事前分布は以下の形式になろう。$$p(\theta) = e^{A{\theta}^2+B\theta + C}$$これを、以下の族としてまとめる。$$p(\theta) \propto \exp{\left(-\frac{1}{2{{\tau}_0}^2}(\theta-\mu_{0})^2\right)}$$すなわち、\(\theta \sim N(\mu_{0}, {\tau_{0}}^2)\)であり、\({\mu}_0, {{\tau}_0}^2\)はハイパーパラメータhyperparameterである。通常ハイパーパラメータは既知であるものと考える。事後分布については、\(\theta\)以外の変数は一定と考えるので、以下の式で与えられる。$$p(\theta \mid y) \propto \exp{\left(-\frac{1}{2}\left(\frac{(y-\theta)^2}{{\sigma}^2} + \frac{(\theta-{\mu}_0)^2}{{{\tau}_0}^2}\right)\right)}$$指数の中身を展開する。$$\begin{eqnarray}\frac{(y-\theta)^2}{{\sigma}^2} + \frac{(\theta-{\mu}_0)^2}{{{\tau}_0}^2} & = & \left(\frac{1}{{\sigma}^2} + \frac{1}{{{\tau}_0}^2}\right){\theta}^2-2\theta\left(\frac{y}{{\sigma}^2} + \frac{{\mu}_0}{{{\tau}_0}^2}\right)+\frac{y^2}{{\sigma}^2}+\frac{{{\mu}_0}^2}{{{\tau}_0}^2}\end{eqnarray}$$である。$$\frac{1}{{{\tau}_1}^2} = \frac{1}{{\sigma}^2} + \frac{1}{{{\tau}_0}^2}, \mu_1 = {{\tau}_1}^2\left(\frac{{\mu}_0}{{{\tau}_0}^2}+\frac{y}{{\sigma}^2}\right) \tag{2.10}\label{eq2.10}$$とすると、これは$$p(\theta \mid y) \propto \exp{\left(-\frac{1}{2{{\tau}_1}^2}(\theta-{\mu}_1)^2\right)} \tag{2.9}\label{eq2.9}$$と書ける。

Precisions of the prior and posterior distributions

正規分布を扱うときに、分散の逆数は重要な役割を果たし、精度precisionと呼ばれる。\eqref{eq2.10}が示すのは、それぞれ精度が予めわかっている正規分布に従う事前分布とデータにおいて、事後分布の精度は事前分布の精度とデータの精度の和になるということである。事後分布の平均\(\mu_1\)の解釈については色々な方法がある。\eqref{eq2.10}において、事後分布の平均は事前分布の平均と観測されたデータ\(y\)の、精度を加味した加重平均になっている。または、次のように観測値\(y\)で補正した表現も可能である。$${\mu}_1 = {\mu}_0 + (y-{\mu}_0)\frac{{{\tau}_0}^2}{{\sigma}^2+{{\tau}_0}^2}$$あるいは、データが事後分布の平均に縮尺していると捉えると、以下のようになる。$${\mu}_1 = y-(y-{\mu}_0)\frac{{\sigma}^2}{{\sigma}^2+{{\tau}_0}^2}$$それぞれの表現は、事後分布の平均が事前分布の平均と観測されたデータの中間点に存在することを示している。極端な例としては、\(y = \mu_0, {{\tau}_0}^2 = 0\)のとき\({\mu}_1 = \mu_0\)となり、\(y = {\mu}_0, {\sigma}^2 = 0\)のとき\(\mu_1 = y\)となる。

\({{\tau}_0}^2 = 0\)のときは、事前分布の精度はデータよりも正確であるということであり、事後分布と事前分布は同一となり、どちらも\({\mu}_0\)に集中する。\({\sigma}^2 = 0\)のときは、データは完全に正確で、事後分布は観察された値に集中する。\(y = {\mu}_0\)のときは、事前分布とデータの平均は一致し、事後分布の平均もこの値と同一になる。

Posterior predictive distribution

将来の観測値\(\tilde{y}\)に対しての予測分布は\(p(\tilde{y}\mid y)\)は、\eqref{eq1.4}を用いると直接計算することができる。すなわち、$$\begin{eqnarray}p(\tilde{y}\mid y) & = & \int{p(\tilde{y}\mid \theta)p(\theta\mid y)d\theta}\\ & \propto & \int{\exp{\left(-\frac{1}{2{\sigma}^2}(\tilde{y}-\theta)^2\right)}\exp{\left(-\frac{1}{2{{\tau}_1}^2}(\theta-{\mu}_1)^2\right)}d\theta}\end{eqnarray}$$一行目の変形は、将来の予測\(\tilde{y}\)と与えられた\(\theta\)の値は観測値\(y\)に依存せず、成立する。二変量正規分布の特性を使用して、\(\tilde{y}\)の分布をより簡単に決定できる。被積分関数の積は、\((\tilde{y}, \theta)\)の\(2\)次関数の指数である。したがって、\(\tilde{y}\)と\(\theta\)は結合した正規事後分布を持ち、\(\tilde{y}\)の周辺事後分布は正規分布となる。

事後の予測分布は、\(E(\tilde{y}\mid \theta) =\theta, \text{var}(\tilde{y}\mid \theta) = {\sigma}^2\)および\eqref{eq2.7}, \eqref{eq2.8}を用いて以下のように決定できる。$$E(\tilde{y}\mid y) = E(E(\tilde{y}\mid \theta, y)\mid y) = E(\theta\mid y) = \mu_1$$および$$\begin{eqnarray}\text{var}(\tilde{y}\mid y) & = & E(\text{var}(\tilde{y}\mid \theta, y)\mid y) + \text{var}(E(\tilde{y}\mid \theta, y)\mid y) \\ & = & E({\sigma}^2\mid y) +\text{var}(\theta\mid y)\\ & = & {\sigma}^2 + {{\tau}_1}^2 \end{eqnarray}$$したがって、\(\tilde{y}\)の事後予測分布は\(\theta\)の事後分布と同一の平均をもち、分散は二つの要素に分解される。一つはモデルの予測分散の\({\sigma}^2\)であり、もう一つは\(\theta\)が不確定であるため生じる事後分散\({\tau}_1\)である。

Normal model with multiple observations

単一の観測モデルは、各々の観測が独立で同一な分布に従う\(y = (y_1, y_2, \cdots, y_n)\)であるという仮定のもと、以下のように拡張される。正式に記述すると、$$\begin{eqnarray}p(\theta \mid y) & \propto & p(\theta)p(y\mid \theta)\\ & = & p(\theta)\prod_{i=1}^{n}{p(y_i\mid \theta)}\\ & \propto & \exp{\left(-\frac{1}{2{{\tau}_0}^2}(\theta-{\mu}_0)^2\right)}\prod_{i=1}^{n}{\exp{\left(-\frac{1}{2{\sigma}^2}(y_i-\theta)^2\right)}} \\ & \propto & \exp{\left(-\frac{1}{2}\left(\frac{1}{{{\tau}_0}^2}(\theta-\mu_0)^2+\frac{1}{{\sigma}^2}\sum_{i = 1}^{n}{(y_i-\theta)^2}\right)\right)}\end{eqnarray}$$である。事後分布はサンプルの平均\(\displaystyle \bar{y} = \frac{1}{n}\sum_{i}{y_i}\)を通じてのみ\(y\)に依存しており、\(\bar{y}\)はこのモデルの十分統計量となる。実際、\(\displaystyle \bar{y}\mid \theta, {\sigma}^2 \sim \text{N}\left(\theta, \frac{{\sigma}^2}{n}\right)\)なので、\(\bar{y}\)を単一の観測のように扱えば、直ちに$$p(\theta\mid y_1, y_2, \cdots, y_n) = p(\theta\mid \bar{y}) = \text{N}(\theta \mid {\mu}_n, {{\tau}_n}^2) \tag{2.11}\label{eq2.11}$$を得る。ここで、$${\mu}_n = \frac{\frac{1}{{{\tau}_0}^2}{\mu}_0+\frac{n}{{\sigma}^2}\bar{y}}{\frac{1}{{{\tau}_0}^2}+\frac{n}{{\sigma}^2}}, \ \ \frac{1}{{{\tau}_n}^2} = \frac{1}{{{\tau}_0}^2} + \frac{n}{{\sigma}^2} \tag{2.12}\label{eq2.12}$$である（証明はそれほど難しくなく、展開して\(\theta\)についてまとめれば良い）。一度に計算せず、順次\(y_1, y_2, \cdots, y_n\)の計算を行っても同様の結果を得る。

事後分布の平均と分散において、事前分布の精度\(\displaystyle \frac{1}{{{\tau}_0}^2}\)とデータの精度\(\displaystyle \frac{n}{{\sigma}^2}\)は同等の役割を果たす。\(n\)が大きいとき、事後分布はサンプル平均\(\bar{y}\)と\({\sigma}^2\)の影響を大きく受ける。例えば、\({{\tau}_0}^2 = {\sigma}^2\)のとき、事前分布は値が\(\mu{}_0\)の\(1\)つの追加観測値と同じ重みを持つ。より詳細に述べると、\(n\)を固定して\({\tau}_0\to \infty\)とするか、あるいは\({{\tau}_0}^2\)を固定して\(n\to \infty\)とすると、$$p(\theta \mid y) \approx \text{N}\left(\theta \mid \bar{y}, \frac{{\sigma}^2}{n}\right) \tag{2.13}\label{eq2.13}$$を得る。これは実用上、事前分布についての信念が、比較的大きい\(\theta\)の範囲で拡散している場合は常に適切な近似値となる。

2.6 Other standard single-parameter models

一般に、事後密度\(p(\theta|y)\)は閉形式を持っていない。正規化定数\(p(y)\)は、積分\eqref{eq1.3}のために計算が特に難しいことが多い。多くの正式なベイズ分析は、閉形式が利用可能な状況に集中している。そのようなモデルは時に非現実的であるが、より現実的なモデルを構築する際に、それらの分析はしばしば有用な出発点を提供する。標準的な分布である二項分布、正規分布、ポアソン分布、指数分布は、単純な確率モデルから自然に派生したものである。すでに説明したように、二項分布は交換可能な結果を数えることから動機づけられ、正規分布は交換可能な・あるいは独立な項の和である確率変数に適用される。我々はまた、すべての正データの対数に正規分布を適用する機会もあり、それは当然多くの独立した乗法因子の積としてモデルされる観察に適用されるであろう。ポアソン分布と指数分布は、すべての時間区間で交換的に発生するとモデル化された事象のカウント数と待ち時間として、それぞれ発生する。一般により複雑な結果に対しては、これらの基本的な分布の組み合わせによって現実的な確率モデルを構築する。例えばセクション22.2では、心理学的実験における分裂病患者の反応時間を、対数スケールの正規分布の二項混合としてモデル化する。これらの標準的なモデルには、それぞれ関連する共役事前分布の系列があり、順番に説明する。

Normal distribution with known mean but unknown variance

\(p(y\mid \theta, {\sigma}^2) = \text{N}(y\mid \theta, {\sigma}^2)\)において、\(\theta\)が既知であり、\({\sigma}^2\)が未知であるとする。このとき、\(n\)個の互いに独立で同一の分布に従う観測データ\(y\)の尤度は以下のように記載される。$$\begin{eqnarray}p(y\mid {\sigma}^2) &\propto& {\sigma}^{-n}\exp{\left(-\frac{1}{2{\sigma}^2}\sum_{i=1}^{n}{(y_i-\theta)^2}\right)}\\ & = & {({\sigma}^2)}^{-\frac{n}{2}}\exp{\left(-\frac{n}{2{\sigma}^2}v\right)}\end{eqnarray}$$ただし、\(v\)は$$v = \frac{1}{n}\sum_{i=1}^{n}{(y_i-\theta)^2}$$で定義される十分統計量である。これに対応する共役事前分布は以下の逆ガンマ分布である。$$p({\sigma}^2)\propto ({\sigma}^2)^{-(\alpha+1)}e^{-\frac{\beta}{{\sigma}^2}}$$ これはhyperparameterとして\((\alpha, \beta)\)を有する。なお逆ガンマ分布とは、$$f_{\alpha, \beta}(x) = \frac{{\beta}^{\alpha}}{\Gamma(\alpha)}x^{-\alpha+1}e^{-\frac{\beta}{x}}\ \ (x > 0)$$で定義される。ちなみに、ガンマ分布は$$f_{\alpha, \beta}(x) = \frac{{\beta}^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x} \ \ (x > 0)$$で定義される。ここで逆ガンマ分布が出てくるのは、平均を固定し分散を推定する場合、逆ガンマ分布が正規分布の共役事前分布となるからである。

パラメータの推定の際に簡便なのは、尺度パラメータ\({\sigma_0}^2\)と自由度\({\nu}_0\)を持つ尺度付き\({\chi}^2\)乗分布を用いることである。これは、$$\text{Inv}{\chi}^2(\theta \mid v, s^2) \sim \frac{\left(\frac{v}{2}\right)^{\frac{v}{2}}}{\Gamma\left(\frac{v}{2}\right)}s^{v}{\theta}^{-\left(\frac{v}{2}+1\right)}\exp{\left(-\frac{vs^2}{2\theta}\right)}$$で表される。つまり、\({\sigma}^2\)の事前分布が\(\displaystyle \frac{{\sigma_0 \nu_0}^2}{X}\)に従うとすればよい。ただし、\(X\sim {{\chi}_{\nu_0}}^2\)である。

ここでは便利であるが、標準的ではない考え方、\({\sigma}^2 \sim \text{Inv}{\chi}^2(\nu_0, {\sigma_0}^2)\)とするものを用いる。この\({\sigma}^2\)の事後分布は$$\begin{eqnarray}p({\sigma}^2 \mid y) & \propto & p({\sigma}^2)p(y\mid {\sigma}^2)\\ & \propto & \left(\frac{{\sigma_0}^2}{{\sigma}^2}\right)^{\frac{{\nu}_0}{2}+1}\exp{\left(-\frac{{\nu}_0{{\sigma_0}^2}}{2{\sigma}^2}\right)}\cdot {({\sigma}^2)}^{-\frac{n}{2}}\exp{\left(-\frac{n}{2}\frac{v}{{\sigma}^2}\right)}\\ & \propto & {({\sigma}^2)}^{-((n+{\nu}_0)/2+1)}\exp{\left(-\frac{1}{2{\sigma}^2}({\nu}_0{\sigma_0}^2+n\nu)\right)} \end{eqnarray}$$となる。したがって、$${\sigma}^2\mid y \sim \text{Inv}-{\chi}^2\left({\nu}_0+n, \frac{{\nu}_0{\sigma_0}^2+n\nu}{{\nu}_0+n}\right)$$である。これは尺度付き逆カイ二乗分布であり、尺度は事前分布とデータの尺度の加重平均であり、自由度は事前分布とデータの自由度の和になっている。事前分布は平均して\({\sigma_0}^2\)の標準偏差を有する自由度\({\nu}_0\)の観測値と等しい情報を提供していると考えることができる。

Poisson model

ポアソン分布はカウントデータを考えるときに自然に生じると言える。例えば、疫学の分野で、疾患の発生数を考えるときなどである。データ\(y\)が発生率\(\theta\)のポアソン分布に従うとき、単一の観測値\(y\)の確率分布は$$p(y\mid \theta) = \frac{{\theta}^{y}e^{-\theta} }{y!}\ \ (y = 0, 1, 2, \cdots)$$となる。独立で同一の分布に従う観測値のベクトル\(y = (y_1, \cdots, y_n)\)に対しては尤度は$$\begin{eqnarray}p(y\mid \theta) & = & \prod_{i=1}^{n}{\frac{1}{y_i!}{\theta}^{y_i}e^{-\theta}}\\ & \propto & {\theta}^{t(y)}e^{-n\theta}\end{eqnarray}$$となる。ただし、\(\displaystyle t(y) = \sum_{i=1}^{n}{y_i}\)であり、\(t(y)\)は十分統計量となる。これは指数型分布族の形式として、$$p(y\mid \theta) \propto e^{-n\theta}e^{t(y)\log{\theta}}$$と書き直せる。自然パラメータは\(\phi(\theta) = \log{\theta}\)であり、自然事前共役分布は$$p(\theta)\propto (e^{-\theta})^{\eta}e^{\nu \log{\theta}}$$となる。\((\eta, \nu)\)はハイパーパラメータとなる。別の言い方をすると、尤度は\({\theta}^a e^{-b\theta}\)の形になり、共役事前分布は\(p(\theta) \propto {\theta}^A e^{-B\theta}\)の形でなければならない。より簡便なパラメータとしては、$$p(\theta) \propto e^{-\beta \theta}{\theta}^{\alpha-1}$$という\((\alpha, \beta)\)をハイパーパラメータとするガンマ分布\(\text{Gamma}(\alpha, \beta)\)の形になる。

\(p(y\mid \theta)\)と\(p(\theta)\)を比較すると、ある意味で事前分布は\(\beta\)の事前観測の中の観測カウント数\(\alpha-1\)と等しいと言える。この共役事前分布を用いると、事後分布は$$\theta \mid y \sim \text{Gamma}(\alpha + b\bar{y}, \beta + n)$$となる。

The negative binomial distribution

共役分布族の場合、周辺分布\(p(y)\)を見つけるのに以下の公式$$p(y) = \frac{p(y\mid \theta)p(\theta)}{p(\theta\mid y)}$$を用いて、よく知られた形の事前分布および事後分布が用いられることがある。例えば、単一の観測値に対する\(y\)に対するポアソン分布の場合、以下の事前予測分布$$\begin{eqnarray}p(y) & = & \frac{\text{Poisson}(y\mid \theta)\text{Gamma}(\theta \mid \alpha, \beta)}{\text{Gamma}(\theta \mid \alpha+y, 1+\beta)} \\ & = & \frac{\Gamma(\alpha+y){\beta}^{\alpha}}{\Gamma(\alpha)y!(1+\beta)^{\alpha+y}}\end{eqnarray}$$となるが、これは$$p(y) = \begin{pmatrix}\alpha+y-1\\ y\end{pmatrix}\left(\frac{\beta}{\beta+1}\right)^{\alpha}\left(\frac{1}{\beta+1}\right)^{y}$$と変換できる。これは、負の二項分布として知られる。$$y\sim \text{Neg-bin}(\alpha, \beta)$$上の式変形は負の二項分布は割合\(\theta\)のポアソン分布とガンマ分布の積であることを示している。つまり、$$\text{Neg-bin}(y\mid \alpha, \beta) = \int{\text{Poisson}(y\mid \theta)\text{Gamma}(\theta\mid \alpha, \beta)d\theta}$$となる。負の二項分布については後の章で触れる。

Poisson model parameterized in terms of rate and exposure

ポアソン分布は時系列のポイント\(y_1, y_2, \cdots, y_n\)に拡張したほうが便利である。$$y_i \sim \text{Poisson}(x_i\theta) \tag{2.14}\label{eq2.14}$$ここで\(x_i\)は説明変数\(x\)の既知の正の値である。\(\theta\)は未知の興味あるパラメータである。疫学ではしばしば\(\theta\)をrateと呼び、\(x_i\)を\(i\)番目の暴露因子という。このモデルでは\(y_i\)に関して交換可能ではないが、\((x, y)_i\)のペアに関しては交換可能である。\(\theta\)の尤度はポアソン分布のモデルを拡張して、$$p(y\mid \theta) \propto \theta^{(\sum_{i=1}^{n}{y_i})}e^{-(\sum_{i=1}^{n}{x_i})\theta}$$となる。ここで、\(\theta\)に依存しない要素は無視している。\(\theta\)に関するガンマ分布が共役となり、$$\theta \sim \text{Gamma}(\alpha, \beta)$$となる。結果として、事後分布は$$\theta \mid y \sim \text{Gamma}\left(\alpha + \sum_{i=1}^{n}{y_i}, \beta + \sum_{i=1}^{n}{x_i}\right) \tag{2.15}\label{eq2.15}$$となる。

Exponential model

指数分布はしばしば待ち時間などの連続で正の実数値を取る確率変数に対して用いられ、しばしば時間スケールで測定される。結果\(y\)のサンプリング分布は与えられたパラメータ\(\theta\)のもとで$$p(y\mid \theta) = \theta \exp(-y\theta)\ \ (y > 0)$$と表される。\(\displaystyle \theta = \frac{1}{E(y\mid \theta)}\)はrateと呼ばれる。数学的には、指数分布はガンマ分布で\((\alpha, \beta) = (1, \theta)\)とした場合である。しかしこの場合、パラメータ\(\theta\)の事前分布としてではなく、結果\(y\)のサンプリング分布として用いられており、ポアソン分布の例とは異なる。

指数分布で特徴的なのは、無記憶性であり、このことが生存時間などの自然なモデルとしてこの分布が採用される理由でもある。つまり、ある被験者が時間\(t\)後に生存している確率は、これまで経過した時間とは無関係であり、どのような\(s, t\) に対しても\(\text{Pr}(y > t+s\mid y > s, \theta) = \text{Pr}(y > t\mid \theta)\)が成り立つ。指数分布のパラメータ\(\theta\)の共役事前分布は、ポアソン分布の平均のときと同様に、\(\text{Gamma}(\theta\mid \alpha, \beta)\)となり、対応する事後分布は\(\text{Gamma}(\theta\mid \alpha+1, \beta + y)\)となる。\(n\)個の独立した指数分布の観測\(y = (y_1, y_2, \cdots, y_n)\)と定数rate\(\theta\)に対するサンプリング分布は$$p(y\mid \theta) = {\theta}^n \exp{(-n\bar{y}\theta)} \ \ \bar{y}\geq 0$$となる。これは\(y\)を固定した\(\theta\)の尤度として見たとき、\(\text{Gamma}(n+1, n\bar{y})\)に比例する。したがって、\(\theta\)に関する事前分布\(\text{Gamma}(\alpha, \beta)\)は、総待ち時間\(\beta\)を有する\(\alpha-1\)の指数分布の観測とみなすことができる。

2.7 Example: informative prior distribution for cancer rates

Bayesian inference for the cancer death rates

詳細は教科書参照。ここでは各郡countyごとの腎臓癌死亡率を以下のモデルで推定する。$$y_j \sim \text{Poisson}(10n_j{\theta}_j) \tag{2.16}\label{eq2.16}$$\(y_j\)は郡\(j\)の1980年から1989年に置ける腎臓癌の死亡者数である。\(n_j\)は郡の人口、\({\theta}_j\)は各郡の単位人数あたりの年死亡率である。このモデルは\({\theta}_j\)が郡ごとに異なるという点で\eqref{eq2.14}と異なる。そこで添字を\(i\)ではなく\(j\)としている。ベイズ推論を行うためには、未知のレート\(\theta_j\)についての事前分布が必要である。ここでは簡単のためポワソン分布の共役分布であるガンマ分布を用いる。後に詳しく論じるように、\(\alpha = 20, \beta = 430,000\)がこの期間の米国における腎臓癌死亡率の妥当な事前分布になる。この事前分布は平均\(\displaystyle \frac{\alpha}{\beta} = 4.65\times 10^{-5}\)と標準偏差\(\displaystyle \frac{\sqrt{\alpha}}{\beta} = 1.04\times 10^{-5}\)を有する。事後分布は$$\theta_j\mid y_j \sim \text{Gamma}(20+y_j , 430,000+10n_j)$$となる。平均と分散は$$\begin{eqnarray}E(\theta_j\mid y_j) & = & \frac{20+y_j}{430,000+10n_j}\\ \text{var}(\theta_j\mid y_j) & = & \frac{20+y_j}{(430,000+10n_j)^2}\end{eqnarray}$$で与えられる。事後平均は実際の死亡率\(\displaystyle \frac{y_j}{10n_j}\)と事前分布の平均\(\displaystyle \frac{\alpha}{\beta} = 4.65\times 10^{-5}\)の重み付き平均になる。

2.8 Noninformative prior distributions

Proper and improper prior distributions

再度、既知の分散\({\sigma}^2\)を有する正規モデル\(\text{N}({\mu}_0, {{\tau}_0}^2)\)の平均\(\theta\)の推定に戻る。事前精度\(\displaystyle \frac{1}{{{\tau}_0}^2}\)がデータの精度\(\displaystyle \frac{n}{{\sigma}^2}\)に比べて小さい場合、事後分布はおおよそ\({{\tau}_0}^2 = \infty\)のように振る舞い、$$p(\theta\mid y) \approx \text{N}\left(\theta\mid \bar{y}, \frac{{\sigma}^2}{n}\right)$$となる。別の言い方をすると、事後分布はおおよそ\(p(\theta)\)が定数\(\theta\in (-\infty, \infty)\)に対して比例すると仮定した場合の分布と同じになる。このような状況は\(p(\theta)\)の積分が無限大になり、確率が\(1\)になるという前提を破壊するので可能ではない。一般に、事前分布\(p(\theta)\)がデータによらず、積分して\(1\)になるとき、これを適切なproperな事前分布という。次の例として、既知の平均を持つが、分散は未知である正規分布を無情報事前分布、その共役事前分布を尺度付き逆\(\chi\)二乗分布として考える。事前の自由度\(\nu_0\)がデータの自由度\(n\)に対して比較的小さい場合、事後分布はおおよそ\(\nu_0 = 0\)としたときのように、$$p({\sigma}^2\mid y) \approx \text{Inv-}{\chi}^2({\sigma}^2\mid n, v)$$となる。この事後分布の限定された形は、\({\sigma}^2\)の事前分布を\(p({\sigma}^2) \approx 1/{\sigma}^2\)と定義することによっても導き出すことができるが、これは範囲\((0, \infty)\)にわたって無限積分を持つという不適切なものである。

Improper prior distributions can lead to proper posterior distributions

二つの例のどちらも事前分布と尤度の組み合わせは適切な結合確率モデル\(p(y, \theta)\)に至らなかった。しかし、ベイズ推定の手法を用いて正規化しない事後密度関数を$$p(\theta\mid y) \propto p(y\mid \theta)p(\theta)$$で定義することができる。不適切な事前分布から得られる事後分布の解釈には十分な注意が必要であり、事後分布が有限積分であることと、意味のある形式であることを常にチェックしなければならない。最も合理的な解釈は尤度が事前密度を支配している場合の近似である。ベイズ解析のこの側面については、第4章でより詳細に述べられる。

Jefferys’ invariance principle

無情報事前分布を定義するときにしばしば用いられる手法は、Jeffryによって導入された。

Harold Jeffreys - Wikipedia

en.wikipedia.org

これはパラメータの一対一変換\(\phi = h(\theta)\)を考慮するものである。$$p(\phi) = p(\theta)\left|\frac{d\theta}{d\phi}\right| = p(\theta)|h^{\prime}(\theta)|^{-1} \tag{2.19}\label{2.19}$$

Jefferys’ general principleとは、事前密度\(p(\theta)\)を決定するためのどんなルールも、変換されたパラメータに適用すれば同等の結果をもたらすはずだということである。つまり\(p(\theta)\)を決定するため\eqref{2.19}を用いて\(p(\phi)\)が計算されたとき、この結果は変換モデル\(p(y, \phi) = p(\phi)p(y\mid \phi)\)を用いて直接得られた\(p(\phi)\)と一致しなければならない。

Jefferys’s principleによると\(p(\theta)\propto [J(\theta)]^{\frac{1}{2}}\)のように無情報事前分布が決定される。ただし、\(J(\theta)\)は\(\theta\)のFisher情報量であり、以下のようである。$$\begin{eqnarray}J(\theta) & = & E\left(\left(\frac{d\log{p(y\mid \theta)}}{d\theta}\right)^2\middle| \theta\right)\\ & = & -E\left(\frac{d^2 \log{p(y\mid \theta)}}{d{\theta}^2}\middle| \theta\right) \tag{2.20}\label{eq2.20}\end{eqnarray}$$