Stan Advent Boot Camp 第2日目 Stanコードの構造と基本的な書き方

もう12月ですね。12月といえばアドカレです。

この記事はStan advent calendar2020の2日目の記事です。

1日目の記事で、Stanの導入はうまくできましたでしょうか。Stanはいろんなプラットホームで動きますが、この記事ではR上で動かすことを前提に書きます。よって、rstanパッケージを使います。

rstanが動く環境を用意できたら、次はStanコードを書いて、RからデータをStanに渡して推定します。この記事では、Stanコードの基本的な構造と書き方、データの渡し方について解説します。

Stanコードの構造

Stanコードは、全部で7種類のブロックがあります。

functionsブロック
dataブロック
transformed dataブロック
parametersブロック
transformed dataブロック
modelブロック
generated quantitiesブロック

の7つです。

しかし、実際は、

dataブロック
parametersブロック
modelブロック

だけを使っても、基本的な推定ができます。

この記事では初心者向けということで、この3つについて解説します。

まずは、Stanコードってどんなのか見ておきましょう。下の例は、二項分布に従う確率変数をデータで得たとき、そのパラメータである成功率pを推定するためのStanコードです。

data{
  int N;
  int Y[N];
  int Trial;
}

parameters{
  real&lt;lower=0,upper=1&gt; p;
}

model{
  for(n in 1:N){
    Y[n] ~ binomial(Trial,p);
  }
}

data{

int N;

int Y[N];

int Trial;

}

parameters{

real<lower=0,upper=1> p;

}

