7.5　回帰分析

データ $x = (1,2,3,4,5,6)$，$y = (1,3,2,4,3,5)$ を $y \sim ax + b$ でフィットする問題です。誤差が正規分布になるモデル $y \sim \mathcal{N}(ax + b, \sigma^2)$ を仮定して，パラメータ $a$，$b$，$\sigma^2$ の事後分布を求めます。

Stanらしいコードは次のようになります（$a$，$b$，$\log \sigma$ について一様な Jeffreys の事前分布を仮定します）：

// ex75.stan

data {
    int n;
    vector[n] x;
    vector[n] y;
}

parameters {
    real a;
    real b;
    real logsigma;
}

model {
    y ~ normal(a * x + b, exp(logsigma));
}

Rには次のように打ち込みます：

xdata = c(1, 2, 3, 4, 5, 6)
ydata = c(1, 3, 2, 4, 3, 5)
data = list(n=length(xdata), x=xdata, y=ydata)
fit = stan("ex75.stan", data=data)

ここで pairs() プロットをしてみましょう：

pairs(fit)

a と b が強く相関していることがわかります。そこで，x と y からそれぞれの平均を引いてみましょう：

data = list(n=length(xdata), x=xdata-mean(xdata), y=ydata-mean(ydata))
fit = stan("ex75.stan", data=data)
pairs(fit)

今度はうまくいきそうです。繰返し数を増やして，結果を見てみます：

fit = stan("ex75.stan", data=data, iter=26000, warmup=1000)

結果は次のようになります：

fit
Inference for Stan model: ex75.
4 chains, each with iter=26000; warmup=1000; thin=1; 
post-warmup draws per chain=25000, total post-warmup draws=1e+05.

          mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
a         0.63    0.00 0.30  0.05  0.47  0.63  0.78  1.21 34931    1
b         0.00    0.00 0.52 -0.99 -0.27  0.00  0.27  1.01 40509    1
logsigma  0.01    0.00 0.40 -0.64 -0.28 -0.04  0.24  0.91 27730    1
lp__     -3.04    0.01 1.67 -7.33 -3.79 -2.61 -1.82 -1.15 19726    1

$a$ の中央値は $0.63$，95%信用区間は $[0.05, 1.21]$ です。

本文の図7.7に相当する詳しい分布も描いてみましょう：

e = extract(fit)
mean(e$a)  # lm(y ~ x) の結果 a = 0.6286 と比較
hist(e$a, breaks=seq(-100,100,0.05), col="gray", xlim=c(-0.5,2), freq=FALSE)
curve(dt((x-0.6286)/0.2100,4) / 0.2100, add=TRUE)  # t分布と一致するか確認

古典的な $t$ 分布とぴったり一致します。

7.5 回帰分析

7.5　回帰分析