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)

Arguments

wl

Vector of nwl wavenumbers at which the spetra are observed.

spc

n_y * nwl Matrix of observed Raman spectra.

peakWL

Vector of locations for each peak (cm^-1)

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

minESS

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

destDir

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

Value

a List containing weighted parameter values, known as particles:

weights

Vector of importance weights for each particle.

beta

npart * npeaks Matrix of regression coefficients for the amplitudes.

scale

npart * npeaks Matrix of scale parameters.

sigma

Vector of npart standard deviations.

alpha

bl.knots * n_y * npart Array of spline coefficients for the baseline.

basis

A dense nwl * bl.knots Matrix containing the values of the basis functions.

expFn

npart * nwl Matrix containing the spectral signature.

ess

Vector containing the effective sample size (ESS) at each SMC iteration.

logEvidence

Vector containing the logarithm of the model evidence (marginal likelihood).

accept

Vector containing the Metropolis-Hastings acceptance rate at each SMC iteration.

sd.mh

niter * 2npeaks Matrix of random walk MH bandwidths at each SMC iteration..

References

Chopin (2002) "A Sequential Particle Filter Method for Static Models," Biometrika 89(3): 539--551, DOI: 10.1093/biomet/89.3.539

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(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.0730000000000075 sec." #> [1] "Step 1: initialization for 2 peaks took 0.00899999999998613 sec." #> [1] "Observation 1 20190429-135642" #> [1] "Baselines took 0.0029999999999859 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.00100000000000477 sec for ESS 124.347094238812" #> [1] "*** Resampling with 141 unique indices took 0 sec ***" #> [1] "MH acceptance rate 0.606 took 0.0220000000000198 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 2 of 41) #> Required 1 iterations to reduce ESS from 500 to 299.811 #> [1] "reweighting took 0.000999999999976353 sec for ESS 299.810748073499" #> [1] "MH acceptance rate 0.318 took 0.0230000000000246 sec." #> previous ESS 299.811 (target: 269.83 for observation 1 of 1; wavenumber 3 of 41) #> Required 1 iterations to reduce ESS from 299.811 to 193.413 #> [1] "reweighting took 0 sec for ESS 193.413068812062" #> [1] "*** Resampling with 200 unique indices took 0.000999999999976353 sec ***" #> [1] "MH acceptance rate 0.514 took 0.0240000000000009 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 4 of 41) #> Required 1 iterations to reduce ESS from 500 to 389.127 #> [1] "reweighting took 0 sec for ESS 389.126886390313" #> [1] "MH acceptance rate 0.262 took 0.0229999999999961 sec." #> previous ESS 389.127 (target: 350.214 for observation 1 of 1; wavenumber 5 of 41) #> Required 1 iterations to reduce ESS from 389.127 to 192.126 #> [1] "reweighting took 0 sec for ESS 192.126190340229" #> [1] "*** Resampling with 199 unique indices took 0.00100000000000477 sec ***" #> [1] "MH acceptance rate 0.44 took 0.0260000000000105 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 6 of 41) #> Required 1 iterations to reduce ESS from 500 to 404.401 #> [1] "reweighting took 0 sec for ESS 404.401469132144" #> [1] "MH acceptance rate 0.234 took 0.0230000000000246 sec." #> previous ESS 404.401 (target: 363.961 for observation 1 of 1; wavenumber 7 of 41) #> Required 1 iterations to reduce ESS from 404.401 to 300.316 #> [1] "reweighting took 0 sec for ESS 300.315815893175" #> [1] "MH acceptance rate 0.248 took 0.0249999999999773 sec." #> previous ESS 300.316 (target: 270.284 for observation 1 of 1; wavenumber 8 of 41) #> Required 1 iterations to reduce ESS from 300.316 to 220.285 #> [1] "reweighting took 0 sec for ESS 220.28527017881" #> [1] "*** Resampling with 213 unique indices took 0.00100000000000477 sec ***" #> [1] "MH acceptance rate 0.462 took 0.0229999999999961 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 9 of 41) #> Required 1 iterations to reduce ESS from 500 to 375.815 #> [1] "reweighting took 0.00100000000000477 sec for ESS 375.815219990133" #> [1] "MH acceptance rate 0.228 took 0.0229999999999961 sec." #> previous ESS 375.815 (target: 338.234 for observation 1 of 1; wavenumber 10 of 41) #> Required 2 iterations to reduce ESS from 375.815 to 303.243 #> [1] "reweighting took 0.00100000000000477 sec for ESS 303.243410856979" #> [1] "MH acceptance rate 0.282 took 0.0219999999999914 sec." #> previous ESS 303.243 (target: 272.919 for observation 1 of 1; wavenumber 12 of 41) #> Required 1 iterations to reduce ESS from 303.243 to 213.711 #> [1] "reweighting took 0.00100000000000477 sec for ESS 213.711277025406" #> [1] "*** Resampling with 207 unique indices took 0.00100000000000477 sec ***" #> [1] "MH acceptance rate 0.364 took 0.0270000000000152 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 13 of 41) #> Required 1 iterations to reduce ESS from 500 to 319.725 #> [1] "reweighting took 0.000999999999976353 sec for ESS 319.72533422697" #> [1] "MH acceptance rate 0.248 took 0.0240000000000009 sec." #> previous ESS 319.725 (target: 287.753 for observation 1 of 1; wavenumber 14 of 41) #> Required 1 iterations to reduce ESS from 319.725 to 262.825 #> [1] "reweighting took 0.00100000000000477 sec for ESS 262.825070707226" #> [1] "MH acceptance rate 0.214 took 0.0270000000000152 sec." #> previous ESS 262.825 (target: 236.543 for observation 1 of 1; wavenumber 15 of 41) #> Required 1 iterations to reduce ESS from 262.825 to 202.944 #> [1] "reweighting took 0 sec for ESS 202.943859981634" #> [1] "*** Resampling with 206 unique indices took 0.000999999999976353 sec ***" #> [1] "MH acceptance rate 0.386 took 0.0260000000000105 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 16 of 41) #> Required 1 iterations to reduce ESS from 500 to 446.262 #> [1] "reweighting took 0.00100000000000477 sec for ESS 446.261994294307" #> [1] "MH acceptance rate 0.222 took 0.0229999999999961 sec." #> previous ESS 446.262 (target: 401.636 for observation 1 of 1; wavenumber 17 of 41) #> Required 2 iterations to reduce ESS from 446.262 to 381.706 #> [1] "reweighting took 0.00100000000000477 sec for ESS 381.706413887462" #> [1] "MH acceptance rate 0.236 took 0.0229999999999961 sec." #> previous ESS 381.706 (target: 343.536 for observation 1 of 1; wavenumber 19 of 41) #> Required 1 iterations to reduce ESS from 381.706 to 288.501 #> [1] "reweighting took 0.00100000000000477 sec for ESS 288.500838557999" #> [1] "MH acceptance rate 0.24 took 0.0229999999999961 sec." #> previous ESS 288.501 (target: 259.651 for observation 1 of 1; wavenumber 20 of 41) #> Required 1 iterations to reduce ESS from 288.501 to 257.735 #> [1] "reweighting took 0 sec for ESS 257.735445628182" #> [1] "MH acceptance rate 0.262 took 0.0279999999999916 sec." #> previous ESS 257.735 (target: 231.962 for observation 1 of 1; wavenumber 21 of 41) #> Required 2 iterations to reduce ESS from 257.735 to 207.139 #> [1] "reweighting took 0.00100000000000477 sec for ESS 207.139355276311" #> [1] "*** Resampling with 209 unique indices took 0.00100000000000477 sec ***" #> [1] "MH acceptance rate 0.486 took 0.0229999999999961 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 23 of 41) #> Required 2 iterations to reduce ESS from 500 to 435.168 #> [1] "reweighting took 0.00100000000000477 sec for ESS 435.168498401019" #> [1] "MH acceptance rate 0.332 took 0.0300000000000011 sec." #> previous ESS 435.168 (target: 391.652 for observation 1 of 1; wavenumber 25 of 41) #> Required 1 iterations to reduce ESS from 435.168 to 388.94 #> [1] "reweighting took 0.00100000000000477 sec for ESS 388.939890819422" #> [1] "MH acceptance rate 0.284 took 0.0279999999999916 sec." #> previous ESS 388.94 (target: 350.046 for observation 1 of 1; wavenumber 26 of 41) #> Required 1 iterations to reduce ESS from 388.94 to 307.664 #> [1] "reweighting took 0.00100000000000477 sec for ESS 307.663501224414" #> [1] "MH acceptance rate 0.338 took 0.0309999999999775 sec." #> previous ESS 307.664 (target: 276.897 for observation 1 of 1; wavenumber 27 of 41) #> Required 1 iterations to reduce ESS from 307.664 to 239.741 #> [1] "reweighting took 0.00100000000000477 sec for ESS 239.74111807824" #> [1] "*** Resampling with 219 unique indices took 0.00100000000000477 sec ***" #> [1] "MH acceptance rate 0.452 took 0.0320000000000107 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 28 of 41) #> Required 2 iterations to reduce ESS from 500 to 440.446 #> [1] "reweighting took 0 sec for ESS 440.445587004102" #> [1] "MH acceptance rate 0.332 took 0.032999999999987 sec." #> previous ESS 440.446 (target: 396.401 for observation 1 of 1; wavenumber 30 of 41) #> Required 1 iterations to reduce ESS from 440.446 to 310.192 #> [1] "reweighting took 0 sec for ESS 310.191801845506" #> [1] "MH acceptance rate 0.302 took 0.0330000000000155 sec." #> previous ESS 310.192 (target: 279.173 for observation 1 of 1; wavenumber 31 of 41) #> Required 2 iterations to reduce ESS from 310.192 to 266.296 #> [1] "reweighting took 0.000999999999976353 sec for ESS 266.296127723442" #> [1] "MH acceptance rate 0.332 took 0.0330000000000155 sec." #> previous ESS 266.296 (target: 239.667 for observation 1 of 1; wavenumber 33 of 41) #> Required 1 iterations to reduce ESS from 266.296 to 238.453 #> [1] "reweighting took 0.00100000000000477 sec for ESS 238.453385931048" #> [1] "*** Resampling with 219 unique indices took 0 sec ***" #> [1] "MH acceptance rate 0.582 took 0.039999999999992 sec." #> previous ESS 500 (target: 450 for observation 1 of 1; wavenumber 34 of 41) #> Required 1 iterations to reduce ESS from 500 to 384.639 #> [1] "reweighting took 0.00100000000000477 sec for ESS 384.638594616592" #> [1] "MH acceptance rate 0.428 took 0.0269999999999868 sec." #> previous ESS 384.639 (target: 346.175 for observation 1 of 1; wavenumber 35 of 41) #> Required 3 iterations to reduce ESS from 384.639 to 304.536 #> [1] "reweighting took 0 sec for ESS 304.53640437768" #> [1] "MH acceptance rate 0.45 took 0.0300000000000011 sec." #> previous ESS 304.536 (target: 274.083 for observation 1 of 1; wavenumber 38 of 41) #> Required 4 iterations to reduce ESS from 304.536 to 244.307 #> [1] "reweighting took 0 sec for ESS 244.306593556333" #> [1] "MH acceptance rate 0.492 took 0.0320000000000107 sec." #> [1] "SMC iter. 1 took 0.980999999999995 sec."