俺俺Stan manual翻訳:5-1章: 線形回帰モデル

本文章はStan Modeling Language
Stan Development Team. 2015. Stan Modeling Language Users Guide and Reference Manual, Version 2.7.0.の5-1章: Linear regressionの翻訳となります。

訳者は勉強がてら翻訳しようと思っただけで、プログラムの専門家でもデータ解析の専門家でもありません。そのため誤訳・意味の取り違えなどあると思いますが、その際はご指摘いただけると幸いです。5章は長いので何回かに分けて進めたいと思います。また、翻訳ペースも気分しだいです。ご了承ください。

5.回帰モデル

Stanでは単純な線形回帰からマルチレベル一般化線形モデルまで、さまざまな種類の回帰モデルを実行することができます。

5.1.線形回帰モデル

ここでは誤差を正規分布と仮定する、1つの説明変数、傾き、切片を持つ、最も単純な回帰モデルについて解説します。このモデルは基本的な回帰モデルとして、以下のように示すことができます。

\begin{aligned}
\mathbf{y}_n = \alpha + \beta \mathbf{x}_n + \mathbf{\epsilon}_n\ where\ \mathbf{\epsilon}_n \sim Normal(0,\sigma).
\end{aligned}

また、以下のように、残差項を省略して書いたり、

\begin{aligned}
\mathbf{y}_n - (\alpha + \beta \mathbf{x}_n)   \sim Normal(0,\sigma).
\end{aligned}

さらに省略して以下のように書くこともあります。

\begin{aligned}
\mathbf{y}_n \sim Normal(\alpha + \beta \mathbf{x}_n, \sigma).
\end{aligned}

Stanにおいて、このモデルは以下のように記述されます。

  data {
    int<lower=0> N;
    vector[N] x;
    vector[N] y;
  }
  parameters {
    real alpha;
    real beta;
    real<lower=0> sigma;
  }
  model {
    y ~ normal(alpha + beta * x, sigma);
  }

data block内のNは観測値の数を示しており、説明変数, 目的変数がそれぞれ持つデータの数に対応しています。alpha, betaはそれぞれ切片のパラメータを表しており、モデルの誤差は分散sigmaの正規分布に従うことがparameters blockに記述されています。このモデルは2つの回帰係数の十分に広い一様分布を持っています。

マトリクスの宣言とベクトル化

前項で示したモデルのsampling statementは以下のようにベクトル化した形でした。

  y ~ normal(alpha + beta * x, sigma);

同じモデルをベクトル化せずに書くと以下のようになります。

  for (n in 1:N)
    y[n] ~ normal(alpha + beta * x[n], sigma);

ベクトル化のメリットとして、簡潔にモデルを記述できること、処理の高速化が挙げられます。

一般的に、Stanでは"normal" のような分布関数の引数にベクトルを入力することができます。他のいずれかの引数がvectors, arraysの場合には、それらのサイズを合わせる必要があります。
他の引数のいずれかがscalarである場合は、各ベクトルの入力が再利用されます。ベクトル化の詳細については、第29章を参照してください。
Stanがこのように振る舞う理由として、Stanの算術演算子は行列の演算演算を実行するには荷が重い事が挙げられます。今回の場合、xがベクトル型, betaは実数型であり、beta*xはベクトル型で表現されます。 Stanではベクトル化データを扱うことができるため、1つ以上の説明変数を持つ多変量モデルであっても、下記のようにmatrix表記を使うことで直接モデルを表記できます。

  data {
    int<lower=0> N; // データの個数
    int<lower=0> K; // 説明変数の数
    matrix[N,K] x; // 説明変数データのmatrix
    vector[N] y; // 目的変数ベクトル
  }
  parameters {
    real alpha; // 切片
    vector[K] beta; // 説明変数の傾き
    real<lower=0> sigma; // 分散の振る舞い
  }
  model {
    y ~ normal(x * beta + alpha, sigma); // 尤度
  }

sigmaの宣言に際し、lower=0が制約として記入されています。これは、モデルの中でsigmaが０位上の値しか取らないことを意味します。モデルブロックに事前分布はなく、効果は非負の実数において、十分に広い一様分布に従います（この文章自信なし）。もちろんより良い事前分布を使っても構いませんが、十分に広い一様分布を用いた場合でも、事後分布が適切であれば許容されます。
上記のモデルでは、xはN × K個の説明変数データであること、K個の説明変数があること、N個のx × betaがN個のデータそれぞれに対応していることが記述されています。これらの目的変数がN個の目的変数ベクトルと対応しているため、行列演算を使ってモデルを記述することができるのです。また、x個の列の中から1つの値を選んだり、alphaのパラメータ（切片のデータ）を取り除いたりすることもできます。
出力される結果は同じですが、上記モデルのsampling statementは、下記のループでモデルをコーディングするアプローチに比べ、より効率的です。

  model {
  for (n in 1:N)
  y[n] ~ normal(x[n] * beta, sigma);
  }

Stanのインデックススキームでは、"x[n]"は行nのmatrixを意味します。これは、betaが列ベクトル、出力される x[n] * betaがreal型のスカラーであるためです。

切片について

以下のモデル表記のように、切片の傾きalphaはどこにも表記されていません。

  y[n] ~ normal(x[n] * beta, sigma);

その代わりに、Stanでは入力した"matrix x"の最初の列は、1の値の列であると仮定しています。このように、出力される"beta[1]"は切片の役割を果たします。これは切片の情報が傾きと異なる事前分布から算出されることを避けるためです。もし切片が傾きと異なる分布から得られるのであれば、その値を使わない方が良いでしょう(自信なし)。しかし、それほど計算時間が大きく変わってくるわけでは無いので、切片の算出は明確な基準に基づいて行われるべきでしょう。