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

投稿日:2018年7月29日 更新日:

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

本記事のポイント

  • Pystanで平均値の差を推定する方法がわかります
  • ぼく自身のPystanの練習をかねた備忘録ですw

Pythonでベイズ統計を勉強する理由

ぼくは普段、Mplusの他に、Rでstan、JAGS、brmsなどいろいろ使用しています。

なのに、Pystan(Pythonでベイズ統計)を練習しはじめた理由は、話題のPythonを勉強したいという動機からです(笑)。

ややミーハーな理由ですいません。

おすすめ
【断言】新しいことに挑戦したい人にオススメの4つのテーマ【2019年版】

きょうごく本記事では「最近、行動の重要性をいろんなところで聞きます。私も行動は重要だと思いますので、新しいことに挑戦したい!けど、何に挑戦したらよいの?」という疑問にお答えします 本記事のポイント こ ...

続きを見る

Pystanの勉強はudemyの動画教材でやっています。

この動画は、ベイズ統計の基礎、MCMCの原理、グラフィカルモデルの基礎、さまざまな統計モデルの基礎からStanコードへの落とし込み方、結果の解釈などなど、めちゃくちゃわかりやすく解説してくれています。

学習効率がとてもよいので、これからベイズ統計を学びたい人や、ぼくのようなエンドユーザーでちょっとかじっているけどももっと腕を上げたい人には特にオススメです。

がっつりレビューを書いたので、ぜひお読みいただけたらと思います。

マジでできることが広がりますので、コスパ最高ですよ。

30日間返金保証があるので、その間なら実質無料ですから、「自分にはハードルが高いかも」という人でもとりあえず試すとよいです。

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

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

続きを見る

Pystanで平均の差の推定

課題:2つの変数の平均に差がある確率を明らかにしたい

ここではかりに「2つの変数の平均に差がある確率を推定する」という課題を設定します。

従来の統計の場合、帰無仮説「2つの変数の平均に差がない」と対立仮説「2つの変数の平均に差がある」を設定します。

そして、帰無仮説のもとで手元にあるデータが得られる確率を計算し、それが5%未満であれば帰無仮説を棄却することになります。

でもこの場合、帰無仮説が間違っているとも言えないし、かりに対立仮説を採用するにしても平均に差がある確率がどの程度あるのかについて何も言えません。

他方、ベイズ統計を活用したら2つの変数の平均に差がある確率を直に推定できるので、「AとBの平均に差がある確率は92%です」などのように知りたいことを知ることができます。

Pystanで平均の差を推定する

2つの変数の平均に差がある確率を推定するために、ここではPystanを用いています。

Pystanの使い方の詳細は動画にゆずりますが、慣れればそんなに難しくないです。

「【PythonとStanで学ぶ】仕組みが分かるベイズ統計学入門」をチェックする

以下がPystanで平均の差を推定するためのコードです。

データはシミュレーションデータを生成しています。

シミュレーションデータはサンプルサイズが50、変数xの平均が4、標準偏差が2、変数yの平均が8、標準偏差が3という設定で生成しました。

StanコードはRstanとまったく同じで、ここでは分散が異なる平均の差を推論しており、 delta_over = delta > 0 ? 1 : 0; で2つの変数の平均に差がある確率を求めています。

#パッケージ読み込み
import pandas as pd
import matplotlib.pyplot as plt
import pystan
%matplotlib inline

#データ読み込み
df = pd.read_excel("./data/data2.xlsx")

#stanコード
stan_model = """
data{
    int<lower=0> N;
    real<lower=0> x[N];
    real<lower=0> y[N];
}
parameters{
    real mu_x;
    real<lower=0> sigma_x;
    real mu_y;
    real<lower=0> sigma_y;
}
model{
    mu_x ~ normal(0,100);
    sigma_x ~ cauchy(0,5);
    mu_y ~ normal(0,100);
    sigma_y ~ cauchy(0,5);
    for(n in 1:N){
        x[n] ~ normal(mu_x, sigma_x);
    }
    for(n in 1:N){
        y[n] ~ normal(mu_y, sigma_y);
    }
}
generated quantities{
    real delta;
    real delta_over;
    delta = mu_x - mu_y;
    delta_over =   delta > 0 ? 1 : 0;
}
"""

#stanコードのコンパイル
sm = pystan.StanModel(model_code=stan_model)

#データの引き渡し
stan_data = {"N":df.shape[0], "x":df["x"], "y":df["y"]}

#Let's stan!
fit = sm.sampling(data = stan_data, iter=2000, warmup=500, chains = 4, seed=123)

#結果
fit

以下は、結果です。

delta_over をみるとわかるように、このシミュレーションデータでは平均に差がある確率は100%ということがわかりました。

Inference for Stan model: anon_model_08f58440e4765f58e92b2739ed38ccd1.
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_x         8.44  5.1e-3    0.4   7.66   8.18   8.44    8.7   9.23   6000    1.0
sigma_x      2.78  3.8e-3   0.29   2.27   2.57   2.75   2.95   3.42   6000    1.0
mu_y         4.06  3.6e-3   0.27   3.52   3.88   4.06   4.24   4.59   5797    1.0
sigma_y       1.9  2.6e-3    0.2   1.56   1.77   1.88   2.02   2.34   5750    1.0
delta        4.38  6.1e-3   0.47   3.48   4.06   4.37    4.7   5.32   6000    1.0
delta_over    1.0     0.0    0.0    1.0    1.0    1.0    1.0    1.0   6000    nan
lp__       -130.6    0.03   1.47 -134.3 -131.3 -130.2 -129.5 -128.8   3248    1.0

Samples were drawn using NUTS at Sun Jul 29 15:22:06 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で学ぶ】仕組みが分かるベイズ統計学入門」をチェックする

また、Stanコードの読み方、書き方は以下の資料がとてもわかりやすいです。

注意点としては、RでStanを使う方法ということと、Stanコードが前のバージョンのままというところですが、これを読み込めば基本的なところはほぼ理解できます。

また、以下の書籍もRでStanを使う方法ですが、Stanの理解には必読なのでPystanを使ってみたい人は入手しましょう。

まとめ

本記事では「はじめてPythonでStanを使いますが、平均の差の推定を行うにはどうしたらよいですか?」という疑問にお答えしました。

これは、医療系でよく使う統計解析技術の1つなので、各自の研究で使ってみてください。

その他、単回帰、重回帰、ロジスティック回帰、状態空間モデル、階層ベイズモデルについては、以下の動画で詳しく学べます。

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

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

続きを見る

30日間返金保証もあるので、試しに見てみるとよいでしょう。

なお、ここではPystanについて紹介しましすが、PyMC3という選択肢もあります。

関心ある人は以下の記事もどうぞです。

おすすめ
【初心者向け】「Pythonでベイズ統計」を学ぶ教材【1冊+1動画+α】

きょうごく本記事では「Pythonでベイズ統計を実行したいけど、難しくてよくわからない。オススメの学習教材はありますか?」という疑問にお答えします。 本記事のポイント 書籍で学ぶなら「Pythonによ ...

続きを見る

おすすめ記事一覧

1

きょうごく本記事では「現在の職場は人間関係が最悪で、労働条件も悪く、日々消耗しています。より好待遇な職場に移りたいです。おすすめの求人転職サイト・エージェントを教えてください」という疑問にお答えします ...

2

きょうごく本記事では「いまの職場は人間関係も労働条件も悪いので、より好待遇な職場で働くために転職したいです。医療福祉系の転職は転職サイトや転職エージェントを使ってもよいですか。おすすめはどこですか」と ...

3

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

4

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

5

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

6

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

7

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

-研究
-

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