library(rstan)
N <- 1001
N_miss <- ceiling(N / 10)
N_obs <- N - N_miss
mu <- 3
sigma <- 2
y_obs <- rnorm(N_obs, mu, sigma)
missing_data_code <-
'
data {
int N_obs;
int N_miss;
real y_obs[N_obs];
}
parameters {
real mu;
real sigma;
real y_miss[N_miss];
}
model {
// add prior on mu and sigma here if you want
y_obs ~ normal(mu,sigma);
y_miss ~ normal(mu,sigma);
}
generated quantities {
real y_diff;
y_diff <- y_miss[101] - y_miss[1];
}
'
results <- stan(model_code = missing_data_code,
data = list(N_obs = N_obs, N_miss = N_miss, y_obs = y_obs))
y_diff <- apply(extract(results, c("y_miss[1]", "y_miss[101]")), 1:2, diff)