RStan 在精确和变分贝叶斯模式下给出不同的结果

Posted 2023-04-19

【中文标题】RStan 在精确和变分贝叶斯模式下给出不同的结果

【英文标题】：RStan gives different results in exact and variational bayes modes

【发布时间】：2019-07-24 15:54:24

【问题描述】：

我正在开发一个依赖于 RStan 的 R package，我似乎在后者中遇到了故障模式。

我使用精确推理 (rstan::stan()) 运行贝叶斯逻辑回归，并使用变分推理 (rstan::vb()) 得到非常不同的结果。以下代码下载 German Statlog Credit 数据并在该数据上运行两个推理：

library("rstan")

seed <- 123
prior_sd <- 10
n_bootstrap <- 1000

# Index of coefficients in the plot and summary statistics
x_index <- 21
y_index <- 22

# Get the dat from online repository
library(data.table)
raw_data <- fread('http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/german/german.data-numeric', data.table = FALSE)

statlog <- list()
statlog$y <- raw_data[, 25] - 1
statlog$x <- cbind(1, scale(raw_data[, 1:24]))

# Bayesian logit in RStan
train_dat <- list(n = length(statlog$y), p = ncol(statlog$x), x = statlog$x, y = statlog$y, beta_sd = prior_sd)

stan_file <- "bayes_logit.stan"
bayes_log_reg <- rstan::stan(stan_file, data = train_dat, seed = seed,
                             iter = n_bootstrap * 2, chains = 1)
stan_bayes_sample <- rstan::extract(bayes_log_reg)$beta

# Variational Bayes in RStan
stan_model <- rstan::stan_model(file = stan_file)
stan_vb <- rstan::vb(object = stan_model, data = train_dat, seed = seed,
                     output_samples = n_bootstrap)
stan_vb_sample <- rstan::extract(stan_vb)$beta

带有模型的斯坦文件bayes_logit.stan是：

// Code for 0-1 loss Bayes Logistic Regression model
data 
  int<lower=0> n; // number of observations
  int<lower=0> p; // number of covariates
  matrix[n,p] x; // Matrix of covariates
  int<lower=0,upper=1> y[n]; // Responses
  real<lower=0> beta_sd; // Stdev of beta

parameters 
  vector[p] beta;

model 
  beta ~ normal(0,beta_sd);
  y ~ bernoulli_logit(x * beta); // Logistic regression

系数 21 和 22 的结果非常不同：

> mean(stan_bayes_sample[, 21])
[1] 0.1316655

> mean(stan_vb_sample[, 21])
[1] 0.3832403

> mean(stan_bayes_sample[, 22])
[1] -0.05473327

> mean(stan_vb_sample[, 22])
[1] 0.1570745

一张图清楚地显示了差异，其中点是精确推断，线条是变分推断的密度：

我在我的机器和 Azure 上得到了相同的结果。我已经注意到，当数据被缩放和居中时，精确推断会给出相同的结果，而变分推断会给出不同的结果，因此我可能会在不知不觉中触发数据处理的不同步骤。

更令人困惑的是，最近在 2019 年 5 月 30 日，使用相同版本的 RStan 的相同代码对两种方法给出了非常相似的结果，如下所示，其中红点大致位于同一位置，但蓝线的位置和比例不同（绿线是我正在实施的方法，我没有包括在最小的可重现示例中）：

有人有提示吗？

情节代码

情节的代码有点长：

requireNamespace("dplyr", quietly = TRUE)
requireNamespace("ggplot2", quietly = TRUE)
requireNamespace("tibble", quietly = TRUE)


#The first argument is required, either NULL or an arbitrary string.
stat_density_2d1_proto <- ggplot2::ggproto(NULL,
                                           ggplot2::Stat,
                                           required_aes = c("x", "y"),

                                           compute_group = function(data, scales, bins, n) 
                                             # Choose the bandwidth of Gaussian kernel estimators and increase it for
                                             # smoother densities in small sample sizes
                                             h <- c(MASS::bandwidth.nrd(data$x) * 1.5,
                                                    MASS::bandwidth.nrd(data$y) * 1.5)

                                             # Estimate two-dimensional density
                                             dens <- MASS::kde2d(
                                               data$x, data$y, h = h, n = n,
                                               lims = c(scales$x$dimension(), scales$y$dimension())
                                             )

                                             # Store in data frame
                                             df <- data.frame(expand.grid(x = dens$x, y = dens$y), z = as.vector(dens$z))

                                             # Add a label of this density for ggplot2
                                             df$group <- data$group[1]

                                             # plot
                                             ggplot2::StatContour$compute_panel(df, scales, bins)
                                           
)

# Wrap that ggproto in a ggplot2 object
stat_density_2d1 <- function(data = NULL,
                             geom = "density_2d",
                             position = "identity",
                             n = 100,
                             ...) 
  ggplot2::layer(
    data = data,
    stat = stat_density_2d1_proto,
    geom = geom,
    position = position,
    params = list(
      n = n,
      ...
    )
  )


append_to_plot <- function(plot_df, sample, method,
                           x_index, y_index) 
  new_plot_df <- rbind(plot_df, tibble::tibble(x = sample[, x_index],
                                               y = sample[, y_index],
                                               Method = method))
  return(new_plot_df)


plot_df <- tibble::tibble()

plot_df  <- append_to_plot(plot_df, sample = stan_bayes_sample,
                           method = "Bayes-Stan",
                           x_index = x_index, y_index = y_index)
plot_df  <- append_to_plot(plot_df, sample = stan_vb_sample,
                           method = "VB-Stan",
                           x_index = x_index, y_index = y_index)

ggplot2::ggplot(ggplot2::aes_string(x = "x", y = "y", colour = "Method"),
                data = dplyr::filter(plot_df, plot_df$Method != "Bayes-Stan")) +
  stat_density_2d1(bins = 5) +
  ggplot2::geom_point(alpha = 0.1, size = 1,
                      data = dplyr::filter(plot_df,
                                           plot_df$Method == "Bayes-Stan")) +
  ggplot2::theme_grey(base_size = 8) +
  ggplot2::xlab(bquote(beta[.(x_index)])) +
  ggplot2::ylab(bquote(beta[.(y_index)])) +
  ggplot2::theme(legend.position = "none",
                 plot.margin = ggplot2::margin(0, 10, 0, 0, "pt"))

【参考方案1】：

变分推理是一种近似算法，我们不希望它提供与通过 MCMC 实现的完整贝叶斯相同的答案。评估变分推理是否接近的最好读物是 Yuling Yao 及其同事的这篇 arXiv 论文，Yes, but does it work? Evaluating variational inference。 Bishop 的机器学习文本中很好地描述了近似值的工作原理。

我认为最近版本之间 Stan 的变分推理算法没有任何变化。与完全贝叶斯相比，变分推理对算法的参数和初始化更加敏感。这就是为什么它在我们所有的界面中仍然被标记为“实验性”。您可以尝试运行旧版本来控制初始化并确保有足够的迭代。变分推理在优化步骤中可能会非常失败，最终会得到次优的近似值。如果最佳变分近似不是很好，它也可能失败。

【讨论】：

我从版本控制确认它是相同的代码运行并给出不同的答案。 RStan 是否可以使用取决于时钟时间的随机初始化，这是两个实验之间唯一明显的变化？

我今天再次运行相同的代码，现在我发现与第一次尝试相同的结果。自从我发布问题以来，出现了一个新版本的 StanHeaders，所以也许新版本解决了这个问题。