◆この記事のテーマ

ベイズファクターを使って平均値の差を検定する方法を整理します。この手法は、例えば下図のようなシーンに適用できます。下図左側では、交通事故対策Aの実施前と実施後の単位時間あたり事故発生件数のデータを蓄積し、対策前後で事故発生件数の平均値に差があるかを検定しています。また、下図右側では交通事故対策Aが実施された交差点と実施されていない交差点の単位時間あたり事故発生件数の平均値に差があるかを検定しています。

※上記の事例は整数データであるからポアソン分布や負の二項分布等を使用するのが理想ですが、ここでは事故発生件数の平均値が十分大きく正規近似が可能であると仮定します。

※左側は対応のある（one-sample）グループ間の検定、右側は対応のない（two-sampe）グループ間の検定です。設定上、後者の設計では対策実施効果以外の要因（異なる交差点を比較に用いることによるランダム効果）が結果に絡む可能性が高くなるので、対策Aの効果を確かめることが分析の目的であれば前者の対応のある検定の設計がより望ましいといえます。

参考文献は以下の通り。ベイズファクターを使って平均値の差の検定が英語なら文献が多い。

effect size $\delta$の導入

ここでは、対応のある場合、対応のない場合に分けて検定の対象となるパラメータ$\delta$を導入します。

対応のある場合

同じ対象に対し、2つの異なる状態で観測された値をそれぞれ確率変数$X_i$、$Y_i$からの観測値$x_i$、$y_i$とおき、独立同分布の確率変数の組$X^n=(X_1,X_2,\ldots,X_n)$、$Y^n=(Y_1,Y_2,\ldots,Y_n)$から観測値$x^n=(x_1,x_2,\ldots,x_n)$、$y^n=(y_1,y_2,\ldots,y_n)$を得たとします。

まずは、簡単に確率変数$X_i$、$Y_i$の差が平均$\mu$の正規分布に従うと記載しちゃいます。

$$ X_i - Y_i \sim \mathrm{Normal}(\mu, \sigma^2), ~~~ i=1,\ldots,n \tag{1} $$

$(1)$式の唐突に現れた分散$\sigma^2$は$(X_i - Y_i)$のスケールを表現するパラメータです。これはそのときどきの検定対象によって異なる値をとります。このとき、平均$\mu$と$\sigma$の間に何の関係も与えない場合、$\mu$の値の大小がどれくらいの意味をもつものなのか把握が困難になってしまいます。そこで、平均$\mu$をeffect sizeと呼ばれる値$\delta$とスケールを指定する値に分け、スケールを指定する値を$\sigma$としちゃいます。

$$ \mu = \delta \sigma \tag{2} $$

これでようやくスケールが単位あたりに調整されたeffect size$\delta$に着目すれば良いことになり、以降の議論を一般的な形で記述できるようになりました。

◆$X_i$及び$Y_i$と$\delta$の関係のモデル化（対応のある場合） $$ X_i - Y_i \sim \mathrm{Normal}(\delta \sigma, \sigma^2), ~~~ i=1,\ldots,n \tag{3} $$

対応のない場合

グループ1、グループ2に所属する対象からの観測値をそれぞれ確率変数$X_i$、$Y_i$からの観測値$x_i$、$y_i$とおき、独立同分布の確率変数の組$X^n=(X_1,X_2,\ldots,X_n)$、$Y^m=(Y_1,Y_2,\ldots,Y_m)$から観測値$x^n=(x_1,x_2,\ldots,x_n)$、$y^m=(y_1,y_2,\ldots,y_m)$を得たとします。

まずは、$X_i$、$Y_j$がともに分散$\sigma$の正規分布に従い、また2つの分布の平均が$\alpha$だけ離れているとおきます。

