Fit the model with Voigt peaks using Sequential Monte Carlo (SMC).

Source: R/fitVoigtPeaksSMC.R

Fit the model with Voigt peaks using Sequential Monte Carlo (SMC).

fitVoigtPeaksSMC(wl, spc, lPriors, conc = rep(1, nrow(spc)),
  npart = 10000, rate = 0.9, mcAR = 0.23, mcSteps = 10,
  minESS = npart/2, destDir = NA)

Arguments

wl

Vector of nwl wavenumbers at which the spetra are observed.
spc

n_y * nwl Matrix of observed Raman spectra.
lPriors

List of hyperparameters for the prior distributions.
conc

Vector of n_y nanomolar (nM) dye concentrations for each observation.
npart

number of SMC particles to use for the importance sampling distribution.
rate

the target rate of reduction in the effective sample size (ESS).
mcAR

target acceptance rate for the MCMC kernel
mcSteps

number of iterations of the MCMC kernel
minESS

minimum effective sample size, below which the particles are resampled.
destDir

destination directory to save intermediate results (for long-running computations)

Examples

wavenumbers <- seq(200,600,by=10)
spectra <- matrix(nrow=1, ncol=length(wavenumbers))
peakLocations <- c(300,500)
peakAmplitude <- c(10000,4000)
peakScale <- c(10, 15)
signature <- weightedLorentzian(peakLocations, peakScale, peakAmplitude, wavenumbers)
baseline <- 1000*cos(wavenumbers/200) + 2*wavenumbers
spectra[1,] <- signature + baseline + rnorm(length(wavenumbers),0,200)
lPriors <- list(scaG.mu=log(11.6) - (0.4^2)/2, scaG.sd=0.4, scaL.mu=log(11.6) - (0.4^2)/2,
   scaL.sd=0.4, bl.smooth=5, bl.knots=20, loc.mu=peakLocations, loc.sd=c(5,5),
   beta.mu=c(5000,5000), beta.sd=c(5000,5000), noise.sd=200, noise.nu=4)
