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.0159999999999911sec." #> [1] "Mean noise parameter sigma is now 237.135625157295" #> [1] "Mean spline penalty lambda is now 5" #> [1] "Step 1: initialization for 2 Voigt peaks took 0.0180000000000007 sec." #> [1] 7592.846 5846.064 #> [1] "Reweighting took 0.00100000000000477sec. for ESS 45.4666155840655 with new kappa 0.01953125." #> [1] "13 M-H proposals accepted. Temp ESS is 45.4666155840655" #> [1] 7759.770 5824.073 #> [1] "Iteration 2 took 0.007000000000005sec. for 1 MCMC loops (acceptance rate 0.26)" #> [1] "Reweighting took 0.00200000000000955sec. for ESS 40.1966003739475 with new kappa 0.0310211181640625." #> [1] "18 M-H proposals accepted. Temp ESS is 40.1966003739475" #> [1] 7258.658 5316.871 #> [1] "Iteration 3 took 0.007000000000005sec. for 1 MCMC loops (acceptance rate 0.36)" #> [1] "Reweighting took 0.0029999999999859sec. for ESS 37.0342443826065 with new kappa 0.0385912656784058." #> [1] "2 M-H proposals accepted. Temp ESS is 37.0342443826065" #> [1] 7359.558 5327.796 #> [1] "Iteration 4 took 0.007000000000005sec. for 1 MCMC loops (acceptance rate 0.04)" #> [1] "Reweighting took 0.00199999999998113sec. for ESS 34.0233554769581 with new kappa 0.0461022714152932." #> [1] "16 M-H proposals accepted. Temp ESS is 34.0233554769581" #> [1] 7078.89 4961.64 #> [1] "Iteration 5 took 0.00600000000000023sec. for 1 MCMC loops (acceptance rate 0.32)" #> [1] "Reweighting took 0.00200000000000955sec. for ESS 31.5095188336404 with new kappa 0.0535545974198612." #> [1] "11 M-H proposals accepted. Temp ESS is 31.5095188336404" #> [1] 6875.868 5055.978 #> [1] "Iteration 6 took 0.00899999999998613sec. for 1 MCMC loops (acceptance rate 0.22)" #> [1] "Reweighting took 0.0040000000000191sec. for ESS 29.1178223709203 with new kappa 0.0609487021275186." #> [1] "9 M-H proposals accepted. Temp ESS is 29.1178223709203" #> [1] 6224.311 4546.760 #> [1] "Iteration 7 took 0.0109999999999957sec. for 1 MCMC loops (acceptance rate 0.18)" #> [1] "Reweighting took 0.00100000000000477sec. for ESS 27.0547716546321 with new kappa 0.0682850403921473." #> [1] "20 M-H proposals accepted. Temp ESS is 27.0547716546321" #> [1] 6064.248 4488.660 #> [1] "Iteration 8 took 0.00600000000000023sec. for 1 MCMC loops (acceptance rate 0.4)" #> [1] "Reweighting took 0.00300000000001432sec. for ESS 25.0371041789325 with new kappa 0.0755640635140837." #> [1] "7 M-H proposals accepted. Temp ESS is 25.0371041789325" #> [1] 6039.314 4409.979 #> [1] "Iteration 9 took 0.00800000000000978sec. for 1 MCMC loops (acceptance rate 0.14)" #> [1] "Reweighting took 0.00100000000000477sec. for ESS 22.9706390766225 with new kappa 0.0827862192678799." #> [1] "*** Resampling with 22 unique indices took 0 sec ***" #> [1] "13 M-H proposals accepted. Temp ESS is 23.1481481481481" #> [1] 7976.527 3272.191 #> [1] "Iteration 10 took 0.00499999999999545sec. for 1 MCMC loops (acceptance rate 0.26)" #> [1] "Reweighting took 0.00100000000000477sec. for ESS 20.3257799062644 with new kappa 0.0971176845918193." #> [1] "*** Resampling with 20 unique indices took 0 sec ***" #> [1] "16 M-H proposals accepted. Temp ESS is 23.5849056603774" #> [1] 8114.949 2820.833 #> [1] "Iteration 11 took 0.00500000000002387sec. for 1 MCMC loops (acceptance rate 0.32)" #> [1] "Reweighting took 0.00100000000000477sec. for ESS 22.0637119499312 with new kappa 0.111225220770072." #> [1] "*** Resampling with 22 unique indices took 0 sec ***" #> [1] "11 M-H proposals accepted. Temp ESS is 18.6567164179104" #> [1] 8348.604 3143.621 #> [1] "Iteration 12 took 0.00600000000000023sec. for 1 MCMC loops (acceptance rate 0.22)" #> [1] "Reweighting took 0.00100000000000477sec. for ESS 16.687051741007 with new kappa 0.138999432621007." #> [1] "*** Resampling with 18 unique indices took 0 sec ***" #> [1] "12 M-H proposals accepted. Temp ESS is 18.6567164179104" #> [1] 8788.414 3487.884 #> [1] "Iteration 13 took 0.00500000000002387sec. for 1 MCMC loops (acceptance rate 0.24)" #> [1] "Reweighting took 0.00199999999998113sec. for ESS 15.841362583786 with new kappa 0.165905700351601." #> [1] "*** Resampling with 18 unique indices took 0 sec ***" #> [1] "17 M-H proposals accepted. Temp ESS is 19.53125" #> [1] 9085.682 3582.569 #> [1] "Iteration 14 took 0.00600000000000023sec. for 1 MCMC loops (acceptance rate 0.34)" #> [1] "Reweighting took 0.00100000000000477sec. for ESS 17.8472366764519 with new kappa 0.178938423783607." #> [1] "*** Resampling with 22 unique indices took 0 sec ***" #> [1] "6 M-H proposals accepted. Temp ESS is 18.3823529411765" #> [1] 9316.776 3549.084 #> [1] "Iteration 15 took 0.00499999999999545sec. for 1 MCMC loops (acceptance rate 0.12)" #> [1] "Reweighting took 0.00100000000000477sec. for ESS 17.0404851957438 with new kappa 0.191767510911988." #> [1] "*** Resampling with 21 unique indices took 0 sec ***" #> [1] "9 M-H proposals accepted. Temp ESS is 20.8333333333333" #> [1] 9516.545 3636.150 #> [1] "Iteration 16 took 0.00499999999999545sec. for 1 MCMC loops (acceptance rate 0.18)" #> [1] "Reweighting took 0.00100000000000477sec. for ESS 18.5162120825855 with new kappa 0.217024776195989." #> [1] "*** Resampling with 18 unique indices took 0 sec ***" #> [1] "8 M-H proposals accepted. Temp ESS is 13.2978723404255" #> [1] 9365.153 3560.630 #> [1] "Iteration 17 took 0.00499999999999545sec. for 1 MCMC loops (acceptance rate 0.16)" #> [1] "Reweighting took 0.00100000000000477sec. for ESS 11.0200628713916 with new kappa 0.241492751939864." #> [1] "*** Resampling with 15 unique indices took 0 sec ***" #> [1] "16 M-H proposals accepted. Temp ESS is 15.2439024390244" #> [1] 9647.994 3643.681 #> [1] "Iteration 18 took 0.00499999999999545sec. for 1 MCMC loops (acceptance rate 0.32)" #> [1] "Reweighting took 0.00100000000000477sec. for ESS 13.3315029217404 with new kappa 0.265196103441743." #> [1] "*** Resampling with 19 unique indices took 0 sec ***" #> [1] "14 M-H proposals accepted. Temp ESS is 11.9047619047619" #> [1] 9745.161 3661.574 #> [1] "Iteration 19 took 0.007000000000005sec. for 1 MCMC loops (acceptance rate 0.28)" #> [1] "Reweighting took 0.00100000000000477sec. for ESS 10.2876198527751 with new kappa 0.288158725209189." #> [1] "*** Resampling with 20 unique indices took 0 sec ***" #> [1] "20 M-H proposals accepted. Temp ESS is 15.4320987654321" #> [1] 9781.120 3587.293 #> [1] "Iteration 20 took 0.00600000000000023sec. for 1 MCMC loops (acceptance rate 0.4)" #> [1] "Reweighting took 0.00100000000000477sec. for ESS 14.0699301827639 with new kappa 0.310403765046402." #> [1] "*** Resampling with 27 unique indices took 0 sec ***" #> [1] "14 M-H proposals accepted. Temp ESS is 17.8571428571429" #> [1] 9921.688 3575.921 #> [1] "Iteration 21 took 0.007000000000005sec. for 1 MCMC loops (acceptance rate 0.28)" #> [1] "Reweighting took 0.00100000000000477sec. for ESS 15.1855866087995 with new kappa 0.353503529731001." #> [1] "*** Resampling with 20 unique indices took 0 sec ***" #> [1] "3 M-H proposals accepted. Temp ESS is 11.6822429906542" #> [1] 9786.734 3609.509 #> [1] "Iteration 22 took 0.00600000000000023sec. for 1 MCMC loops (acceptance rate 0.06)" #> [1] "Reweighting took 0.000999999999976353sec. for ESS 9.87502166026076 with new kappa 0.434315588514626." #> [1] "*** Resampling with 13 unique indices took 0 sec ***" #> [1] "11 M-H proposals accepted. Temp ESS is 9.12408759124088" #> [1] 9805.618 3628.044 #> [1] "Iteration 23 took 0.00599999999997181sec. for 1 MCMC loops (acceptance rate 0.22)" #> [1] "Reweighting took 0sec. for ESS 8.94276542291468 with new kappa 0.57573669138597." #> [1] "*** Resampling with 14 unique indices took 0 sec ***" #> [1] "10 M-H proposals accepted. Temp ESS is 8.62068965517241" #> [1] 9742.735 3586.530 #> [1] "Iteration 24 took 0.00499999999999545sec. for 1 MCMC loops (acceptance rate 0.2)" #> [1] "Reweighting took 0sec. for ESS 7.82298623208552 with new kappa 0.681802518539477." #> [1] "*** Resampling with 18 unique indices took 0 sec ***" #> [1] "15 M-H proposals accepted. Temp ESS is 13.7362637362637" #> [1] 9770.341 3602.253 #> [1] "Iteration 25 took 0.00499999999999545sec. for 1 MCMC loops (acceptance rate 0.3)" #> [1] "Reweighting took 0sec. for ESS 12.3958510074088 with new kappa 0.721577203722042." #> [1] "*** Resampling with 22 unique indices took 0 sec ***" #> [1] "15 M-H proposals accepted. Temp ESS is 17.3611111111111" #> [1] 9739.542 3568.772 #> [1] "Iteration 26 took 0.00699999999997658sec. for 1 MCMC loops (acceptance rate 0.3)" #> [1] "Reweighting took 0sec. for ESS 14.8984795713768 with new kappa 0.860788601861021." #> [1] "*** Resampling with 19 unique indices took 0 sec ***" #> [1] "9 M-H proposals accepted. Temp ESS is 15.4320987654321" #> [1] 9749.091 3519.827 #> [1] "Iteration 27 took 0.00499999999999545sec. for 1 MCMC loops (acceptance rate 0.18)" #> [1] "Reweighting took 0sec. for ESS 15.6206834105581 with new kappa 1." #> [1] "*** Resampling with 19 unique indices took 0 sec ***" #> [1] "5 M-H proposals accepted. Temp ESS is 13.1578947368421" #> [1] 9796.170 3474.857 #> [1] "Iteration 28 took 0.00499999999999545sec. for 1 MCMC loops (acceptance rate 0.1)"