$$ \begin{cases} X_i \sim \mathrm{Normal}\left(\mu - \cfrac{\alpha}{2}, \sigma^2\right),~~~i=1,\ldots,n \\\\ Y_j \sim \mathrm{Normal}\left(\mu + \cfrac{\alpha}{2}, \sigma^2\right),~~~j=1,\ldots,m \end{cases} \tag{4} $$

この場合もやはり同様に、2つの分布の平均値の差$\alpha$と分散$\sigma$の間に何の関係も与えない場合、$\alpha$の値の大小がどれくらいの意味をもつものなのか把握が困難になってしまいます。そこで、やはり同様に

$$ \alpha=\delta \sigma \tag{5} $$

とおいてやれば、ようやくスケールが単位あたりに調整されたeffect size$\delta$に着目すれば検定が可能になり、こちらも以降の議論を一般的な形で記述できるようになります。

◆$X_i$及び$Y_i$と$\delta$の関係のモデル化（対応のない場合） $$ \begin{cases} X_i \sim \mathrm{Normal}\left(\mu - \cfrac{\delta \sigma}{2}, \sigma^2\right),~~~i=1,\ldots,n \\
Y_j \sim \mathrm{Normal}\left(\mu + \cfrac{\delta \sigma}{2}, \sigma^2\right),~~~j=1,\ldots,m \end{cases} \tag{6} $$

仮説の設定

$(3)$、$(6)$式の$\delta$に着目し、従来通りに帰無仮説と対立仮説を設定します。

◆帰無仮説と対立仮説 従来の統計的仮説検定に対応する仮説は、 $\delta$の正負が明らかでない両側検定の場合、 $$ \mathrm{Unrestricted~~Model:} \begin{cases} H_0:\delta=0 \\\\ H_1:\delta\neq0 \end{cases} $$ $\delta$が0以上の値をとることが明らかな片側検定の場合、 $$ \mathrm{Order-restricted~~Model:} \begin{cases} H_0:\delta=0 \\\\ H_2:\delta > 0 \end{cases} $$

以降では説明の簡略化のため、両側検定の場合の仮説を前提にします。

事前分布の設定

さて、ベイズの枠組みではすべての母数に事前分布を設定してやる必要があります。この事前分布の役割は、分析者が前もって知りうる情報を反映させるものである…と思っていませんか？
この考え方は主観ベイズ（subjective Bayes）と呼ばれており、間違ってはいません。しかし今回の目的は検定です。それもベイズファクターを用いた検定であり、事前分布の設定が検定結果に大きく左右されることからも、なるべく事前分布に主観を与えたくありません。このような態度を、客観ベイズ（objective Bayes）と呼び、ベイズファクターを用いた仮説検定の場合は客観ベイズの立場で事を考えることが重要になってきます（Rouder et al.(2009)）。

古今東西、客観ベイズの分野（それもなぜか心理学研究分野）では、研究者の間で合意のとれた検定のための汎用性のある事前分布の設定方法が盛んに研究されてきました。平均値の差の検定のための汎用的な事前分布もいくつか開発されています。ここでは天下り的に、JZSの事前分布（JZS proir）と単位情報事前分布（unit information proir）について整理します。

◆JZSの事前分布

JZSの事前分布は検定対象のパラメータ$\delta$まわりの事前分布を以下のように設定する。

$$ \begin{cases} \delta \sim \mathrm{Normal}(0,g) \\
g \sim \mathrm{Inverse}~\mathrm{Gamma}\left(\cfrac{1}{2}, \cfrac{r^2}{2}\right) \\
\pi\left(\sigma^2\right) ∝ \cfrac{1}{\sigma^2} \\
\end{cases} \tag{7} $$

ここで、$r$はscale factorとよばれ、effect sizeのスケールを設定するために分析者が指定するパラメータである。また$\pi(x)$はパラメータ$x$の事前分布である。

これは、下記設定に等しい。

$$ \begin{cases} \delta \sim \mathrm{Cauchy}(0,r) \\
\pi\left(\sigma^2\right) ∝ \cfrac{1}{\sigma^2} \\
\end{cases} \tag{8} $$

