<p>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, <a href="https://arxiv.org/abs/1802.02538" rel="nofollow noreferrer">Yes, but does it work? Evaluating variational inference</a>.  There's a good description of how the approximations work in Bishop's machine learning text.</p>

<p>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.  </p>
RStan gives different results in exact and variational bayes modes
<p>I am working on an <a href="https://github.com/alan-turing-institute/PosteriorBootstrap/" rel="nofollow noreferrer">R package</a> that depends on RStan and I seem to have hit a failure mode in the latter.</p>

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

<pre class="lang-r prettyprint-override"><code>library("rstan")

seed &lt;- 123
prior_sd &lt;- 10
n_bootstrap &lt;- 1000

# Index of coefficients in the plot and summary statistics
x_index &lt;- 21
y_index &lt;- 22

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

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

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

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

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

<p>The Stan file <code>bayes_logit.stan</code> with the model is:</p>

<pre><code>// Code for 0-1 loss Bayes Logistic Regression model
data {
  int&lt;lower=0&gt; n; // number of observations
  int&lt;lower=0&gt; p; // number of covariates
  matrix[n,p] x; // Matrix of covariates
  int&lt;lower=0,upper=1&gt; y[n]; // Responses
  real&lt;lower=0&gt; beta_sd; // Stdev of beta
}
parameters {
  vector[p] beta;
}
model {
  beta ~ normal(0,beta_sd);
  y ~ bernoulli_logit(x * beta); // Logistic regression
}
</code></pre>

<p>The results for coefficients 21 and 22 are very different:</p>

<pre><code>&gt; mean(stan_bayes_sample[, 21])
[1] 0.1316655

&gt; mean(stan_vb_sample[, 21])
[1] 0.3832403

&gt; mean(stan_bayes_sample[, 22])
[1] -0.05473327

&gt; mean(stan_vb_sample[, 22])
[1] 0.1570745
</code></pre>

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

<p><a href="https://i.stack.imgur.com/tZerD.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/tZerD.png" alt="Plot of exact and variational inferences"></a></p>

<p>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.</p>

<p>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):</p>

<p><a href="https://i.stack.imgur.com/Rn7SU.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/Rn7SU.png" alt="Plot of exact and variational inference in a previous version"></a></p>

<p>Does anyone have a hint?</p>

<h1>Code for the plot</h1>

<p>The code for the plot is a bit long:</p>

<pre class="lang-r prettyprint-override"><code>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 &lt;- 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 &lt;- c(MASS::bandwidth.nrd(data$x) * 1.5,
                                                    MASS::bandwidth.nrd(data$y) * 1.5)

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

                                             # Store in data frame
                                             df &lt;- 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 &lt;- data$group[1]

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

# Wrap that ggproto in a ggplot2 object
stat_density_2d1 &lt;- 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 &lt;- function(plot_df, sample, method,
                           x_index, y_index) {
  new_plot_df &lt;- rbind(plot_df, tibble::tibble(x = sample[, x_index],
                                               y = sample[, y_index],
                                               Method = method))
  return(new_plot_df)
}

plot_df &lt;- tibble::tibble()

plot_df  &lt;- append_to_plot(plot_df, sample = stan_bayes_sample,
                           method = "Bayes-Stan",
                           x_index = x_index, y_index = y_index)
plot_df  &lt;- 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"))

</code></pre>