R
RStan

Stan で相関係数を推定する

More than 3 years have passed since last update.

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