Stanによる外れ値のモデリング – 伊東生態統計研究室

『ベイズモデリングの世界』（伊庭幸人編, 岩波書店, 2018）の講義3にある外れ値のモデルをStanで実装してみました。

準備

グラフ作成用にggplot2と、Stanのモデリングとあてはめのためにrstanパッケージを読み込みます。

library(ggplot2)
library(rstan)
## 要求されたパッケージ StanHeaders をロード中です
## 
## rstan version 2.32.7 (Stan version 2.32.2)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)

データ

模擬データを生成します。100個の観測値のうち、最後の4個が外れ値になっています。

set.seed(1234)

alpha <- 1    # intercept
beta <- 0.8   # slope
sigma <- 1.2  # standard deviation

n1 <- 96
n2 <- 4
x1 <- runif(n1, 0, 20)
y1 <- rnorm(n1, alpha + beta * x1, sigma)
x2 <- runif(n2, 0, 10)
y2 <- runif(n2, 10, 20)
df <- data.frame(x = c(x1, x2), y = c(y1, y2),
                 outlier = c(rep("FALSE", n1), rep("TRUE", n2)))

図示します。

ggplot(df, aes(x, y, colour = outlier)) +
  geom_point(size = 3)

モデル

Stanによるモデルです。コンパイルしてstan_modelというオブジェクトに格納しておきます。

外れ値である確率の事前確率をデータとしてQに与えます。観測値Yは、外れ値でないときはNormal(alpha + beta * X, sigma)から、外れ値の場合はUniform(-100, 100)からそれぞれ発生するとします。

generated quantitiesブロックで、各観測値が外れ値である事後確率pを計算します。

Stan User’s GuideのFinite Mixturesを参考にしました。

data {
  int<lower=0> N;
  vector[N] X;
  vector[N] Y;
  real<lower=0, upper=1> Q; // prior outlier prob. 
}

transformed data {
  real B = 100; // range for uniform dist.
}

parameters {
  real alpha;                 // intercept
  real beta;                  // slope
  real<lower=0> sigma;        // std. dev.
}

transformed parameters {
  matrix[N, 2] lp;

  // bernoulli_lpmf(1 | Q) + uniform_lpdf(Y[n] | -B, B)
  lp[, 1] = rep_vector(log(Q) + log(1 / (2 * B)), N);
  for (n in 1:N) {
    lp[n, 2] = log1m(Q)
               + normal_lpdf(Y[n] | alpha + beta * X[n], sigma);
  }
}

model {
  for (n in 1:N) {
    target += log_sum_exp(lp[n]);
  }
  
  //priors
  alpha ~ normal(0, 10);
  beta ~ normal(0, 10);
  sigma ~ normal(0, 5);
}

generated quantities {
  vector[N] p;   // prob. outlier for each observation

  for (n in 1:N) {
    p[n] = exp(lp[n, 1] - log_sum_exp(lp[n]));
  }
}

あてはめ

MCMCを実行して、モデルをあてはめます。

Qの事前確率は0.05に設定しました。

stan_data <- list(N = nrow(df), X = df$x, Y = df$y, Q = 0.05)
fit <- sampling(stan_model, data = stan_data,
                refresh = 0)

結果

alpha, beta, sigmaの事後分布の要約を表示します。それぞれデータの生成に用いたものと近い値が得られました。

print(fit, pars = c("alpha", "beta", "sigma"))
## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##       mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
## alpha 1.18    0.01 0.23 0.73 1.03 1.18 1.33  1.63  1569    1
## beta  0.79    0.00 0.02 0.74 0.77 0.79 0.80  0.83  1547    1
## sigma 1.20    0.00 0.10 1.01 1.12 1.19 1.26  1.42  2019    1
## 
## Samples were drawn using NUTS(diag_e) at Mon Mar 17 19:20:20 2025.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

pの最後の6個の事後分布の要約を表示します。97〜100番目が外れ値です。

print(fit, pars = c(paste0("p[", 95:100, "]")))
## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##        mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
## p[95]  0.00    0.00 0.00 0.00 0.00 0.00 0.00  0.00  1886    1
## p[96]  0.00    0.00 0.00 0.00 0.00 0.00 0.00  0.00  1959    1
## p[97]  1.00    0.00 0.00 1.00 1.00 1.00 1.00  1.00  2117    1
## p[98]  0.68    0.01 0.23 0.19 0.51 0.73 0.88  0.98  2021    1
## p[99]  1.00    0.00 0.00 1.00 1.00 1.00 1.00  1.00  1987    1
## p[100] 1.00    0.00 0.00 1.00 1.00 1.00 1.00  1.00  2136    1
## 
## Samples were drawn using NUTS(diag_e) at Mon Mar 17 19:20:20 2025.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

98番目は外れ値でも、もっとも外れ値でないところ近いものなので、外れ値である事後確率もやや小さくなっています。