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

投稿日:2018年8月10日 更新日:

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

本記事のポイント

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

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

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

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

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

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

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

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

続きを見る

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

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

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

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

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

続きを見る

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

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

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

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

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

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

けども、「相関がある(相関≠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 ...

続きを見る

おすすめ記事一覧

1

きょうごく本記事では「コミュニケーション能力を高めたい。けど、そもそもコミュニケーション能力って何だろう?コミュニケーション能力の本質を理解したい」という疑問にお答えします こんな方におすすめ コミュ ...

2

きょうごく本記事では研究を学べるおすすめ本を紹介した記事のまとめます こんな方におすすめ 研究を学べるおすすめ本を知りたい 選りすぐりのおすすめ研究本で学びたい 本記事を書いているぼくは作業療法士であ ...

3

大学院に進学したい作業療法士向けのまとめ記事。 本記事では、大学院進学に興味のある作業療法士向けに書いた記事をまとめました。 今後も記事は増える予定ですが、これから大学院進学を検討している人はぜひ参考 ...

4

きょうごく本記事では、このブログで解説した作業療法士を目指す高校生のための記事をまとめています 本記事の内容 『作業療法士が何なのかよく分からない』という人から、『作業療法士になりたい』といった方向け ...

5

きょうごく本記事では学会発表・学会参加で失敗したくない人向けの記事をまとめています こんな方におすすめ 学会発表・学会参加のコツを知りたい 本記事では学会発表・学会参加で失敗したくない 学会発表のレベ ...

-研究
-

Copyright© 京極真の研究室 , 2019 All Rights Reserved.