はじめに
前回扱った有限混合モデルでは、クラスター数が恣意的に決定されたり、パラメータの識別性が弱かったりといったデメリットがあった。そこで、今回は特にクラスター数の恣意性を克服すべく、ディリクレ過程モデルを扱う。
一般にディリクレ過程によるモデルはノンパラメトリックな手法として知られるが、まずノンパラメトリックな手法を概観する。
以下のパラメトリックな設定を考える。
パラメトリックな設定では、$\mathcal F^*$は$\mathcal F=\lbrace \mathcal Y上のすべての分布\rbrace$に比べて小さいことがわかる。
したがって、ノンパラメトリックな設定では$\mathcal F$の部分集合のうち、より大きな（柔軟な）ものを考える。そのアプローチの一つが、確率過程によるものである。
出力のベクトルを$y=(y_1, ,\cdots, y_n)$とし、以下の階層的なディリクレ過程混合モデルを考える。
ただし、$k(\theta_i)$はパラメトリックな密度関数。$\alpha$は集中度パラメータで、$P_0$が$P$に「収縮」する程度を示す。
①,②,③の表記は$y_i|P$を$f(y|P)$からサンプリングするのと同義である。
$P(\cdot)$は$P_0(\theta)$から生成されたアトム$\theta_h$の重み付き和として表現でき、これを棒折り過程という。
ただし、$\pi_h = v_h\Pi_{l<h}(1-v_h)$, $v_h = Beta(1,\alpha)$, $\theta_h \sim P_0.$
一般に、$\pi_1 = v_1$, $\pi_i = v_i\Pi_{j=1}^{l-1}(1-v_j)$である。
データの生成
平均がそれぞれ2,5,10のポアソン分布に従うサンプルを200個ずつ発生させ、分布を確認する。
set.seed(12345)
lambdas <- c(2, 5, 10)
n_per_group <- 200
df <- tibble::tibble(
group = rep(paste0("lambda=", lambdas), each = n_per_group),
y = c(
rpois(n_per_group, lambdas[1]),
rpois(n_per_group, lambdas[2]),
rpois(n_per_group, lambdas[3])
)
)
df |>
ggplot2::ggplot(ggplot2::aes(x = y, fill = group)) +
ggplot2::geom_histogram(binwidth = 1, alpha = 0.6, position = "identity") +
ggplot2::facet_wrap(~ group, scales = "free_y") +
ggplot2::theme_minimal() +
ggplot2::labs(x = "y", y = "Count")
データの生成過程
$y_i \in \mathbb{N}, i=1,\cdots, n$でサンプルを表す。
棒折り過程によって、$y_i$が従う分布を表現する。
また、以下で棒折り過程を表現する。
Stanファイルの実装
data {
int<lower=0> sample_type;
int<lower=1> cluster_type;
array[sample_type] int<lower=0> outcome;
int<lower=0> outcome_max;
}
parameters {
real<lower=0> intensity;
vector[cluster_type] phi;
vector<lower=0, upper=1>[cluster_type - 1] breaks;
real<lower=0> sigma_phi;
}
transformed parameters {
simplex[cluster_type] pi;
{
real sum = 0;
pi[1] = breaks[1];
sum = pi[1];
for (c in 2:(cluster_type - 1)) {
pi[c] = (1 - sum) * breaks[c];
sum += pi[c];
}
pi[cluster_type] = 1 - sum;
}
}
model {
for (n in 1:sample_type) {
vector[cluster_type] case_vector;
for (c in 1:cluster_type) {
case_vector[c] = bernoulli_lpmf(1 | pi[c])
+ poisson_log_lpmf(outcome[n] | phi[c]);
}
target += log_sum_exp(case_vector);
}
phi ~ normal(0, sigma_phi);
breaks ~ beta(1, intensity);
intensity ~ gamma(1, 1);
sigma_phi ~ gamma(0.1, 0.1);
}
generated quantities {
vector[outcome_max + 1] density;
matrix[outcome_max + 1, cluster_type] cluster_level_density;
for (y in 0:outcome_max) {
for (c in 1:cluster_type) {
real lambda = exp(phi[c]);
cluster_level_density[y + 1, c] = pi[c] * lambda^y * exp(-lambda) / tgamma(y + 1);
}
density[y + 1] = sum(cluster_level_density[y + 1, ]);
}
}
パラメータの推定
クラスターの最大数$C$は10とした。
stan_data <- list(
sample_type = nrow(df),
cluster_type = 10,
outcome = df$y,
outcome_max = max(df$y)
)
stan_model <- rstan::stan_model("dirichlet_process_mixture.stan")
vb_fit <- rstan::vb(stan_model,
data = stan_data,
output_samples = 1000,
iter = 10000,
tol_rel_obj = 0.01)
各クラスターが選ばれる確率（つまり、$\pi_c$）を可視化する。
pi_summary <- rstan::summary(vb_fit, pars = "pi")$summary
pi_est <- pi_summary |>
tibble::as_tibble() |>
tibble::rownames_to_column("variable") |>
dplyr::mutate(index = stringr::str_extract(variable, "[0-9]+") |>
as.numeric() |>
as.factor())
ggplot2::ggplot(pi_est) +
ggplot2::geom_point(ggplot2::aes(x = index, y = mean)) +
ggplot2::geom_segment(ggplot2::aes(x = index, xend = index, y = `2.5%`, yend = `97.5%`)) +
ggplot2::labs(x = "Cluster", y = "Probability")
最初の3つのクラスターでデータのほとんどを説明していることがわかる。
また、生成したサンプルの密度と事後分布による密度を比較する。
density_summary <- rstan::summary(vb_fit, pars = "density", probs = c(0.05, 0.95))$summary
density <- density_summary |>
tibble::as.tibble() |>
tibble::rownames_to_column("variable") |>
dplyr::mutate(
Y = stringr::str_extract(variable, "[0-9]+") |>
as.numeric() - 1,
mean = `mean`,
q5 = `5%`,
q95 = `95%`
)
df_data <- tibble::tibble(Y = df$y)
ggplot2::ggplot() +
ggplot2::geom_bar(data = df_data, ggplot2::aes(x = Y, y = ggplot2::after_stat(prop)),
fill = "grey70", colour = "black") +
ggplot2::geom_point(data = density, ggplot2::aes(x = Y, y = mean),
colour = "red", size = 1.5) +
ggplot2::geom_segment(data = density,
ggplot2::aes(x = Y, xend = Y, y = q5, yend = q95),
colour = "red") +
ggplot2::labs(x = "Y", y = "Density / Proportion") +
ggplot2::theme_minimal()
最後に、サンプルの分布と実際にクラスタリングに使われた分布を可視化する。
density_draws <- as.data.frame(vb_fit, pars = "density")
cluster_level_density_draws <- as.data.frame(vb_fit, pars = "cluster_level_density")
density_mean <- colMeans(density_draws)
cluster_level_density_mean <- colMeans(cluster_level_density_draws)
outcome_max <- length(dens_mean) - 1
C <- ncol(dens_c_draws) / (outcome_max + 1)
density_long <- tibble::tibble(
y = rep(0:outcome_max, C),
component = rep(1:C, each = outcome_max + 1),
prob = cluster_level_density_mean
) |>
dplyr::mutate(component = factor(component, labels = paste0("group", 1:C)))
density_df <- tibble::tibble(
y = 0:outcome_max,
prob = density_mean
)
ggplot2::ggplot() +
ggplot2::geom_col(data = density_df, ggplot2::aes(x = y, y = prob),
fill = "grey70", alpha = 0.6, width = 0.9) +
ggplot2::geom_line(data = density_long, ggplot2::aes(x = y, y = prob, color = component), size = 1) +
ggplot2::theme_minimal() +
ggplot2::labs(x = "Count",
y = "Probability"
)
平均がそれぞれ2,5,10のポアソン分布に従うサンプルを発生させていたが、事後密度の推定においては平均が1,3,9あたりのポアソン分布が使われているように見える。
注意：ディリクレ過程では、データ数が大きくなるにつれて推定されるクラスター数が増え、ひとつのデータポイントからなるクラスターがたくさん生成されるとともに、各クラスターのサイズがどんどん小さくなるとされる。
真のデータ生成過程がクラスター数の有限混合モデルであるときには、ディリクレ過程混合モデルで得られた推定値は一致性を持たない。（こちらのqiitaより引用。）
参考文献
https://ito4303.sakura.ne.jp/posts/2024-02-15-Poisson-Dirichlet-Stan/
https://qiita.com/gen_nospare/items/15c6a94d8e43121a686f