model{

for(n in 1:N){

Y[n] ~ binomial(Trial,p);

}

上から順番に3つのブロックがありますね。それがそれぞれ、dataブロック、parametersブロック、modelブロックです。

では、それぞれについて解説します。

dataブロック

dataブロックは、RからStanにデータを渡すためのブロックです。このブロックでR上から読み込む変数を宣言します。宣言するときの型については後述します。ここは宣言だけするところなので、文法などはあまり関係ないです。

parametersブロック

parametersブロックは、Stanで推定するパラメータを宣言します。ここも宣言するだけのところです。上の例にあるように、パラメータの取りうる値の範囲もここで明示します。

modelブロック

modelブロックは、確率モデルを記述するところで、ここでプログラミングっぽい文章を書きます。データについての尤度関数や、パラメータの事前分布をここで書きます。

コードの構造はこんな感じです。

次からは具体的な書き方について書きます。

変数の型の宣言

Stanでは変数の型の扱いが厳格です。最初は、ここに躓く人がいるかもしれません。ただ、パターンを覚えればそんなに難しくはありません。

型の種類

まず、変数が連続か離散か、の区別があります。

int型
real型

この2つの区別はとても重要です。int型は整数型ともいわれ、離散的な変数につかう型です。二項分布とかポアソン分布に従うデータなどはint型を使います。0以上の整数しか扱えません。

real型は、連続型の変数に使われ、少数とか負の値とかをとりえます。正規分布とかに従うデータはこちらを使います。

宣言の例として、連続型の変数Yを宣言するときは、

real Y;

と書きます。

さて、上の解説を見ると全部real型にしておけば大丈夫では？と思うかもしれませんが、そういうわけではないのです。二項分布やポアソン分布のような離散確率分布は、データとしてint型しか受け付けてくれません。このように、Stanは型の扱いが非常に厳格なのです。

次に、ベクトルやマトリックスに関する型の違いについてです。int型とreal型で宣言する場合、それはスカラー（ベクトルではなく、1つだけの数字）の型になります。そして、スカラー以外に、

vector型
matrix型

の2つがあります。これらの型は、線形代数の演算が定義されているという点がポイントです。あとで解説する配列とvector, matrixの違いは、線形代数の演算が可能か否か、という点で異なります。

最後に配列について。Stanでは配列というのはそういう型があるのではなくて、int, real ,vector, matrixそれぞれの型について、配列が定義できます。

たとえば、

vector[N] Y;

と宣言したら、N次元のベクトルを宣言したことになります。そして、

real Y[N];

と宣言したとき、real型の変数YがN個の要素を持つ配列になっている、ということを意味します。どちらもN個の数がはいった変数を宣言している点は同じですが、前者は線形代数の計算ができますが、後者はできません。

また、ベクトルやマトリックス自体を配列にできます。たとえば、

vector[P] Y[N];
matrix[P,P] Y[N];

という感じです。上の宣言はP次元のベクトルYがN個まとめた配列、ということです。下の宣言はP×Pの正方行列YがN個ある配列ということです。これは少し高度な分析のときに使う宣言方法なので覚えておきましょう。

宣言するときのパターン

個人的にオススメな宣言方法を書いておきます。

まず、最初にサンプルサイズ、変数の数といった基本的な定数を宣言します。配列の要素になるこれらの値はすべてint型で宣言しなければならないので、

int N;
int P;

といった感じで宣言して、そのあと変数を宣言していきます。サンプルサイズNの目的変数Yと説明変数Xがある場合、個人的には目的変数は配列で、説明変数はベクトルで宣言することをオススメします。

real Y[N];
vector[N] X;

こんな感じで。理由は、目的変数は離散型の場合はvector型が使えないので、どのみち配列で宣言しないといけないこと、そして目的変数に対して線形代数的な処理を行うことがほぼないこと、の2つが理由です。数少ない例として、多変量正規分布を使う場合は、vector型で宣言しないといけませんが、初心者の人がそれを使うことはないと思うので、大丈夫です。とりあえず目的変数は配列で宣言しておくことにしておけば、離散型で宣言するときも

int Y[N];

というようにrealをintに変えるだけでいいので、楽です。

一方、説明変数は線形モデルの場合は回帰係数と説明変数を行列の掛け算をつかって予測値y~を計算することがあります。また、説明変数に離散型を使うことはないので、vector型かmatrix型で宣言するといいでしょう。

また、Stanでは離散型のパラメータは推定できないので、parametersブロックでint型を使うことはありません。なので、重回帰分析などのように複数の回帰係数のパラメータを宣言するときは、

vector[P] beta;

というように、vector型で宣言するほうが、使い勝手はよいです。

パラメータの下限と上限の宣言

続いて、parametersブロックで重要なのは、推定するパラメータの取りうる値の範囲を宣言に含めておく必要があるということです。たとえば、上の例にあるように、二項分布の成功率パラメータpは、確率ですから、0~1の範囲しかとりえません。そのような場合は

real<lower = 0, upper = 1> p;

というように、lowerとupperを<>で囲って書きます。下限しか制約がない場合は、

real<lower=0> q;

というようにlowerだけ書きます。

この制約は、型の直後に書きます。たとえばvector型なら、

vector<lower=0>[P] alpha;

というように、次元数を表す大カッコの前に書きます。

また、特殊な制約を与えたい場合は、Stanにはそれぞれに便利な型が用意されています。たとえばベクトルだけど大きさの順序が決まっているようなパラメータを推定したい場合、

ordered[P] beta;

と宣言します。配列が1番めが一番小さく、P番目が一番大きくなります。順序性があり、かつすべて非負の値を取らせたい場合は

positive_ordered[P] beta;

とします。

合計して1になる、正の値をとるベクトルは

simplex[P] theta;

というsimplex型を使います。

これら以外にもいろんな型があるので、Stanのユーザーズガイドを見てもらえるといいかと思います。

modelブロックの書き方

さて、ようやくmodelブロックまできました。ここでは確率モデルを記述します。

今回の例は、二項分布に従う確率変数Yをデータで得て、そこから成功率パラメータpを推定するというストーリーで解説します。

いま、コインを投げるという試行を30回独立に行い（Trial=30）、表が出た回数を記録します。それを1セットとして、10セット（N=10）行ったとします。このデータを次のように乱数で発生させました。なお、成功率pの真値(p_true)は0.5としました。

N &lt;- 10
Trial &lt;- 30
p_true &lt;- 0.5

Y &lt;- rbinom(N,Trial,p_true)

N <- 10

Trial <- 30

p_true <- 0.5

Y <- rbinom(N,Trial,p_true)

このデータを使って、成功率pを推定してみます。

もう一度Stanコードを見てみましょう。

data{
  int N;
  int Y[N];
  int Trial;
}

parameters{
  real&lt;lower=0,upper=1&gt; p;
}

model{
  for(n in 1:N){
    Y[n] ~ binomial(Trial,p);
  }
}

data{

int N;

int Y[N];

int Trial;

}

parameters{

real<lower=0,upper=1> p;

}

model{

for(n in 1:N){

Y[n] ~ binomial(Trial,p);

}

まずdataブロックを見て分かるように、変数Yは10個の要素を持つ配列で宣言されています。これについて二項分布を仮定する場合、for文を用いて、上のように書きます。

~（チルダ）を使って、データYのn番目の要素は、binomial分布に従うよ、ということを書いています。Trialは二項分布の試行数のパラメータですが、これは殆どの場合外在的に与えます。よって、推定するのは成功率pだけです。10セットそれぞれ同じコインを使っているので、成功率は同じであると仮定しています。

じつは上のコードは、もう少し簡単にかけて、

data{
  int N;
  int Y[N];
  int Trial;
}

parameters{
  real&lt;lower=0,upper=1&gt; p;
}

model{
  Y ~ binomial(Trial,p);
}

data{

int N;

int Y[N];

int Trial;

}

parameters{

real<lower=0,upper=1> p;

}

model{

Y ~ binomial(Trial,p);

}

とすることもできます。for文を省略して、いっきにpを推定しています。これは、binomialという関数が、配列やベクトルに対しても適用できるようになっているからです。ただ、最初のうちは、上の例のように、for文を書くようにするほうがいいと思います。配列などの扱い方をきちんと身につけてから、省略形を覚えましょう。

Rからの実行

Stanコードが書けたら、さっそくRからrstanを実行して推定してみましょう。

次のようなコードをRで書きます。

library(rstan)

stan_data &lt;- list(N=N,Trial=Trial,Y=Y)

model &lt;- stan_model("prob.stan")

fit &lt;- sampling(model,stan_data)

fit

library(rstan)

stan_data <- list(N=N,Trial=Trial,Y=Y)

model <- stan_model("prob.stan")

fit <- sampling(model,stan_data)

fit

まず、rstanパッケージを読み込みます。

続いて、Stanに渡すためのデータをリスト化します。リスト内では、左がStanで定義した変数名、右がRのオブジェクト名です。

続いて、Stanコードをコンパイルします。コンパイルとは、StanコードをC++言語に変換し、機械語に翻訳することをいいます。このコンパイルが曲者で、Stan初心者が最初に躓くポイントです。

というのは、Stanコードに文法的な間違いがあると、ここでエラーメッセージがでてくるからです。そこで、初心者の人がエラーメッセージで心が折れないように、Stanコードを書くときの格言をお教えしましょう。

文末にはセミコロン
半角文字しか使っちゃダメよ
Stanは間違えない、常に君が間違えている
人間は、lを1 と打ち間違えるときがある
エラーメッセージは友達
文末にはセミコロン
In readLines(file, warn = TRUE) :はエラーじゃないよ、最後の1行をブランクにするんだよ

これらを胸に、コードを修正してみてください。上にも書いてますが、ほとんどのエラーは、文末にセミコロンがないか、スペルミスかのどちらかです。人間はびっくりするような打ち間違いをするので、よーく確認しましょう。探しものと同じです。こんなところにあるはずがないと思ってるところにあります。

コンパイルがうまくできたら、ようやくサンプリングです。MCMCの設定については3日目の記事に委ねます（疲れた）。上のコードを走らせると、一瞬で分析が終わります。

<span class="GD15MCFCKUB ace_keyword">&gt; </span><span class="GD15MCFCKTB ace_keyword">fit
</span><span class="GD15MCFCEUB">Inference for Stan model: prob.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
p       0.50    0.00 0.03    0.44    0.48    0.50    0.52    0.56  1489    1
lp__ -209.87    0.02 0.78 -212.13 -210.04 -209.56 -209.38 -209.33  1442    1

Samples were drawn using NUTS(diag_e) at Tue Dec  1 21:29:41 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
</span>

> fit

Inference for Stan model: prob.

4 chains, each with iter=2000; warmup=1000; thin=1;

post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

p 0.50 0.00 0.03 0.44 0.48 0.50 0.52 0.56 1489 1

lp__ -209.87 0.02 0.78 -212.13 -210.04 -209.56 -209.38 -209.33 1442 1

Samples were drawn using NUTS(diag_e) at Tue Dec 1 21:29:41 2020.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at

convergence, Rhat=1).

みたところ、うまく真値を復元できているようです。

いかがでしたか？（一度使ってみたかった）　Stanコードの書き方について簡単ですがまとめてみました。

Stanを習熟するには、まず写経です。コピペ、ダメ、絶対。自分で書き写してこそ、身につくものがあります。がんばって、エラーメッセージと友だちになってください。