以上が対応のある平均値の差の検定における事前分布の設定である。対応のない平均値の差の検定ではさらに$\mu$に以下のとおり事前分布を設定する。

$$ \pi(\mu) ∝ 1 \tag{9} $$

memo

コーシー分布は下記式で定義される。 $$ \mathrm{Cauthy}(y | \mu, \sigma) = \cfrac{1}{\pi \sigma}\cfrac{1}{1 + ((y - \mu)/\sigma)^2} $$ $\mu=0$、$\sigma=1$としたコーシー分布を標準コーシー分布と呼び、下記式で定義される。 $$ \mathrm{Cauthy}(y | \mu=1, \sigma=1) = \cfrac{1}{\pi}\cfrac{1}{1 + y^2} $$ 逆ガンマ分布は下記式で定義される。 $$ \mathrm{Inverse}~\mathrm{Gamma}\left( y | \alpha, \beta \right) = \cfrac{\beta^{\alpha}}{\Gamma(\alpha)}y^{-\alpha-1}\exp\left( -\cfrac{\beta}{y} \right) ~~~(y>0) $$

◆単位情報事前分布

単位情報事前分布は検定対象のパラメータ$\alpha$まわりの事前分布を以下のように設定する。

$$ \begin{cases} \delta \sim \mathrm{Normal}(0,r^2) \\
\sigma_{\alpha}^2 ∝ \cfrac{1}{\sigma^2} \\
\end{cases} \tag{10} $$

$$ \pi(\mu) ∝ 1 \tag{11} $$

JZSの事前分布と単位情報事前分布の唯一の違いはeffect size $\delta$の事前分布（正規分布）の分散の設定です。JZSの事前分布は分散$g$が逆カイ二乗分布に従うとした一方、単位情報事前分布は分散固定です。その意味では、単位情報事前分布の方がモデルに与える事前情報が多いことになります。また$r=1$のとき、2つの設定は$\delta$の事前分布を標準正規分布にするか、標準コーシー分布にするかの違いに帰着することが分かります。

標準正規分布と標準コーシー分布を比較した図が以下です。

図からも、標準コーシー分布は標準正規分布と比較し裾が重く、0から外れた値もとりやすい分布になっていることから、$r=1$の場合、JZSの事前分布は単位情報事前分布より事前情報が少ないことが確認できます。

また、両設定に共通のパラメータ$\sigma$及び、対応のある検定の時のみ必要な$\mu$については、ともにJeffeysの事前分布と呼ばれる事前分布を用いています。Jeffreysの事前分布はフィッシャー情報量から導かれるもので、事前分布と事後分布のカルバック・ライブラー情報量を最大化する効果がある、という点で無情報事前分布として使えるとされているようです。このあたりは詳しくないですが、Appendixに定義から正規分布の平均及び分散に対するJeffreysの事前分布を導出する方法を載せています。

ちなみに、$\sigma$のように検定の対象となるパラメータとは別に存在し、かつ両モデルに共通の役割をもつパラメータはnuisance parameterと称され、こうしたパラメータはベイズファクターの値の決定にクリティカルでないことが証明されています。これに関しては以前の記事を参照してくだせえ。

ベイズファクターの解法

前節でモデルの設定を整理したので、本節ではベイズファクターの求め方を整理します。

解析的に求める

先に説明したJZSの事前分布と単位情報事前分布は解析的に計算することが可能です。ここでは導出や計算式までは深堀りしません。あまり面白くないので…。式を確認したい場合は参考文献を参照のこと。この文献の著者らはJZSの事前分布や単位情報事前分布の解析解を求めることができるWebページやRパッケージを開発しているので、実践編で解析解を求める際はこのあたりを活用させてもらいます。

MCMC法で近似的に求める

