Fit the model using Sequential Monte Carlo (SMC).
fitSpectraSMC(wl, spc, peakWL, lPriors, conc = rep(1, nrow(spc)), npart = 10000, rate = 0.9, minESS = npart/2, destDir = NA)
| wl | Vector of |
|---|---|
| spc |
|
| peakWL | Vector of locations for each peak (cm^-1) |
| 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). |
| minESS | minimum effective sample size, below which the particles are resampled. |
| destDir | destination directory to save intermediate results (for long-running computations) |
a List containing weighted parameter values, known as particles:
weightsVector of importance weights for each particle.
betanpart * npeaks Matrix of regression coefficients for the amplitudes.
scalenpart * npeaks Matrix of scale parameters.
sigmaVector of npart standard deviations.
alphabl.knots * n_y * npart Array of spline coefficients for the baseline.
basisA dense nwl * bl.knots Matrix containing the values of the basis functions.
expFnnpart * nwl Matrix containing the spectral signature.
essVector containing the effective sample size (ESS) at each SMC iteration.
logEvidenceVector containing the logarithm of the model evidence (marginal likelihood).
acceptVector containing the Metropolis-Hastings acceptance rate at each SMC iteration.
sd.mhniter * 2npeaks Matrix of random walk MH bandwidths at each SMC iteration..
Chopin (2002) "A Sequential Particle Filter Method for Static Models," Biometrika 89(3): 539--551, DOI: 10.1093/biomet/89.3.539
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(scale.mu=log(11.6) - (0.4^2)/2, scale.sd=0.4, bl.smooth=10^11, bl.knots=20, beta.mu=5000, beta.sd=5000, noise.sd=200, noise.nu=4) result <- fitSpectraSMC(wavenumbers, spectra, peakLocations, lPriors, npart=500)#> [1] "SMC with 1 observations at 1 unique concentrations, 500 particles, and 41 wavenumbers." #> [1] "Step 0: computing 20 B-spline basis functions took 0 sec." #> [1] "Step 1: initialization for 2 peaks took 0.0199999999999996 sec." #> [1] "Observation 1 20200206-171236" #> [1] "Baselines took 0 sec" #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 1 of 41) #> Required 1 iterations to reduce ESS from 500 to 124.347 #> [1] "reweighting took 0 sec for ESS 124.347094238812" #> [1] "*** Resampling with 141 unique indices took 0 sec ***" #> [1] "MH acceptance rate 0.628 took 0.0100000000000007 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 2 of 41) #> Required 1 iterations to reduce ESS from 500 to 320.74 #> [1] "reweighting took 0 sec for ESS 320.740377088985" #> [1] "MH acceptance rate 0.326 took 0 sec." #> previous ESS 320.74 (target: 288.666 for observation 1 of 1; wavenumber 3 of 41) #> Required 1 iterations to reduce ESS from 320.74 to 214.813 #> [1] "reweighting took 0 sec for ESS 214.812562402227" #> [1] "*** Resampling with 213 unique indices took 0.0199999999999996 sec ***" #> [1] "MH acceptance rate 0.51 took 0 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 4 of 41) #> Required 1 iterations to reduce ESS from 500 to 378.378 #> [1] "reweighting took 0 sec for ESS 378.378221512996" #> [1] "MH acceptance rate 0.262 took 0.00999999999999979 sec." #> previous ESS 378.378 (target: 340.54 for observation 1 of 1; wavenumber 5 of 41) #> Required 1 iterations to reduce ESS from 378.378 to 212.431 #> [1] "reweighting took 0 sec for ESS 212.43115567311" #> [1] "*** Resampling with 209 unique indices took 0 sec ***" #> [1] "MH acceptance rate 0.52 took 0.0200000000000005 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 6 of 41) #> Required 1 iterations to reduce ESS from 500 to 418.27 #> [1] "reweighting took 0 sec for ESS 418.270434952814" #> [1] "MH acceptance rate 0.266 took 0.00999999999999979 sec." #> previous ESS 418.27 (target: 376.443 for observation 1 of 1; wavenumber 7 of 41) #> Required 1 iterations to reduce ESS from 418.27 to 313.419 #> [1] "reweighting took 0 sec for ESS 313.419097710268" #> [1] "MH acceptance rate 0.27 took 0 sec." #> previous ESS 313.419 (target: 282.077 for observation 1 of 1; wavenumber 8 of 41) #> Required 1 iterations to reduce ESS from 313.419 to 234.918 #> [1] "reweighting took 0 sec for ESS 234.918384156081" #> [1] "*** Resampling with 221 unique indices took 0.0200000000000005 sec ***" #> [1] "MH acceptance rate 0.53 took 0.0199999999999996 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 9 of 41) #> Required 1 iterations to reduce ESS from 500 to 382.101 #> [1] "reweighting took 0 sec for ESS 382.10113787029" #> [1] "MH acceptance rate 0.252 took 0 sec." #> previous ESS 382.101 (target: 343.891 for observation 1 of 1; wavenumber 10 of 41) #> Required 2 iterations to reduce ESS from 382.101 to 340.752 #> [1] "reweighting took 0 sec for ESS 340.75175185415" #> [1] "MH acceptance rate 0.24 took 0.00999999999999979 sec." #> previous ESS 340.752 (target: 306.677 for observation 1 of 1; wavenumber 12 of 41) #> Required 1 iterations to reduce ESS from 340.752 to 257.469 #> [1] "reweighting took 0 sec for ESS 257.468567934346" #> [1] "MH acceptance rate 0.248 took 0.0199999999999996 sec." #> previous ESS 257.469 (target: 231.722 for observation 1 of 1; wavenumber 13 of 41) #> Required 1 iterations to reduce ESS from 257.469 to 163.637 #> [1] "reweighting took 0 sec for ESS 163.636924827613" #> [1] "*** Resampling with 196 unique indices took 0 sec ***" #> [1] "MH acceptance rate 0.38 took 0 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 14 of 41) #> Required 2 iterations to reduce ESS from 500 to 383.54 #> [1] "reweighting took 0 sec for ESS 383.539733218961" #> [1] "MH acceptance rate 0.236 took 0.00999999999999979 sec." #> previous ESS 383.54 (target: 345.186 for observation 1 of 1; wavenumber 16 of 41) #> Required 1 iterations to reduce ESS from 383.54 to 329.147 #> [1] "reweighting took 0 sec for ESS 329.146760390885" #> [1] "MH acceptance rate 0.246 took 0.0200000000000014 sec." #> previous ESS 329.147 (target: 296.232 for observation 1 of 1; wavenumber 17 of 41) #> Required 2 iterations to reduce ESS from 329.147 to 262.715 #> [1] "reweighting took 0 sec for ESS 262.71519124484" #> [1] "MH acceptance rate 0.266 took 0 sec." #> previous ESS 262.715 (target: 236.444 for observation 1 of 1; wavenumber 19 of 41) #> Required 1 iterations to reduce ESS from 262.715 to 172.819 #> [1] "reweighting took 0 sec for ESS 172.819040418002" #> [1] "*** Resampling with 182 unique indices took 0 sec ***" #> [1] "MH acceptance rate 0.546 took 0.00999999999999979 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 20 of 41) #> Required 2 iterations to reduce ESS from 500 to 445.054 #> [1] "reweighting took 0 sec for ESS 445.054009605662" #> [1] "MH acceptance rate 0.31 took 0.0199999999999996 sec." #> previous ESS 445.054 (target: 400.549 for observation 1 of 1; wavenumber 22 of 41) #> Required 2 iterations to reduce ESS from 445.054 to 343.64 #> [1] "reweighting took 0 sec for ESS 343.640127656283" #> [1] "MH acceptance rate 0.31 took 0.00999999999999979 sec." #> previous ESS 343.64 (target: 309.276 for observation 1 of 1; wavenumber 24 of 41) #> Required 2 iterations to reduce ESS from 343.64 to 305.199 #> [1] "reweighting took 0 sec for ESS 305.198702121203" #> [1] "MH acceptance rate 0.302 took 0.0199999999999996 sec." #> previous ESS 305.199 (target: 274.679 for observation 1 of 1; wavenumber 26 of 41) #> Required 1 iterations to reduce ESS from 305.199 to 259.667 #> [1] "reweighting took 0 sec for ESS 259.667238819587" #> [1] "MH acceptance rate 0.302 took 0 sec." #> previous ESS 259.667 (target: 233.701 for observation 1 of 1; wavenumber 27 of 41) #> Required 1 iterations to reduce ESS from 259.667 to 210.444 #> [1] "reweighting took 0.0200000000000014 sec for ESS 210.444067988262" #> [1] "*** Resampling with 205 unique indices took 0 sec ***" #> [1] "MH acceptance rate 0.484 took 0 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 28 of 41) #> Required 2 iterations to reduce ESS from 500 to 445.484 #> [1] "reweighting took 0 sec for ESS 445.48425386648" #> [1] "MH acceptance rate 0.3 took 0.00999999999999979 sec." #> previous ESS 445.484 (target: 400.936 for observation 1 of 1; wavenumber 30 of 41) #> Required 1 iterations to reduce ESS from 445.484 to 321.032 #> [1] "reweighting took 0 sec for ESS 321.031879468013" #> [1] "MH acceptance rate 0.336 took 0 sec." #> previous ESS 321.032 (target: 288.929 for observation 1 of 1; wavenumber 31 of 41) #> Required 2 iterations to reduce ESS from 321.032 to 287.012 #> [1] "reweighting took 0.0199999999999996 sec for ESS 287.012114827201" #> [1] "MH acceptance rate 0.324 took 0 sec." #> previous ESS 287.012 (target: 258.311 for observation 1 of 1; wavenumber 33 of 41) #> Required 2 iterations to reduce ESS from 287.012 to 175.12 #> [1] "reweighting took 0 sec for ESS 175.120181697452" #> [1] "*** Resampling with 166 unique indices took 0 sec ***" #> [1] "MH acceptance rate 0.57 took 0.00999999999999979 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 35 of 41) #> Required 3 iterations to reduce ESS from 500 to 391.914 #> [1] "reweighting took 0 sec for ESS 391.914312402872" #> [1] "MH acceptance rate 0.376 took 0.0199999999999996 sec." #> previous ESS 391.914 (target: 352.723 for observation 1 of 1; wavenumber 38 of 41) #> Required 1 iterations to reduce ESS from 391.914 to 349.783 #> [1] "reweighting took 0 sec for ESS 349.782811791971" #> [1] "MH acceptance rate 0.338 took 0 sec." #> previous ESS 349.783 (target: 314.805 for observation 1 of 1; wavenumber 39 of 41) #> Required 1 iterations to reduce ESS from 349.783 to 301.755 #> [1] "reweighting took 0 sec for ESS 301.755196114894" #> [1] "MH acceptance rate 0.41 took 0.0100000000000016 sec." #> previous ESS 301.755 (target: 271.58 for observation 1 of 1; wavenumber 40 of 41) #> Required 2 iterations to reduce ESS from 301.755 to 268.804 #> [1] "reweighting took 0 sec for ESS 268.803973327769" #> [1] "MH acceptance rate 0.444 took 0.0199999999999996 sec." #> [1] "SMC iter. 1 took 0.46 sec."