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

Community
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)
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) は架空のデータ分析者であり、日本の若手のデータ分析者集団のペンネームである。当初このデータ分析者集団は秘密結社として活動し、ホクソエムを一個人として活動させ続けた。
Machine Learning and Data Analysis Company for Your Smiles :)