Stan (PyStan) の試用

PyStan（Python用Stan）の試用

注意（2023.09）：PyStan 3（PyStanの最新版）は、WindowsではWindowsにUbuntuをインストールすることにより使用できる。UbuntuのWindowsへのインストールについては

http://y-okamoto-psy1949.la.coocan.jp/Python/tips/ubuntu_wsl/

を参照されたい。PyStan 3の使用例については

http://y-okamoto-psy1949.la.coocan.jp/Python/tips/TryPyStan3/

などに用意した。

注意：　現在（2022.07.29）、Windows上では、anacondaでのpystanはMVSCコンパイラが選ばれて使えない。しかし、pystanとPythonのバージョンを選べば、使えることが確認できた。Python 3.7の仮想環境においてpystan 2.17.1.0をインストールしてみると問題なくpystanが使えた。サンプリングの実行時に、パラメータ「n_jobs=1」の設定が必要であることに注意。

Pystan 2.17.1.0はコマンド

conda install pystan==2.17.1.0

でインストールできる。等号=が２つ==であることに注意。

仮想環境の作成は、ここで説明している。

Linux (Ubuntu 22, CentOS 8 Stream)の場合、AnacondaのPython 3.9仮想環境においてPyStan 2.19.1.1でエラーが出た（2023.03.21）。PyStan 2.19.1.1をPython 3.8の仮想環境で試みると問題は無いようであった。

まず、PyStanのインストールについて説明する。その後、Stan (PyStan)の使い方について説明する。

PyStan (Stan)のインストール

Python用StanとしてPyStanがある。以下に、PyStanの基礎的なことについて説明するが、PyStanについては＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞でも説明している。

PythonにおいてStanを使うためには、ライブラリPyStanをインストールする必要がある。PyStanの実行時にC++コンパイラが使用されるが、PyStanとC++コンパイラの相性の問題がある。この相性の問題は、AnacondaにおいてPyStanをインストールすると、自動的にインストールするPyStanと相性の良いC++コンパイラの組み合わせでインストールされる。

Anacondaのインストールと使用法はこのウエブサイトで説明している。

Anacondaがインストールされたら、自動的に作成されている仮想環境baseにおいて次のコマンドで更新しておく。

conda update anaconda

続いて次のコマンドによる更新も行っておく。

conda update --all

なお、Windowsの場合、標準のPythonもhttps://www.python.org/からインストールしておくと、Python IDLEのエディタが使えるので便利である。Python IDLEのエディタでスクリプトを作成して、実行はAnacondaの仮想環境で行えばよい。Linuxの場合は、geditなど適当なエディタでスクリプトを作成すればよい。Windows用Pythonのインストールについては、このウェブサイトで説明している。

仮想環境の作成

Anacondaのインストール後、使用するPyStan用の仮想環境を構築した方が良い。PyStanが更新されて互換性の問題が出たときに、PyStanの旧版用と新版用の仮想環境を用意しておくと旧版を使用し続けながら新版も試すことができる。仮想環境は、例えば、次のコマンドで「py39」という名前の「Python 3.9の最新版」がインストールされた仮想環境が作成される。

conda create –n py39 python=3.9

インストール後、コマンド

conda activate py39

を実行すると、仮想環境py39がアクティブになる。詳しくは、このウェブサイトを参照されたい。

仮想環境py39がアクティブな状態での端末（コンソール）のコマンドプロンプトは以下のようになっている。

(base) PS C:\Users\Yasuharu> conda activate py39

(py39) PS C:\Users\Yasuharu>

コマンドプロンプトの先頭に表示されている（base）や（py39）はアクティブな仮想環境を表している。Anaconda起動時は仮想環境baseがアクティブである。この仮想環境baseは、Anacondaのインストール時に自動的に作成されている。

仮想環境py39がアクティブな状態で、pystanをインストールするが、データ分析でよく使うライブラリnumpy、pandas、matplotlib、seaborn、scipyなどもインストールしてくと便利である。これらのインストールは次のコマンドでまとめてインストールできる。

(py39) PS C:\Users\Yasuharu> conda install pystan numpy pandas matplotlib seaborn scipy

