【Pythonで学ぶデータ分析】相関があるかどうかをベイズ統計で調べる 〜 年齢と原付事故死傷者数に関係はあるのか？：やさしい推測統計（ベイズ統計編）
バイクに乗る若い人の運転はどうも危なっかしく感じられることがあります。そこで、手軽に乗れる原付一種（かつての排気量50ccまでのバイク）を対象に、年齢と交通事故死傷者数との間に関係があるかどうかを調べてみたいと思います。つまり、相関のベイズ検定に取り組もうというわけです。『社会人1年生から学ぶやさしいデータ分析』ベイズ統計編の第8回です。
相関の検定をベイズ統計で行うには？ 連載『社会人1年生から学ぶ、やさしいデータ分析』のベイズ統計編、第8回となる今回は、古典的な「無相関の検定」に代わるベイズ相関検定の方法を解説します。事前分布として、台を−1〜1に拡張したベータ分布（以降、拡張ベータ分布と呼びます）を使い、尤度（ゆうど）関数として多変量正規分布を使います。また、LKJ事前分布を使って、拡張ベータ分布のコードを簡潔に書く例も紹介します。といっても、分析の進め方はこれまでとほとんど同じです。パラメーターとして指定する分散・共分散行列の形式や意味などを確認しながら、ゆっくりと見ていきましょう。
この連載では、簡単な事例を通してベイズ統計の考え方と分析の進め方を解説します。新しい用語や考え方が幾つも出てきますが、全てを理解しなくても大丈夫です。登場するたびに、分かりやすく丁寧に説明するので、その時点で分からないことがあっても気にせず、先に進んでください。回を追うごとに、少しずつ理解が深まります。どんな考え方で、どんな手順で分析を進めるのか、といった大きな流れを捉えてください。
連載：
この連載では、データをさまざまな角度から分析し、その背後にある有益な情報を取り出す方法を学ぶ『社会人1年生から学ぶ、やさしいデータ分析』シリーズの「記述統計と回帰分析編」「確率分布編」「推測統計（区間推定編・仮説検定編）」に続く「ベイズ統計編」です。
これまでの推測統計を土台に、近年活用が広がっているベイズ統計の考え方と分析手順を、古典的な手法との違いを整理しながら解説します。初めての方でも無理なく理解できるよう具体例を通して進めるとともに、ベイズ的なアプローチの特徴やメリットを実感できる構成とし、「どのように考え、どう使い分けるのか」に重点を置いて解説していきます。
筆者紹介： IT系ライターの傍ら、かつて、非常勤講師として東大で情報・プログラミング関連の授業を、一橋大でAI関連の授業を担当。まだ発展途上のようだが、生成AIへの質問と回答を分野ごとにまとめ、体系的に整理、補足しつつ、さらに個別の学びを促す記述を加えた資料を一貫して作成できるツールが実用レベルで本格的に使えるようになると、書籍の概念が変わるのではないかと思ったりしている。一斉授業だったものが、ユーザーに合った個別指導のできる書籍といったイメージ。そうなると、出版社のビジネスモデルも大きく変わり、私の仕事もなくなってしまいそうだけど、まあ、そのときは、縁側でお茶でもすすりながら、庭の景色を眺めて暮らすことにしようと思う。ウチには縁側も庭もないけど。
相関係数をベイズ相関検定によって分析する
これまでは、平均の差のベイズt検定など、主に平均についての分析を行ってきました。今回からは関係に注目して分析を進めていきたいと思います。図1に示したように、年齢と原付事故死傷者数に関係があるかどうかを調べます。若い人ほど事故を起こしやすいという先入観が正しいかどうか、データを集めて調べてみようというわけです。年齢が上がるほど、事故が少なくなるという関係であれば、負の相関となるはずです。
図1 相関があるかどうかをベイズ相関検定により確かめる
年齢層によって人口が異なるだけでなく、原付を運転する人の割合も異なるので、単純に年齢と死傷者数の関係を調べても適切に分析できない。そこで、年齢別の原付購入者数当たりの死傷者数との関係を調べることにする。結果を見ると、やはり若い人の方が事故で死傷しやすいことが分かる。その結果を得るまでの分析方法を見ていく。なお、年齢と死傷者数の相関は疑似相関であるという可能性もあることに注意（例えば、走行距離と死傷者数に正の相関があり、年齢とともに走行距離が減るのかもしれない）。今回の趣旨は手法の理解なので、このまま進める。
年齢別原付事故死傷者数は、警察庁の交通事故統計で公表されている2024年のデータを基にしています。このデータはe-Statで提供されている「一般原付自転車乗車中の年齢層別死傷者数の推移」のリンクからダウンロードできます。また、年齢別の原付購入者数は、日本自動車工業会の2025年度二輪車市場動向調査報告書を基にしています。
分析に使うデータを用意する
年齢別原付事故死傷者数については、年齢が5歳刻みの階級に分けられているのですが、原付の購入者数については10歳刻みの階級になっています。そこで、線形補間を行い、5歳刻みで予測した購入者数を求めたデータを使うことにします。ただ、そのための前処理については、本筋の話から外れるので、ここでは分析に使う前処理済みのデータだけを以下に示しておきます。分析には標準化された値を使います。
| 年齢 | 死傷者数 （原付購入者数当たり） |
標準化された年齢 | 標準化された死傷者数 （原付購入者数当たり） |
|---|---|---|---|
| 17.0 | 38.0 | -1.62 | 1.86 |
| 22.0 | 44.9 | -1.39 | 2.39 |
| 27.0 | 29.0 | -1.16 | 1.18 |
| 32.0 | 19.3 | -0.93 | 0.45 |
| 37.0 | 15.0 | -0.69 | 0.12 |
| 42.0 | 10.3 | -0.46 | -0.24 |
| 47.0 | 9.6 | -0.23 | -0.29 |
| 52.0 | 10.3 | 0.00 | -0.24 |
| 57.0 | 7.2 | 0.23 | -0.47 |
| 62.0 | 4.2 | 0.46 | -0.70 |
| 67.0 | 2.8 | 0.69 | -0.80 |
| 72.0 | 3.5 | 0.93 | -0.75 |
| 77.0 | 3.7 | 1.16 | -0.74 |
| 82.0 | 2.6 | 1.39 | -0.82 |
| 87.0 | 1.0 | 1.62 | -0.94 |
年齢は15〜19、20〜24、……といった階級に分けられているので、その中央値を階級の代表値として記してある。85歳以上については、階級の代表値を便宜的に87.0とした。この表では小数点以下を1桁または2桁で丸めた値を示してある。
サンプルファイルには、元のデータに対して前処理を行い、標準化されたデータ（上の表の右2列）を返すmakedata関数を作成してあります。まず、分析の準備としてmakedata関数を実行してデータを作っておいてください。リンクを開くと、説明やPythonのプログラムが表示されるので、最初のコードセルをクリックし、［Shift］＋［Enter］キーを押せば実行できます。ただし、最初のコードセルではmakedata関数の定義を行うだけなので、実行しても結果は何も表示されません。
サンプルファイルには前処理の方法も記してあるので、興味のある方はご覧ください。なお、原付購入者数のデータはアンケートにより得られたものです。また、80歳以上のデータがなく、外挿（範囲外の値を推測）しているので、やや信頼性に欠けることをご承知おきください。
事前分布と尤度関数を定義して、サンプリングを行う
今回も分析の全体像から見ていきましょう。図2のようになります。相関係数の範囲は−1〜1なので、その確率密度をうまく表現できる事前分布が必要になります。形としてはベータ分布が適しているのですが、ベータ分布は台が0〜1です。そこで、台を−1〜1になるようにベータ分布を変換して使います。
また、それぞれの変数は独立したものではなく、相関がある（その可能性がある）ので、尤度関数としては個別の正規分布ではなく、それぞれの相関も併せて表現できる多変量正規分布を使います。
図2 事前分布と尤度関数の考え方
相関係数ρの分布は、台を−1〜1に拡張したベータ分布（拡張ベータ分布）とする。図中のExtBetaという関数はPyMCに存在しないので、ベータ分布を変換して使う。尤度関数は多変量正規分布（この場合は2変数）。多変量正規分布を定義するためのMvNormalクラスには、各変数の平均μ0,μ1と、分散・共分散行列Σを指定する。Σは合計を表す記号ではなく、分散・共分散行列を表す変数の名前であることに注意。
今回は、拡張ベータ分布と多変量正規分布という2つの分布が新しく登場しました。それらの分布がどのようなものであるかを確認してから、事前分布と尤度関数の定義を見ていきましょう。
以下の図3は、α＝2, β＝2のベータ分布と拡張ベータ分布の確率密度関数を可視化したものです。
図3 ベータ分布と拡張ベータ分布（α＝2, β＝2）
拡張ベータ分布はベータ分布の台を−1〜1に引き延ばし、全体の面積を1にするため高さを半分にしたもの（幅が2倍で高さが半分のオレンジの曲線）。相関係数の分布を表すのに使う。事前分布としてこのような分布を利用する。
拡張ベータ分布への変換方法が少し分かりにくいかもしれませんが、図3の横軸を引き延ばし、縦軸を半分にするイメージです。手順は以下の通りです。
- 元の分布の0〜1という範囲を2倍する → 範囲は0〜2になる
- そこから1を引く → 横軸の範囲が−1〜1になる
- 確率密度を2で割る → 縦軸の高さが半分になる
ちなみに、拡張ベータ分布を、元のベータ分布に戻すには、上とは逆の計算を行います。つまり、
- 拡張ベータ分布の−1〜1という範囲に1を足す → 範囲は0〜2になる
- それを2で割る → 横軸の範囲が0〜1になる
- 確率密度に2を掛ける → 縦軸の高さが倍になる
というわけです。
図4は、多変量正規分布の確率密度関数を可視化したものです。この例では、それぞれの変数の平均はμ0＝0,μ1＝0とし、分散・共分散行列Σは、標準偏差σ0＝1, σ1＝1、と相関係数ρ＝0.8を基に
という式で定義された値を使っています。具体的には、
の場合のグラフです。
図4 多変量正規分布（μ0＝0, μ1＝0, Σ＝[[1, 0.8],[0.8, 1]]）
この例は、2変数の多変量正規分布（3変数以上になると可視化は難しいが考え方は同じ）。立体的な正規分布といったイメージ。尤度関数としてこのような分布を利用し、事前分布を事後分布に更新する。
図3と図4のグラフを描くためのコードも参考として、サンプルファイルに含めてあります。
事前分布を定義する
新たに登場した分布のイメージが確認できたので、事前分布の定義から見ていきましょう。年齢と（購入者数当たりの）死傷者数は標準化されているので、それぞれの平均μ0, μ1については、平均が0、標準偏差が1の正規分布とします。また、標準偏差σ0, σ1については、いずれも元の正規分布の標準偏差を1とした半正規分布とします。ここまでは、前回までのお話で見た通りです。
相関係数ρの事前分布は、台が−1〜1の拡張ベータ分布とします。分布の形としては、0を中心として左右対称なものとしましょう。
ベータ分布にはαとβという2つのパラメーターがありますが、これらはどちらも分布の形を決めるもので、その大小関係で山の位置が決まります。αとβが等しい場合は左右対称な分布になります。一般に、拡張ベータ分布のパラメーターとしては、αやβの逆数であるκ（カッパと読みます）という値もよく使われます（κ ＝ 1/α、κ ＝ 1/βです）。ここでは、αとβをそのまま使い、いずれも値を2とします（図3で見たような0付近に山のある分布）。
ここまでをまとめると事前分布は以下のようになります。
- 平均μ0、μ1： μ＝0、σ＝1の正規分布
- 標準偏差σ0、σ1： 元の正規分布の標準偏差がσ＝1である半正規分布
- 相関係数ρ： α＝2、β＝2の拡張ベータ分布
では、これらの定義を書いてみましょう。平均や標準偏差の事前分布は2つずつありますが、shape引数に2を指定すれば、1つのコードでまとめて2つ定義できます（リスト1）。
# 平均の事前分布
mu = pm.Normal("mu", mu=0, sigma=1, shape=2) # shape=2により、[mu0, mu1]の形式になる
# 標準偏差の事前分布 (正の値のみをとるHalfNormal分布)
sigma = pm.HalfNormal("sigma", sigma=1, shape=2) # shape=2により、[sigma0, sigma1]の形式になる
# 相関係数の事前分布 (拡張ベータ分布)
rho_orig = pm.Beta("beta_raw", alpha=2, beta=2) # kappa=0.5
rho = pm.Deterministic("rho", 2 * rho_orig - 1) # 0〜1を2倍して0〜2とし、1を引いて-1〜1とする
ここまで見てきた事前分布の定義を素直に書いたもの。shape引数により、複数の事前分布をまとめて定義している。相関係数の事前分布については、rho_origがベータ分布に従うものとしている（rho_origの範囲は0〜1となる）ので、rho_origを2倍して1を引き、範囲を−1〜1に拡張したrhoを使う。
pymc.Normalが返す確率変数は2つの要素を持つリストの形になっており、muに代入されています。この場合、mu[0]がμ0に、mu[1]がμ1に当たります。また、sigmaに代入されているpymc.HalfNormal関数の返り値については、sigma[0]がσ0に、sigma[1]がσ1に当たります。
前回もお話ししましたが、上に示した「事前分布の定義」とは、正確には「ある事前分布に従う確率変数の定義」ということです。例えば、mu = pm.Normal("mu", ...)というコードであれば、正規分布に従う確率変数を定義していることになります。モデル内で利用される確率変数の名前が第1引数に指定した"mu"で、その確率変数を表すオブジェクトが返り値として返され、代入されたmuという変数名で参照できるようになります。さらに混乱させるような話ですが、確率変数は「変数」と名乗っているものの、数学的には「試行の結果を数値に対応させる関数」です（実用的には気にしなくてもほとんど問題がないので、そんな話もあったよなという程度にさらっと流していただいて結構なのですが）。
尤度関数を定義する
続いて尤度関数の定義です。尤度関数は多変量正規分布で、そのパラメーターには平均μ0とμ1のリスト（配列）と、分散・共分散行列Σを指定します。Σについては、それぞれの標準偏差σ0,σ1と相関係数ρを基に、（1）式で見たような配列として定義します。念のため（1）式を再掲しておきます。
以下のリスト2が分散・共分散行列Σを定義して、尤度関数（多変量正規分布）の確率変数を定義するpymc.MvNormalクラスに引数として指定した例です。
# 分散・共分散行列（多変量正規分布のパラメーターとして指定）
cov_mat = pm.Deterministic(
"cov_matrix",
pm.math.stack(
[
[sigma[0]**2, rho * sigma[0] * sigma[1]],
[rho * sigma[0] * sigma[1], sigma[1]**2]
]
)
)
# 尤度関数（多変量正規分布）
pm.MvNormal("obs", mu=mu, cov=cov_mat, observed=data)
分散・共分散行列の確率変数は、sigmaとrhoを使って計算で求めるので、pymc.Deterministic関数で記述する。その中ではパラメーターを2×2の行列にまとめるためにpymc.math.stack関数を使う。単なる値の行列ではなく、計算の手順をまとめたデータ（計算グラフと呼ばれる）を作成していることに注意。
ここでは、分散・共分散行列を定義するために、pymc.math.stack関数を使っていることに注目してください。引数としては、（1）式で見たΣを構成するパラメーターを指定しているだけです。後は、pymc.MvNormalの引数に、平均のリストmuと分散・共分散行列の定義であるcov_matを引数として指定するだけです。
後で紹介しますが、LKJ事前分布を利用すると、上で見た拡張ベータ分布や分散・共分散行列を記述しなくても、同様の分析ができてしまいます。それなら、最初から楽な方法でやってくださいよ、と言われそうですが、自分で記述した方が仕組みがよく理解できるので、あえて一歩ずつ進めています。
サンプリングを行う
サンプリングのためのコードは前回までに見た例とほぼ同じです。つまり、pymcモジュールのsample関数を使うだけです。例によって、これまで断片的に記したコードも含め、結果の要約を表示するところまで、コードの全体を掲載します。サンプルファイルの4番目のコードセルに入力されています。
import numpy as np
import pymc as pm
import arviz as az
# データを取得する
data = makedata()
with pm.Model() as model:
# 平均の事前分布
mu = pm.Normal("mu", mu=0, sigma=1, shape=2) # shape=2により、[mu0, mu1]の形式になる
# 標準偏差の事前分布 (正の値のみをとるHalfNormal分布)
sigma = pm.HalfNormal("sigma", sigma=1, shape=2) # shape=2により、[sigma0, sigma1]の形式になる
# 相関係数の事前分布 (拡張ベータ分布)
rho_orig = pm.Beta("beta_raw", alpha=2, beta=2) # kappa=0.5
rho = pm.Deterministic("rho", 2 * rho_orig - 1) # 0〜1を2倍して0〜2とし、1を引いて-1〜1とする
# 分散・共分散行列（多変量正規分布のパラメーターとして指定）
cov_mat = pm.Deterministic(
"cov_matrix",
pm.math.stack(
[
[sigma[0]**2, rho * sigma[0] * sigma[1]],
[rho * sigma[0] * sigma[1], sigma[1]**2]
]
)
)
# 尤度関数（多変量正規分布）
pm.MvNormal("obs", mu=mu, cov=cov_mat, observed=data)
# サンプリング
trace = pm.sample(draws=8000, chains=2, random_seed=42)
# 結果の要約を表示
az.summary(trace, var_names=["rho"])
# サンプリングというコメント以降が、これまで断片的に記したコードの続き。後でベイズ因子を求めるために、やや多めに（8000個×2チェーンの）サンプリングを行った。事後分布の要約として、相関係数rhoのみの推定値などを表示する。
pymc.sample関数によってサンプリングが実行されると、各変数のサンプルが自動的に作られます。上のコードでは、それらをarvizモジュールのsummary関数で要約しています。実行すると、図5のような結果が表示されます。
図5 事後分布の要約
相関係数ρの点推定値は、rhoの行のmean列の値、つまり−0.776となる。94％信用区間は−0.943〜−0.575。負の部分のみとなっているので、負の相関があるものと示唆される。なお、データから直接計算した相関係数は−0.882だが、今回は事前分布が0を中心とした分布であるため、その影響で点推定値が標本相関係数より少し0寄りの値になっている。r_hatが1.0である（＝収束している）ことも確認しておこう。
ベイズ因子を求める 〜 ベイズ相関検定を行う
続けて、Savage-Dickey法を使ってベイズ因子を求めてみます。事前分布（拡張ベータ分布）での、0に対する確率密度と、事後分布をスムーズな曲線にした場合の0に対する確率密度との比を取って、ベイズ因子BF10の値を求めます。改めて帰無仮説と対立仮説を確認しておきましょう。
- 帰無仮説H0 : ρ＝0
- 対立仮説H1 : ρ ≠ 0
以下のコードはサンプルファイルの5番目のコードセルに入力されています（リスト4）。
from scipy.stats import gaussian_kde, beta
# 事後サンプル
rho_samples = trace.posterior["rho"].values.flatten()
# KDEによる事後分布の確率密度
kde = gaussian_kde(rho_samples) # 滑らかな曲線にする
posterior = kde.pdf(0)[0] # 0の位置での確率密度
# 事前分布の確率密度
rho = 0
prior = beta.pdf((rho + 1) / 2, a=2, b=2) * 0.5 # 0.75となる
# ベイズ因子
bf10= prior / posterior
print(f"BF10: {bf10:.3f}")
# 出力例：
# BF10: 783.822
サンプリングで得られた事後分布を、scipy.statsモジュールのgaussian_kdeクラスを使って滑らかな曲線にする。事後分布と事前分布のρ＝0に対する確率密度関数の値を求め、前者を後者で割ってベイズ因子を求める。結果は783.822となる。
事前分布のρ＝0に対する確率密度関数の値を求めるためにbeta.pdf((rho + 1) / 2, a=2, b=2) * 0.5というコードが書かれています。これは、以下のような意味の計算です。
- 拡張ベータ分布の横軸を、元のベータ分布の横軸に戻す（図3の下で述べたように、1を足して2で割ればよい）
- 戻した点で確率密度関数の値を求め、2で割る（高さを半分にする）
結果はBF10＝783.822となりました。この値はJeffreys（1961）の評価尺度に照らし合わせると、H1を支持する「決定的な証拠となっている」ことが分かります。従って、対立仮説ρ ≠ 0（相関係数は0ではない）が支持されます。相関係数が負なので、実質的に負の相関があると言えます。つまり、若い人ほど原付事故の死傷者数が多く、年齢が上がると死傷者数が少なくなるということになります。
既にお気付きかと思いますが、この連載を続けて読んでくださっている方には、今回の記事についても、幾度となく既視感を覚えたのではないでしょうか。というのも、これまでに見てきたさまざまな分析の手順はどれもこれも同じだからです（もう飽きちゃったよ、とは言わないでくださいね）。異なるのはどのようなモデルを立てるか、ということだけです。なお、対立仮説が正しい確率を求める方法や、SMCサンプリングにより周辺尤度を求めてベイズ因子を計算する方法もこれまでと同じなので、ここには掲載せず、サンプルファイルに含めておくにとどめます。
LKJ事前分布を利用する
上で予告したように、LKJ事前分布を利用すると、上で見た拡張ベータ分布や分散・共分散行列を記述しなくても、同様の分析ができます。いきなりですが、コードを見てみましょう（リスト5）。説明はコードの後に記します。
import numpy as np
import pymc as pm
import arviz as az
# データを取得する
data = makedata()
with pm.Model() as model:
# 平均の事前分布
mu = pm.Normal("mu", mu=0, sigma=1, shape=2) # shape=2により、[mu0, mu1]の形式になる
# 標準偏差の分布
sd_dist = pm.HalfNormal.dist(sigma=1, shape=2) # 確率変数ではなく分布を求める
# 分散・共分散行列（多変量正規分布のパラメーターとして指定）
chol, corr, sigmas = pm.LKJCholeskyCov(
"chol_cov", n=2, eta=2, sd_dist=sd_dist
) # コレスキー分解された分散・共分散行列、相関行列、標準偏差が自動的に作られる
# 相関係数の計算
rho = pm.Deterministic("rho", corr[0, 1])
# 尤度関数（多変量正規分布）
obs = pm.MvNormal("obs", mu=mu, chol=chol, observed=data)
# サンプリング
trace = pm.sample(draws=8000, chains=2, random_seed=42)
# 結果の要約
az.summary(trace, var_names=["rho"]) # 点推定値はrho=-0.779、94％信用区間は-0.938〜-0.575
pymc.LKJCholeskyCovクラスを利用する。引数のnは変数の個数。etaは分布の形を決めるパラメーターで、ベータ分布のα, βと同じ値を指定すればよい。sd_distにはpymc.HalfNormal.distメソッドなどで得られた標準偏差の分布を指定する。返り値はcholがコレスキー分解された分散・共分散行列、corrが相関行列、sigmasが、各変数の標準偏差。corr[0, 1]が相関係数に当たり、sigmas[0]とsigmas[1]が、それぞれσ0とσ1に当たる。後の処理はベータ分布を利用した場合のコードと全く同じ。
理屈は後回しにして、コードを実行してみてください。相関係数について、−0.779という点推定値と−0.938〜−0.575という94％信用区間が得られます。
さて、ここで事前分布として使われるLKJ分布は、Lewandowski, D., Kurowicka, D., & Joe, H.（2009）で提案された分布で、論文の著者の頭文字から取られた名前です。詳細についてはかなり難しくなるので、これ以上は触れませんが、興味のある方は元の論文に当たってみてください（こちらで読むこともできます）。
pymc.LKJCholeskyCovクラスに指定する引数や返り値については、リスト5のキャプションにも記しましたが、引数sd_distの指定と、返り値、特にcholのコレスキー分解の意味について補足しておきます。
引数sd_distにはpymc.HalfNormal.distメソッドなどで得られた標準偏差の分布を指定します。これは、pymc.HalfNormalクラスで定義される確率変数とは異なるものです。正確には「計算グラフに含まれる分布オブジェクト」とでも言うべきものですが、とりあえずは「分布そのもの」を渡す、といったイメージで捉えていただいてもいいかと思います。
返り値のcholはコレスキー分解された分散・共分散行列を表す確率変数です。コレスキー分解とは、行列を下三角行列（対角成分の右上の成分が全て0の行列）とその転置の積に分解するような計算です。コレスキー分解を行うと、逆行列の計算が簡単になり、計算の結果が安定するといったメリットがあります。……が、とりあえずは、分散・共分散行列を変換して使いやすいように表現したもの、と考えておいてもらってけっこうです。コレスキー分解がどのようなものであるかは、後のコラムで簡単に解説します。
相関行列（相関係数）のLKJ事前分布や標準偏差の事前分布はpymc.LKJCholeskyCovの中で作られます。それらの確率変数が返り値のcorrとsigmasで参照されるというわけです。なお、2変数の場合、LKJ分布は拡張ベータ分布と同じものになります。
コラム コレスキー分解とは
コレスキー分解とは上でも触れたように、行列を下三角行列に分解することです。
具体的な計算方法（アルゴリズム）は割愛しますが、どのように分解されるかという簡単な例を紹介しておきます。以下のような行列があるものとします。
コレスキー分解した下三角行列Lは、以下の通りです。右上の成分が0になっていますね。
Lを転置した上三角行列LTは、以下の行列です。
このとき、A＝LLTとなります。このことを確認するためのコードもサンプルファイルに含めてあります。
少しばかり昔話を。40年以上前のことです。大学2回生の時に、心理学専攻の学生のために計算機実習の授業があったのですが、そのときの期末課題がPascalでコレスキー分解により逆行列を求めるプログラムを作れ、というものでした。文学部の学生で、数学は苦手、プログラミングも初心者だった私は全く歯が立たず、結局、コレスキー分解ができずに、逆行列を求める関数を使った数行のコードを提出することしかできませんでした。今でこそ分かったようなことを書いていたりしますが、最初はそんなものです。ほぼ諦めていた単位は、なぜか「優」でした（平常点のおかげかもです）。苦い思い出です。故N先生ごめんなさい。
サンプルファイルのリスト5の後には、Savage-Dickey法によってベイズ因子を求めるコードも入力してあります。コードはリスト4と全く同じですが、結果はBF10＝378.957となります。2変数の場合、LKJ事前分布を利用した場合と、拡張ベータ分布を利用した場合は同じ処理になりますが、リスト4で求めたBF10＝783.822という値とずいぶんと異なります。これは、内部のアルゴリズムなどが異なるためと考えられます。大きな違いがあるように思われるかもしれませんが、実際にはそれほどではなく、「決定的な証拠」であるという傾向は変わりません。
拡張ベータ分布の例とLKJ事前分布の例でベイズ因子に大きな差が出るのは、2つの要因が重なった結果です。1つは、上で触れたように、拡張ベータ分布の例とLKJ事前分布の例で、内部のアルゴリズムが異なるためです。そのため、サンプリングされた事後分布にわずかな違いが生じます。加えて、Savage-Dickey法により、特定の点（この場合ρ＝0）に対する確率密度を比較することが差を大きくしています。今回の例のように、事後分布の山がρ＝0から大きく外れた位置にある場合、ρ＝0近辺のサンプルが少なくなり、確率密度関数の推定値に誤差が生じます。その値が分母となるので、わずかな違いによりベイズ因子の値が大きく変わります。それらの2つの要因によりベイズ因子に大きな差が出たというわけです。なお、観測データの相関係数の絶対値が極端に大きな値でなければ、Savege-Dickey法で求めた拡張ベータ分布とLKJ事前分布でのベイズ因子はほぼ等しくなります。
ちなみに、SMCサンプリングにより周辺尤度の比を求めてベイズ因子を計算するコードも参考としてサンプルファイルに含めてあります。その場合、BF10＝752.815となります。
ただし、この値と、リスト4で得られた拡張ベータ分布でのベイズ因子の値BF10＝783.822とがほぼ一致するのは「たまたま」です。例えば、乱数の初期設定を変えてrandom_seed=0とすると、拡張ベータ分布を利用した場合のベイズ因子はBF10＝2778.753となり、LKJ事前分布の場合はBF10＝568.046となります。このように（今回の例のように）、事後分布の中心が事前分布の中心から極端に外れている場合、Savage-Dickey法での結果は不安定になることがあります。しかし、上でも述べたように、全体的な傾向は変わりません（これらの値は指数関数的に変わるので、10を底とする対数を取ってみると、それぞれ3.44、2.75となり、それほど大きな差ではないことが分かります。Jeffreysの評価尺度でも階級の幅が10を底とする指数関数の値になっています）。ベイズ因子のより正確な値を求めたい場合には、SMCサンブリングによって周辺尤度の比を求める方法や、この後のコラムで紹介する解析解を求める方法が適しています。
というわけで、今回のメインのお話はここまでです。が、解析解を求める方法も次のコラムで紹介しておきましょう。解析解はBF10＝752.226です。
コラム 相関係数のベイズ因子を解析的に求める方法
実は、pingouinモジュールのbayesfactor_pearson関数や、統計パッケージJASPなどでは、相関係数のベイズ因子を解析的に求める方法を使っています。これは、Ly, Verhagen & Wagenmakers（2015）による以下の式を利用したものです（両側検定の場合の式）。
ただし、α＝1/κ、rは相関係数、Bはベータ関数（scipy.special.beta）、Γはガンマ関数（scipy.special.gamma）、2F1はガウス超幾何関数（scipy.special.hyp2f1）と呼ばれるものです。かっこ内はそれぞれの値を求めるために利用できる関数の名前です。いずれも、かなり専門的な関数ですが、答えは単に引数を指定するだけで得られます（最後の「今回使った新出のクラス・関数・メソッド」のところにも説明を追加してありますが、実際のところ、これらの関数を直接使う機会はあまりないかと思います。わけが分からないよ、となっても、あまり気にしないでください）。
今回の例であれば、α＝2, n＝15なので、以下のようなコードで解析解が求められます。上の（2）式をそのままコードとして表しただけです。
import numpy as np
from math import pi, sqrt
from scipy.special import beta, gamma, hyp2f1
# データを取得する
data = makedata()
# パラメーター
kappa = 0.5
a = 1 / kappa
r = np.corrcoef(data[:, 0], data[:, 1])[0, 1] # 相関係数
n=15
# 解析解を求める
numerator = 2**(1-2*a) * sqrt(pi) * gamma((n + 2 * a - 1) / 2) # 分子
denominator = beta(a, a) * gamma((n + 2 * a) / 2) # 分母
f = hyp2f1((n-1)/2, (n-1)/2, (n+2*a)/2, r**2) # ガウス超幾何関数
# ベイズ因子
bf10 = numerator / denominator * f
print(f"BF10: {bf10:.3f}")
# 出力例：
# BF10: 752.226
上のコードは両側検定のみに対応した（2）式を素直に記述したもの。元の論文には片側検定のための補正を行う式も書かれているが、かなり複雑なので、ここでは省略した。結果はBF10＝752.226となる。
さらに、pingouinモジュールのbayesfactor_pearson関数を利用すれば、自分で式を書かなくても、値を与えるだけでベイズ因子が求められます。しかも、上のような両側検定だけでなく、片側検定のための補正も行えます。その方法の説明も含めて、サンプルファイルにコードを入力しておきました。なお、bayesfactor_pearson関数の内部では、リスト6のように各項をそのまま計算するのではなく、誤差の拡大を防ぐために対数を取って計算しています（例えばbeta関数の代わりにbetaln関数を使うなど）。
今回は、「関係」に注目して、相関係数のベイズ推定とベイズ検定を行いました。事前分布として拡張ベータ分布やLKJ分布を利用したことと、尤度関数として多変量正規分布を利用したことが、これまでと異なる主要なポイントでした。
さて、次回も「関係」に注目します。いわゆるカイ二乗検定（独立性の検定）をベイズ統計によって行ってみます。「雨男や雨女は本当に存在するのか？」という事例を取り上げます。次回もお楽しみに！
今回使った新出のクラス・関数・メソッド（主なもの）
引数も主なものだけを掲載します。引数に「=値」と書かれたもの（例：axis=0）は既定値を表します。
Beta： ベータ分布に従う確率変数を定義するクラス
- モジュール： pymc
- 形式： Beta(name, alpha, beta)
- 引数：
- name： 確率変数に付ける名前を文字列で指定する。この名前でModelに登録される
- alpha： 分布の形を決めるパラメーター。betaと合わせて平均α/(α＋β)や分散が決まる
- beta： 分布の形を決めるパラメーター。alphaと合わせて平均α/(α＋β)や分散が決まる
stack： 確率変数をベクトルや行列にまとめる関数
- モジュール： pymc.math
- 形式： stack(tensors, axis=0)
- 引数：
- tensors： ひとまとめにしたい確率変数などを指定する
- axis： 結合の方向
- 返り値： 計算グラフのノード（計算の方法を定義したものと考えるとよい）
※axisについては、既定値の0であれば行方向に並び、1であれば列方向に並ぶイメージ。例えば、話を簡単にするために、tensorsを
a = [1, 2, 3]
b = [4, 5, 6]
のような定数とすると、pm.math.stack([a, b], axis=0).eval()の場合、返り値は、
[[1, 2, 3],
[4, 5, 6]]
となり、pm.math.stack([a, b], axis=1).eval()の場合、返り値は
[[1, 4],
[2, 5],
[3, 6]]
となる。
LKJCholeskyCov： LKJ分布に従う確率変数を定義するクラス
- モジュール： pymc
- 形式： LKJCholeskyCov(name, eta, n, sd_dist)
- 引数：
- name： 確率変数に付ける名前を文字列で指定する。この名前でModelに登録される
- eta： 分布の形状を表すパラメーター。2変数の場合、拡張ベータ分布の1/κと同じ値を指定する
- n： 変数の個数
- sd_dist： 標準偏差の分布。HalfNormal.distメソッドなどで作成した分布を指定する
MvNormal： 多変量正規分布に従う確率変数を定義するクラス
- モジュール： pymc
- 形式： MvNormal(name, mu, cov, observed)
- 引数：
- name： 確率変数に付ける名前を文字列で指定する。この名前でModelに登録される。
- mu： 平均を表すパラメーターのリスト
- cov： 分散・共分散行列を表すパラメーターのリスト（配列）
- covの代わりにコレスキー分解された分散・共分散行列を引数cholに指定してもよい
- observed： 観測データ。この引数の指定がない場合は事前分布として扱われ、この引数の指定がある場合は尤度関数として扱われる
corrcoef： 相関行列を求める関数
- モジュール： numpy
- 形式： corrcoef(x, y)
- 引数：
- x： 値のリストや配列
- y： 値のリストや配列（省略可能）
- 返り値： xとyの相関行列。2変数の場合、相関係数は返り値の[0, 1]や[1, 0]の位置にある値
※xやyに二次元のリストや配列を指定すると、それぞれの行同士の相関行列が求められる。
beta： ベータ関数の値を求める関数
- モジュール：scipy.special
- 形式： beta(a, b)
- 引数：
- a, b： ベータ関数に与える引数の値
- 返り値： ベータ関数の値
※ベータ関数B(a, b)は次に示すガンマ関数Γ(z)を使うと以下のように表される。
○gamma： ガンマ関数の値を求める関数
- モジュール： scipy.special
- 形式： gamma(z)
- 引数：
- z： ガンマ関数に与える引数の値
- 返り値： ガンマ関数の値
※ガンマ関数は階乗の考え方を複素数の範囲に広げた関数。考え方としては以下のような関係。
ガンマ関数の値は、以下の式で求められる。
○hyp2f1： ガウス超幾何関数の値を求める関数
- モジュール： scipy.special
- 形式： hyp2f1(a, b, c, z)
- 引数：
- a, b, c, z： 超幾何関数に与える引数の値
- 返り値： 超幾何関数の値
※|z|<1について、ガウス超幾何関数2F1は以下のような式で表される
ただし、(a)nは、
とする（この(a)nはポッホハマー記号と呼ばれる）。
参考文献
Lewandowski, D., Kurowicka, D., & Joe, H. (2009). Generating random correlation matrices based on vines and extended onion method. Journal of Multivariate Analysis, 100(9), 1989-2001.
Ly, Alexander & Verhagen, A J & Wagenmakers, Eric-Jan. (2015). Harold Jeffreys’s default Bayes factor hypothesis tests: Explanation, extension, and application in psychology. Journal of Mathematical Psychology.
Copyright© Digital Advantage Corp. All Rights Reserved.