2

I am working on an R package that depends on RStan and I seem to have hit a failure mode in the latter.

I run a Bayesian logistic regression with exact inference (rstan::stan()) and get very different results with variational inference (rstan::vb()). The following code downloads the German Statlog Credit data and runs both inferences on that data:

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

The Stan file bayes_logit.stan with the model is:

// 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
}

The results for coefficients 21 and 22 are very different:

> 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

And a plot clearly shows the difference, where the dots are exact inference and the lines are the density for variational inference:

Plot of exact and variational inferences

I get the same results on my machine and on Azure. I have noted that exact inference gives the same results when the data is scaled and centered and variational inference gives different results, so I may unwittingly trigger a different step of data processing.

Even more confusing is that the same code with the same version of RStan, as recently as May 30th 2019, was giving very similar results for the two methods, as shown below, where the red dots are roughly in the same place but the blue lines are different in location and scale (and the green lines are for the method I am implementing, which I did not include in the minimal reproducible example):

Plot of exact and variational inference in a previous version

Does anyone have a hint?

Code for the plot

The code for the plot is a bit long:

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 Answer 1

4

Variational inference is an approximate algorithm and we don't expect it to provide the same answer as full Bayes implemented through MCMC. The best thing to read on evaluating whether variational inference even gets close is this arXiv paper by Yuling Yao and colleagues, Yes, but does it work? Evaluating variational inference. There's a good description of how the approximations work in Bishop's machine learning text.

I don't think anything has changed in Stan's variational inference algorithm between versions recently. Variational inference can be much more sensitive to the parameters of the algorithm and to initializations than full Bayes. That's why it's still marked as "experimental" in all of our interfaces. You might try running old versions controlling for initialization and making sure there are enough iterations. Variational inference can fail pretty badly on the optimization step, winding up with suboptimal approximations. It can also fail if the best variational approximation is not very good.

2
  • I confirm from version control that it's the same code running and giving different answers. Could RStan use a random initialization that depends on clock time, which is the only obvious change in between the two experiments? Jul 26, 2019 at 14:13
  • I ran the same code again today and now I find the same results as in my first attempt. A new version of StanHeaders appeared since I posted the question, so maybe new versions fixed the problem. Sep 11, 2019 at 21:54

Your Answer

Reminder: Answers generated by Artificial Intelligence tools are not allowed on Stack Overflow. Learn more

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.