LoginSignup
4

More than 5 years have passed since last update.

Stan で相関係数を推定する

Posted at

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

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
4