Up

Stanによるベイズ分析の準備

 

データは、確率モデルで表されることが多い。データをで、確率モデルのパラメータをで表すと、確率モデルは、で表すことができる。

パラメータの確率分布を考え、で表す。

パラメータとデータの同時確率を考え、次式で与える。

 

 

ベイズ分析は、式（１）に基づいて分析が行われる（Gelman et al., 2014）。岡本（2019）には、ベイズ分析の入門レベルの解説がある。

 

ベイズ分析では、データが与えられたときのパラメータの条件付確率が求められる。ここで、式（１）より、

 

 

である。

 

データ が与えられたときの条件付確率は、Stanを用いればシミュレーションで生成した乱数により推定できる。図E1に、Stanスクリプトの例を示す。

Stanは、CmdStanPyモジュールを用いればPythonで利用できる。CmdStanPyは、WindowsではWSL（Windows Subsystem for Linux）にインストールしたLinuxにおいてMiniconda(Anaconda)を利用すると簡単にインストールすることができる。インストール後の使用も、基本的には簡単である。以下において、WindowsへのLinuxのインストール、Linuxの簡単な使用例、LinuxへのMiniconda(Anaconda)のインストール、CmdStanPyのインストール、CmdStanPyによるベイズ分析について説明する。これらは、他のウェブサイトで説明したものであるが、ここにまとめて編集したものである。

 

 

WindowsへのLinuxのインストール

 

Windowsの場合、WSL（Windows Subsystem for Linux）を用いて、WindowsにインストールされたLinux上にcondaコマンドでCmdStanPyを簡単にインストールすることができる。WSLによるWindowsへのLinuxのインストールについて説明する。

端末ソフトWindows PowerShellをスタートメニュから選んで起動する（図A1）。

 

図A1

 

端末PowerShellが起動して、端末の画面が表示されたら、コマンド

 

wsl --list --online

 

を実行してみる。「--」は「-」が２つであることに注意。端末入力において、パラメータの前の「-」は、１個の場合と２個の場合があるので注意。

上のコマンドを実行すると、現在インターネット経由でダウンロード・インストールできるLinuxのリストが表示される（図A2）。

 

図A2

 

上のリストにある「Ubuntu-24.04」をインストールするときは、コマンド

 

wsl --install Ubuntu-24.04

 

を実行する。「install」の前の「--」は「-」が2個である。

上のコマンドを実行すると、Ubuntu-24.04がダウンロードされ、インストールが始まる。

端末に表示される指示に従って、入力していけばよい。

最後に、パスワードを設定したら、コマンド

 

exit

 

を入力して、端末を閉じる（終了する）。コマンドexitは、１回目はインストール処理からのexit、２回目はPowershellからのexitと、２回必要である。

 

スタートメニュにおいてインストールしたUbuntu-24.04を選び、起動する（図A3）。立ち上げたスタートメニュに表示されないときは、「すべて」を選ぶなどして探す。

 

図A3

 

Ubuntuの端末が起動したら、まず更新コマンドを実行する（図A4）。

起動してしばらく待っても端末内にキーボードによる入力ができないときは、端末内の画面をマウスでクリックする。何回かクリックすると、端末内でカーソルが点滅して、点滅位置へのキーボード入力が可能になる。

 

図A4

 

更新後は、コマンドexitで閉じて、再起動する。

 

コマンド

 

sudo apt update && sudo apt upgrade

 

によるUbuntuの更新は、時々行っておく。例えば、1週間に１回ぐらい。

 

 

 

Linuxの簡単な使用例

 

WSLのUbuntuで扱うファイルやフォルダの操作は、Windowsのエクスプローラで出来る。WSLのUbuntuの端末でのカレントディレクトリ操作などは、端末でのコマンドで行う。ディレクトリの操作とドライブの操作について説明する。

 

 

 

フォルダ（ディレクトリ）の移動

 

 

カレントディレクトリ（フォルダ）の移動は、コマンドcdで行うことができる。コマンド

 

cd 移動先フォルダ

 

を実行すれば、移動先フォルダに移動できる。

一番上のディレクトリは「/」で表されるので、

 

cd  /

 

を実行すると、一番上のディレクトリにカレントディレクトリが移動する。移動後、カレントディレクトリの内容をコマンドlsで表示すると図B1のようになる。

 