本記事と実践編の本題はこちらです。これにはネストされたモデル同士の比較で有効なSavage-Dickey法を使ってベイズファクターを計算します。Savage-Dickey法の概要を以下に示します。詳細は以前の記事を参照のこと。

◆Savage-Dickey法

モデル$M_0$と$M_1$に共通のパラメータ$\boldsymbol{\delta}=(\delta_1, \ldots, \delta_m)$があるとき、$M_0$と$M_1$が

$$ p(\mathcal{D} | M_0) = p(\mathcal{D} | \boldsymbol{\delta} = \boldsymbol{\delta}_0, M_1) \tag{12} $$

とネストされた関係にある場合、ベイズファクター$BF_{01}$は下記のとおり計算できる。

$$ BF_{01} = \cfrac{p(\boldsymbol{\delta}=\boldsymbol{\delta}_0 | \mathcal{D}, M_1)}{p(\boldsymbol{\delta}=\boldsymbol{\delta}_0 |M_1)} \tag{13} $$

ここで、$\mathcal{D}$は全データセットである。

$(11)$、$(12)$式において、$H_0$、$H_1$にそれぞれ対応するモデルを$M_0$、$M_1$と読み替えれば、$H_1$に対応するモデル$M_1$における対象パラメータの事前分布と事後分布の一点比較をするだけでベイズファクター$BF_{01}$を求められることが分かります。

事前分布$p(\boldsymbol{\delta}=\boldsymbol{\delta}_0 |M_1)$の方は、MCMCの結果からではなく、設定したモデルから値を算出します。たとえば、従来通りの仮説検定の枠組みでは、JZSの事前分布、単位情報事前分布の双方ともに$\delta=0$が帰無仮説に対応するパラメータ空間上の点ですから、

$$ \begin{cases} JZSの事前分布の場合： p(\delta=0 |M_1) = \cfrac{1}{\pi r}~~ &\therefore&~~ BF_{01} = \pi r \times p(\delta=0 | \mathcal{D}, M_1)　\\\\ 単位情報事前分布の場合： p(\delta=0 |M_1) = \cfrac{1}{\sqrt{2\pi r^2}}~~ &\therefore&~~ BF_{01} = \sqrt{2\pi r^2} \times p(\delta=0 | \mathcal{D}, M_1)　\tag{14} \end{cases} $$

これが求めたいベイズファクターになります。

あとは$p(\delta=0 | \mathcal{D}, M_1)$が分かればベイズファクターを求められるので、以降はこれを求めることを考えます。

$\delta$の事後分布の求め方

条件の無い周辺事後分布を使う

こちらは全パラメータの同時分布から$\delta$以外のパラメータ$\boldsymbol{\psi}$を周辺化消去することで$\delta$の事後分布を求めています。JZSの事前分布の場合、$\boldsymbol{\psi}=(\sigma,g,(\mu))$、単位情報事前分布の場合、$\boldsymbol{\psi}=\sigma$ですね。$r$は分析者が事前に設定するものなのでこれに含まれません。

$$ p(\delta=0 | \mathcal{D}, M_1) = \int_{\boldsymbol{\psi}}p(\delta=0,\boldsymbol{\psi} | \mathcal{D}, M_1) d\boldsymbol{\psi} = \mathrm{E}_{\boldsymbol{\psi}}[p(\delta=0, \boldsymbol{\psi} |\mathcal{D},M_1)] \tag{15} $$

$(14)$式はMCMCのサンプルを用いて近似的に求めるなら、$\boldsymbol{\psi}$を無視し、$\delta$のサンプルの分布をカーネル密度推定、もしくは以前の記事のように対数スプライン推定をするだけでイケます。簡単で便利です。

条件付き周辺事後分布を使う

こちらは$\boldsymbol{\psi}$に関する条件付きの$\delta$の事後分布を使って$p(\delta=0 | x^n, M_1)$を求める方法です。参考文献に倣ってConditional MArginal Density Estimator を略し、CMDEと呼びます。

