Collapsed Gibbs Sampler with GSS Slice Sampling

1 Description

Uses a collapsed (marginalizes over a, b, amp numerically) gibbs sampler which samples from the conditionals with slice sampling. Also re-parameterizes density with (sqrt(gwidth^2+lwidth^2), gwidth-lwidth) instead of (gwidth, lwitdh) if min(gwidth/lwidth, lwidth/gwidth) > 0.05 (otherwise transformation becomes unstable). Finds global maxima in posterior using population based golden section search and uses it as the starting point for the slice sampler. About x10 times faster than the previous implementation. Using Zell Prior could be another x10.

2 Data

library('plotly')

Loading required package: ggplot2

Attaching package: 'plotly'

The following object is masked from 'package:ggplot2':

    last_plot

The following object is masked from 'package:stats':

    filter

The following object is masked from 'package:graphics':

    layout

library(spectralmc)

#toluene data
cdata <- get_toluene_measured_data()
x <- cdata$x
y <- cdata$y

3 Algorithm

peak_spawner <- function() {
  params <- list(
          amp = rnorm(1, 1, 1), # doesnt matter
          lwidth = rnorm(1, 0.1, 0.001), # matters!!!
          gwidth = rnorm(1, 5, 0.001), # matters!!!
          pos = runif(1, 1, 1), # doesnt matter
          a = 0, # doesnt matter
          b = 0 # doesnt matter
  )
  return(params)
}

fit_kp <- 30
max_a <- 1e-2
lower <- c(-Inf, -max_a, 0)
upper <- c(Inf, max_a, Inf)
noise_sigma <- 100 # estimated by eye; intentionally to big such that tiny peaks are overnoised
signal_model <- voigt.model

res <- chains.rjmcmc_like_collapsed_slicer_gibbs(
        x,
        y,
        signal_model,
        peak_spawner,
        lower,
        upper,
        noise_sigma,
        fit_kp = fit_kp,
        add_peak_every = 10,
        iter = 20000,
        print_bar = TRUE,
        half_steps = TRUE,
        decorr_trafo = TRUE,
        save_loc = save_loc,
        checkpoint_every = 50,
)

4 Results

4.1 Fit

4.1.1 All

plt.fit.chain.interactive(x, y, res$samples, signal_model, step_width = 200)

4.1.2 Burn-in

plt.fit.chain.interactive(x, y, res$samples[seq(1, 500)], signal_model, step_width = 5)

4.1.3 Traceplots

plt.traceplot(res$samples, "pos")

plt.traceplot(res$samples, "amp")

plt.traceplot(res$samples, "lwidth")

plt.traceplot(res$samples, "gwidth")

plt.traceplot(res$samples, "a")

plt.traceplot(res$samples, "b")

4.1.4 Posterior plot

plt.posterior(res$post_vals, steps = NULL)

4.1.5 Fit vs True

plt.samples_vs_true(samples, 'gwidth')

plt.samples_vs_true(samples, 'lwidth')

5 Todo

• Critical Bug amp, a, b are not sampled instead they are set to roughly their means, easy to fix.
• Not a bug: Uses a safeguard which prevents messing up the fit if the slice sampler fails, this however destroys the gibbs sampler property and causes rejection plateaus.

---

1. RWTH Aachen, philipp.meissner@rwth-aachen.de↩︎

2. This work was developed under a Seed Fund Project (2021) of the RWTH Aachen University, funded under "the Excellence Strategy of the Federal Government and the Länder". Project: "Spectra-Bayes: A Bayesian statistical machine learning model for spectral reconstruction" (Project Leaders: Prof. Dr. Maria Kateri, Prof. Dr.-Ing. Hans-Jürgen Koß) ↩︎