図B1

 

Windowsのディレクトリは、ディレクトリmntにある。

ディレクトリmntに移動して、中を見てみる。

 

図B2

 

カレントディレクトリは、コマンドラインの先頭の「＄」の左側に表示されていることに注意。

ディレクトリc、すなわち、WindowsのCドライブに移動して、コマンドlsで中を見てみる（図B3）。Windo0wのCドライブは、Linuxでは小文字cで表されている。

 

図B3

 

WindowsのドライブCの内容が表示されている。

 

なお、カレントディレクトリは「./」で、カレントディレクトリの１つ上のディレクトリは「../」で表される。

 

ディレクトリ名やファイル名を入力するとき、途中まで入力した後、Tabキーを押すと、残りが補完される。

 

図B4

 

上の例では、「Fuj」と入力されている。この状態で、Tabキーを押すと、図B5のように「Fujitsu/」と補完されている。

 

図B5

 

「Fijitsu/」の最後の「/」は、ディレクトリなので付けられているもので無視してよい。

 

このTabキーによるファイル名の補完は、途中まで同じ文字列のものがあれば、その途中までが表示される。その後はキー入力して補えばよい。途中まで補った後、その後の部分に同じものが無ければ、Tabキーを押せば、残りが補完される。

 

また、コマンド入力時に、上向き矢印キーを押すと、１回押すことに１つ前に入力したコマンド、下向き矢印キーを押すと現在表示されているコマンドの次に入力したコマンドが表示される。

 

 

 

WSL上のLinuxのドライブ

 

 

WindowsのディレクトリへのUbuntuのカレントディレクトリの移動は、コマンドcd /mnt/で行うことができる。

 

図B6

 

コマンドlsで表示されたリストに小文字cで表されているのが、WindowsのCドライブである。

使用PCには、USBメモリが２台接続されているが、表示されていない（図B7）。すなわち、Ubuntuからは見えていない。

 

図B7

 

USBメモリをUbuntuから使用できるようにするために、マウントする。

マウントは次の形式でmountコマンドを実行する。

 

sudo  mount  –t  drvfs  Windowsのドライブ名  ubuntuのマウント先

 

先頭の「sudo」は、以下のコマンドを管理者権限で実行することを意味する。オプション –t  drvfsが必要である。Windowsのドライブ名は、例えば、E:である。E:を指定すると、ドライブE:がマウントされるが、ドライブ名は小文字を用いる。E:ドライブであればe:である。

Ubuntuのマウント先は、例えば、ディレクトリ/mntにusb_memと作成しておく。

実行例を、図B8に示す。

 

図B8

 

アンマウントする（Ubuntuから切り離す）ときは、umountコマンドを用いる。図９に例を示す。

 

図B9

 

アンマウントしても、ubuntuのディレクトリusb_memは残っている。このディレクトリも削除するときは、rmdirコマンドを用いる。図B10 に、例を示す。

 

図B10

 

 

 

 

LinuxへのMiniconda(Anaconda)のインストール

 

LinuxにMiniconda（Anaconda）をインストールしておくとCmdStanPyなどのモジュール（ライブラリ）の準備（インストール）が簡単になる。MinicondaはAnacondの簡易版である。以下では、Minicondaについて説明するが、Anacondaの場合も同様である。

 

Minicondaのファイルをダウンロードするために

 

https://www.anaconda.com/docs/getting-started/miniconda/main

 

あるいは

 

https://www.anaconda.com/download

 

からMiniconda（Anaconda）のウェブサイトに入る。ウェブサイトの動作は、アクセスした時期によって変わっていくので、ここでの説明は参考程度と了解して頂きたい。

 

ウェブサイトの指示に従って、Minicondaのダウンロードのウェブサイトに進む。

 

https://www.anaconda.com/download/success

 

下図のような画面になる。

 

図C1

 

「Linux」タブを選ぶ。

 

図C2

 

インストールするPCのCPUに合わせてダウンロードするファイルを選ぶ。上は、IntelのCPUの場合である。

図C3には、ダウンロードしたMinicondaのファイルがWindowsのフォルダC:\WorkSpace\Download\tempにあることが示されている。

 

図C3

 

Ubuntuのカレントディレクトリを、ダウンロードしたファイルの置かれているフォルダ、図C3の場合はWindowsのフォルダC:\WorkSpace\Download\tempに移す。