$$ p(\delta=0 | \mathcal{D}, M_1) = \int_{\boldsymbol{\psi}}p(\delta=0 |\boldsymbol{\psi}, \mathcal{D}, M_1) d\boldsymbol{\psi} = \mathrm{E}_{\boldsymbol{\psi}}[p(\delta=0, |\boldsymbol{\psi},\mathcal{D},M_1)] \tag{16} $$

$(14)$式の計算では切り捨てていた$\boldsymbol{\psi}$について、その事後分布も$x^n$についての情報を含んでいるのだから、$\boldsymbol{\psi}$も考慮して$p(\delta=0 | x^n, M_1)$を求めてやるべきだ、というのがCMDEの発想ですが、実際、$(16)$式による計算は$(15)$式よりもはるかに$p(\delta=0 | \mathcal{D}, M_1)$の推定精度が良いことが参考文献による実験で解明されています。その理由はラオ・ブラックウェルの定理により明らかとのこと。勉強できたらAppendixに載せたいと思います。

$(16)$式は、MCMCサンプルを用いて次のように近似的に計算されます。

$$ \mathrm{E}_ {\boldsymbol{\psi}}[p(\delta=0, |\boldsymbol{\psi}, \mathcal{D},M_1)] \eqsim \cfrac{1}{T} \sum_{t=1}^{T} p(\delta=0 | \boldsymbol{\psi}^{(t)}, \mathcal{D},M_1) \tag{17} $$

$(17)$式は、条件付き事後分布$p(\delta|\boldsymbol{\psi},\mathcal{D},M_1)$の関数形が明らかであれば解析的に計算可能です。そこで、JZSの事前分布、単位情報事前分布それぞれの場合の条件付き事後分布$p(\delta|\boldsymbol{\psi},\mathcal{D},M_1)$及びその導出過程を、対応のある場合とない場合に分けて提示します。

◆対応のある平均値の差の検定におけるCMDE

JZSの事前分布を設定したときの条件付き事後分布$p(\delta | \boldsymbol{\psi}=(\sigma,g), \mathcal{D}=(x^n,y^n), M_1)$は以下のとおりである。

$$ p(\delta | \boldsymbol{\psi},\mathcal{D}, M_1) = \mathrm{Normal}\left(\cfrac{\sum_{i=1}^{n}(x_i-y_i)}{\sigma\left(n + \cfrac{1}{g}\right)}, \cfrac{1}{n + \cfrac{1}{g}}\right) \tag{18} $$

単位情報事前分布を設定したときの条件付き事後分布$p(\delta | \boldsymbol{\psi}=\sigma, \mathcal{D}=(x^n,y^n), M_1)$は、

$$ p(\delta | \boldsymbol{\psi},\mathcal{D}, M_1) = \mathrm{Normal}\left(\cfrac{\sum_{i=1}^{n}(x_i-y_i)}{\sigma\left(n + \cfrac{1}{r^2}\right)}, \cfrac{1}{n + \cfrac{1}{r^2}}\right) \tag{19} $$

