CmdStanPyの試運転
先日CmdStanPyをインストールしてみました。（詳細は下記参照）
publicjournal.hatenablog.com
RのStanやPyStanと結構違うので、どういうものかと一通り使ってみました。ひとまず“Hello, World”を参照にしながら動作確認をしました。
同じような例題が書籍「Pythonによるベイズ統計モデリング」のP48～57にあるので、同書と同じことをCmdStanPyを使ってやってみたいと思います。
例題
こんな例題を考えます。
コイン投げ問題を考える。コインを10回投げて裏と表を記録する。表を１、裏を０として記録したところ、[1,0,1,0,0,0,0,0,1,0]となった。このデータをもとに、コインが偏っているかを判断する。コインの偏りを調べるために、次のモデルを用いる。
プログラムの全貌
プログラムは次の通り。
Stanのファイル
Stanのファイルは次の通り
解説
統計モデルの作成
上記のStanファイルにそって統計モデルを作成します。
モデルオブジェクトの生成
先ほどのStanファイルを読み込みます。
PyStanならばpystan.stan(model_code=XXX, data=XXXX, iter=1000, chains=4)などと書いていましたが、CmdStanPyならばStanファイルを指定して次の通りになります。
### Stanファイルを読み込んでオブジェクトを生成する stan_file_path = 'sample_model.stan' model = CmdStanModel(stan_file=stan_file_path)
model_codeでモデルをベタ書きできるのかは調査中です。
サンプリングの実行
sample()でハミルトニアンモンテカルロ法のノーUターンサンプラが走るそうです。
### MCMCを実行して事後分布をサンプリングする iter_sampling = 1000 # サンプリングの数.デフォルト1000 chains = 4 # 並列サンプリングの数．デフォルト4 # データの作成 front_and_back_list = [1,0,1,0,0,0,0,0,1,0] data = { "N" : len(front_and_back_list), "y" : front_and_back_list } # 実行 fit_sm = model.sample(data=data, iter_sampling=iter_sampling, iter_warmup = iter_warmup, chains = chains, seed=1234)
sample()の引数として、よく使うのは次のものかと思います。
data：データの辞書を渡す。これは他のStanと同じ。iter_sampling：サンプル数を指定する。他のStanのiterと同じっぽい。iter_warmup：バーンインの数を指定する。他のStanのwarmupと同じっぽい。chains：並列して走るチェーンの数。これは他のStanと同じ。seed：seed値を設定。これは他のStanと同じ。
下記を参考にしました。
またノーUターンサンプラは下記に詳細が載っています。
診断と結果の解釈
ここら辺はもっと便利なものがあるのかどうか分かりませんが、ひとまずPyMC3の収束診断っぽいものを作ってみました。
### 事後分布の診断 la = fit_sm.stan_variables() plt.figure(figsize=(10,3)) bins=15 for i, k in enumerate(la.keys()): boxplot_list =[] for j in range(chains): x = la[k][j*iter_sampling:(j+1)*iter_sampling] hist, bin_edges = np.histogram(x, bins=bins) # 度数分布に変換 #print(len(hist),len(bin_edges)) plt.subplot(len(la.keys()), 3, 3*i+1) plt.plot(bin_edges[:bins], hist, alpha=0.6, lw=0.6) plt.subplot(len(la.keys()), 3, 3*i+2) plt.plot(x, alpha=0.5, lw=0.8) boxplot_list.append(x) plt.subplot(len(la.keys()), 3, 3*i+3) plt.boxplot(boxplot_list) plt.show() ### 結果を解釈する samples = az.from_cmdstanpy(posterior = fit_sm) print(az.summary(samples))
また結果の要約についてもarvizというライブラリを使うと比較的簡単に出来るようです。
Data variables:
theta float64 1.002
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
theta 0.335 0.133 0.095 0.577 0.003 0.002 1462.0 1934.0 1.0
CmdStanPyのインストール
MCMCのライブラリとして昔はPyMC3を使っていましたが、仕様がコロコロ変わったり、TensorFlowベースになる予定だったPyMC4を開発中止したりと、何かと悪い印象がありました。そんなわけで、半年前くらいからStanに切り替えました。RのStanかcolaboratory上からPyStanを使っています。「Stanは使いやすくていいな」なんて思いながらしばらく使っていましたが、最近CmdStanPyというのを見つけて、インストールしようとしました。CmdStanPyというのは下記のものです。
www.wenyanet.com
公式ドキュメントに沿ってインストールしたつもりですが、苦戦しました。
StanをPythonで呼び出そうとすると何かと苦労する印象です。そういえば、何年か前にMCMCの勉強会がありました。PythonユーザはPyStanを使おうと頑張りましたが、結局インストールできた人は居なかった気が・・・。Pythonユーザは諦めてreticulateでRのStanを呼び出していました。（その時は私はPyMC3を使っていました）
例に倣ってCmdStanPyのインストールもすんなり行きませんでした。いろいろ試行錯誤した結果よさそうな方法を模索したので、メモしておきます。
インストール
デスクトップにインストール
私のデスクトップ環境は下記のとおりです
- Windows11
- python 3.7.11
手順としては最初にcmdstanpyをインストールしてから、次にcmdstanをインストールするようです。
まずはターミナルからconda-forgeでCmdStanPyをインストールしました。というかconda-forgeからはインストール出来ますが、pipからではインストールできませんでした。コマンドは下記のとおりです。
conda install -c conda-forge cmdstanpy
その次に、下記のようにCオプションを使ってCmdStanをインストールしました。
python -c 'import cmdstanpy; cmdstanpy.install_cmdstan()'
もしくは、
python -m cmdstanpy.install_cmdstan
これでなんとか行けました。
試したこと
すんなりインストール出来たわけではなく、色々と試行錯誤しました。
pipからcmdstanpyをインストールすると、次段階のcmdstanがインストール出来ません。原因不明です。
なので、公式ドキュメントに書いてある、
pip install --upgrade cmdstanpy
とか
pip install --upgrade cmdstanpy[all]
からインストールすると、次段階のcmdstanをインストールする際に、「ファイルがありません」というようなエラーが出力されます。
下記も同じです。
install_cmdstan
同じくPythonから下記を呼び出してもダメでした。
import cmdstanpy cmdstanpy.install_cmdstan()
なぜ駄目なのかソースコードを見てみましたが、段々面倒になってきてやめました。
Dockerを利用してインストールする
('23/1/5追記)
Cmdstanpyを再インストールしたところ、ローカルの環境がぐちゃぐちゃになってしまいました・・・。こういうのに時間を取られるのが嫌なので、Docker上にインストールします。というかDocker上にインストールするのが一番楽だと思います。
巻末の参考３を参照してDockerfileを下記のように書いてインストールしました。参考元はubuntuのイメージを使っていますが、私はminicondaを使いました。
FROM continuumio/miniconda3 # conda create RUN conda create -n chemodel python==3.9 # install conda package RUN conda install -c conda-forge cmdstanpy=1.0.4 RUN conda install -c conda-forge cmdstan=2.30.0
Pythonでよく使うコード：その２ 描画編
何度も引き直すのでメモ記事を作成しました。下記の続きです。
publicjournal.hatenablog.com
publicjournal.hatenablog.com
描画系
基礎
qiita.com
時間の表記
テクニック
描画の保存
その他の描画ライブラリ
Gitのアクセストークン
結構面倒で，ちょっとハマってしまいました．．．
【git】パスワードを受け付けなくなりアクセストークンを取得後に行った設定
docs.github.com
Google ColaboratoryでNotebookをHTMLに変換する
Google Colaboratoryはどこでも同じ環境で作業ができるので便利ですが，jupyter labにあるようなHTMLダウンロードの機能がないようです．
そこで下記を参考にしました．
qiita.com
ただ，ファイルの場所を聞かれたりとそのまま使うことが出来なかったので，すこし書き換えました．
何度も書き換えるのが面倒なので，このブログに書き下して，コピペして使えるようにします．
from google.colab import files import re import shutil # そのままHTML化すると上書きされてしまうので，一旦コピーする dir = '/content/drive/MyDrive/・・・コピーしたいファイルが格納されているディレクトリ・・・/' file_name = '・・・今使っているファイルの名前・・・.ipynb' fn = dir+'・・・コピー後の名前・・・.ipynb' shutil.copyfile(dir+file_name , fn) # HTML化 fn_s = re.escape(fn) output_fn = fn.split('.', 1)[0]+'.html' output_fn_s = re.escape(output_fn) !jupyter nbconvert --to html $fn_s files.download(output_fn) !rm $fn_s
