本記事のポイント
- Pystanで平均値の差を推定する方法がわかります
- ぼく自身のPystanの練習をかねた備忘録ですｗ
Pythonでベイズ統計を勉強する理由
ぼくは普段、Mplusの他に、Rでstan、JAGS、brmsなどいろいろ使用しています。
なのに、Pystan（Pythonでベイズ統計）を練習しはじめた理由は、話題のPythonを勉強したいという動機からです（笑）。
ややミーハーな理由ですいません。
-
-
【断言】新しいことに挑戦したい人にオススメの4つのテーマ【2020年】
きょうごく本記事では「最近、行動の重要性をいろんなところで聞きます。私も行動は重要だと思いますので、新しいことに挑戦したい！けど、何に挑戦したらよいの？」という疑問にお答えします 本記事のポイント こ ...
続きを見る
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; で２つの変数の平均に差がある確率を求めています。
#パッケージ読み込み
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によ ...
続きを見る