Windowsのドライブは、Linuxではディレクトリ/mnt（/は、一番上のディレクトリであることを示す）にあり、WindowsのドライブCは小文字ｃで表されるので、次のコマンドを実行する。

 

cd /mnt/c/WorkSpace/Download/temp/

 

コマンドlsで移動後のディレクトリの内容を示すことができる。

 

yasuh@pcyo:/mnt/c/WorkSpace/Download/temp$ ls

Miniconda3-latest-Linux-x86_64.sh

 

ダウンロードしたファイルMiniconda3-latest-Linux-x86_64.shをコマンドbashで実行する。

 

bash Miniconda3-latest-Linux-x86_64.sh

 

上記の一連の端末上での操作を図C4に示す。

 

図C4

 

コマンドbashを実行するとインストールが始まり、Enterキーを押すことが求められる（図C5）。

 

図C5

 

Enterキーを押すと、表示が進み、［yes/no］の選択が求められるので（図Ｃ６）、「yes」と入力する。

 

図C6

 

「yes」と答えて、画面が進むと、インストール先の選択が求められる（図C7）。

 

図C7

 

Enterキーを押すと、インストールが進む。

インストールの最後に図C8のような問い合わせが表示される。

 

図C8

 

「yes」と答えておくと、後の作業が楽である。

 

図C9

 

インストールが終了したら、コマンドexitを実行して、一旦、Ubuntuを終了する。

 

Ubuntuを再起動する。

 

図C10

 

コマンドラインの先頭に「(base)」とあるのは、Minicondaの仮想環境baseの端末であることを表している。Minicondaのコマンドcondaが使えるが、Ubuntu（Linux）のコマンドも使える。

 

再起動後、まずMinicondaのアップデートをコマンド

 

conda update –n base conda

 

を実行して行う。

 

図C11

 

パラメータ「-n base」は、baseという名前の仮想環境のアップデートであることを示し、続くパラメータ「conda」はcondaのアップデートであることを示す。

 

続いて、コマンド

 

conda update  --all

 

を実行してすべてのパッケージのアップデートを行う（図C12）。

 

図C12

 

パラメータ「--all」は、現在のアクティブな仮想環境、コマンドラインの先頭に示されている仮想環境で図C12の場合base、内のすべてのパッケージをアップデートすることを表す。

 

 

 

CmdStanPyのインストール

 

 

CmdStanPyのインストールの前に、g++が使えることを確認しておく。

Ubuntuの端末において、コマンドg++を入力したときに、入力ファイル名が指定されていませんというようなメッセージが返ってきたらOKである。コマンドg++はありませんというようなメッセージが返ってきたときは、メッセージとともにインストールコマンドも示されているので、それを実行する。簡単なコマンドでg++のインストールができる。

例えば、コマンドは

 

sudo apt install g++

 

である。先頭の「sudo」は、コマンドを管理者権限で実行することを表す。

 

CmdStanPyを、MiniCondaの仮想環境としてmystanという名前の仮想環境を用意して、mystanにCmdStanPyをインストールするときは、次のコマンドを実行すればよい。

 

conda create -n mystan -c conda-forge cmdstanpy

 

図D1

 

パラメータ「-n mystan」は、仮想環境を「mystan」という名前で作成することを示す。「mystan」は、他の名前でもよい。パラメータ「-c conda-forge」は、仮想環境の作成に必要な情報を得る場所を示す。パラメータ「cmdstanpy」は、cmdstanpyの仮想環境を作成することを示す。

 

仮想環境を作成したら、更新を行う。

まず、作成した仮想環境mystanをアクティブにする。

 

conda activate mystan

 

アクティブになった仮想環境mystanにおいて、パッケージの更新を行う。

 

conda update --all

 

図D2に実行例を示す。コマンドラインの先頭にアクティブな仮想環境名が表示されていることに注意。

 

図Ｄ２

 

CmdStanPyについては、このウェブサイトも参照されたい。

 

 

 

CmdStanPyによるベイズ分析

 

 

まず、簡単な例を示し、次に、分析結果のグラフ表示例を示す。なお、サンプルスクリプトは、このウェブサイトにも用意した。

 

Stanスクリプト（binomial.stan）（図E1）

 

