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:
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):
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"))