◆導出 (18)式を導出する。ベイズルールより、 $$ p(\delta | \mathcal{D},\boldsymbol{\psi}, M_1) \propto p(\mathcal{D} | \delta, \boldsymbol{\psi}, M_1 )p(\delta | \boldsymbol{\psi}) $$ $(3)$式、$(7)$式を用いるため、 $$ \begin{cases} p(\mathcal{D} | \delta, \boldsymbol{\psi}, M_1 ) = ∏_{i=1}^{n} p(x_i-y_i | \delta, \sigma) = ∏_{i=1}^{n} \mathrm{Normal} (x_i-y_i | \delta\sigma, \sigma^2 ) \\\\ p(\delta | \boldsymbol{\psi}) = p(\delta | g) = \mathrm{Normal}(\delta |0,g) \end{cases} $$ よって、 $$ \begin{split} p(\delta | \mathcal{D},\boldsymbol{\psi}, M_1) &=& ∏_{i=1}^{n} \mathrm{Normal} (x_i-y_i | \delta\sigma, \sigma^2 ) \times \mathrm{Normal}(\delta | 0,g) \\\\ &=& ∏_{i=1}^{n} \cfrac{1}{\sqrt{2\pi}\sigma}\exp\left(-\cfrac{((x_i - y_i) - \delta\sigma)^2}{2\sigma^2}\right) \times \cfrac{1}{\sqrt{2\pi g}}\exp\left(-\cfrac{\delta^2}{2g}\right) \\\\ &\propto& \exp \left( -\cfrac{n\sigma^2\delta^2 - 2\sum_{i=1}^{n}(x_i - y_i)\delta\sigma + \sum_{i=1}^{n}(x_i - y_i)^2} {2\sigma^2} - \cfrac{\delta^2}{2g} \right) \\\\ &\propto& \exp\left(-\cfrac{1}{2\sigma^2}n\sigma^2\left(\delta - \cfrac{1}{n\sigma}\sum_{i=1}^{n}(x_i - y_i)\right)^2 - \cfrac{\delta^2}{2g}\right) \\\\ &\propto& \exp\left( -\cfrac{1}{2} \left( \left(n + \cfrac{1}{g}\right)\delta^2 -2\cfrac{1}{\sigma}\sum_{i=1}^{n}(x_i - y_i)\delta \right) \right) \\\\ &\propto& \exp\left(-\cfrac{1}{2\cfrac{1}{n + \cfrac{1}{g}}}\left(\delta - \cfrac{\cfrac{1}{\sigma}\sum_{i=1}^{n}(x_i - y_i)}{n + \cfrac{1}{g}}\right)^2\right) \\\\ &\propto& \mathrm{Normal}\left(\cfrac{\sum_{i=1}^{n}(x_i - y_i)}{\sigma\left(n + \cfrac{1}{g}\right)}, \cfrac{1}{n + \cfrac{1}{g}}\right) \end{split} $$ $(19)$式は$(3)$、$(10)$式を用いるので、上述の導出で$g$を$r^2$に置き換えることで同様に導出できる。

◆対応のない平均値の差の検定におけるCMDE

JZSの事前分布を設定したときの条件付き事後分布$p(\delta | \boldsymbol{\psi}=(\sigma,g,\mu), \mathcal{D}=(x^n,y^m), M_1)$は以下のとおりである。

$$ p(\delta | \boldsymbol{\psi},\mathcal{D}, M_1) = \mathrm{Normal}\left(\cfrac{2g\left( \sum_{j=1}^{m}y_j - \sum_{i=1}^{n} x_i + (n-m)\mu \right)}{\left(4 + (n+m)g \right)\sigma }, \cfrac{4g}{4+(n+m)g}\right) \tag{20} $$

単位情報事前分布を設定したときの条件付き事後分布$p(\delta | \boldsymbol{\psi}=(\sigma, \mu), \mathcal{D}=(x^n,y^m), M_1)$は、

$$ p(\delta | \boldsymbol{\psi},\mathcal{D}, M_1) = \mathrm{Normal}\left(\cfrac{2r^2\left( \sum_{j=1}^{m}y_j - \sum_{i=1}^{n} x_i + (n-m)\mu \right)}{\left(4 + (n+m)r^2 \right)\sigma }, \cfrac{4r^2}{4+(n+m)r^2}\right) \tag{21} $$

