【Stan】Stanの構造を勉強中

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形式でデータを格納していることを指定しています。



stan解析用データ - Google スプレッドシート

使ったデータが気になる方はGoogleスプレッドシートをご覧ください。

ちなみにこのデータ以下の式でεがそれぞれ平均0,分散1の正規分布になるように乱数を生成しています。

y=exp(ln(10+\varepsilon_b)\times{x}+\varepsilon_y)

 

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を計算により変換出来ます。

今回はb2=ln(b1)とだけ指定してみます。

 

modelブロック

model {
  log_Y ~ normal( b2 * X, sigma);
  b1 ~normal(0,1000);
  sigma ~ normal(0,1000);
}

 

ここが一番大事で他のブロックと違ってmodelブロックは必須らしい。

dataとparametersの関係式を指定したり、parametersの事前分布を指定することが出来る。

 

今回はln(Y)=b2\times{X}としたが、transformedされていないdataとparametersを使って記載するとY=exp(ln(b1)\times{X})となるのでややシンプルになっていることが分かる。

 

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

 

実行するとこういう結果になりました。