「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。データ分析・語学に力点を置いています。 →現在、コンサルタントの雛になるべく、少しずつ勉強中です(※2024年1月21日改訂)。

MENU

ベイズ統計学への入門(その03/X)

はじめに

 さまざまなテキスト

などを参照しながらベイズ統計学について学んでいきます。
 また理論だけでなく、可能な限りシミュレーションを含めていくこととし、それも\mathrm{R},\mathrm{Stan},\mathrm{Python}\mathrm{Julia}など幅広い言語で実装していきたい。

各種バージョン情報

  • OS

     Windows 11 Home 22H2
  • R

     R-4.1.3
  • RStudio

     RStudio 2022.02.2+485 "Prairie Trillium" Release (8acbd38b0d4ca3c86c570cf4112a8180c48cc6fb, 2022-04-19) for Windows
  • Python

     3.11.0
  • Jupyter Notebook

     6.4.12
  • Julia

     1.8.0

今回のまとめ

  • 分布に関する\mathrm{Bayes}の定理により、
    \begin{aligned}f(\theta|x)=\displaystyle{\frac{f(x,\theta)}{f(x)}}=\displaystyle{\frac{f(x|\theta)f(\theta)}{\displaystyle{\int_{-\infty}^{\infty}f(x|\theta)f(\theta)d\theta}}}\end{aligned}
    が成り立つ。
  • 確率分布において母数および変数を含む部分をカーネルという。これに対して正規化定数(係数)を全事象の確率が1になるように規格化するための定数部分を指す。
  • \mathrm{Bayes}推定においては、具体的に推定値を得るのにたとえば①EAP定量、②MAP推定量、③MED推定量を用いる。

3. Bayes推定

3.3 Bayes推定の一般的な枠組み

 真の分布G(\boldsymbol{x})に従う確率変数を\boldsymbol{X}とし、標本\boldsymbol{X}_n=\{\boldsymbol{x}_1,\cdots,\boldsymbol{x}_n\},\boldsymbol{x}_i\in\mathbb{R}^m,i=1,\cdots,nが得られたとする(すなわち各\boldsymbol{x}_iが各試行により得た1つ1つの標本を指す。)。また各標本は独立に真の分布に従うとする。
 目標としてはこの標本を基に真の分布を推測するであるが、そのために統計モデル\{f(\boldsymbol{x}|\boldsymbol{\theta}; \boldsymbol{\theta}\in\Theta)\}を導入する。すなわちG(\boldsymbol{x})を推測するという問題を、統計モデル\{f(\boldsymbol{x}|\boldsymbol{\theta}; \boldsymbol{\theta}\in\Theta)\}用いて予測分布を推定する問題に置き換えるのである*1
 定数0\lt\beta\lt\infty(逆温度という)を導入し、逆温度\betaの母数\thetaの事後分布を


\begin{aligned}
f(\boldsymbol{\theta}|\boldsymbol{X}_n)=\displaystyle{\frac{1}{Z_n(\beta)}\pi(\boldsymbol{\theta})}f(\boldsymbol{X}_n|\boldsymbol{\theta})^{\beta}
\end{aligned}

で定義する。ここでパラメータの確率分布\pi(\boldsymbol{\theta})を事前分布として設定している。またZ_n(\beta)は分配関数と呼ばれ、


\begin{aligned}
Z_n(\beta)=\displaystyle{\int_{\Theta}\pi(\boldsymbol{\theta})f(\boldsymbol{X}_n|\boldsymbol{\theta})^{\beta}d\boldsymbol{\theta}}
\end{aligned}

で定義される規格化のための定数である。\beta=1のときZ_n(1)を周辺尤度という。
 事後分布f(\boldsymbol{\theta}|\boldsymbol{X}_n)を用いて統計モデルf(\boldsymbol{x}|\boldsymbol{\theta}; \boldsymbol{\theta}\in\Theta)の期待値を取ったもの


