A demo using the R packages 'serrsBayes' and 'hyperSpec' to read and analyse a Raman spectrum of ethanol.

EthanolDemo.R

library(serrsBayes)

### read the data ----
library(hyperSpec)
etoh <- read.spc("EtOH.spc")
plot(etoh)

wv <- wl(etoh)
i1 <- which.min(wv < 300)
i2 <- which.min(wv < 1800)
abline(v=wv[i1], col=2, lty=2)
abline(v=wv[i2], col=2, lty=2)
wv <- wv[i1:i2]
nW <- length(wv)
y <- etoh$spc[i1:i2]
dim(y) <- c(1, nW)

peakLocations <- c(430, 880, 1055, 1090, 1280, 1460)
nPK <- length(peakLocations)
pkIdx <- numeric(nPK)
for (i in 1:nPK) {
  pkIdx[i] <- which.min(wv < peakLocations[i])
}
plot(wv, y[1,], type='l', col=4,
     xlab=expression(paste("Raman shift (cm"^{-1}, ")")), ylab="Intensity (a.u.)", main="Observed Raman spectrum for ethanol (EtOH)")
points(peakLocations, y[1,pkIdx], cex=2, col=2)
text(peakLocations + c(100,20,40,0), y[1,pkIdx] + c(0,700,400,700), labels=peakLocations)
rug(peakLocations, col=2)

### fit the model ----
set.seed(1234)
lPriors <- list(loc.mu=peakLocations, loc.sd=rep(50,nPK), scaG.mu=log(16.47) - (0.34^2)/2,
                scaG.sd=0.34, scaL.mu=log(25.27) - (0.4^2)/2, scaL.sd=0.4, noise.nu=5,
                noise.sd=50, bl.smooth=1, bl.knots=50)
tm <- Sys.time() # approx. 10 min for 3000 particles with max. 50 MCMC steps
res1 <- fitVoigtPeaksSMC(wv, y, lPriors, npart=3000, rate=0.499, mcSteps=50)
res1$time <- Sys.time() - tm
print(res1$time)
save(res1, file="etohVoigt.rda")

### results ----
print(paste(format(res1$time), "for",length(res1$ess),"SMC iterations."))
samp.idx <- sample.int(nrow(res1$location), 200, replace = FALSE)
plot(wv, y[1,], type='l', lwd=3,
     xlab=expression(paste("Raman shift (cm"^{-1}, ")")), ylab="Intensity (a.u.)", main="Fitted model with Voigt peaks")
samp.mat <- resid.mat <- matrix(0,nrow=length(samp.idx), ncol=nW)
samp.sigi <- samp.lambda <- numeric(length=nrow(samp.mat))
for (pt in 1:length(samp.idx)) {
  k <- samp.idx[pt]
  samp.mat[pt,] <- mixedVoigt(res1$location[k,], res1$scale_G[k,],
                              res1$scale_L[k,], res1$beta[k,], wv)
  samp.sigi[pt] <- res1$sigma[k]
  samp.lambda[pt] <- res1$lambda[k]
  
  Obsi <- y[1,] - samp.mat[pt,]
  g0_Cal <- length(Obsi) * samp.lambda[pt] * res1$priors$bl.precision
  gi_Cal <- crossprod(res1$priors$bl.basis) + g0_Cal
  mi_Cal <- as.vector(solve(gi_Cal, crossprod(res1$priors$bl.basis, Obsi)))
  
  bl.est <- res1$priors$bl.basis %*% mi_Cal # smoothed residuals = estimated basline
  lines(wv, bl.est, col="#C3000009")
  lines(wv, bl.est + samp.mat[pt,], col="#00C30009")
  resid.mat[pt,] <- Obsi - bl.est[,1]
}

for (j in 1:nPK) {
  hist(res1$location[,j], main=paste("Peak",j),
       xlab=expression(paste("Peak location (cm"^{-1}, ")")),
       freq=FALSE, xlim=range(res1$location[,j], peakLocations[j]))
  lines(density(res1$location[,j]), col=4, lty=2, lwd=3)
  abline(v=peakLocations[j], col=2, lty=3, lwd=3)
}
newLocs <- colMeans(res1$location)

for (j in 1:nPK) {
  hist(res1$beta[,j], main=paste("Peak",j),
       xlab="Peak amplitude (a.u.)", freq=FALSE)
  lines(density(res1$beta[,j]), col=4, lty=2, lwd=3)
}