オンライン教育プラットフォーム|Thriver Project
研究教育
PR

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

京極真
記事内にプロモーションを含む場合があります
京極真
京極真

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

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

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

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

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

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

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

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

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

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

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

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

https://www.kyougokumakoto.com/2018/08/python-stan-udemy-review.html

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

コードの読み方

各コードの読み方は上記に # でコメントを記しましたが、より詳細については紹介した動画でどうぞです。

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

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

Stan超初心者入門 from Hiroshi Shimizu

Stanコードの書き方 中級編 from Hiroshi Shimizu

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

created by Rinker
¥2,970 (2024/05/21 22:44:57時点 Amazon調べ-詳細)

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

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

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

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

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

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

https://www.kyougokumakoto.com/2018/08/python-pystan-correlation.html
https://www.kyougokumakoto.com/2018/07/python-bayes-material.html
著者紹介
京極 真
1976年大阪府生まれ。Ph.D、OT。Thriver Project代表。吉備国際大学ならびに同大学大学院・教授(役職:保健科学研究科長、(通信制)保健科学研究科長、人間科学部長、他)。首都大学東京大学院人間健康科学研究科博士後期課程・終了。『医療関係者のための信念対立解明アプローチ』『OCP・OFP・OBPで学ぶ作業療法実践の教科書』『作業で創るエビデンス』など著書・論文多数。
記事URLをコピーしました