これでpystanが使用できるようになった。

実行したいプログラムのあるディレクトリに移ってpythonコマンドを実行すればよい。

例えば、以下のように行えばよい。

(base) PS C:\Users\Yasuharu> conda activate py39

(py39) PS C:\Users\Yasuharu> cd s:\Samples

(py39) PS S:\Samples> ls

ディレクトリ: S:\Samples

Mode LastWriteTime Length Name

---- ------------- ------ ----

-a---- 2021/07/30 10:17 234 StanSample.py

(py39) PS S:\Samples> cat StanSample.py

import pystan

StanScript = """

data {int N; int k;}

parameters {real<lower=0.0, upper=1.0> p;}

model {k ~ binomial(N, p);}

"""

fit = pystan.stan(model_code=StanScript, data={'N':10, 'k':9}, n_jobs=1)

print(fit)

(py39) PS S:\Samples> python StanSample.py

上のコマンドpython StanSample.pyの実行で、以下の出力を得る。

Inference for Stan model: anon_model_3e3224cef22eeb31257b8976794ce02d.

4 chains, each with iter=2000; warmup=1000; thin=1;

post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

p 0.83 2.7e-3 0.1 0.59 0.77 0.85 0.91 0.98 1474 1.0

lp__ -5.96 0.03 0.8 -8.3 -6.14 -5.65 -5.46 -5.41 920 1.0

Samples were drawn using NUTS at Fri Jul 30 10:21:51 2021.

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

(py39) PS S:\Samples>

Windowsの場合、ディレクトリを移るときは端末（コンソール）の選択に気を付ける必要がある。

WindowsにおいてAnacondaをインストールすると端末がAnaconda Powershell PromptとAnaconda Promptの２つ用意される。Anaconda Promptでは、Cドライブ以外のドライブには移れない。コマンドcdでCドライブ以外のドライブにも移りたいときは、Anaconda Powershell Promptを選ぶとよい。

端末の終了はexitコマンドで出来る。

(py39) PS S:\Samples> exit

とコマンドを実行すれば端末は閉じられる。

Stanの使用

PyStanはPython用のStanインターフェースである。PyStanのドキュメントは

https://pystan2.readthedocs.io/en/latest/

から得ることができる。Stanのドキュメントは

https://mc-stan.org/users/documentation/

から得ることができる。

PyStanとStanの説明は＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。

Stanはベイズ分析のためのソフトであるが、ベイズ分析についてはこのウェブサイトで説明している。

Stanスクリプトは、関数宣言部、データ宣言部、データ変換部、パラメータ宣言部、パラメータ変換部、モデル宣言部、サンプル生成部から構成されるが、これらの中から必要なものを選べばよい。Stanスクリプトの文法は、C/C++のものと似ているので、C/C++を知っていると使い易い。

ベイズ分析の基本モデルは、データとパラメータの同時分布

であるが、これをStanスクリプトで表すと、事後分布からのサンプルを得ることができる。

具体的に次のデータについて考える。

data = [57, 30, 34, 59, 63, 26, 60, 58, 34, 56]

上のデータは、1次元の長さNの配列として表されるので、Stanスクリプトのデータ宣言部を以下のように用意する。

data {

int N;

real Data[N];

}

NとData[N]の具体的なデータ値は、Stanスクリプト実行時に指定する。

上のデータに対して1変量正規分布モデルによる分析を行う。

１変量正規分布は次式

で表され、平均値パラメータと標準偏差パラメータをもつ。すなわち

（２）　　　　

である。

基本モデル（１）を構成するためにはモデル（２）の他にパラメータとの確率分布の設定が必要である。設定を行わなかった場合は、一様分布が設定されるが、設定した方がサンプリングが安定すると説明されている。パラメータに対する制約を出来るだけ弱くするためにパラメータには分散の大きい正規分布、パラメータには分散の大きい指数分布を個別に設定すると、パラメータ全体に対する確率分布はそれらの積として与えられる。すなわち、

となる。

設定（２）と（３）をStanスクリプトで次のように表す。

parameters {

real mu;

real<lower = 0.0> sgm;

}

