2018.08.10 2023.10.14

【初心者向け】Pythonでベイズ統計学入門【相関がある確率を推定する】

京極真

記事内にプロモーションを含む場合があります

京極真

本記事では「PythonでStanを使って、相関がある確率を推論するにはどうしたらよいですか？」という疑問にお答えします

本記事のポイント

PyStanで相関係数と「相関がある」という仮説の確からしさを推論する方法がわかります
ぼく自身のPystanの練習をかねた備忘録ですｗ

この記事を書いているぼくは普段、主にRとStanでベイズ統計を活用しています。

現在、Pythonを使えるようになりたくて、PyStanを勉強しているところです。

PyStanはumdemyの動画で学習中でして、以下の記事でがっつりレビューしています。

数十時間しっかりやればできることが増えるので、費用対効果は高めだと思います。

また30日間返金保証がありますので、その間は実質無料ですから「自分には理解できないかも」という人は、ひとまずお試しすることを強くオススメします。

ベイズ統計で相関分析すると実現できる未来

ベイズ統計で相関分析すると「相関がある」といえる確率を直に推定できます。

従来の統計でできることは、帰無仮説（相関がない（相関＝0））と対立仮説（相関がある（相関≠0））を立てて、帰無仮説のもとで手元のデータが得られる確率を計算し、5％未満あるいは1％未満であれば帰無仮説を棄却し、対立仮説を採用するというものでした。

これは、帰無仮説（相関がない（相関＝0））は真であるという仮定のもとで統計量を調べているだけであって、研究者が本当に知りたい対立仮説（相関がある（相関≠0））の確からしさについては何も検討していません。

他方、ベイズ統計は「相関がある（相関≠0）」という仮定の確からしさを直に推論することができます。

基礎からのベイズ統計学: ハミルトニアンモンテカルロ法による実践的入門

created by Rinker

¥3,520 (2023/12/01 12:19:44時点 Amazon調べ-詳細)

もちろん、これは得られたデータや統計モデルに依存しています。

けども、「相関がある（相関≠0）」という研究仮説が正しい確率について直接言及できるのは、ベイズ統計の大きな利点です。

多くの人たちは、忙しい仕事の合間をぬって研究する以上、やっぱ設定した研究仮説がどの程度確からしいのかを知りたいはずです。

ベイズ統計はそれを実現できる可能性のツールというわけです。

Pythonによるベイズ統計で相関がある確率の推定

課題：2つの変数間に相関がある確率を明らかにしたい

ここでは「2つの変数間に相関があるという研究仮説が正しい確率」を明らかにするという課題を設定します。

Pythonで乱数を発生させて、相関のある2つのシミュレーションデータを作成しました。

#シミュレーションデータの生成

np.random.seed(1) x = np.random.randint(0, 50, 100)

y = x + np.random.normal(0, 10, 100)

np.corrcoef(x, y)

以下の通り、約0.78の相関係数を示すシミュレーションデータができました。

array([[1. , 0.78510866],

[0.78510866, 1. ]])

PyStanで相関がある確率を推定する

次に、PyStanを使って2変数の「相関係数」と「相関があるという研究仮説が正しい確率」を計算します。

まずはパッケージを読み込みます。

# パッケージの読み込み

import numpy as np

import pandas as pd

import pystan

次に、Stanにデータを渡すために、np/c_ という関数を使いました。#Stanにデータを渡す X = np.c_[x, y]

Stanコードは以下の通りでして、これはPearsonの相関分析のベイズ版です。

Stanコードは以下の資料を参考にしています。

Stan超初心者入門 from Hiroshi Shimizu

すこし補足しておくと、相関があるという仮説の確からしさは以下のコードが対応しています。

delta_over = rho > 0 ? 1 : 0; は相関が0以上の確率を推論するように指示しています。

これは、研究仮説に応じて、相関が0.4以上 delta_over = rho > 0.4 ? 1 : 0; とか、0.8以上delta_over = rho > 0.8 ? 1 : 0; などのように柔軟に変えることができます。

#Stanコード stan_model = “””

data{

int N; #サンプルサイズ

vector[2] X[N]; #データ

}

parameters{

vector[2] mu; #平均

vector<lower=0>[2] sigma; #標準偏差

real<lower=-1,upper=1> rho; #相関

}

transformed parameters{

cov_matrix[2] Sig;

Sig[1,1] = sigma[1]^2;

Sig[2,2] = sigma[2]^2;

Sig[1,2] = rho*sigma[1]*sigma[2];

Sig[2,1] = rho*sigma[1]*sigma[2];

}

model{

mu ~ normal(0,100);

sigma ~ cauchy(0,5);

rho ~ uniform(-1,1);

X ~ multi_normal(mu,Sig);

}

generated quantities{

real delta_over;

delta_over = rho > 0 ? 1 : 0; #相関がある（0以上）という仮説が正しい確率

}

“””

# Stanコードの読み込み

sm = pystan.StanModel(model_code=stan_model)

# Stanにデータを渡す

stan_data = {‘N’ : X.shape[0], ‘X’ : X}

# Let’s Stan! fit = sm.sampling(data = stan_data, iter=2000, warmup=500, chains = 4)

# 結果の出力

print(fit)

結果は以下の通りです。

Rhatをみると、全部1.0なのできちんと収束していることがわかります。

相関係数は0.77であり、最初に確認した結果とほぼ同じです（当たり前）。

また相関があるという研究仮説が正しい確率は1.0、すなわち100％の確率で確かそうだということがわかりました。

Inference for Stan model: anon_model_7d53975e06a8e1e50479d9a4b80eaf60. 4 chains, each with iter=2000; warmup=500; thin=1; post-warmup draws per chain=1500, total post-warmup draws=6000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat mu[0] 22.07 0.02 1.43 19.29 21.1 22.07 23.03 24.88 3399 1.0 mu[1] 21.77 0.03 1.63 18.53 20.65 21.8 22.89 24.94 3410 1.0 sigma[0] 14.2 0.02 0.99 12.46 13.5 14.14 14.82 16.31 3281 1.0 sigma[1] 16.0 0.02 1.13 13.96 15.22 15.94 16.7 18.37 3343 1.0 rho 0.77 7.0e-4 0.04 0.68 0.75 0.77 0.8 0.84 3558 1.0 Sig[0,0] 202.57 0.5 28.35 155.16 182.36 199.95 219.69 265.91 3245 1.0 Sig[1,0] 176.4 0.56 28.82 127.31 156.1 173.9 193.81 240.22 2693 1.0 Sig[0,1] 176.4 0.56 28.82 127.31 156.1 173.9 193.81 240.22 2693 1.0 Sig[1,1] 257.26 0.64 36.49 194.79 231.62 254.04 278.9 337.28 3277 1.0 delta_over 1.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 6000 nan lp__ -597.8 0.03 1.59 -601.7 -598.6 -597.4 -596.6 -595.7 2219 1.0

Samples were drawn using NUTS at Tue Aug 7 15:30:22 2018. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).