「Stan で欠測データの相関係数を推定してみた」補足資料 #R

実行したプログラムの全体です。

model.stan

data{
  int<lower=0> N;
  vector[2] y[N];
  int<lower=0> Nmiss;
  real x[Nmiss];
}
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);
  for(i in 1:Nmiss)
    x[i] ~ normal(mu[1], sqrt(sigma[1,1]));
#    increment_log_prob(-0.5 * log(2 * pi() * sigma[1,1]) - 0.5 * ((x[i] - mu[1])^2)/sigma[1,1]);
}
generated quantities {
  real cor;
  cor <- cov / sqrt(var_x * var_y);
}

main.R

library(rstan) 

set.seed(123)
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)
data1 <- subset(data, x >= 60)
data2 <- subset(data, x < 60)

set.seed(314)
datastan <- list(N=nrow(data1), y=data1, x=data2$x, Nmiss=nrow(data2))
fit <- stan(file="model.stan", 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.17    0.01  0.31    49.57    49.96    50.16    50.37    50.78   709 1.01
mu[2]         47.91    0.11  2.02    43.73    46.63    47.91    49.19    51.87   320 1.01
var_x         99.02    0.16  4.66    90.17    95.93    98.97   102.08   108.24   872 1.00
var_y        153.43    1.57 26.75   109.73   134.96   151.41   168.48   211.08   292 1.02
cov           99.73    0.73 13.26    74.21    91.13    99.59   108.24   126.37   331 1.02
sigma[1,1]    99.02    0.16  4.66    90.17    95.93    98.97   102.08   108.24   872 1.00
sigma[1,2]    99.73    0.73 13.26    74.21    91.13    99.59   108.24   126.37   331 1.02
sigma[2,1]    99.73    0.73 13.26    74.21    91.13    99.59   108.24   126.37   331 1.02
sigma[2,2]   153.43    1.57 26.75   109.73   134.96   151.41   168.48   211.08   292 1.02
cor            0.81    0.00  0.04     0.71     0.79     0.81     0.84     0.87   785 1.01
lp__       -3171.17    0.07  1.73 -3175.50 -3172.02 -3170.86 -3169.92 -3168.93   577 1.00