【結論】文章力がない原因とその対策【研究者が語る】
京極真
京極真のWebマガジン
【初心者向け】Pythonでベイズ統計学入門【平均に差がある確率を推定する】
京極真
記事内にプロモーションを含む場合があります
本記事では「はじめてPythonでStanを使いますが、平均の差の推定を行うにはどうしたらよいですか？」という疑問にお答えします
- Pystanで平均値の差を推定する方法がわかります
- ぼく自身のPystanの練習をかねた備忘録ですｗ
Contents
Pythonでベイズ統計を勉強する理由
ぼくは普段、Mplusの他に、Rでstan、JAGS、brmsなどいろいろ使用しています。
なのに、Pystan（Pythonでベイズ統計）を練習しはじめた理由は、話題のPythonを勉強したいという動機からです（笑）。
ややミーハーな理由ですいません。
Pystanの勉強はudemyの動画教材でやっています。
この動画は、ベイズ統計の基礎、MCMCの原理、グラフィカルモデルの基礎、さまざまな統計モデルの基礎からStanコードへの落とし込み方、結果の解釈などなど、めちゃくちゃわかりやすく解説してくれています。
学習効率がとてもよいので、これからベイズ統計を学びたい人や、ぼくのようなエンドユーザーでちょっとかじっているけどももっと腕を上げたい人には特にオススメです。
がっつりレビューを書いたので、ぜひお読みいただけたらと思います。
マジでできることが広がりますので、コスパ最高ですよ。
30日間返金保証があるので、その間なら実質無料ですから、「自分にはハードルが高いかも」という人でもとりあえず試すとよいです。
Pystanで平均の差の推定
課題：2つの変数の平均に差がある確率を明らかにしたい
ここではかりに「2つの変数の平均に差がある確率を推定する」という課題を設定します。
従来の統計の場合、帰無仮説「2つの変数の平均に差がない」と対立仮説「2つの変数の平均に差がある」を設定します。
そして、帰無仮説のもとで手元にあるデータが得られる確率を計算し、それが5％未満であれば帰無仮説を棄却することになります。
でもこの場合、帰無仮説が間違っているとも言えないし、かりに対立仮説を採用するにしても平均に差がある確率がどの程度あるのかについて何も言えません。
他方、ベイズ統計を活用したら2つの変数の平均に差がある確率を直に推定できるので、「AとBの平均に差がある確率は92％です」などのように知りたいことを知ることができます。
Pystanで平均の差を推定する
2つの変数の平均に差がある確率を推定するために、ここではPystanを用いています。
Pystanの使い方の詳細は動画にゆずりますが、慣れればそんなに難しくないです。
以下がPystanで平均の差を推定するためのコードです。
データはシミュレーションデータを生成しています。
シミュレーションデータはサンプルサイズが50、変数xの平均が4、標準偏差が2、変数yの平均が8、標準偏差が3という設定で生成しました。
StanコードはRstanとまったく同じで、ここでは分散が異なる平均の差を推論しており、 delta_over = delta > 0 ? 1 : 0; で２つの変数の平均に差がある確率を求めています。
#パッケージ読み込み
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).
コードの読み方
各コードの読み方は上記に # でコメントを記しましたが、より詳細については紹介した動画でどうぞです。
また、Stanコードの読み方、書き方は以下の資料がとてもわかりやすいです。
注意点としては、RでStanを使う方法ということと、Stanコードが前のバージョンのままというところですが、これを読み込めば基本的なところはほぼ理解できます。
Stan超初心者入門 from Hiroshi Shimizu
Stanコードの書き方 中級編 from Hiroshi Shimizu
また、以下の書籍もRでStanを使う方法ですが、Stanの理解には必読なのでPystanを使ってみたい人は入手しましょう。
created by Rinker
¥3,234
(2023/12/06 19:12:44時点 Amazon調べ-詳細)
まとめ：【初心者向け】Pythonでベイズ統計学入門【平均に差がある確率を推定する】
本記事では「はじめてPythonでStanを使いますが、平均の差の推定を行うにはどうしたらよいですか？」という疑問にお答えしました。
これは、医療系でよく使う統計解析技術の1つなので、各自の研究で使ってみてください。
その他、単回帰、重回帰、ロジスティック回帰、状態空間モデル、階層ベイズモデルについては、以下の動画で詳しく学べます。
30日間返金保証もあるので、試しに見てみるとよいでしょう。
なお、ここではPystanについて紹介しましすが、PyMC3という選択肢もあります。
著者紹介
1976年大阪府生まれ。Ph.D、OT。Thriver Project代表。吉備国際大学ならびに同大学大学院・教授。作業療法学科長、保健科学研究科長、（通信制）保健科学研究科長。首都大学東京大学院人間健康科学研究科博士後期課程・終了。『医療関係者のための信念対立解明アプローチ』『作業療法リーズニングの教科書』『作業で創るエビデンス』など著書・論文多数。
