DockerにNumPyroをインストールする

先日までCmdStanPyを使っていました。が、Anaconda をアンインストールしてPipやconda forgeから再インストールしたら、何故かうまく動作しませんでした。。。やっぱりStan系はインストールが大変です、、。
「諦めてPyStanを使うか、PyCM3に戻ろうかなぁ」なんて思っていた矢先、こういう記事を発見しました。

NumPyroというパッケージがあって、これがなかなか興味深いです。ノーUターンサンプラでMCMCの計算が出来るので、CmdStanPyの代わりになりそうです。NumPyroの公式ドキュメントを以下に貼り付けておきます。

さて、CmdStanPyのインストールを試しているうちに、ローカル環境がぐちゃぐちゃになってしまいました、、。これを反省して、今回からはDocker上に環境構築をします。
(「Numpyro インストール」でググってみるとGoogle Colaboratoryを利用した方法が多く紹介されていましたが、練習もかねてDockerを使いました)

インストール

Pipを利用するといろいろと上手くいかなかったので、miniconda上に構築しました。Docker fileを次のように書きました。割とシンプルな例ですが、結構試行錯誤しました。

FROM continuumio/miniconda3  

RUN conda create -n chemodel python==3.10 
  
RUN conda install numpy 
RUN conda install scipy 
RUN conda install -c conda-forge jax 
RUN conda install -c conda-forge numpyro

CmdStanPyの試運転

先日CmdStanPyをインストールしてみました。(詳細は下記参照)
publicjournal.hatenablog.com
RのStanやPyStanと結構違うので、どういうものかと一通り使ってみました。ひとまず“Hello, World”を参照にしながら動作確認をしました。

同じような例題が書籍「Pythonによるベイズ統計モデリング」のP48~57にあるので、同書と同じことをCmdStanPyを使ってやってみたいと思います。

例題

こんな例題を考えます。

コイン投げ問題を考える。コインを10回投げて裏と表を記録する。表を1、裏を0として記録したところ、[1,0,1,0,0,0,0,0,1,0]となった。このデータをもとに、コインが偏っているかを判断する。コインの偏りを調べるために、次のモデルを用いる。


\begin{aligned}
\theta  \sim Beta \big(1,1\big) \\
y  \sim Bin  \big(n=10, p=\theta\big) 
\end{aligned}

プログラムの全貌

プログラムは次の通り。

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

f:id:shu10038:20220227182108j:plain

また結果の要約についても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

arviz-devs.github.io

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のインストールもすんなり行きませんでした。いろいろ試行錯誤した結果よさそうな方法を模索したので、メモしておきます。


インストール

デスクトップにインストール

私のデスクトップ環境は下記のとおりです

手順としては最初に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上にインストールするのが一番楽だと思います。
巻末の参考3を参照して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でよく使うコード:その2 描画編

何度も引き直すのでメモ記事を作成しました。下記の続きです。
publicjournal.hatenablog.com
publicjournal.hatenablog.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