rstanで打ち切りデータがあるときのパラメータ推定をする

これに、打ち切りデータがあるときの平均値の推定問題がある。

StanとRでベイズ統計モデリング (Wonderful R)

データが正規分布に従うのだろうが、25を下回るデータは<25となっているので、このデータを無視して平均値を推定すると、本当にデータがあったときより過大評価するだろう、ということ。
打ち切りの場合に<25を無視するのではなく、<25となったときにの区間までの（累積）確率分布を考慮すれば、少しはマシな推定になる、という話。
オリジナルのコードは

data {
  int N_obs;
  int N_cens;
  real Y_obs[N_obs];
  real L;
}
parameters{
  real<lower=0> m;
  real<lower=0> s_Y;
}
model{
  for(n in 1:N_obs){
    Y_obs[n] ~ normal(m, s_Y);
  }
  target += N_cens * normal_lcdf(L | m, s_Y);
}

となっているが、打ち切りのデータにインデックスをつけてforを回すパターンにした。また、サンプリングはY ~ normal() でもよいが、target 記法で合わせてみようと思って
サンプリングの部分はtarget += _lpdf(data | parameters)（そのデータをサンプリングする確率密度lpdf）
打ち切りの部分はtarget += _lcdf(censor | parameters)（ までの累積確率分布lcdf）
となる。

rstan で推定した、正規分布のパラメータ の事後分布である。この事後分布の最頻値は、シミュレーションで設定した真の値 50 より（たまたま）小さい。 は40 に設定したが、 を下回った値は<40 と記録されてしまうとして、このデータを除外して平均を推定した場合はOmit の縦線となり、Censor の縦線より大きくなっている。

さて、rstanで打ち切りデータを考慮した場合と、除外して平均を推定した場合とで、どれくらい推定値が変わるかをシミュレーションした。
除外した場合は、真の値よりほとんどの場合で大きい値を推定してしまう。
打ち切りを考慮した場合は、だいたい真の値を推定しているようだった。

library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())


# 偏差値っぽい感じで
m_true <- 50
s_true <- 10
n <- 30
y <- rnorm(n, m_true, s_true)

L <- 40
idx <- which(y < L)
y0 <- replace(y, idx, L)


code <- "
data {
  int N;
  real Y[N];
  int I[N];
  real L;
}
parameters{
  real<lower=0> m;
  real<lower=0> s;
}
model{
  for(i in 1:N){
    if(I[i] == 0){ // データがある場合
      target += normal_lpdf(Y[i] | m, s);
    } else {       // 打ち切りの場合
      target += normal_lcdf(L | m, s);
    }
  }
}
"

m0 <- stan_model(model_code=code)
standata <- list(N=length(y), Y=y, I=replace(rep(0, length(y)), y < L, 1), L=L)
fit <- sampling(m0, standata, iter=1000, warmup=300, chain=4)
ex <- extract(fit, pars=head(fit@model_pars, -1))


library(vioplot)
par(mar=c(3, 2, 3, 2))
plot(y, rep(0.99, n), pch=16, col=2, xlab="", ylab="", ylim=c(0.95, 1.5), yaxt="n")
vioplot(ex$m, horizontal=TRUE, side="right", ylim=range(y), add=TRUE)
pa <- par()$usr
abline(v=m_true); text(m_true, pa[4], "True", xpd=TRUE, srt=0, pos=3)
abline(v=median(ex$m)); text(median(ex$m), pa[4], "Censor", xpd=TRUE, srt=0, pos=3, offset=1.5)
abline(v=mean(y[y>L])); text(mean(y[y>L]), pa[4], "Omit", xpd=TRUE, srt=0, pos=3)
abline(v=L); text(L, pa[4], "L", xpd=TRUE, srt=0, pos=3)


# 除外した平均と、rstanで推定した平均の違いをシミュレーション
iter <- 1000
M <- matrix(0, iter, 2)
for(i in 1:iter){
  y <- rnorm(n, m_true, s_true)
  standata <- list(N=length(y), Y=y, I=replace(rep(0, length(y)), y < L, 1), L=L)
  fit <- sampling(m0, standata, iter=1000, warmup=300, chain=4)
  ex <- extract(fit, pars=head(fit@model_pars, -1))
  M[i, ] <- c(mean(y[y > L]), median(ex$m))
}


xylim <- range(M)
par(mar=c(5, 5, 3, 3), cex.lab=1.5, cex.axis=1.2)
plot(M, xlim=xylim, ylim=xylim, xlab="L 以下を除外した平均", ylab="累積確率を考慮した平均", pch=16)
abline(v=m_true, h=m_true, lty=3, col=4)
abline(0, 1, lty=3, col=4)
vioplot(M[,1], at=par()$usr[4], add=TRUE, side="right", xpd=TRUE, horizontal=TRUE, wex=3)
vioplot(M[,2], at=par()$usr[2], add=TRUE, side="right", xpd=TRUE, wex=3)