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

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

本記事のポイント

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

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

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

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

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

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

おすすめ
【レビュー】【PythonとStanで学ぶ】仕組みが分かるベイズ統計学入門

Udemyの「【PythonとStanで学ぶ】仕組みが分かるベイズ統計学入門」に興味あるにゃーぶーたん きょうごくそれなら、ぼくは実際に購入したので、わかりやすくレビューするよ 本記事のポイント 【P ...

続きを見る

さて、本記事ではぼくの勉強をかねて、PyStanで相関分析をさくっと紹介します。

たぶん、これからPyStanをやってみたいエンドユーザーの参考になるはず。

PyStanで差がある確率を推論する方法は以下の記事で紹介しましたので、関心がある人はあわせてお読みください。

おすすめ
【初心者向け】Pythonでベイズ統計学入門【平均に差がある確率を推定する】

きょうごく本記事では「はじめてPythonでStanを使いますが、平均の差の推定を行うにはどうしたらよいですか？」という疑問にお答えします。 本記事のポイント Pystanで平均値の差を推定する方法が ...

続きを見る

研究計画書作成・研究論文執筆チェックリスト

今すぐ無料ダウンロードする

※効率的・効果的に書けるようになろう！

IMRaDを使った研究論文の書き方講座

今すぐ無料Webセミナーに参加する

※参加者募集！

研究計画書の書き方講座

今すぐ無料Webセミナーに参加する

※参加者募集！

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

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

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

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

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

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

posted with ヨメレバ
豊田 秀樹 朝倉書店 2015-06-25
Amazon
Kindle
楽天ブックス

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

けども、「相関がある（相関≠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コードは以下の資料を参考にしています。

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

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).

まとめ

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

本記事で解説したように、ベイズ統計が使えると研究者が知りたい「研究仮説が正しい確率（仮説が確からしい程度）」を直に推論できます。

これは、実践領域で研究する人たちにとって、とても魅力的な特徴だと思います。

本記事ではPythonでベイズ統計をさくっと紹介しましたが、ぼく自身は動画でまだまだ学習中です。

「こういうことができるようになりたいなぁ」という人は、以下の記事でがっつりレビューしていますので、ぜひご覧になってください。

おすすめ
【レビュー】【PythonとStanで学ぶ】仕組みが分かるベイズ統計学入門

Udemyの「【PythonとStanで学ぶ】仕組みが分かるベイズ統計学入門」に興味あるにゃーぶーたん きょうごくそれなら、ぼくは実際に購入したので、わかりやすくレビューするよ 本記事のポイント 【P ...

続きを見る

-研究教育