result <- fitVoigtPeaksSMC(wavenumbers, spectra, lPriors, npart=50, mcSteps=1)
#> [1] "SMC with 1 observations at 1 unique concentrations, 50 particles, and 41 wavenumbers."
#> [1] "Step 0: computing 24 B-spline basis functions (r=21.0526315789474) took 0.289999999999999sec."
#> [1] "Mean noise parameter sigma is now 251.707012172913"
#> [1] "Mean spline penalty lambda is now 5"
#> [1] "Step 1: initialization for 2 Voigt peaks took 0.0199999999999996 sec."
#> [1] 6635.002 6180.201
#> [1] "Reweighting took 0sec. for ESS 45.0219890215427 with new kappa 0.02734375."
#> [1] "17 M-H proposals accepted. Temp ESS is 45.0219890215427"
#> [1] 6263.616 6076.730
#> [1] "Iteration 2 took 0.0199999999999996sec. for 1 MCMC loops (acceptance rate 0.34)"
#> [1] "Reweighting took 0sec. for ESS 40.6913327506059 with new kappa 0.04254150390625."
#> [1] "11 M-H proposals accepted. Temp ESS is 40.6913327506059"
#> [1] 5998.973 5780.276
#> [1] "Iteration 3 took 0sec. for 1 MCMC loops (acceptance rate 0.22)"
#> [1] "Reweighting took 0sec. for ESS 36.3819571500594 with new kappa 0.0575017929077148."
#> [1] "10 M-H proposals accepted. Temp ESS is 36.3819571500594"
#> [1] 5914.771 5548.921
#> [1] "Iteration 4 took 0.0199999999999996sec. for 1 MCMC loops (acceptance rate 0.2)"
#> [1] "Reweighting took 0sec. for ESS 32.3199226209689 with new kappa 0.0722283273935318."
#> [1] "11 M-H proposals accepted. Temp ESS is 32.3199226209689"
#> [1] 5468.343 5197.696
#> [1] "Iteration 5 took 0.0199999999999996sec. for 1 MCMC loops (acceptance rate 0.22)"
#> [1] "Reweighting took 0sec. for ESS 28.58510501772 with new kappa 0.0867247597780079."
#> [1] "11 M-H proposals accepted. Temp ESS is 28.58510501772"
#> [1] 5629.045 5066.258
#> [1] "Iteration 6 took 0sec. for 1 MCMC loops (acceptance rate 0.22)"
#> [1] "Reweighting took 0.0199999999999996sec. for ESS 24.8852320535079 with new kappa 0.100994685406476."
#> [1] "*** Resampling with 21 unique indices took 0 sec ***"
#> [1] "16 M-H proposals accepted. Temp ESS is 22.7272727272727"
#> [1] 4953.456 4185.926
#> [1] "Iteration 7 took 0.0199999999999996sec. for 1 MCMC loops (acceptance rate 0.32)"
#> [1] "Reweighting took 0sec. for ESS 20.1822388001957 with new kappa 0.115041643447."
#> [1] "*** Resampling with 23 unique indices took 0 sec ***"
#> [1] "11 M-H proposals accepted. Temp ESS is 18.6567164179104"
#> [1] 4454.265 3934.373
#> [1] "Iteration 8 took 0.0200000000000014sec. for 1 MCMC loops (acceptance rate 0.22)"
#> [1] "Reweighting took 0sec. for ESS 17.0638534353983 with new kappa 0.121955380607571."
#> [1] "*** Resampling with 27 unique indices took 0 sec ***"
#> [1] "9 M-H proposals accepted. Temp ESS is 25.5102040816327"
#> [1] 4390.723 3725.357
#> [1] "Iteration 9 took 0.0199999999999996sec. for 1 MCMC loops (acceptance rate 0.18)"
#> [1] "Reweighting took 0sec. for ESS 22.461867856145 with new kappa 0.135674827785577."
#> [1] "*** Resampling with 22 unique indices took 0 sec ***"
#> [1] "12 M-H proposals accepted. Temp ESS is 21.551724137931"
#> [1] 5103.714 3876.222
#> [1] "Iteration 10 took 0sec. for 1 MCMC loops (acceptance rate 0.24)"
#> [1] "Reweighting took 0sec. for ESS 19.9746123938649 with new kappa 0.149179908601428."
#> [1] "*** Resampling with 18 unique indices took 0 sec ***"
#> [1] "12 M-H proposals accepted. Temp ESS is 19.53125"
#> [1] 5553.547 3681.095
#> [1] "Iteration 11 took 0sec. for 1 MCMC loops (acceptance rate 0.24)"
#> [1] "Reweighting took 0sec. for ESS 16.7016852074047 with new kappa 0.175768036457633."
#> [1] "*** Resampling with 22 unique indices took 0 sec ***"
#> [1] "8 M-H proposals accepted. Temp ESS is 19.2307692307692"
#> [1] 5543.164 3722.722
#> [1] "Iteration 12 took 0.00999999999999979sec. for 1 MCMC loops (acceptance rate 0.16)"
#> [1] "Reweighting took 0sec. for ESS 17.6771492494723 with new kappa 0.188646660887983."
#> [1] "*** Resampling with 17 unique indices took 0 sec ***"
#> [1] "16 M-H proposals accepted. Temp ESS is 25.5102040816327"
#> [1] 5720.547 3180.812
#> [1] "Iteration 13 took 0sec. for 1 MCMC loops (acceptance rate 0.32)"
#> [1] "Reweighting took 0sec. for ESS 23.1239742072579 with new kappa 0.239356244582484."
#> [1] "*** Resampling with 25 unique indices took 0 sec ***"
#> [1] "9 M-H proposals accepted. Temp ESS is 20.1612903225806"
#> [1] 6298.792 3004.667
#> [1] "Iteration 14 took 0sec. for 1 MCMC loops (acceptance rate 0.18)"
#> [1] "Reweighting took 0sec. for ESS 18.336896992616 with new kappa 0.286896479296078."
#> [1] "*** Resampling with 22 unique indices took 0 sec ***"
#> [1] "14 M-H proposals accepted. Temp ESS is 17.8571428571429"
#> [1] 6449.332 3051.928
#> [1] "Iteration 15 took 0sec. for 1 MCMC loops (acceptance rate 0.28)"
#> [1] "Reweighting took 0.0199999999999996sec. for ESS 16.1983228556296 with new kappa 0.331465449340073."
#> [1] "*** Resampling with 19 unique indices took 0 sec ***"
#> [1] "13 M-H proposals accepted. Temp ESS is 20.1612903225806"
#> [1] 6597.955 2861.793
#> [1] "Iteration 16 took 0.0199999999999996sec. for 1 MCMC loops (acceptance rate 0.26)"
#> [1] "Reweighting took 0sec. for ESS 17.9813650806255 with new kappa 0.352357154048196."
#> [1] "*** Resampling with 26 unique indices took 0 sec ***"
#> [1] "15 M-H proposals accepted. Temp ESS is 21.1864406779661"
#> [1] 6765.555 2978.587
#> [1] "Iteration 17 took 0.0200000000000014sec. for 1 MCMC loops (acceptance rate 0.3)"
#> [1] "Reweighting took 0.00999999999999979sec. for ESS 19.3603429285685 with new kappa 0.362476573516193."
#> [1] "*** Resampling with 26 unique indices took 0 sec ***"
#> [1] "5 M-H proposals accepted. Temp ESS is 18.1159420289855"
#> [1] 6848.542 2928.699
#> [1] "Iteration 18 took 0.00999999999999979sec. for 1 MCMC loops (acceptance rate 0.1)"
#> [1] "Reweighting took 0sec. for ESS 16.3145035262298 with new kappa 0.382399180593812."
#> [1] "*** Resampling with 24 unique indices took 0 sec ***"
#> [1] "15 M-H proposals accepted. Temp ESS is 26.0416666666667"
#> [1] 6901.877 2881.806
#> [1] "Iteration 19 took 0sec. for 1 MCMC loops (acceptance rate 0.3)"
#> [1] "Reweighting took 0sec. for ESS 23.9521036768412 with new kappa 0.401699206200256."
#> [1] "*** Resampling with 24 unique indices took 0 sec ***"
#> [1] "14 M-H proposals accepted. Temp ESS is 26.0416666666667"
#> [1] 6733.589 3028.415
#> [1] "Iteration 20 took 0.0199999999999996sec. for 1 MCMC loops (acceptance rate 0.28)"
#> [1] "Reweighting took 0.00999999999999979sec. for ESS 22.7662892950956 with new kappa 0.43909300581274."
#> [1] "*** Resampling with 22 unique indices took 0 sec ***"
#> [1] "8 M-H proposals accepted. Temp ESS is 17.6056338028169"
#> [1] 7097.211 3013.068
#> [1] "Iteration 21 took 0.00999999999999979sec. for 1 MCMC loops (acceptance rate 0.16)"
#> [1] "Reweighting took 0sec. for ESS 14.9787087456132 with new kappa 0.474149692949443."
#> [1] "*** Resampling with 20 unique indices took 0 sec ***"
#> [1] "12 M-H proposals accepted. Temp ESS is 13.1578947368421"
#> [1] 7377.151 3067.840
#> [1] "Iteration 22 took 0sec. for 1 MCMC loops (acceptance rate 0.24)"
#> [1] "Reweighting took 0sec. for ESS 11.374196803767 with new kappa 0.507015337140103."
#> [1] "*** Resampling with 17 unique indices took 0 sec ***"
#> [1] "12 M-H proposals accepted. Temp ESS is 10.2459016393443"
#> [1] 7678.229 3230.539
#> [1] "Iteration 23 took 0sec. for 1 MCMC loops (acceptance rate 0.24)"
#> [1] "Reweighting took 0sec. for ESS 9.09202415093272 with new kappa 0.537826878568847."
#> [1] "*** Resampling with 16 unique indices took 0 sec ***"
#> [1] "11 M-H proposals accepted. Temp ESS is 11.3636363636364"
#> [1] 7857.739 3300.104
#> [1] "Iteration 24 took 0sec. for 1 MCMC loops (acceptance rate 0.22)"
#> [1] "Reweighting took 0sec. for ESS 9.96441413979713 with new kappa 0.595598518747741."
#> [1] "*** Resampling with 16 unique indices took 0 sec ***"
#> [1] "16 M-H proposals accepted. Temp ESS is 12.7551020408163"
#> [1] 8050.136 3148.604
#> [1] "Iteration 25 took 0sec. for 1 MCMC loops (acceptance rate 0.32)"
#> [1] "Reweighting took 0sec. for ESS 11.3063638020647 with new kappa 0.79779925937387."
#> [1] "*** Resampling with 14 unique indices took 0 sec ***"
#> [1] "13 M-H proposals accepted. Temp ESS is 15.625"
#> [1] 8743.733 3133.301
#> [1] "Iteration 26 took 0.0199999999999996sec. for 1 MCMC loops (acceptance rate 0.26)"
#> [1] "Reweighting took 0sec. for ESS 13.1415203588948 with new kappa 0.898899629686935."
#> [1] "*** Resampling with 20 unique indices took 0 sec ***"
#> [1] "7 M-H proposals accepted. Temp ESS is 13.8888888888889"
#> [1] 8747.184 3105.890
#> [1] "Iteration 27 took 0.0199999999999996sec. for 1 MCMC loops (acceptance rate 0.14)"
#> [1] "Reweighting took 0sec. for ESS 12.4626912816888 with new kappa 0.949449814843468."
#> [1] "*** Resampling with 19 unique indices took 0 sec ***"
#> [1] "11 M-H proposals accepted. Temp ESS is 15.4320987654321"
#> [1] 8879.895 3051.112
#> [1] "Iteration 28 took 0.0200000000000014sec. for 1 MCMC loops (acceptance rate 0.22)"
#> [1] "Reweighting took 0sec. for ESS 14.060139388755 with new kappa 1."
#> [1] "*** Resampling with 20 unique indices took 0 sec ***"
#> [1] "21 M-H proposals accepted. Temp ESS is 18.9393939393939"
#> [1] 8818.267 3008.041
#> [1] "Iteration 29 took 0.0199999999999996sec. for 1 MCMC loops (acceptance rate 0.42)"