Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
4
Help us understand the problem. What is going on with this article?
@hoxo_m

Stan で相関係数を推定する

More than 5 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
4
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
hoxo_m
ホクソエム (hoxo_m) は架空のデータ分析者であり、日本の若手のデータ分析者集団のペンネームである。当初このデータ分析者集団は秘密結社として活動し、ホクソエムを一個人として活動させ続けた。
hoxom
Machine Learning and Data Analysis Company for Your Smiles :)

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
4
Help us understand the problem. What is going on with this article?