【Pythonで学ぶデータ分析】対応のあるデータの母平均に差があるかどうかをベイズt検定で調べる 〜 ホラー映画を観ると握力は上がるのか？：やさしい推測統計（ベイズ統計編）
手に汗握るホラー映画を観た後では、観る前よりも握力が強くなったような気がしませんか？ 同じ人の2回の測定値の差を求め、ベイズ統計により検定します。事前分布のパラメーターを変えても結果が安定するかどうかを調べる「感度分析」にも触れます。『社会人1年生から学ぶ、やさしいデータ分析』ベイズ統計編の第7回です。
母平均の差の検定をベイズ統計で行うには？ 連載『社会人1年生から学ぶ、やさしいデータ分析』のベイズ統計編、第7回となる今回は、古典的な「対応のあるデータのt検定」に代わるベイズt検定の方法を解説します。ホラー映画を観た後の握力と観る前の握力の差を求め、その値を基に分析します。従って、第3回で見た、1群の母平均のベイズ推定やベイズ検定と同じ方法が使えます。これまで、ベイズ検定では、主にSavage-Dickey法を利用してベイズ因子を求めていましたが、今回は周辺尤度の比によりベイズ因子を求める方法についても詳しく解説します。また、感度分析により、事前分布のパラメーターを変えても安定した結果が得られるかどうかを確かめます。
この連載では、簡単な事例を通してベイズ統計の考え方と分析の進め方を解説します。新しい用語や考え方が幾つも出てきますが、全てを理解しなくても大丈夫です。登場するたびに、分かりやすく丁寧に説明するので、その時点で分からないことがあっても気にせず、先に進んでください。回を追うごとに、少しずつ理解が深まります。どんな考え方で、どんな手順で分析を進めるのか、といった大きな流れを捉えてください。
連載：
この連載では、データをさまざまな角度から分析し、その背後にある有益な情報を取り出す方法を学ぶ『社会人1年生から学ぶ、やさしいデータ分析』シリーズの「記述統計と回帰分析編」「確率分布編」「推測統計（区間推定編・仮説検定編）」に続く「ベイズ統計編」です。
これまでの推測統計を土台に、近年活用が広がっているベイズ統計の考え方と分析手順を、古典的な手法との違いを整理しながら解説します。初めての方でも無理なく理解できるよう具体例を通して進めるとともに、ベイズ的なアプローチの特徴やメリットを実感できる構成とし、「どのように考え、どう使い分けるのか」に重点を置いて解説していきます。
筆者紹介： IT系ライターの傍ら、かつて、非常勤講師として東大で情報・プログラミング関連の授業を、一橋大でAI関連の授業を担当。まだ発展途上のようだが、生成AIへの質問と回答を分野ごとにまとめ、体系的に整理、補足しつつ、さらに個別の学びを促す記述を加えた資料を一貫して作成できるツールが実用レベルで本格的に使えるようになると、書籍の概念が変わるのではないかと思ったりしている。一斉授業だったものが、ユーザーに合った個別指導のできる書籍といったイメージ。そうなると、出版社のビジネスモデルも大きく変わり、私の仕事もなくなってしまいそうだけど、まあ、そのときは、縁側でお茶でもすすりながら、庭の景色を眺めて暮らすことにしようと思う。ウチには縁側も庭もないけど。
対応のあるデータの母平均の差をベイズt検定によって分析する
前回取り上げたベイズt検定は、独立した2群の母平均について、差があるかどうかを調べるというものでした。今回は、対応のある2群のデータを基に、母平均の差を分析します。対応のあるデータとは、同じ人がテストを2回受けたときの、1回目の値と2回目の値のようなデータです。図1をご覧ください。ホラー映画を観る前と観た後の握力のデータなので、この例も対応のあるデータであることが分かりますね。ただし、データは架空のものです。実際のところ、ホラー映画を観ると握力が上がるという研究はないようなので、図1の結果は「虚構」でしかありません。あくまで、手法を学ぶことが趣旨です（念のため）。
図1 対応のある2群の母平均の差をベイズt検定により確かめる
例えば、先頭の値の43.5と44.5は同じ人（1人目）の握力で、次の44.5と48.4も同じ人（2人目）の握力。つまり、これらのデータは対応のある2群のデータとなる。このような場合、2回の測定値の差を求め、その差について、1群の母平均の分析を行えばよい。ただし、今回は効果量に注目して分析を行う。さらに感度分析により、事前分布のパラメーターが異なっても安定した結果が得られるかどうかを調べる。
独立した2群の母平均の差を分析する場合は、文字通り母平均の差を分析しました。一方、対応のある2群の母平均の差の場合は、差の母平均を分析することになります。つまり、図1の中央右の「1.0, 3.9, -1.3, ... ,8.9」という1群について、母平均の推定や検定を行えばいいということです。なお、前提として、2群の差は正規分布しているものとします。
1群の母平均の分析については、すでにこの連載の第3回で取り扱いましたが、今回は効果量δに注目して分析を行うことにします。前回お話ししたように、効果量を使うと、検定の対象となる母数が1つだけで済む（正規分布の場合、母数がμとσの2つ）といったメリットや、データの単位が違っても統一的に分析できるというメリットがあります。さらに、今回の感度分析では、検定の対象となる母数が1つだけである方が簡単です。この例の場合、効果量δは以下の式で求められます。
つまり、標準偏差σと比較して、差の平均μが0からどれだけ離れているか、ということですね。では、コードを見ていきましょう。なお、ここから、ベイズ因子を求めるところまでは、ほぼ第3回と第6回で学んだ知識がそのまま使えます。前回までの記事を既に読まれた方には、同じような話でやや退屈かもしれませんが、どの部分が違っているかに注目して読み進めてください。その後で、今回取り上げる新たな内容（「周辺尤度の比によりベイズ因子を求める」「感度分析を行う」）をお話しします。前回までの話はよく分かった、という方は、以下を斜め読みで済ませて、「周辺尤度の比によりベイズ因子を求める」まで一気に進んでいただいても構いません。
事前分布と尤度関数を定義して、サンプリングを行う
最初にデータを確認しておきます。ホラー映画を観る前の握力をx1、観た後の握力をx2とします。これらの差をそのまま使っても構わないのですが、標準偏差が1になるようにスケールを調整しておきましょう。そのためのコードが以下のリスト1です。コードはこちらのサンプルファイルに含めてあります。リンクを開くと、説明やPythonのプログラムが表示されるので、最初のコードセルをクリックし、［Shift］＋［Enter］キーを押して、実行しておいてください。
import numpy as np
x1 = [43.5, 44.5, 50.7, 39.0, 38.4, 38.4, 51.1, 45.4, 36.7, 43.8,
36.8, 36.7, 41.7, 26.6, 27.9, 36.1, 32.9, 42.2, 33.6, 30.1,
50.3, 38.4, 40.5, 30.0, 36.2, 40.8, 31.9, 42.6, 35.8, 38.0]
x2 = [44.5, 48.4, 49.4, 52.3, 46.5, 36.3, 56.1, 39.6, 34.1, 48.8,
44.5, 41.6, 45.1, 29.1, 24.5, 36.5, 34.6, 51.5, 39.3, 25.3,
55.9, 40.5, 41.1, 37.1, 45.4, 49.5, 31.7, 45.1, 41.5, 46.9]
# スケールを調整する
d = np.subtract(x2, x1) # 差を求める
x = d / d.std(ddof=1) # 標準偏差で割る
x2とx1の差を、stdメソッドで求めた標準偏差で割れば、標準偏差が1にそろえられる。引数のddof=1は不偏標準偏差を求めるという指定。単にスケールをそろえるだけであれば、標本標準偏差（ddof=1を指定しない）を使ってもいいが、ここでは、Cohen's dとの対応を重視して、母標準偏差の推定値である不偏標準偏差を使っている。以降の分析はこのデータ（x）を使って行うことになる。なお、ここでのスケール調整は厳密な「標準化」とは異なり、差の平均（d.mean()）を引いて、平均を0にそろえてはいけない（基準値である0との差を調べることができなくなる）。
続いて、事前分布と尤度関数を定義していきましょう。全体像は図2のようになります。
図2 事前分布と尤度関数の考え方
注目しているパラメーターは効果量δ。その事前分布についてはコーシー分布とする。また、差の標準偏差σの分布についても定義しておく。これらは、尤度関数のパラメーターμを求めるために使う。μは、効果量δと差の標準偏差σから求めることができる。計算方法については後述する。
ちょっとしたおさらいですが、尤度関数とは「データを観測したときにパラメーターのさまざまな値がどの程度尤（もっと）もらしいかを表す関数」ですね。ここでは（図2）では、その関数が正規分布Normal(μ,σ)の形になっているというわけです。
事前分布を定義する
まず、事前分布の定義からです。効果量δの事前分布としては、中心αが0、幅を表すパラメーターβが0.707のコーシー分布とします（一般的に使われる設定です）。差の標準偏差σについては、元の正規分布の標準偏差を1とした半正規分布とします。δとσは上でも少し触れたように、尤度関数のパラメーターμを求めるために使います。ここまでをまとめると、以下のようになります。
- 効果量δ： α＝0、β＝1/√2＝0.707のコーシー分布
- 差の標準偏差σ： 元の正規分布の標準偏差が1である半正規分布
これらの定義を書くと以下のようになります（リスト2）。
import pymc as pm
with pm.Model() as model:
# 事前分布
delta = pm.Cauchy("delta", alpha=0, beta=0.707) # 効果量
sigma = pm.HalfNormal("sigma", sigma=1) # 標準偏差
ここまで見てきた事前分布の定義を素直に書いたもの。コーシー分布はpymcモジュールのCauchyクラスを利用して定義する。引数alphaに分布の中心を、引数betaに分布の幅を表す値を指定すればよい。引数betaには0.707という値が使われることが多い。標準偏差は正の値しか取らないので半正規分布（HalfNormalクラス）とした。
「事前分布を定義する」と何気なく表現していますが、正確には「ある事前分布に従う確率変数を定義する」ということです。例えば、delta = pm.Cauchy("delta", ...)というコードであれば、コーシーｰ分布に従う確率変数を定義していることになります。モデル内で利用される確率変数の名前が第1引数に指定した"delta"で、その確率変数を表すオブジェクトが返り値として返され、代入されたdeltaという変数名で参照できるようになります。
尤度関数を定義する
続いて、尤度関数の定義です。データの母集団が正規分布しているものと想定される場合、尤度関数は正規分布となります（その根拠はちょっと難しいので、ここでは割愛しますが、数学的に導き出されます）。パラメーターは差の平均μと差の標準偏差σです。μについては、効果量δと差の標準偏差σから計算して求めます。図3に示したような方法です。
新たな変数を計算によって定義するので、pymcモジュールのDeterministic関数を使って記述します（リスト3）。
mu = pm.Deterministic("mu", delta * sigma)
図3の右下の式をそのまま記述すればよい。muは効果量（delta）×標準偏差（sigma）で求められる。
尤度関数は簡単です。上で作成したmuを使って、以下のように記述します。
pm.Normal("obs", mu=mu, sigma=sigma, observed=x)
握力の差の分布は、平均mu、標準偏差sigmaの正規分布に従う。観測されたデータ（observed引数）には、スケールを調整したデータxを指定する。
サンプリングを行う
サンプリングのためのコードは前回までに見た例とほぼ同じです。つまり、pymcモジュールのsample関数を使うだけです。これまで断片的に記したコードも含め、結果の要約を表示するところまで、コードの全体を掲載します。サンプルファイルの2番目のコードセルに入力されています。
import pymc as pm
import arviz as az
# データは最初のセルで入力されたものを使う
with pm.Model() as model:
# 事前分布
delta = pm.Cauchy("delta", alpha=0, beta=0.707) # 効果量
sigma = pm.HalfNormal("sigma", sigma=1) # 標準偏差
mu = pm.Deterministic("mu", delta * sigma) # 平均（効果量と標準偏差から逆算。基準値を0と仮定）
# 尤度関数：平均mu、標準偏差sigmaの正規分布に従うものとする
pm.Normal("obs", mu=mu, sigma=sigma, observed=x)
# サンプリング（後でベイズ因子を求めるので、少し多めにサンプリングする）
trace = pm.sample(draws=8000, chains=2, random_seed=42)
# 事後分布の要約
az.summary(trace, var_names=["delta", "mu", "sigma"])
# サンプリングというコメント以降が、これまで断片的に記したコードの続き。後でベイズ因子を求めるために、やや多めに（8000個の）サンプリングを行った。事後分布の要約として、効果量deltaと握力の差の平均muとその標準偏差sigmaの推定値などを表示する。なお、この平均はスケール調整された値であることに注意。
pymc.sample関数によってサンプリングが実行されると、各変数のサンプルが自動的に作られます。上のコードでは、それらをarvizモジュールのsummary関数で要約しています。実行すると、図4のような結果が表示されます。
図4 事後分布の要約
効果量δの点推定値は、deltaの行のmean列の値、つまり0.682となる。94％信用区間はhdi_3％〜hdi_97％なので、0.314〜1.076。0以下の値が含まれないので、効果量は正、つまりx2>x1であることが示唆される。r_hatが1.0である（＝収束している）ことも確認しておこう。
ベイズ因子を求める 〜 対応のある2群のベイズt検定を行う
続けましょう。まず、Savage-Dickey法を使ってベイズ因子を求めてみます。事前分布（コーシー分布）での、0に対する確率密度と、事後分布をスムーズな曲線にした場合の0に対する確率密度との比を取って、ベイズ因子BF10の値を求めればいいですね。帰無仮説と対立仮説を確認した後、コードを見ていきます。
- 帰無仮説H0 : δ＝0
- 対立仮説H1 : δ ≠ 0
以下のコードはサンプルファイルの3番目のコードセルに入力されています。これも前回までのコードそのままです。
from scipy.stats import gaussian_kde, cauchy
# 事後分布
delta_posterior = trace.posterior['delta'].values.flatten()
kde = gaussian_kde(delta_posterior) # 滑らかな曲線にする
posterior = kde.evaluate(0)[0] # δ=0の位置の確率密度関数の値
# 事前分布（コーシー分布）
prior = cauchy.pdf(0, loc=0, scale=0.707) # δ=0の位置の確率密度関数の値
# ベイズ因子
bf10 = prior / posterior
print(f"BF10: {bf10:.3f}")
# 出力例：
# BF10: 70.702
事後分布では、scipy.statsモジュールのgaussian_kdeクラスを使って、サンプリングされたデータのヒストグラムを滑らかな曲線にする。事後分布と事前分布のδ＝0に対する確率密度関数の値を求め、その比を取ればベイズ因子が求められる。結果は70.702となる。
上の結果をJeffreys（1961）の評価尺度に照らし合わせると、H1を支持する「非常に強い証拠となっている」ことが分かります。従って、対立仮説δ ≠ 0（効果量は0ではない）が支持されます。効果量が正なので、実質的にμは0よりも大きいと言えます。つまり、ホラー映画を観た後では握力は上がるということになります（くどいようですが、この結果は架空のデータによるものなので、実際にそうなるわけではありません）。
対立仮説が正しい確率を求める
さらに、以下の式を使って、対立仮説H1が正しい確率を求めましょう。この式を導き出す方法は前回紹介しました。
となります。リスト6で求めたBF10＝70.702を代入すると、
という極めて高い確率になります。一方、帰無仮説H0が正しい確率は、
となります。
周辺尤度の比によりベイズ因子を求める
前回までは、ベイズ因子を主にSavage-Dickey法で求めてきました。Savage-Dickey法では、ある点における事前分布の確率密度関数の値と事後分布の確率密度関数の値の比を取ってベイズ因子としました。一方、前回のコラムでは、ベイズ因子が周辺尤度の比であることを紹介しました。実は、Savage-Dickey法によって求めた値は周辺尤度の比と等しくなります。そこで、周辺尤度の比を求めることによって、ベイズ因子を計算してみましょう。前回までは、サンプルプログラム中に参考としてそのためのコードを含めていただけですが、モヤモヤしていた方もおられるでしょうから、ここで解説することにします。
といっても、考え方そのものは簡単です。H0によるモデルとH1によるモデルでサンプリングを行い、周辺尤度の比を求めるだけです。しかし、残念なことに、pymc.sample関数では周辺尤度を求めることができません。そこで、逐次モンテカルロ法（SMC：Sequential Monte Carlo）と呼ばれる方法によりサンプリングを行うpymc.sample_smc関数を使います。pymc.sample_smc関数では、周辺尤度を求めることができます（ただし、周辺尤度の対数を取った対数周辺尤度がチェーンごとに得られます）。
マルコフ連鎖モンテカルロ法（MCMC）は、1つの値を動かしながら、山を上り下りするかのように事後分布を求めるイメージです。
一方の、逐次モンテカルロ法（SMC）は、多数の値を「ばらまいて」、確率密度の高い値を残しながら山の形に近づけていく（事後分布に近づけていく）イメージです。
引数の指定方法はほぼ同じなので、早速コードを見てみましょう。帰無仮説H0のモデル（model0）では効果量deltaを0としています。一方、対立仮説H1のモデル（model1）では、効果量deltaがコーシー分布に従うものとしています。後はリスト5とほとんど同じです。その点をまず確認しておき、以下のコードをざっと眺めてみてください。詳細については、コードの後で説明します。
import pymc as pm
import numpy as np
from scipy.special import logsumexp
with pm.Model() as model0: # 帰無仮説（delta = 0）
# 事前分布
delta = 0 # 効果量は0 ← ここに注目
sigma = pm.HalfNormal("sigma", sigma=1)
mu = pm.Deterministic("mu", delta * sigma) # 結局mu=0となるので、尤度関数の引数にmu=0を指定しても同じ
# 尤度関数
pm.Normal("x_obs", mu=mu, sigma=sigma, observed=x) # 平均は0とする
# SMCサンプリング ↓ ここに注目
idata0 = pm.sample_smc(draws=8000, chains=2, threshold=0.8, random_seed=42, progressbar=False)
with pm.Model() as model1: # 対立仮説（delta != 0）
# 事前分布
delta = pm.Cauchy("delta", alpha=0, beta=0.707) # 効果量はコーシー分布に従う ← ここに注目
sigma = pm.HalfNormal("sigma", sigma=1)
mu = pm.Deterministic("mu", delta * sigma)
# 尤度関数
pm.Normal("x_obs", mu=mu, sigma=sigma, observed=x)
# SMCサンプリング ↓ ここに注目
idata1 = pm.sample_smc(draws=8000, chains=2, threshold=0.8, random_seed=42, progressbar=False)
# H0の対数周辺尤度とH1の対数周辺尤度（チェーン数の個数だけ求められる）
log_m0_array = [v[-1] for v in idata0.sample_stats.log_marginal_likelihood.values] # H0
log_m1_array = [v[-1] for v in idata1.sample_stats.log_marginal_likelihood.values] # H1
# ベイズ因子は以下のコードでも求められるが、オーバーフローなどのトラブルを避けるために
# 後のコードを使う（コメントを外して実行すれば同じ結果となる）
# m0 = np.mean(np.exp(log_m0_array)) # 対数周辺尤度を指数関数に指定し、平均を求める
# m1 = np.mean(np.exp(log_m1_array))
# print(f"BF10: {m1 / m0:.3f}")
# 周辺尤度の平均の対数
log_m0 = logsumexp(log_m0_array) - np.log(len(log_m0_array))
log_m1 = logsumexp(log_m1_array) - np.log(len(log_m1_array))
# ベイズ因子
log_bf10 = log_m1 - log_m0
bf10 = np.exp(log_bf10)
print(f"BF10: {bf10:.3f}")
# 出力例：
# BF10: 78.417
pymc.sample_smc関数を利用してサンプリングを行い、H0の場合の周辺尤度とH1の場合の周辺尤度を求める。それらの比を取ればベイズ因子BF10の値が得られる。リスト6とはBF10の値が異なるが、いずれも乱数を使ったシミュレーションのため（傾向に大きな違いはない）。
まずは、リスト7（4番目のコードセル）を実行して、結果が求められることを確認しておいてください。続いて、コードについて説明します。
これまでと異なる点は、もちろん、pymc.sample関数ではなく、pymc.sample_smc関数を使ってサンプリングを行っている点です。帰無仮説のモデル（model0）について、その部分を抜き出してみます。
idata0 = pm.sample_smc(draws=8000, chains=2, threshold=0.8, random_seed=42, progressbar=False)
引数の指定方法はリスト5とほぼ同じです。違いはthreshold=0.8とprogressbar=Falseという指定を追加してあるところだけです。thresholdは、事後分布に近づけていく「歩幅」をどれぐらい小さくするかの基準といったイメージの引数です。thresholdの値が大きいほど、歩幅が小さくなり、精度が上がります。実は、仕様なのかバグなのか、原因は不明ですが、thresholdの値を既定値の0.5のままにしたり、極端な値を指定したりすると、後述する対数周辺尤度の配列が期待した形にならない場合があります。そこで、対症療法ではありますが、既定値よりも精度を上げた設定にしてあります。progressbar=Falseについては、サンプリングの進行状況を表示しないという指定です。処理そのものとは関係ありませんが、できるだけ速く結果を表示したいので省略しただけです。
pymc.sample_smc関数にtune引数がないことにも気付いたでしょうか。SMCではウォームアップ（バーンイン）が不要なので、tune引数がありません（指定するとエラーになります）。
この、pymc.sample_smc関数の返り値から、各チェーンで求められた対数周辺尤度が得られます。それらを基に平均の周辺尤度を求めればいいのですが、データの形式が複雑なので、以降のコードもやや長くなってしまいます。ただ、コードの書き方は本質的なお話ではないので、Pythonのプログラミングにあまり慣れていない方は考え方だけ理解していただければ十分です。
得られた対数周辺尤度の値を取り出すに当たって、どのような値が返されているかを確認用のコードを使って見ておきましょう。こちらも帰無仮説のモデルから得られた結果（idata0）の例で見てみます。リスト7の次（5番目のコードセル）に以下のようなコードが入力されているので、実行してみてください。
print("idata0の対数周辺尤度")
print(idata0.sample_stats.log_marginal_likelihood.values)
# 出力例：
# idata0の対数周辺尤度
# [[nan nan nan nan np.float64(-50.56121127475562)]
# [nan nan nan nan np.float64(-50.57304795830348)]]
idata0.sample_stats.log_marginal_likelihood.valuesに各チェーンで求められた対数周辺尤度が記録されている。各行の最後の列がその値となっている。
出力例を見ると、チェーン数と同じ行数のリストとなっています。対数周辺尤度は逐次更新されていくのですが、最終的な結果だけが最後の列に記録されており、それまでの部分（途中経過）にはnanという値が入れられています。このデータから、全てのチェーンの最後にある値を取り出し、それらの対数周辺尤度の値を基に周辺尤度の平均を求めればいい、ということです。
まず、各チェーンの対数周辺尤度の値だけを取り出します。それぞれのチェーンの最後の要素だけ取り出して並べるだけです。最後の要素はインデックスとして-1を指定すれば取り出せます。
# 最後の要素だけを取り出す
log_m0_array = [v[-1] for v in idata0.sample_stats.log_marginal_likelihood.values]
print(log_m0_array) # 取り出された対数周辺尤度のリストを表示する
# 出力例：
# [np.float64(-50.56121127475562), np.float64(-50.57304795830348)]
繰り返し処理によって、idata0.sample_stats.log_marginal_likelihood.valuesに含まれる各行の最後の値を取り出せばよい。コードを簡潔にするため、リスト内包表記を利用して記述している。
これで、各チェーンの対数周辺尤度が求められました。ここから、周辺尤度の平均を求めるのですが、単純に対数周辺尤度を平均して指数関数に指定してもうまくいきません。対数周辺尤度を指数関数に指定してから平均する必要があります。例えば、得られた対数周辺尤度が2つで、それらをlog A, log Bとしたとき、
ではなく、
を求める必要があります。numpyのexp関数を使えば指数関数の値が求められるので、
import numpy as np
m0 = np.mean(np.exp(log_m0_array))
のnp.exp(log_m0_array)で、（2）式の分子の各項、つまり、exp(log A)とexp(log B)に当たる値が求められます。それをnp.mean関数で平均すれば、（2）式のような平均が求められるというわけです。
idata1についても同様に表せるので、ベイズ因子BF10は以下のコードで求められます。
m0 = np.mean(np.exp(log_m0_array))
m1 = np.mean(np.exp(log_m1_array))
print(f"BF10: {m1 / m0:.3f}")
ただ、このまま計算すると、np.exp関数の引数に大きな値や小さな値が与えられた場合にオーバーフローやアンダーフローなどのトラブルが起こる可能性があるので、リスト7では（2）式を変形したコードを使っています。それについては、ますます本筋から外れるので、（2）式で計算しなければならない理由と併せて、コラムにまとめておきます。ここまでのお話を踏まえて、もう一度リスト7を一通り見直していただくといいかと思います。
コラム 対数周辺尤度の配列を基に、周辺尤度の平均を求める方法
まず、（1）式ではなく（2）式で周辺尤度の平均を求めなければならない理由から見ていきます。厳密な証明ではありませんが、簡単な例を使って確認しましょう。
例えば、周辺尤度がA＝2とB＝3であるものとしましょう。その平均はもちろん2.5です。対数周辺尤度は、log(周辺尤度)という値です。この例であれば、対数周辺尤度はlog A ＝ log 2≈ 0.69とlog B ＝log 3≈ 1.10となります。もし、（1）式に従って計算してしまうと、
となり、正しい結果になりません。（2）式だと、
となり、正しい結果が得られます。要するに、元のAとBの値を求めて平均するということです。
続いて、（2）式を変形して平均の周辺尤度を求める方法です。リスト7で使っているlog_m0 = logsumexp(log_m0_array) - np.log(len(log_m0_array))というコードがちょっと謎ですね。このコードは周辺尤度の平均の対数を求めるものです。2つのチェーンで求められた対数周辺尤度がlog A, log Bであるものとします。
まず、周辺尤度の平均を表す（2）式の対数を取ります。
この式は、周辺尤度の平均の対数を表しています。ところで、対数の場合、割り算は引き算にできます。つまり、以下の公式が使えます。
この公式を（3）式に適用すると、
となります。式の下に書いたコードは、リスト7との対応を示したものです。scipy.special.logsumexp関数は、指数関数の値を求め、それらを合計し、その対数を取る関数です。引数には対数の配列を指定します。つまり、枠で囲んだ値の配列を指定するということです。log 2の「2」はチェーンの個数です。というわけで、（5）式に沿って書いた以下のコードで、帰無仮説H0の場合の周辺尤度の平均の対数が求められます。
log_m0 = logsumexp(log_m0_array) - np.log(len(log_m0_array))
さらに、log_m0を指数関数に指定し、np.exp(log_m0)とすれば、周辺尤度の平均が求められますが、まだ対数のままにしておきます。対立仮説H1についても同様に、以下のコードで周辺尤度の平均の対数が求められます。
log_m1 = logsumexp(log_m1_array) - np.log(len(log_m1_array))
ここから、ベイズ因子を求める計算も見ていきましょう。H0に対する周辺尤度をM0とし、H1に対する周辺尤度をM1として、（4）の公式をもう一度利用します。XにM1を当てはめ、YにM0を当てはめて、周辺尤度の比の対数を取ってみます。
両辺を指数関数に指定してみましょう。
exp(log X) ＝ Xなので、
となります。左辺は周辺尤度の比なのでベイズ因子ですね。というわけで、左辺の式をコードとして書けば、以下のようにベイズ因子が求められます。これがリスト7で使ったコードです。log M1がlog_m1に当たり、log M0がlog_m0に当たります（念のため再確認です）。
log_bf10 = log_m1 - log_m0
bf10 = np.exp(log_bf10)
print(f"BF10: {bf10:.3f}")
なお、サンプルファイルには、σの事前分布としてJeffreys事前分布p(σ) ∝ 1/σを使い、数値積分により周辺尤度を求めてベイズ因子を計算した例を含めてあります。これについては、第5回で少し解説しました。サンプルファイルにも、数値積分の方法と併せて簡単な説明を掲載しています。ここで詳しく解説したいのはヤマヤマなのですが、今のところPyMCを使って分析することを中心にお話ししているので、今回は残念ながら割愛します（いずれ機会を見て紹介したいと思います）。結果は74.42で、第5回で紹介した神ツールJASPで求められる結果と同じ値になります。なお、σの事前分布が異なるモデルなので、リスト6やリスト7の結果と直接比較できる値ではありません。
事前分布のパラメーターを変えても安定した結果が得られるかを調べる 〜 感度分析を行う
事前分布については、何らかの根拠を基にパラメーターを決めていました。しかし、その根拠があいまいな場合もあります。今回の例では、効果量の事前分布としてコーシー分布を利用し、広がりを表すパラメーターβの値をよく使われる0.707としました。しかし、そのように決め打ちするのではなく、より幅の狭い値（例えば0.5）や幅の広い値（例えば1.5）を指定した場合にも一貫してベイズ因子BF10の値が大きければ、対立仮説を支持する結果の信頼性が高まります。このように、パラメーターの値を変えてどのような結果が得られるかを調べることを感度分析と呼びます。
感度分析の方法は、パラメーターの値を変えてベイズ因子BF10の値を幾つか求めるだけです。ただし、サンプリングには時間がかかるので、パラメーターの値を数多く指定するのはあまり現実的ではありません。リスト5をベースに、コーシー分布のパラメーターをbeta=0.5、beta=0.707、beta=1.5に変えてBF10を求めた結果は以下のようになります（これらはサンプリングによる推定値なので、図5の数値積分による値とは異なりますが、傾向は同じです）。
- beta=0.5の場合： BF10＝63.204
- beta=0.707の場合： BF10＝70.702
- beta=1.5の場合： BF10＝38.071
Jeffreysの評価尺度によると、いずれも対立仮説を支持する「非常に強い証拠となっている」ことが分かります。なお、上の値は手作業でパラメーターの値を変えて求めても大した手間ではありませんが、指定されたコーシー分布のパラメーターを用いてサンプリングを行い、事後分布の（0に対する）確率密度関数の値を求めるような関数を作成しておけば、一気に結果が得られます。サンプルファイルにはそのためのコードも含めておきました。
さらに、コーシー分布のパラメーターを変えながら、数値積分による方法でベイズ因子を求め、感度分析を行った例もサンプルファイルに含めてあります。数値積分の場合はサンプリングを行わないので、高速かつ正確に計算ができます。そこで、パラメーターの値を0.01から1.5まで50個に刻んでベイズ因子BF10を求め、グラフ化してみました。結果が見やすいので以下に掲載しておきます。
「えーっ！ そのコードこそぜひ説明してくださいよ」という声が聞こえてきそうですが、数値積分により周辺尤度を求めてベイズ因子を計算するコードを関数にして、コーシー分布のパラメーターを変えながらベイズ因子を順に求めてグラフにしているだけです。上でもお話ししましたが、Jeffreys事前分布を利用し、ベイズ因子を数値積分によって求める方法については、数値積分の方法も含めていずれ詳しく解説したいと思います。いましばらくお待ちください。
図5 数値積分による方法で感度分析を行った場合のグラフ
横軸rはコーシー分布の幅を表すパラメーター（betaに当たるもの）の値を0.01から1.5まで50個に分けた値。縦軸はそれらの値に対するベイズ因子。縦軸が対数目盛りになっていることに注意。対数目盛りとしたのは、Jeffreys（1961）の評価尺度の各レベルが等しい間隔になるようにするため。なお、JASPでもこれとほぼ同じグラフが作成される。
図5を見ると、コーシー分布の幅が変わっても、BF10の値が広い範囲にわたって10以上の値になっていることが分かります。つまり、一貫して「強い証拠となっている」という評価（以上）であることが視覚的に確認できます。
今回は、効果量に注目して、対応のある2群の平均に差があるかどうかをベイズt検定により分析しました。これまで見てきたSavage-Dickey法だけでなく、逐次モンテカルロ法（SMC）を利用して周辺尤度を求め、その比を求めてベイズ因子を計算する方法も解説しました。また、感度分析の考え方と簡単な例も紹介しました。
さて、これまで、平均についての分析を中心に見てきましたが、次回は相関係数のベイズ検定に取り組みます。事例としては、年齢と原付事故の死傷者数に関係があるかどうかをベイズ的に調べてみます。次回もお楽しみに！
今回使った新出のクラス・関数・メソッド（主なもの）
引数も主なものだけを掲載します。引数に「=値」と書かれたもの（例：=True）は既定値を表します。
sample_smc： SMCサンプリングを行う関数
- モジュール： pymc（定義はpymc.smc.sample_smcだが、pymc.sample_smcで呼び出せる）
- 形式： sample_smc(draws=2000, chains=None, cores=None, threshold=0.5, random_seed=None, progressbar=True)
- 引数：
- draws： 最終的に取り出すサンプルの個数
- chains： サンプルのチェーン数。None（既定値）の場合は、引数coresと2の大きい方が自動的に設定される
- cores： 並行して実行されるチェーンの数。None（既定値）の場合はCPUのコア数が指定されたものと見なされる
- threshold： 事後分布を求めるための「歩幅」に当たるβ値（コーシー分布のbetaとは別のもの）をどの程度変化させるかというパラメーター。値が大きくなると少しずつ変化させる（精度が上がる）。0〜1の範囲で指定する
- random_seed： 乱数の初期設定。値を指定しない場合は毎回異なる乱数となる
- progressbar： サンプリングの進行状況を表示する場合はTrueを、表示しない場合はFalseを指定する
- 返り値： サンプリングの履歴（サンプリングされた値や診断結果など）。通常、arvizによって分析したり要約したりする
※SMCサンプリングでは、ウォームアップ（バーンイン）が不要なので、tune引数はありません。
logsumexp： 指数関数の合計の対数を求める関数
- モジュール： scipy.special
- 形式： logsumexp(a)
- 引数：
- a： 合計を求めたい値。通常、リストや配列を指定する
- 返り値： np.log(np.sum(np.exp(a)))で求められる値
Copyright© Digital Advantage Corp. All Rights Reserved.