\begin{aligned}
f^{*}(\boldsymbol{x})=\displaystyle{\int_{\Theta}f(\boldsymbol{x}|\boldsymbol{\theta}) f(\boldsymbol{\theta}|\boldsymbol{X}_n)d\boldsymbol{\theta}}
\end{aligned}

を予測分布という。\mathrm{Bayes}推定とは、「真の分布はだいたいのところ予測分布であろう」と推測することに他ならない。
 \beta=1の場合、\mathrm{Bayes}の定理に帰着でき、


\begin{aligned}
\pi(\boldsymbol{\theta}|\boldsymbol{X}_n)=\displaystyle{\frac{f(\boldsymbol{X}_n|\boldsymbol{\theta})\pi(\boldsymbol{\theta})}{\displaystyle{\int f(\boldsymbol{X}_n|\boldsymbol{\theta})\pi(\boldsymbol{\theta}) d\boldsymbol{\theta}}}} \propto f(\boldsymbol{X}_n|\boldsymbol{\theta})\pi(\boldsymbol{\theta})
\end{aligned}

が得られる。
 母数の事前分布から得た


\begin{aligned}
f_{\mathrm{anterior}}(\boldsymbol{x})=\displaystyle{\int_{\Theta}f(\boldsymbol{x}|\boldsymbol{\theta}) f(\boldsymbol{X}_n|\boldsymbol{\theta})d\boldsymbol{\theta}}
\end{aligned}


を事前予測分布といい、これに対して事後分布を用いて評価した


\begin{aligned}
f_{\mathrm{posterior}}(\boldsymbol{x})=\displaystyle{\int_{\Theta}f(\boldsymbol{x}|\boldsymbol{\theta}) f(\boldsymbol{\theta}|\boldsymbol{X}_n)d\boldsymbol{\theta}}
\end{aligned}

を事後予測分布という。
 確率変数としてX,\thetaがあり、\thetaXの確率分布の母数だとする。このときX\thetaで条件づけられた確率をF(\theta)とすれば、X|\theta\sim F(\theta)と書く。

3.4 事前分布の設定

 事前分布は、\mathrm{Bayes}統計学のアイディアから言えば自由に選べばよい。しかし

  ①科学的な客観性を担保するのに問題が生じ得る

  ②個人の想定する確率を数理的に表現しきることが困難である

という問題がある。そのため、どのような事前分布を選ぶかが重要な問題になる。

3.4.1 自然共役事前分布

 \mathrm{Bayes}統計学では事後分布が常に計算可能だとは限らないという問題がある。そこで事後分布が求まるような事前分布を選ぶという方針によりこれを解消することができる。自然共役な事前分布は、計算機がそこまで発展していなかった時代に用いられることが多かった。今では数値的に計算することが容易になってきたことから、それほど使われなくなってきている。


図表1 自然共役事前分布と尤度の組み合わせ例

      
尤度
事前分布
事後分布
     (1) Bernoulli分布 Beta分布 Beta分布
     (2) 二項分布 Beta分布 Beta分布
     (3) Poisson分布 Gamma分布 Gamma分布
     (4) 正規分布の平均 正規分布 正規分布
     (5) 正規分布の分散 逆ガンマ分布 逆ガンマ分布

出典:豊田(2015)*2

3.4.2 無情報事前分布

 事前に何の情報もないことを事前分布として与えた場合、これを無情報事前分布と呼ぶ。
 具体的な分布の与え方にはいくつかの考え方がある。1つは理由不充分の原則と呼ばれるもので、これはどのようなパラメータ値が出るのかは全く未知、すなわち確率が一様に分布していると考え、その状況を上手く表すような分布を事前分布として仮定するものである。


\begin{aligned}
\pi(\boldsymbol{\theta})=C,\ \boldsymbol{\theta}\in A,\ C\in\mathbb{R}
\end{aligned}

 後述するように、無情報事前分布に対しては批判が存在するものの、

 ①現実的なデータ解析が採用した尺度上で利用されるのが普通で、(後述する)尺度変換を考慮する必要性が薄い場合が少なくない、

 ②一様分布を用いると\mathrm{MAP}定量最尤推定量に一致するため便利である、

