Adaptive non-parametric learning
Sampling
We assign an independent normal prior with standard deviation 10 to each covariate, and generate 1000 posterior samples for each method: Bayesian logistic regression, variational Bayes, and Adaptive Non-Parametric Learning.
Note: variational Bayes in RStan is still experimental. We have had situations in the past where the same code with the same versions produces very different results. This thread on StackOverflow mentions this problem. In this code, variational Bayes provides the samples for our method, so the results may vary in the future because of the results from variational Bayes in RStan. Whether variational Bayes gives the expected answers or not, our approach will always give an answer that is in-between the results from variational Bayes and Bayesian logistic regression.
We tune the settings for the sampling with the number of draws, setting the seed, and getting the data:
seed <- 123
prior_sd <- 10
n_bootstrap <- 1000
german <- PosteriorBootstrap::get_german_credit_dataset()For Bayesian logistic regression, we draw samples using the
rstan package and the Bayesian logistic regression model
that ships with this package:
if ("Windows" != Sys.info()["sysname"]) {
train_dat <- list(n = length(german$y), p = ncol(german$x), x = german$x, y = german$y, beta_sd = prior_sd)
stan_file <- PosteriorBootstrap::get_stan_file()
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
}
#>
#> SAMPLING FOR MODEL 'bayes_logit' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.000198 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.98 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 1.17893 seconds (Warm-up)
#> Chain 1: 1.3607 seconds (Sampling)
#> Chain 1: 2.53963 seconds (Total)
#> Chain 1:For variational Bayes, we obtain a mean-field variational Bayes
sample using automatic differentiation variational inference using the
rstan package and the same logistic regression model as
above. The number of samples is n_bootstrap, as these
samples serve for the Adaptive Non-Parametric Learning algorithm:
if ("Windows" != Sys.info()["sysname"]) {
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
}
#> Chain 1: ------------------------------------------------------------
#> Chain 1: EXPERIMENTAL ALGORITHM:
#> Chain 1: This procedure has not been thoroughly tested and may be unstable
#> Chain 1: or buggy. The interface is subject to change.
#> Chain 1: ------------------------------------------------------------
#> Chain 1:
#> Chain 1:
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.00012 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.2 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Begin eta adaptation.
#> Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
#> Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
#> Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
#> Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
#> Chain 1: Iteration: 200 / 250 [ 80%] (Adaptation)
#> Chain 1: Success! Found best value [eta = 1] earlier than expected.
#> Chain 1:
#> Chain 1: Begin stochastic gradient ascent.
#> Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
#> Chain 1: 100 -593.338 1.000 1.000
#> Chain 1: 200 -591.578 0.501 1.000
#> Chain 1: 300 -593.057 0.335 0.003 MEDIAN ELBO CONVERGED
#> Chain 1:
#> Chain 1: Drawing a sample of size 1000 from the approximate posterior...
#> Chain 1: COMPLETED.
#> Warning: Pareto k diagnostic value is 1.23. Resampling is disabled. Decreasing
#> tol_rel_obj may help if variational algorithm has terminated prematurely.
#> Otherwise consider using sampling instead.We now run the Adaptive Non-Parametric Learning algorithm, using the
samples from variational Bayes, for two different values of the
concentration hyper-parameter. Because our sampling method with
c = 1000 takes longer than the timeout allowed on CRAN, we
use only c = 500 for this vignette, but please see the GitHub
page for the same figure as in the paper (with c = 1,
c = 1,000, and c = 20,000).
if ("Windows" != Sys.info()["sysname"]) {
if (identical(Sys.getenv("NOT_CRAN", unset = "true"), "true")) {
concentrations <- c(1, 500)
anpl_samples <- list()
for (concentration in concentrations) {
anpl_sample <- PosteriorBootstrap::draw_logit_samples(x = german$x, y = german$y,
concentration = concentration,
n_bootstrap = n_bootstrap,
posterior_sample = stan_vb_sample,
threshold = 1e-8,
show_progress = TRUE)
anpl_samples[[toString(concentration)]] <- anpl_sample
}
}
}
#> ================================================================================================================================================================Preparing the plot
We now prepare a dataframe with all the samples ready for plotting.
if ("Windows" != Sys.info()["sysname"]) {
if (identical(Sys.getenv("NOT_CRAN", unset = "true"), "true")) {
# Initialise
plot_df <- tibble::tibble()
# Index of coefficients in the plot
x_index <- 21
y_index <- 22
# Create a plot data frame with all the samples
for (concentration in concentrations) {
plot_df <- append_to_plot(plot_df, sample = anpl_samples[[toString(concentration)]],
method = "PosteriorBootstrap-ANPL",
concentration = concentration,
x_index = x_index, y_index = y_index)
plot_df <- append_to_plot(plot_df, sample = stan_bayes_sample,
method = "Bayes-Stan",
concentration = concentration,
x_index = x_index, y_index = y_index)
plot_df <- append_to_plot(plot_df, sample = stan_vb_sample,
method = "VB-Stan",
concentration = concentration,
x_index = x_index, y_index = y_index)
}
}
}Plotting the results
And now we plot the result of the algorithm:
if ("Windows" != Sys.info()["sysname"]) {
if (identical(Sys.getenv("NOT_CRAN", unset = "true"), "true")) {
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::facet_wrap(~concentration, nrow = 1,
scales = "fixed",
labeller = ggplot2::label_bquote(c ~" = "~
.(concentration))
) +
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"))
}
}
#> Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
#> ℹ Please use tidy evaluation ideoms with `aes()`
The Non-parametric update effectively corrects the variational Bayes
approximation for small values of the concentration c. As
we increase the concentration parameter, the samples drawn with our
method move closer to those from the variational Bayes approximation,
showing that the concentration parameter is an effective sample size
tuning of the weight to give to the data compared to the approximate
model.
Of course, local variational methods can provide more accurate posterior approximations to Bayesian logistic posteriors, though these too are approximations, that our algorithm can correct.
