相関係数のベイズ推定

相関係数のベイズ推定を、2変量正規分布を適用して行うスクリプトをPyMC3で書いてみた。

データは、

Data = [[179, 121], [185, 89], [181, 106], [210, 100], [202, 104],

[197, 91], [227, 107], [208, 101], [188, 99], [189, 112],

[240, 111], [203, 103], [202, 85], [190, 95], [221, 96],

[150, 83], [205, 103], [177, 90], [221, 92], [172, 86],

[190, 87], [218, 101], [205, 100], [208, 120], [224, 98]]

である。

上のデータに対するPyMC3を用いたMCMCサンプリングのPythonスクリプトを以下のように用意した。完全なスクリプトはこのウェブサイトの後ろにリスト１として示した。

with pm.Model() as two_v_n:

μ = pm.Normal('μ', mu = 0.0, sigma = 1000, shape = 2)

σ1 = pm.HalfNormal('σ1', sigma = 1000)

σ2 = pm.HalfNormal('σ2', sigma = 1000)

ρ = pm.Uniform('ρ', lower = -1.0, upper = 1.0)

cov_ = pm.math.stack([[σ1**2, ρ*σ1*σ2], [ρ*σ1*σ2, σ2**2]])

cov = pm.Deterministic('cov', cov_)

Y = pm.MvNormal('Y', mu = μ, cov = cov, observed = Data)

trace = pm.sample(3000, chains = 4)

リスト１を実行すると、PyMC3によるMCMCサンプリングが行われ、結果が以下のように図示される。実行は、Ubuntu上で行った。

サンプリングのトレースは下図（図１）のようである。

図１

図１は、スクリプト

az.plot_trace(trace)

によって描かれたものである。サンプリングに問題はないと思われる。

相関係数ρの事後分布を、次のスクリプトで描いた。

az.plot_posterior(trace['ρ'], ref_val = 0.0)

plt.title('ρ', fontsize = 20)

図２に示すグラフが描かれる。相関係数の平均値はρ=0.31で、94% HPDは［-0.045, 0.62］であることが示されている。ref_val = 0.0で指定された値より相関係数ρの値が小さい、すなわち負の値の確率は5.4%である。

図２

なお、PyMC3の説明は、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。

リスト１　2変量正規分布による相関係数のベイズ推定

import pymc3 as pm

import arviz as az

import numpy as np

import matplotlib.pyplot as plt