Stan で相関係数を推定する

Stan で相関係数を推定してみます。

library(rstan) 

N <- 1000
x <- rnorm(N, mean = 50, sd = 10)
y <- 10 + 0.8 * x + rnorm(N, mean =0, sd = 7)
data <- data.frame(x, y)

stancode <- '
  data{
    int<lower=0> N;
    vector[2] y[N];
  }
  parameters {
    vector[2] mu;
    real<lower=0> var_x;
    real<lower=0> var_y;
    real cov;
  }
  transformed parameters {
    matrix[2,2] sigma;
    sigma[1,1] <- var_x;
    sigma[2,1] <- cov;
    sigma[1,2] <- cov;
    sigma[2,2] <- var_y;
  }
  model{
    for(i in 1:N)
      y[i] ~ multi_normal(mu, sigma);
  }
  generated quantities {
    real cor;
    cor <- cov / sqrt(var_x * var_y);
  }
'
datastan <- list(N=N, y=data)
fit <- stan(model_code=stancode, data=datastan, iter=1000, chain=4)
traceplot(fit, ask=FALSE)
print(fit, digit=2)

結果

               mean se_mean   sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
mu[1]         50.04    0.01 0.33    49.41    49.82    50.04    50.26    50.69   658 1.01
mu[2]         50.51    0.01 0.34    49.83    50.28    50.50    50.73    51.18   693 1.01
var_x        106.88    0.18 4.69    98.30   103.69   106.76   109.84   116.71   652 1.00
var_y        119.10    0.21 5.53   109.02   115.21   118.86   122.81   130.27   679 1.00
cov           87.69    0.19 4.55    79.13    84.55    87.59    90.67    96.63   587 1.00
sigma[1,1]   106.88    0.18 4.69    98.30   103.69   106.76   109.84   116.71   652 1.00
sigma[1,2]    87.69    0.19 4.55    79.13    84.55    87.59    90.67    96.63   587 1.00
sigma[2,1]    87.69    0.19 4.55    79.13    84.55    87.59    90.67    96.63   587 1.00
sigma[2,2]   119.10    0.21 5.53   109.02   115.21   118.86   122.81   130.27   679 1.00
cor            0.78    0.00 0.01     0.75     0.77     0.78     0.79     0.80   887 1.00
lp__       -5248.75    0.05 1.48 -5252.37 -5249.52 -5248.47 -5247.65 -5246.71   824 1.00

普通に計算した値と一致することを確かめます。

cor(x, y)

結果

[1] 0.7771426