という利点が存在する。よく一様分布が用いられるが、他に分散が充分に大きい正規分布を与える場合がある。
 一様分布を採用した場合、実数全体で定義すれば積分したときに1にならない場合があり得る。その場合を非正則といい、1になる場合を正則という。非正則な事前分布を選んだとしても、事後分布が非正則になるとは限らない。したがって非正則な事前分布を選択することもあり得る。
 他方で、無情報事前分布を用いることにはいくつかの問題がある:

  • 常に漠然としている事前分布を探すことは常に誤解を招くように見える。もし尤度が与えられた問題において真に優勢であるならば、相対的に平坦な事前確率密度の範囲から選ぶのは重要ではない。無情報事前分布として特定の特徴を持たせることは自動的にそして場合によっては不適切な利用を助長するように見える。
  • 多くの問題において、特に情報の無い事前分布の明確な選択肢がない。なぜならば、あるパラメータ化で平坦ないし一様な確率密度は別のパラメータ化ではそうでないからである。すなわち母数変換すると平坦さや一様性が失われるのである。
  • 不適切な事前分布を持つ競合するモデルのモデル平均化法において更なる困難が生じる。
3.4.3 局所一様事前分布

 無情報事前分布として一様分布を用いることは、実際には無情報を表し得ない側面がある。それは尺度変換に対して一様性が保てない点を指している。これに対処するのに、\mathrm{Fisher}情報量の平方根に比例する分布を導入でき、これを局所一様事前分布もしくは\mathrm{Jeffreys}の事前分布という。
 具体的には、\mathrm{Jeffreys}(1961)は事前情報が無い場合に設定すべき事前分布として、


\begin{aligned}
\pi(\boldsymbol{\theta})\propto\left|J(\boldsymbol{\theta})\right|^{\frac{1}{2}}
\end{aligned}

を提案した。ここでJ(\boldsymbol{\theta})\mathrm{Fisher}情報行列


\begin{aligned}
J(\boldsymbol{\theta})=\displaystyle{-\int\left[\frac{\partial^2f(\boldsymbol{x}|\boldsymbol{\theta})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\top}}\right]f(\boldsymbol{x}|\boldsymbol{\theta}) d \boldsymbol{x}}
\end{aligned}

