-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathREADME.Rmd
119 lines (98 loc) · 5.06 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
---
output: github_document
---
# serrsBayes
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, echo = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
warning = FALSE,
message = FALSE,
comment = "#>",
fig.path = "inst/image/README-"
)
```
[](https://cran.r-project.org/package=serrsBayes)
[](https://zenodo.org/badge/latestdoi/121410558)
[](/~https://github.com/r-hub/cranlogs.app)
```{r hexSticker, echo=FALSE, message=FALSE, warning=FALSE, fig.align='right', fig.asp=1, out.width='100px',out.height='100px'}
tmp <- tempfile()
download.file("https://pure.strath.ac.uk/portal/files/43595106/Figure_2.zip", tmp)
tmp2 <- unzip(tmp, "Figure 2/T20 SERS spectra/T20_1_ REP1 Well_A1.SPC")
library(hyperSpec)
spcT20 <- read.spc (tmp2)
spectra <- spcT20[1,,600~1800]
wavenumbers <- wl(spectra)
nWL <- length(wavenumbers)
library(serrsBayes)
data("result", package = "serrsBayes")
samp.idx <- sample.int(length(result$weights), 50, prob=result$weights)
samp.mat <- resid.mat <- matrix(0,nrow=length(samp.idx), ncol=nWL)
samp.sigi <- samp.lambda <- numeric(length=nrow(samp.mat))
spectra <- as.matrix(spectra)
library(hexSticker)
sticker(expression({plot(wavenumbers, spectra[1,], type='l', ann=FALSE, xaxt='n', yaxt='n', bty='n', lwd=2, col="#000000C0", xlim=c(1100,1800))
for (pt in 1:length(samp.idx)) {
k <- samp.idx[pt]
samp.mat[pt,] <- mixedVoigt(result$location[k,], result$scale_G[k,],
result$scale_L[k,], result$beta[k,], wavenumbers)
samp.sigi[pt] <- result$sigma
samp.lambda[pt] <- result$lambda
Obsi <- spectra[1,] - samp.mat[pt,]
g0_Cal <- length(Obsi) * samp.lambda[pt] * result$priors$bl.precision
gi_Cal <- crossprod(result$priors$bl.basis) + g0_Cal
mi_Cal <- as.vector(solve(gi_Cal, crossprod(result$priors$bl.basis, Obsi)))
bl.est <- result$priors$bl.basis %*% mi_Cal # smoothed residuals = estimated basline
lines(wavenumbers, bl.est, col="#0000FF40")
lines(wavenumbers, bl.est + samp.mat[pt,], col="#FF00002F", lwd=0.5)
resid.mat[pt,] <- Obsi - bl.est[,1]
}}),package="serrsBayes",p_size=8, s_x=0.75, s_y=.75, s_width=2.7, s_height=2.5, p_color = 2, h_color = "#0000FFC0", h_fill ="#FFFFFFFF", filename = "inst/image/README-logo.png", p_y=0.7,
url="https://mooresm.github.io/serrsBayes/", u_x=0.92, u_y=0.08)
knitr::include_graphics("inst/image/README-logo.png", dpi=600)
```
`serrsBayes` provides model-based quantification of surface-enhanced resonance Raman spectroscopy (SERRS) using sequential Monte Carlo (SMC) algorithms. The details of the Bayesian model and informative priors are provided in the arXiv preprint, Moores et al. (2016; v2 2018) "[Bayesian modelling and quantification of Raman spectroscopy.](https://arxiv.org/abs/1604.07299)" Development of this software was supported by the UK
Engineering & Physical Sciences Research Council (EPSRC) programme grant "[In Situ Nanoparticle Assemblies for Healthcare Diagnostics and Therapy](https://gtr.ukri.org/projects?ref=EP%2FL014165%2F1)" (ref: EP/L014165/1).
# Installation Instructions
Stable releases, including binary packages for Windows & Mac OS, are available from CRAN:
* <https://CRAN.R-project.org/package=serrsBayes>
```{r eval=FALSE}
install.packages("serrsBayes")
```
The current development version can be installed from GitHub:
```{r eval=FALSE}
devtools::install_github("mooresm/serrsBayes")
```
# Example Usage
To simulate a synthetic Raman spectrum with known parameters:
```{r example}
set.seed(1234)
library(serrsBayes)
wavenumbers <- seq(700,1400,by=2)
spectra <- matrix(nrow=1, ncol=length(wavenumbers))
peakLocations <- c(840, 960, 1140, 1220, 1290)
peakAmplitude <- c(11500, 2500, 4000, 3000, 2500)
peakScale <- c(10, 15, 20, 10, 12)
signature <- weightedLorentzian(peakLocations, peakScale, peakAmplitude, wavenumbers)
baseline <- 1000*cos(wavenumbers/200) + 2*wavenumbers
spectra[1,] <- signature + baseline + rnorm(length(wavenumbers),0,200)
plot(wavenumbers, spectra[1,], type='l', xlab=expression(paste("Raman shift (cm"^{-1}, ")")), ylab="Intensity (a.u.)")
lines(wavenumbers, baseline, col=2, lty=4)
lines(wavenumbers, baseline + signature, col=4, lty=2, lwd=2)
```
Fit the model using SMC:
```{r SMC, results='hide'}
lPriors <- list(scale.mu=log(11.6) - (0.4^2)/2, scale.sd=0.4, bl.smooth=10^11, bl.knots=50,
beta.mu=5000, beta.sd=5000, noise.sd=200, noise.nu=4)
tm <- system.time(result <- fitSpectraSMC(wavenumbers, spectra, peakLocations, lPriors))
```
Sample 200 particles from the posterior distribution:
```{r plotting}
print(tm)
samp.idx <- sample.int(length(result$weights), 200, prob=result$weights)
plot(wavenumbers, spectra[1,], type='l', xlab=expression(paste("Raman shift (cm"^{-1}, ")")), ylab="Intensity (a.u.)")
for (pt in samp.idx) {
bl.est <- result$basis %*% result$alpha[,1,pt]
lines(wavenumbers, bl.est, col="#C3000009")
lines(wavenumbers, bl.est + result$expFn[pt,], col="#0000C309")
}
```