◆導出 (20)式を導出する。ベイズルールより、 $$ p(\delta | \mathcal{D},\boldsymbol{\psi}, M_1) \propto p(\mathcal{D} | \delta, \boldsymbol{\psi}, M_1 )p(\delta | \boldsymbol{\psi}) $$ $(6)$式、$(7)$式を用いるため、 $$ \begin{cases} p(\mathcal{D} | \delta, \boldsymbol{\psi}, M_1 ) = ∏_{i=1}^{n} \mathrm{Normal}(x_i | \mu - \cfrac{\delta\sigma}{2}, \sigma^2) \times ∏_{j=1}^{m} \mathrm{Normal}(y_i | \mu + \cfrac{\delta\sigma}{2}, \sigma^2) \\\\ p(\delta | \boldsymbol{\psi}) = p(\delta | g) = \mathrm{Normal}(\delta |0,g) \end{cases} $$ よって、 $$ \begin{split} p(\delta | \mathcal{D},\boldsymbol{\psi}, M_1) &=& ∏_{i=1}^{n} \mathrm{Normal}(x_i | \mu - \cfrac{\delta\sigma}{2}, \sigma^2) \times ∏_{j=1}^{m} \mathrm{Normal}(y_i | \mu + \cfrac{\delta\sigma}{2}, \sigma^2) \times \mathrm{Normal}(\delta |0,g) \\\\ &=& \prod_{i=1}^{n}\cfrac{1}{\sqrt{2\pi}\sigma}\exp\left\{-\cfrac{\left( x_i - \left( \mu - \cfrac{\delta\sigma}{2} \right) \right)^2}{2\sigma^2}\right\} \times \prod_{j=1}^{m}\cfrac{1}{\sqrt{2\pi}\sigma}\exp\left\{-\cfrac{\left( y_j - \left( \mu + \cfrac{\delta\sigma}{2} \right) \right)^2}{2\sigma^2}\right\} \times \cfrac{1}{\sqrt{2\pi g}}\exp\left({\cfrac{\delta^2}{2g}}\right) \\\\ &\propto& \exp\left\{ -\cfrac{1}{2\sigma^2}\left( \delta\sigma\left( \sum_{i=1}^{n} x_i - \sum_{j=1}^{m} y_j \right) - (n-m)\mu\sigma\delta + \cfrac{n+m}{4}\sigma^2\delta^2 \right) - \cfrac{1}{2g}\delta^2 \right\} \\\\ &\propto& \exp\left\{ -\cfrac{1}{2} \left( \left( \cfrac{1}{g} + \cfrac{n+m}{4} \right) \delta^2 + \cfrac{\sum_{i=1}^{n}x_i - \sum_{j=1}^{m}y_j -(n-m)\mu}{\sigma}\delta \right) \right\} \\\\ &\propto& \exp \left\{ -\cfrac{1}{2\cfrac{4g}{4+(n+m)g}} \left( \delta - \cfrac{2g\left( \sum_{j=1}^{m}y_j - \sum_{i=1}^{n} x_i + (n-m)\mu \right)}{\left(4 + (n+m)g \right)\sigma } \right)^2 \right\} \\\\ &\propto& \mathrm{Normal}\left(\cfrac{2g\left( \sum_{j=1}^{m}y_j - \sum_{i=1}^{n} x_i + (n-m)\mu \right)}{\left(4 + (n+m)g \right)\sigma }, \cfrac{4g}{4+(n+m)g}\right) \end{split} $$ $(21)$式は$(6)$、$(10)$式を用いるので、上述の導出で$g$を$r^2$に置き換えることで同様に導出できる。

これでベイズファクターを使って平均値の差の検定をする準備が整いました。実践編は実践編で！

Appendix

Jaffrayの事前分布

◆Jeffreysの事前分布

Jeffreysの事前分布$p_J(\theta)$は、下記のとおり定義される。

$$ P_J(\theta) = \sqrt{I(\theta)} \tag{22} $$

ここで、$I(\theta)$はフィッシャー情報量とよばれ、$\theta$は（通常確率変数の）パラメータである。

$$ I(\theta) = \mathbb{E}\left[ \left. \left( \cfrac{\partial}{\partial \theta} \mathrm{ln} ~ L(\theta | x) \right)^2 \right| \theta \right] \tag{23} $$

