Simple Example

library("spectralmc")

logprior <- LogPrior(
  amp = prior.Positive(),
  lwidth = prior.Gamma(1.01, 0.01), #shape rate
  gwidth = prior.Gamma(1.01, 0.01), #shape rate
  pos = prior.Unif(0, 1000), # min, max
  a = prior.Norm(0, 0.00001), # mean, sd
  b = prior.Unif(-1000, 1000)
)

proposal <- Proposal(
  amp = prop.Norm(0.05), #sd
  lwidth = prop.Norm(0.2),
  gwidth = prop.Norm(0.2), # can also use prop.Unif(radius) for uniform proposal
  pos = prop.Norm(1),
  a = prop.Norm(0.0000005),
  b = prop.Norm(0.02) #b = prop.Unif(0.02)
)

true_kp <- 3
fit_kp <- 3
noise_sigma <- 0.001
signal_model <- voigt.model

data <- data.synthetic.stormyclouds(seed = 2021, kp = true_kp, noise = 0.001)
x <- data$x
y <- data$y

posterior <- Posterior(Loglikeli(signal_model, noise_sigma = noise_sigma), logprior)
start_params <- prop.voigt.rnd_start(seed = 33, kp = fit_kp)

res <- chains.mhmc(x, y, start_params, posterior, proposal, iter = 20000, print_bar = FALSE)

  |================================================================================| 100% elap: 26s eta:  0s

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

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

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

# can also do normal coda plots
library(coda)
samples_matrix <- utils.samples_to_matrix(res$samples) #can also filter variables with which=c("pos", "amp")
chain <- mcmc(samples_matrix)
summary(chain)

Iterations = 1:20000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 20000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                Mean        SD  Naive SE Time-series SE
amp[1]     4.195e+00 2.055e+00 1.453e-02      2.078e+00
amp[2]     3.445e+00 1.667e+00 1.179e-02      1.464e+00
amp[3]     5.293e+00 2.536e+00 1.793e-02      2.395e+00
lwidth[1]  1.726e+01 9.138e+00 6.462e-02      7.861e+00
lwidth[2]  2.155e+01 2.177e+00 1.539e-02      8.107e-01
lwidth[3]  2.957e+01 4.973e+00 3.517e-02      4.575e+00
gwidth[1]  2.676e+01 1.297e+01 9.172e-02      1.281e+01
gwidth[2]  2.866e+01 2.826e+00 1.999e-02      1.827e+00
gwidth[3]  2.414e+01 7.403e+00 5.235e-02      8.336e+00
pos[1]     2.208e+02 1.389e+01 9.824e-02      8.791e+00
pos[2]     8.580e+02 1.268e+01 8.965e-02      5.444e+00
pos[3]     5.355e+02 3.743e+00 2.647e-02      7.081e-01
a         -7.130e-07 7.864e-06 5.561e-08      4.981e-06
b          1.820e+01 3.718e+00 2.629e-02      1.514e+00

2. Quantiles for each variable:

                2.5%        25%        50%       75%     97.5%
amp[1]     7.782e-01  1.961e+00  4.824e+00 6.072e+00 6.760e+00
amp[2]     5.500e-01  1.996e+00  3.907e+00 4.868e+00 5.534e+00
amp[3]     1.862e-01  3.621e+00  6.179e+00 7.456e+00 7.867e+00
lwidth[1]  1.927e+00  9.695e+00  1.985e+01 2.475e+01 2.934e+01
lwidth[2]  1.860e+01  1.957e+01  2.149e+01 2.303e+01 2.643e+01
lwidth[3]  2.241e+01  2.366e+01  3.066e+01 3.345e+01 3.640e+01
gwidth[1]  2.284e+00  1.875e+01  3.215e+01 3.809e+01 4.081e+01
gwidth[2]  2.397e+01  2.668e+01  2.732e+01 3.151e+01 3.395e+01
gwidth[3]  1.385e+01  1.606e+01  2.516e+01 3.128e+01 3.447e+01
pos[1]     2.090e+02  2.116e+02  2.129e+02 2.294e+02 2.560e+02
pos[2]     8.126e+02  8.524e+02  8.623e+02 8.651e+02 8.713e+02
pos[3]     5.282e+02  5.334e+02  5.355e+02 5.377e+02 5.444e+02
a         -1.255e-05 -3.602e-06 -2.044e-06 4.001e-07 2.042e-05
b          4.084e+00  1.938e+01  1.938e+01 1.939e+01 1.939e+01

plot(chain)

autocorr.plot(chain)

---

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ß) ↩︎