Stanでは、5つのブロックがあるらしい。
data: 観測データ
transformed data: dataに基づいて計算される変数
parameters: モデルのパラメータ
transformed parameters: prametersに基づいて計算される変数
model: モデル
generated quantities: サンプリングされたパラメータを用いて計算される変数
それぞれどう使うのか、今一度まとめておこうかと思います。
Stanコード全体
以下が今回使ったStanコード全体。それぞれのブロックについて説明していきます。
data { int<lower=0> N; vector[N] X; vector[N] Y; } transformed data { vector[N] log_Y = log(Y); } parameters { real b1; real<lower=0> sigma; } transformed parameters { real b2 = log(b1); } model { log_Y ~ normal( b2 * X, sigma); b2 ~normal(0,1000); sigma ~ normal(0,1000); } generated quantities { vector[N] y_sim; y_sim = exp(b2 * X); }
dataブロック
data { int<lower=0> N; vector[N] X; vector[N] Y; }
まずは実際のデータを入力するdataブロック。ここではデータの個数がN個である事。XとYがvector形式でデータを格納していることを指定しています。
使ったデータが気になる方はGoogleスプレッドシートをご覧ください。
ちなみにこのデータ以下の式でεがそれぞれ平均0,分散1の正規分布になるように乱数を生成しています。
transformed dataブロック
続いて観測dataを変換できるtransformed dataブロック。
transformed data {
vector[N] log_Y = log(Y);
}
正直使ったことはなかったのだが、後のmodelやparametersを簡単にして、計算を速くすることも出来る(のだと思う)。
parametersブロック
parameters { real b1; real<lower=0> sigma; }
parametersブロックでは計算したいパラメーターの値を設定します。
今回はXとYの関係を表すb1と誤差のsigmaを指定します。
transformed parametersブロック
transformed parameters { real b2 = log(b1); }
このブロックではparametersを計算により変換出来ます。
今回はとだけ指定してみます。
modelブロック
model {
log_Y ~ normal( b2 * X, sigma);
b1 ~normal(0,1000);
sigma ~ normal(0,1000);
}
ここが一番大事で他のブロックと違ってmodelブロックは必須らしい。
dataとparametersの関係式を指定したり、parametersの事前分布を指定することが出来る。
今回はとしたが、transformedされていないdataとparametersを使って記載するととなるのでややシンプルになっていることが分かる。
generated quantitiesブロック
generated quantities {
vector[N] y_sim;
y_sim = exp(b2 * X);
}
最後にgenerated quantitiesブロック。こちらも使ったことがなかったのですが、サンプリングの後の計算結果を使って、新しくデータを生成することが出来るらしい。
計算後のparametersを使って予測値を出したり、計算後のparametersをさらに特定の計算式で変換したり、といったことが出来るそう。
今回はサンプリング後のb2とXの値を使って、Yの予測値を返してみた。
ちなみに調べているうちに知ったのだが、functionsというブロックで独自の関数を設定することも出来る模様。これはそのうちという事で。
参考:python部分
import os import time import numpy as np import arviz as az from matplotlib import pyplot as plt from cmdstanpy import CmdStanModel import pandas as pd #stanファイルのコンパイル stan_file = "./test.stan" exe_file = "./test" if not os.path.exists(exe_file): model = CmdStanModel(stan_file=stan_file) else: model = CmdStanModel(exe_file=exe_file) #観測データの読み込み csv_file = "./test.csv" df = pd.read_csv(csv_file) #実行 data = { "N": 91, "X": df['X'].values, "Y": df['Y'].values } import multiprocessing num_cpu = multiprocessing.cpu_count() fit = model.sample( data=data, chains=4, # chain数 seed=1, # seed固定 iter_warmup=1000, # warmupの数 iter_sampling=2000, # samplingの数 parallel_chains=num_cpu, # 並列数 save_warmup=True, # warmupもCSVに保存 thin=1, # サンプリング間隔 #output_dir=output_dir, # 出力先 show_console=False, # 標準出力 show_progress=True # progress出力 ) az.style.use("arviz-darkgrid") cmdstanpy_data = az.from_cmdstanpy( posterior=fit, log_likelihood="lp__", ) ll_data = cmdstanpy_data.log_likelihood az.plot_trace(ll_data, compact=False); fit.summary()
実行するとこういう結果になりました。