この定義に沿って正規分布の平均・分散に対するJeffreysの事前分布を導出する。

正規分布の平均$\mu$に対するJeffreysの事前分布の導出 正規分布に従う確率変数$X$について考える。 $$ L(\mu | x) = \cfrac{1}{\sqrt{2\pi}\sigma}\exp\left\{ - \cfrac{(x - \mu)^2}{2\sigma^2} \right\} $$ より、 $$ \mathrm{ln} ~ L(\mu | x) = - \log(\sqrt{2\pi}\sigma) - \cfrac{(x - \mu)^2}{2\sigma^2} $$ だから、 $$ \cfrac{\partial}{\partial \mu} \mathrm{ln} ~ L(\mu | x) = \cfrac{(x - \mu)^2}{\sigma^2} $$ フィッシャー情報量$I(\mu)$は、 $$ I(\mu) = \mathbb{E}\left[ \left( \cfrac{x - \mu}{\sigma^2} \right)^2 \right] = \cfrac{1}{\sigma^2} \mathbb{E} \left[ \left( \cfrac{x -\mu}{\sigma} \right)^2 \right] = \cfrac{\alpha_{X,2}}{\sigma^2} $$ ここで、$\alpha_{X,2}$は確率変数$X$の2次の標準化積率である。標準化積率は $$ \alpha_{X,2} = \mathbb{E}\left[ \cfrac{X - \mu_{X}}{\sigma_X}^2 \right] = \cfrac{\mu_{X,2}}{\sigma_X^2} $$ と表現されること(ここで$\mu_{X,k}$は$X$の中心積率、$\sigma_X$は$X$の標準偏差)、$\mu_{X,2}$は$X$の分散であることから、 $$ I(\mu) = \cfrac{\alpha_{X,2}}{\sigma^2} = \cfrac{1}{\sigma^2} $$ よって、 $$ P_J(\mu) = \sqrt{I(\mu)} = \cfrac{1}{\sigma} \propto 1 $$

正規分布の分散$\sigma^2$に対するJeffreysの事前分布の導出 正規分布に従う確率変数$X$について考える。 $$ L(\sigma^2 | x) = \cfrac{1}{\sqrt{2\pi}\sigma}\exp\left\{ - \cfrac{(x - \mu)^2}{2\sigma^2} \right\} $$ ここで$\sigma^2 = g$とおいて、 $$ L(g | x) = \cfrac{1}{\sqrt{2\pi g}}\exp\left\{ - \cfrac{(x - \mu)^2}{2g} \right\} $$ このとき $$ \cfrac{\partial}{\partial g} \mathrm{ln} ~ L(g | x) = \cfrac{(x-\mu)^2 - g }{2g^2} $$ より、フィッシャー情報量$I(\sigma)$は、 $$ \begin{split} I(g) &=& \mathbb{E}\left[ \left( \cfrac{(x-\mu)^2 - g }{2g^2} \right)^2 \right] \\\\ &=& \cfrac{1}{4g^2}\mathbb{E}\left[ \cfrac{(x - \mu)^4}{g^2} \right] - 2\cfrac{1}{g^2}\mathbb{E} \left[ \cfrac{(x-u)^2}{g} \right] + \cfrac{1}{4g} \\\\ &=& \cfrac{1}{4g^2} \alpha_{X,4} - 2\cfrac{1}{g^2}\alpha_{X,2} + \cfrac{1}{4g} \\\\ &=& \cfrac{1}{2g^2} \end{split} $$ 途中、$\alpha_{X,4}=3$、$\alpha_{X,2}=1$を用いた。よって $$ P_J(\sigma^2) = P_J(g) = \sqrt{I(g)} = \cfrac{1}{\sqrt{2}\sigma} \propto \cfrac{1}{\sigma} $$

2021.06.21 00:00

ベイズファクターを用いた仮説検定～平均値の差の検定～