model {

mu ~ normal(0.0, 1000.0);

sgm ~ exponential(0.001);

for (i in 1:N) {

Data[i] ~ normal(mu, sgm);

}

配列の添え字が１から始まっていることに注意。

以上の設定を１つのStanスクリプトとしてまとめるとリスト１となる。

リスト１　１変量正規分布モデルのStanスクリプト（ファイル名UniNormStan.stan）

data {

int N;

real Data[N];

}

parameters {

real mu;

real<lower = 0.0> sgm;

}

model {

mu ~ normal(0.0, 1000.0);

sgm ~ exponential(0.001);

for (i in 1:N) {

Data[i] ~ normal(mu, sgm);

}

Stanスクリプトはコンパイルされてから実行される。リスト１のStanスクリプトのファイルUniNormStan.stanのコンパイルは次のコードで行うことができる。

import pystan

sm = pystan.StanModel(file = 'UniNormStan.stan')

コンパイルの結果は変数smで指し示されている。

このコンパイルされたものを使って、上のデータに対してMCMCサンプリングによる事後分布のサンプルを次のコードにより求めることができる。

data = [57, 30, 34, 59, 63, 26, 60, 58, 34, 56]

fit = sm.sampling(data = {'Data':data, 'N':len(data)}, iter = 4000,

n_jobs = 1) # The parameter n_jobs = 1 is for Windows

関数samplingにおけるパラメータn_jobs=1はWindowsでの実行のときは必要である。Linuxではこのパラメータは不要である。Linuxにおいてパラメータn_jobs=1が無いときは、使用しているパソコンのCPUのコア数に応じてMCMCサンプリングの複数のチェインが分割並行して実行され、n_jobs=1が指定されたときより計算が速くなる。

リスト１のStanスクリプトを使ったPythonスクリプトの完全なリストをリスト２に示す。

リスト２　Pythonスクリプト（ファイル名UniNormalTwoFiles.py）

import pystan

sm = pystan.StanModel(file = 'UniNormStan.stan')

data = [57, 30, 34, 59, 63, 26, 60, 58, 34, 56]

fit = sm.sampling(data = {'Data':data, 'N':len(data)}, iter = 4000,

n_jobs = 1) # The parameter n_jobs = 1 is for Windows

print(fit)

リスト２のPythonスクリプトは次のコマンドで実行できる。

(py39) PS D:\StanSamples> python UniNormalTwoFiles.py

実行結果は以下のようである。

Inference for Stan model: anon_model_1ca69dd893129cb3ba4fd44d6085fc80.

4 chains, each with iter=4000; warmup=2000; thin=1;

post-warmup draws per chain=2000, total post-warmup draws=8000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

mu 47.84 0.09 5.61 36.58 44.33 47.78 51.29 59.5 3922 1.0

sgm 17.1 0.08 4.85 10.51 13.75 16.18 19.32 29.2 3585 1.0

lp__ -29.78 0.02 1.08 -32.68 -30.21 -29.44 -29.0 -28.71 2687 1.0

Samples were drawn using NUTS at Fri Jul 30 14:27:17 2021.

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

(py39) PS D:\StanSamples>

リスト２のPythonスクリプトではStanスクリプトは別のファイル（リスト１）として用意されている。StanスクリプトをPythonスクリプト内に文字列として用意することができる。例をリスト３に示す。

リスト３　Stanスクリプトを文字列として含むPythonスクリプト

import pystan

UniVarNormStan = """

data {

int N;

real Data[N];

}

parameters {

real mu;

real<lower = 0.0> sgm;

}

model {

mu ~ normal(0.0, 1000.0);

sgm ~ exponential(0.001);

for (i in 1:N) {

Data[i] ~ normal(mu, sgm);

}

"""

sm = pystan.StanModel(model_code = UniVarNormStan)

data = [57, 30, 34, 59, 63, 26, 60, 58, 34, 56]

fit = sm.sampling(data = {'Data':data, 'N':len(data)}, iter = 4000,

n_jobs = 1) # The parameter n_jobs = 1 is for Windows

print(fit)

統計分析をStanを用いて行っている例をこのウェブサイトの後半に用意している。

Home