R/fitVoigtPeaksSMC.R
fitVoigtPeaksSMC.Rd
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)
wl | Vector of |
---|---|
spc |
|
lPriors | List of hyperparameters for the prior distributions. |
conc | Vector of |
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) |
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)"