図E1　２項分布モデルのStanスクリプト

 

を用いて、ベイズ分析を行うPythonスクリプト（binomial.py）（図E2）を実行する。

 

図E2

 

これらのスクリプトファイルは、標準Pythonのテキストエディタを用いて作成したが、テキストエディタであれば何でもよい。標準Pythonのインストールについては、このウェブサイトで説明している。

 

StanスクリプトとPythonスクリプトは同じフォルダ内に置いておく。図E1と図E2の場合、フォルダC:\WorkSpace\Samples内に置かれている（図E3）。

 

図E3

 

仮想環境mystanをアクティブにして、カレントディレクトリをC:\WorkSpace\Samplesに移す。

これらの操作は、図E4に示す通りである。

 

図E4

 

Pythonスクリプトをpythonコマンドで実行する。

 

python binomial.py

 

上のコマンドを実行すると、次のコード

 

model = CmdStanModel(stan_file='binomial.stan')

 

により、Stanスクリプト（binomial.stan）のコンパイルが始まる。コンパイルには、多少時間が掛かる。

コンパイルの結果は、modelで表される。

このmodelで示されるモデルによるMCMCサンプリングは、次のコードで行われ、サンプリングの結果がfitで示されている。

 

fit = model.sample(data={'k':7, 'n':10})

 

分析対象のデータは、パラメータdata={'k':7, 'n':10})で与えられている。

MCMCサンプリングの統計が、次のコードで出力されている。

 

print(fit.summary())

 

上記の操作と端末画面を図E5に示す。

 

図E5

 

 

 

MCMCサンプルのグラフ表示を行ってみる。

 

 

サンプルスクリプトは、次のリストE1とリストE2である。

 

 

リストE1　Stanスクリプト（xtable.stan）

 

functions {

    vector f(vector p_fm, vector pf, vector pm) {

        vector[4] v;

        v[1] = p_fm[1]*pf[1];

        v[2] = p_fm[1]*pf[2];

        v[3] = p_fm[2]*pm[1];

        v[4] = p_fm[2]*pm[2];

        return v;

    }

}

data {

    int o11;   int o12;

    int o21;   int o22;

}

transformed data {

    array[4] int y;

    y[1] = o11;

    y[2] = o21;

    y[3] = o12;

    y[4] = o22;

}

parameters {

    simplex[2] p_fm;

    simplex[2] pf;

    simplex[2] pm;

}

model {

    p_fm ~ dirichlet([1,1]);

    pf ~ dirichlet([1,1]);

    pm ~ dirichlet([1,1]);

    y ~ multinomial(f(p_fm, pf, pm));

}

 

 

 

リストE2　Pythonスクリプト（xtabel.py）

 

from cmdstanpy import CmdStanModel

import numpy as np

import pandas as pd

import arviz as az

import matplotlib.pyplot as plt

import seaborn as sb

#

#    Spiegelhalter(2019), Table 10.2

#        Female  Male

data = [[5,      17],   #   Left arm on top

        [9,      23]]   #   Right arm on top

 

Data = {'o11':data[0][0], 'o12':data[0][1],

        'o21':data[1][0], 'o22':data[1][1]}

 

model = CmdStanModel(stan_file='xtable.stan')

fit = model.sample(data=Data)

 

infdata = az.from_cmdstanpy(fit)  # Arviz InfereceData型へ変換

summary = az.summary(infdata)

print('\nSummary...\n', summary)

 

az.plot_trace(infdata) 

plt.tight_layout()

plt.show()

 

dframe = fit.draws_pd()           # Pandas DataFrame型へ変換

print(dframe.keys())

 

sb.kdeplot(dframe['p_fm[1]'], label='p_fm[1]')

sb.kdeplot(dframe['p_fm[2]'], label='p_fm[2]')

plt.legend()

plt.show()

 

sb.kdeplot(dframe['pf[1]'] - dframe['pm[1]'])

plt.xlabel('pf[1]-pm[1]')

prop = np.mean(dframe['pf[1]'] - dframe['pm[1]'] > 0)

plt.title(f'P(pf[1] > pm[1])={prop:.3f}', fontsize=14)

plt.show()

 

 

 

上のスクリプトは、次のクロス表を分析するものである。

 

 

女性

男性

左腕が上

o11

o12

右腕が上

o21