である。局所一様事前分布は母数の1対1変換\boldsymbol{\psi}=(r_1(\boldsymbol{\theta},\cdots,r_p(\boldsymbol{\theta}))^{\top}に対して不偏、すなわち\pi(\boldsymbol{\psi})\propto\left|J(\boldsymbol{\psi})\right|^{\frac{1}{2}}になる点がメリットである。

3.4.4  弱情報事前分布

 事前分布が、妥当ではあるものの実際に事前知識として用いることのできるあらゆるものよりもそれが持つ情報が意図的に少なくなっているのであれば、その事前分布は弱情報を持つという。一般論として、あらゆる問題において弱情報分布によるモデルを許容し得る何らかの制約条件が存在する。完全な無情報をモデル化しようとするよりも、事後分布が意味を成すことが充分に保障されるが実世界の情報を少量だけ含む弱情報事前分布を用いる方が大半の問題において好ましい。

 あらゆる統計モデルが実質的に弱情報を有する。入力やそれらを結合する関数形の選び方により、統計モデルは常に何らかの情報を有している。しかし確率空間に関する事柄に関する事前の信条をすべて符号化することは不可能だし、場合によってはむしろそうしないことが望ましい。これを踏まえた上で、弱情報事前分布を設定するのに2つの原則を提案できる。

  • 無情報事前分布を変形させたものから始め、推測が合理的であるように制約づけるのに十分な情報を加える。
  • 強情報を持つ事前分布から始め事前の信条として存在する、また任意の経歴に基盤を持つ事前情報を新しいデータへ適用できるような不確実を説明するように拡張する。

3.5 事後分布の評価

 すべての\mathrm{Bayes}推測は事後分布を用いて行う。特定の値を推定によって得た値として用いて母数を評価する方法を点推定という。
 点推定は、その推定した値(推定値という)が過大ないし過少に評価された場合に損失が増大すると見なし、その損失がなるべく小さくなるような値を推定値とするという考え方の下で構成される。より具体的には、母数\thetaおよびその推定量*3\hat{\theta}の損失関数L(\theta,\hat{\theta})を定義し、その損失関数の大小から得るものとする。
 \mathrm{Bayes}推測では以下のような損失関数を用いるのが代表的である:


\begin{aligned}
二乗誤差損失:L(\theta,\hat{\theta})&=(\theta-\hat{\theta})^2,\\
絶対誤差損失:L(\theta,\hat{\theta})&=\left|\theta-\hat{\theta}\right|,\\
0-1損失:L(\theta,\hat{\theta})&=1-\boldsymbol{1}_{\theta}\left(\hat{\theta}\right)\\
&=\begin{cases}0,\hat{\theta}=\theta\\1,\hat{\theta}\neq\theta\end{cases}
\end{aligned}

上述の通り、損失関数を具体的に計算するには\theta,\hat{\theta}それぞれの値が分からないとできない。しかし\thetaの値が不明であるからこそ、このような理論を構築しているのであるから、具体的にその値を計算することはできないことになる。そのため代わりに母数の事後分布による損失関数の期待値で定義されるリスク関数


\begin{aligned}
R\left(\theta|\boldsymbol{x}\right)=E\left[L(\theta,\hat{\theta})|\boldsymbol{x}\right]=\displaystyle{\int_{0}^{1}L(\theta,\hat{\theta})f(\theta|\boldsymbol{x})d\theta}
\end{aligned}

を最小化するような推定値\hat{\theta}を求めることとする。

3.5.1 期待事後推定

 1つ目は事後分布を用いた期待値をその母数の推定値として用いる方法であり、これを事後平均(\mathrm{EAP}: \mathrm{Expected\ A\ Posteriori})という。すなわち


\begin{aligned}
\hat{\theta}_{\mathrm{EAP}}=E\left[\theta|\boldsymbol{x}\right]=\displaystyle{\int\theta f(\theta|\boldsymbol{x})d\theta}=\displaystyle{\int\theta\frac{f(\boldsymbol{x}|\theta)f(\theta)}{f(\boldsymbol{x})}d\theta}
\end{aligned}

である。これは平均二乗誤差、すなわち二乗誤差損失関数の期待値を最小にするような推定量である。
 実際、二乗誤差損失関数の期待値を計算すると、


\begin{aligned}
R\left(\theta|\boldsymbol{x}\right)=&E\left[(\theta-\hat{\theta})^2|\boldsymbol{x}\right]\\
=&E\left[\left(\theta-E[\theta|\boldsymbol{x}]+E[\theta|\boldsymbol{x}]-\hat{\theta}\right)^2|\boldsymbol{x}\right]\\
=&E\left[(\theta-E[\theta|\boldsymbol{x}])^2|\boldsymbol{x}\right]-2\left(\theta-E[\theta|\boldsymbol{x}]\right)\left(E[\theta|\boldsymbol{x}]-\hat{\theta}\right)\\
&+E\left[\left(E[\theta|\boldsymbol{x}]-\hat{\theta}\right)^2|\boldsymbol{x}\right]\\
=&E\left[(\theta-E[\theta|\boldsymbol{x}])^2|\boldsymbol{x}\right]+E\left[\left(E[\theta|\boldsymbol{x}]-\hat{\theta}\right)^2\left|\right.\boldsymbol{x}\right]
\end{aligned}

である。最右辺の第1項は事後分布の分散で\hat{\theta}とは無関係である。したがって最右辺第2項を最小にするような\hat{\theta}を求めればよく、その値\hat{\theta}_{\mathrm{EAP}}


\begin{aligned}
\hat{\theta}_{\mathrm{EAP}}&=E\left[\theta|\boldsymbol{x}\right]=\displaystyle{\int\theta f(\theta|\boldsymbol{x})d\theta}\\
&=\displaystyle{\int\theta\frac{f(\boldsymbol{x}|\theta)f(\theta)}{f(\boldsymbol{x})}d\theta}
\end{aligned}

に他ならない。

3.5.2 事後中央値

 2つ目は事後分布に基づいた中央値をパラメータの推定値として用いる方法であり、これを事後中央値(\mathrm{MED}: \mathrm{Median\ A\ Posteriori})という。すなわち


\begin{aligned}
\hat{\theta}_{\mathrm{MED}}=\theta_0\ \mathrm{s.t.}\ \displaystyle{\int_{-\infty}^{\theta_0}f(\theta|\boldsymbol{x})d\theta=\frac{1}{2}}
\end{aligned}

である。これは平均絶対誤差、すなわち絶対誤差損失関数の期待値を最小にするような推定量である。
 実際、


\begin{aligned}
R\left(\theta|\boldsymbol{x}\right)=&E\left|\theta-\hat{\theta}|\boldsymbol{x}\right|\\
=&\displaystyle{\int_{0}^{1}\left|\theta-\hat{\theta}\right|f(\theta|\boldsymbol{x})d\theta}\\
=&\displaystyle{\int_{0}^{\hat{\theta}}\left(\hat{\theta}-\theta\right)f(\theta|\boldsymbol{x})d\theta}+\displaystyle{\int_{\hat{\theta}}^{1}\left(\theta-\hat{\theta}\right)f(\theta|\boldsymbol{x})d\theta}\\
=&\hat{\theta}F\left(\hat{\theta}|\boldsymbol{x}\right)-\displaystyle{\int_{0}^{\hat{\theta}}\theta f(\theta|\boldsymbol{x})d\theta}\\
&+\displaystyle{\int_{\hat{\theta}}^{1}\theta f(\theta|\boldsymbol{x})d\theta}-\hat{\theta}\left(1-F\left(\hat{\theta}|\boldsymbol{x}\right)\right)\\
=&2\hat{\theta}F\left(\hat{\theta}|\boldsymbol{x}\right)-\hat{\theta}-2\displaystyle{\int_{0}^{\hat{\theta}}\theta f(\theta|\boldsymbol{x})d\theta}+\displaystyle{\int_{0}^{1}\theta f(\theta|\boldsymbol{x})d\theta}\\
\end{aligned}

である。
 ここで最右辺第3項について部分積分から


\begin{aligned}
\displaystyle{\int_{0}^{\hat{\theta}}\theta f(\theta|\boldsymbol{x})d\theta}&=\left[\theta F(\theta|\boldsymbol{x})\right]_{0}^{\hat{\theta}}-\displaystyle{\int_{0}^{\hat{\theta}}F(\theta|\boldsymbol{x})d\theta}
\end{aligned}

であること、また最右辺第4項が事後分布による\thetaの期待値E\left[\theta|\boldsymbol{x}\right]に他ならないことに注意すれば、


\begin{aligned}
R\left(\theta|\boldsymbol{x}\right)&=2\hat{\theta}P\left(\hat{\theta}|\boldsymbol{x}\right)-\hat{\theta}-2\displaystyle{\int_{0}^{\hat{\theta}}\theta f(\theta|\boldsymbol{x})d\theta}+\displaystyle{\int_{0}^{1}\theta f(\theta|\boldsymbol{x})d\theta}\\
&=2\hat{\theta}P\left(\hat{\theta}|\boldsymbol{x}\right)-\hat{\theta}-2\left\{\left[\theta F(\theta|\boldsymbol{x})\right]_{0}^{\hat{\theta}}-\displaystyle{\int_{0}^{\hat{\theta}}F(\theta|\boldsymbol{x})d\theta}\right\}+E\left[\theta|\boldsymbol{x}\right]\\
&=2\displaystyle{\int_{0}^{\hat{\theta}}F(\theta|\boldsymbol{x})d\theta}-\hat{\theta}+E\left[\theta|\boldsymbol{x}\right]
\end{aligned}

を得る。これを最小にするような値を\hat{\theta}_{\mathrm{MED}}とおけば、それが満たすべき1階と2階の条件は


\begin{aligned}
\nabla_{\hat{\theta}}R(\hat{\theta}_{\mathrm{MED}}|\boldsymbol{x})&=2F\left(\hat{\theta}_{\mathrm{MED}}|\boldsymbol{x}\right)-1=0,\\
\nabla_{\hat{\theta}}^2R(\hat{\theta}_{\mathrm{MED}}|\boldsymbol{x})&=2f(\hat{\theta}_{\mathrm{MED}}|\boldsymbol{x})
\end{aligned}

である。2階条件は確率密度関数の性質から\hat{\theta}_{\mathrm{MED}}の値に関わらず必ず満たされるから、求めるべき\hat{\theta}_{\mathrm{MED}}


\begin{aligned}
F\left(\hat{\theta}_{\mathrm{MED}}|\boldsymbol{x}\right)=\displaystyle{\frac{1}{2}}
\end{aligned}

を満たすような推定量であるが、これは\hat{\theta}_{\mathrm{MED}}が事後分布の中央値であることに他ならない。

3.5.3 事後確率最大推定

 3つ目は事後確率を最大にするパラメータ値を推定値として用いる方法であり、これを事後確率最大(\mathrm{MAP}: \mathrm{Maximum\ A\ Posteriori})という。すなわち


\begin{aligned}
\hat{\theta}_{\mathrm{MAP}}=\displaystyle{\arg\max_{\theta} f(\theta|\boldsymbol{x})}
\end{aligned}

である。これは0-1誤差の平均、すなわち0-1誤差損失関数の期待値を最小にするような推定量である。
 実際、


\begin{aligned}
R\left(\theta|\boldsymbol{x}\right)&=E\left[1-\boldsymbol{1}_{\theta}\left(\hat{\theta}\right)|\boldsymbol{x}\right]\\
&=1-E\left[\boldsymbol{1}_{\theta}\left(\hat{\theta}\right)|\boldsymbol{x}\right]\\
&=1-f\left(\hat{\theta}|\boldsymbol{x}\right)
\end{aligned}

であるから、求めたい推定量\hat{\theta}_{\mathrm{MAP}}とおけば、\hat{\theta}_{\mathrm{MAP}}


\begin{aligned}
R\left(\theta|\boldsymbol{x}\right)=1-f\left(\hat{\theta}|\boldsymbol{x}\right)
\end{aligned}

を最小にするような、すなわち事後確率密度関数を最大にするような\thetaが求めたい\hat{\theta}_{\mathrm{MAP}}であり、これは事後分布のモードに他ならない。
 なお実務上は、最尤推定と同様に、事後分布の対数を取り、それを母数で微分したものを0とおいて母数について解くことで事後確率最大値を得る。

3.5.4 事後分散・事後標準偏差

 これまでで述べた母数の推定値の“良さ”を評価するために、その散らばり具合を表す尺度である事後分散・事後標準偏差を定義する。
 事後分散は


\begin{aligned}
V[\theta]=V[\theta|\boldsymbol{x}]=E[(\theta-\hat{\theta}_{\mathrm{EAP}})^2|\boldsymbol{x}]=\displaystyle{\int(\theta-\hat{\theta}_{\mathrm{EAP}})^2f(\theta|\boldsymbol{x})d\theta}
\end{aligned}

で定義され、\mathrm{EAP}定量の散らばりの尺度とする。

参考文献

*1:真の分布は未知である場合が殆どであるから、統計モデルの想定する分布が真の分布と同一の分布族に属するとは限らない。

*2:豊田秀樹・編(2015)「基礎からのベイズ統計学」朝倉書店 P.50参照。

*3:定量は母数を標本の関数としてみたもので、標本が具体的に得られたときに具体的に得られた推定量の値を推定値という。

プライバシーポリシー お問い合わせ