Stanで「丸め」の実装 - Ichijoji Sagarimatsu

Stan reference guideにも記載がある「丸め」について検証します。
丸められたデータとは、本当は連続値なのに、四捨五入（あるいは切り捨て）などで離散値のようになっているデータを指すようです。
打ち切りと似ている気がしますが、レファレンスでは別項目として記載されています。

library(rstan)
set.seed(123)

Y <- rnorm(50, 8, 2) # 平均8, 標準偏差2の正規分布
Y <- round(Y, 0) # 四捨五入

f:id:IchijojiSagarimatsu:20201127124602p:plain
こんなデータの平均と標準偏差を求めます。

Stan reference guideの1つ目のやり方

data{
  int N;
  int Y[N];
}
parameters{
  real mu;
  real<lower=0> sigma;
}
model{
  for(n in 1:N)
    target += log(
      Phi((Y[n] + 0.5 - mu) / sigma) - Phi((Y[n] - 0.5 - mu) / sigma)
      );
}

データY[n]の"真の値"は Y[n]-0.5 ~ Y[n]+0.5 の区間にあるので、その確率を求めて、対数尤度を計算しています。
正規分布の標準化（平均0標準偏差1に変換）を行いPhiで標準正規累積分布関数として実装します。

modelブロックは以下のように書くこともできます。

model{
  for(n in 1:N)
    target += log(normal_cdf(Y[n]+0.5, mu, sigma) - normal_cdf(Y[n]-0.5, mu, sigma));
}

行っていることは同じです。
(Y[n]+0.5 | mu, sigma)と書くとエラーが出ます。~cdfはコンマでつなげるんですね。

Stan reference guideの2つ目のやり方

整数に丸め込まれたデータYに誤差項Y_err(-0.5 ~ 0.5) をパラメータとして宣言し、潜在変数 z[n] につなげます。

data{
  int N;
  int Y[N];
}
parameters{
  real mu;
  real<lower=0> sigma;
  vector<lower=-0.5, upper=0.5>[N] Y_err;
}
transformed parameters{
  vector[N] z;
  for(n in 1:N)
    z[n] = Y[n] + Y_err[n];
}
model{
  z ~ normal(mu, sigma);
}

キックするRコードと、事後分布を確認します（ほぼ同じ推定をしてくれるので方法2の結果だけ載せます）。
また方法2では、丸め誤差 Y_err も推定してくれるわけなので、分布の「端」...ここではYの最大となるY[k]と付随する誤差項Y_err[k]の事後分布も確認してみましょう。
Y[k]は真の平均値から離れているので、Y_err[k]の事後確率密度は -0.5周辺 → 0.5周辺に向かうにつれて小さくなることが予想されます。

data <- list(N=length(Y), Y=Y)
fit <- stan("marume1.stan", data=data, seed=123)
fit2 <- stan("marume2.stan", data=data, seed=123)

max <- which(Y==max(Y)) # Yが最大値をとる行を調べる
stan_hist(fit2, pars=c("mu", "sigma", "z[16]", "Y_err[16]")) # Y[16]でYが最大値を取る

f:id:IchijojiSagarimatsu:20201127130938p:plain