o22

 

受講生を女性と男性に分けて、それぞれの腕組みの腕がどちらが上であったかが示されている（Spiegelhalter, 2019, Table 10.2）。

ベイズ分析のための確率モデルは、以下のように設定されている。

 

受講者が女性である確率を、p_fm[1]で、男性である確率をp_fm[2]で表す。

女性の受講生において、左腕が上である確率をpf[1]で、右腕が上である確率をpf[2]で表す。

男性の受講生において、左腕が上である確率をpm[1]で、右腕が上である確率をpm[2]で表す。

 

このとき、クロス表の度数を多項分布に従うとしてモデルを構成している。

 

上のスクリプトの実行において、いくつかのライブラリが追加で必要であるが、現在（2025.11）のCmdStanPyのPythonとの相性には問題がある。

Pythonを起動してバージョンを確認すると、図Ｂ６のようにPython 3.14.0と表示される。

 

図E６

 

このPythonを１つ前のバージョンであるPython 3.13のものに置き換えるために、図E７に示すPythonのインストールを行う。

 

図E７

 

この後、ライブラリArvizをバージョン0.21.0を指定してインストールする（図E８）。

 

図E８

 

続いて、ライブラリmatplotlibとseabornをインストールする（図E９）。

 

図E９

 

 

コマンドlsを実行して、カレントディレクトリにリストE1のファイル（xtable.stan）とリストE2のファイル（xtable.py）があることを確認して、コマンド

 

python xtable.py

 

を実行する（図E10）。

 

図E10

 

実行すると、Stanスクリプトのコンパイルが始まるが、コンパイルには多少時間が掛かる。コンパイル後、MCMCサンプリングが行われ、終了するとトレース図が表示される（図E11）。

 

図E11

 

これは、次のスクリプトによるものである。

 

infdata = az.from_cmdstanpy(fit)  # Arviz InfereceData型へ変換

summary = az.summary(infdata)

print('\nSummary...\n', summary)

 

az.plot_trace(infdata) 

plt.tight_layout()

plt.show()

 

 

MCMCサンプルがaz.from_cmdstanpy(fit)によりArvizのInferenceData型に変換され、az.plot_trace(infdata)によってトレース図が描かれている。

図E11のグラフより、pfとpmの分布はほとんど同じであることが分かる。

p_fmの事後分布は、p_fm[1]とp_fm[2]が明確に分かれている。

図E11のWindowを閉じると、図E12のグラフが表示される。

 

図E12

 

図E12のグラフは、次のスクリプトにより描かれている。

 

 

dframe = fit.draws_pd()           # Pandas DataFrame型へ変換

print(dframe.keys())

 

sb.kdeplot(dframe['p_fm[1]'], label='p_fm[1]')

sb.kdeplot(dframe['p_fm[2]'], label='p_fm[2]')

plt.legend()

plt.show()

 

 

MCMCサンプルを、dframe = fit.draws_pd()によりPandasのDataFrame型に変換している。

DataFrameのキーは、print(dframe.keys())により端末に以下の様に書き出されている。

 

 

Index(['chain__', 'iter__', 'draw__', 'lp__', 'accept_stat__', 'stepsize__',

       'treedepth__', 'n_leapfrog__', 'divergent__', 'energy__', 'p_fm[1]',

       'p_fm[2]', 'pf[1]', 'pf[2]', 'pm[1]', 'pm[2]'],

      dtype='object')

 

 

キーは、変数名に対応する文字列で表されていることが分かる。

図E12のWindowを閉じると、図E13のグラフが表示される。

 

図E13

 

確率の差pf[1]-pm[1]の事後分布のグラフと、P(pf[1]>pm[1])の事後確率が表示されている。

事後確率は、次のコードで算出している。

 

prop = np.mean(dframe['pf[1]'] - dframe['pm[1]'] > 0)

 

事後確率＝0.350で、ほぼ原点０を中心に跨った分布である。２つの確率pf[1]とpm[1]は、ほぼ同じであると考えられる。

 

図E13のWindowを閉じると、スクリプトの実行終了である。

 

 

 

参考文献

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, 3rd. Ed. CRC Press.  https://sites.stat.columbia.edu/gelman/book/

岡本安晴（2019）　いまさら聞けないPythonでデータ分析．丸善出版

 

 

Up