R/RcppExports.R

In serrsBayes: Bayesian Modelling of Raman Spectroscopy

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' Compute the spectral signature using Lorentzian peaks.
#' 
#' Calculates the value of the Lorentzian function at the given wavelengths,
#' given the parameters of the peaks. This function is thread-safe.
#' 
#' @param location Vector of location parameters of the peaks.
#' @param scale Vector of scale parameters of the peaks.
#' @param amplitude Vector of amplitudes of the peaks.
#' @param wavelengths Vector of wavenumbers at which to compute the function.
#' @return The value of the Lorentian function at the given wavelengths.
#' @examples
#'   Cal_V <- seq(300,400,by=5)
#'   loc <- c(320,350,375)
#'   sca <- c(10,5,18)
#'   amp <- c(1000,5000,2000)
#'   weightedLorentzian(loc,sca,amp,Cal_V)
weightedLorentzian <- function(location, scale, amplitude, wavelengths) {
    .Call('_serrsBayes_weightedLorentzian', PACKAGE = 'serrsBayes', location, scale, amplitude, wavelengths)
}

#' Compute the effective sample size (ESS) of the particles.
#' 
#' The ESS is a "rule of thumb" for assessing the degeneracy of
#' the importance distribution:
#' \deqn{ESS = \frac{(\sum_{q=1}^Q w_q)^2}{\sum_{q=1}^Q w_q^2}}
#' 
#' @param log_weights logarithms of the importance weights of each particle.
#' @return the effective sample size, a scalar between 0 and Q
#' @examples
#' x <- runif(100)
#' effectiveSampleSize(log(x))
#' @references
#' Liu, JS (2001) "Monte Carlo Strategies in Scientific Computing." Springer, NY, pp. 34--36.
effectiveSampleSize <- function(log_weights) {
    .Call('_serrsBayes_effectiveSampleSize', PACKAGE = 'serrsBayes', log_weights)
}

#' Update the importance weights of each particle.
#' 
#' @param spectra \code{n_y * nwl} Matrix of observed Raman spectra.
#' @param peaks \code{nwl * npart} Matrix containing the spectral signatures for each observation.
#' @param baselines \code{nwl * npart} Matrix containing the current values of the baselines.
#' @param i index of the current observation to use in calculating the likelihood
#' @param start index of the next wavelength to use in calculating the likelihood, permuted by \code{idx}
#' @param sigma Vector of \code{npart} standard deviations for each particle.
#' @param old_weights logarithms of the importance weights of each particle.
#' @param alpha the target learning rate for the reduction in effective sample size (ESS).
#' @param idx permutation of the indices of the wavelengths.
#' @return a List containing:
#' \describe{
#'   \item{\code{ess}}{The effective sample size, after reweighting.}
#'   \item{\code{weights}}{Vector of updated importance weights.}
#'   \item{\code{index}}{index of the last wavelength used.}
#'   \item{\code{evidence}}{SMC estimate of the logarithm of the model evidence.}
#' }
#' @references
#' Pitt, dos Santos Silva, Giordani & Kohn (2012)
#' "On some properties of Markov chain Monte Carlo simulation methods based on the particle filter"
#' J. Econometrics 171(2): 134--151,
#' DOI: \doi{10.1016/j.jeconom.2012.06.004}
#' 
#' Zhou, Johansen & Aston (2015) "Towards Automatic Model Comparison: An Adaptive Sequential Monte Carlo Approach"
#' \href{https://arxiv.org/abs/1303.3123}{arXiv:1303.3123} [stat.ME]
reWeightParticles <- function(spectra, peaks, baselines, i, start, sigma, old_weights, alpha, idx) {
    .Call('_serrsBayes_reWeightParticles', PACKAGE = 'serrsBayes', spectra, peaks, baselines, i, start, sigma, old_weights, alpha, idx)
}

#' Compute an ancestry vector for residual resampling of the SMC particles.
#' 
#' @param log_wt logarithms of the importance weights of each particle.
#' @return Vector of indices to the particles that will be propagated forward to the next generation (i.e. the parents)
#' @references
#' Liu & Chen (1998) "Sequential Monte Carlo methods for dynamic systems," JASA 93(443): 1032-1044,
#' \doi{10.1080/01621459.1998.10473765}
#' 
#' Douc, Cappe & Moulines (2005) "Comparison of resampling schemes for particle filtering"
#' In Proc. 4th IEEE Int. Symp. ISPA, pp. 64-69, \doi{10.1109/ISPA.2005.195385}
residualResampling <- function(log_wt) {
    .Call('_serrsBayes_residualResampling', PACKAGE = 'serrsBayes', log_wt)
}

#' Resample in place to avoid expensive copying of data structures, using a permutation
#' of the ancestry vector.
#' 
#' @param log_weights logarithms of the importance weights of each particle
#' @param ampMx \code{npeaks x npart} Matrix of amplitudes for each particle.
#' @param scaleMx \code{npeaks x npart} Matrix of scale parameters for each particle.
#' @param peaks \code{nwl x npart} Matrix containing the expectation of the Lorentzian mixture.
#' @param baselines \code{nwl x n_y x npart} Array of smoothing splines.
#' @param n_y number of observations
#' @param nwl number of wavenumbers
#' @return Vector of indices to the parents of the resampled particles.
#' @references
#' Murray, L.M., Lee, A. & Jacob, P.E. (2015) "Parallel resampling in the particle filter" \href{https://arxiv.org/abs/1301.4019}{arXiv:1301.4019v3}
#' @seealso \code{\link{residualResampling}}
resampleParticles <- function(log_weights, ampMx, scaleMx, peaks, baselines, n_y, nwl) {
    .Call('_serrsBayes_resampleParticles', PACKAGE = 'serrsBayes', log_weights, ampMx, scaleMx, peaks, baselines, n_y, nwl)
}

#' Compute the weighted arithmetic means of the particles.
#' 
#' This SMC estimate of the means can be used to centre independent Metropolis-Hastings proposals.
#' 
#' @param particles \code{npeaks * npart} Matrix of parameter values for each particle.
#' @param log_weights logarithms of the importance weights of each particle.
#' @return A vector of means, one for each row.
#' @seealso \code{\link[stats:weighted.mean]{weighted.mean}}
weightedMean <- function(particles, log_weights) {
    .Call('_serrsBayes_weightedMean', PACKAGE = 'serrsBayes', particles, log_weights)
}

#' Compute the weighted variance of the particles.
#' 
#' This SMC estimate of the variance can be used to scale the bandwidth of adaptive,
#' Gaussian random walk Metropolis-Hastings proposals.
#' 
#' @param particles \code{npeaks * npart} Matrix of parameter values for each particle.
#' @param log_weights logarithms of the importance weights of each particle.
#' @param mean Vector of weighted means of each particle.
#' @return A vector of variances, one for each row.
#' @seealso \code{\link[Hmisc:wtd.stats]{wtd.var}}
weightedVariance <- function(particles, log_weights, mean) {
    .Call('_serrsBayes_weightedVariance', PACKAGE = 'serrsBayes', particles, log_weights, mean)
}

#' Compute the spectral signature using Gaussian peaks.
#' 
#' Calculates the value of the squared exponential radial basis function at the given wavelengths,
#' given the parameters of the peaks. This function is thread-safe.
#' 
#' @param location Vector of location parameters of the peaks (mean).
#' @param scale Vector of scale parameters of the peaks (standard deviation).
#' @param amplitude Vector of amplitudes of the peaks.
#' @param wavelengths Vector of wavenumbers at which to compute the function.
#' @return The value of the Gaussian function at the given wavelengths.
#' @examples
#'   Cal_V <- seq(300,400,by=5)
#'   loc <- c(320,350,375)
#'   sca <- c(10,5,18)
#'   amp <- c(1000,5000,2000)
#'   weightedGaussian(loc,sca,amp,Cal_V)
weightedGaussian <- function(location, scale, amplitude, wavelengths) {
    .Call('_serrsBayes_weightedGaussian', PACKAGE = 'serrsBayes', location, scale, amplitude, wavelengths)
}

#' Update all of the parameters using a single Metropolis-Hastings step.
#' 
#' Updates all of the parameters using a single Metropolis-Hastings step, such that the
#' baseline cancels out in the MH ratio, using the marginalisation identity of Chib (1995).
#' If \code{npart > 1}, then multiple MCMC chains will be executed independently in parallel using OpenMP.
#' This means that all functions used for the proposal distributions and to evaluate the MH ratio
#' need to be thread-safe. Specifically, no calls to \code{R::rnorm}, \code{R::dnorm}, nor their
#' Rcpp equivalents, can be made from within the parallel portion of the code.
#' 
#' @param spectra \code{n_y * nwl} Matrix of observed Raman spectra.
#' @param n number of observations to use in calculating the likelihood
#' @param conc Vector of \code{n} nanomolar (nM) dye concentrations
#' @param wavelengths Vector of \code{nwl} wavenumbers at which the spetra are observed.
#' @param peakWL Vector of locations for each peak (cm^-1)
#' @param betaMx \code{npeaks * npart} Matrix of regression coefficients to update.
#' @param scaleMx \code{npeaks * npart} Matrix of scale parameters to update.
#' @param sigma Vector of \code{npart} standard deviations to update.
#' @param expMx \code{nwl * npart} Matrix of expectations of the Lorentzian or Gaussian function.
#' @param baselines \code{nKnots * n_y * npart} Array of smoothing splines.
#' @param sd_mh Vector of \code{2 * npeaks} bandwidths for the random walk proposals.
#' @param priors List of hyperparameters for the prior distributions.
#' @return The number of RWMH proposals that were accepted.
#' @references
#' Chib (1995) "Marginal Likelihood from the Gibbs Output," JASA 90(432): 1313--1321,
#' \doi{10.1080/01621459.1995.10476635}
#' 
#' Rosenthal (2000) "Parallel computing and Monte Carlo algorithms" Far East J. Theor. Stat. 4(2): 207--236,
#' URL: \href{https://www.pphmj.com/abstract/1961.htm}{https://www.pphmj.com/abstract/1961.htm}
marginalMetropolisUpdate <- function(spectra, n, conc, wavelengths, peakWL, betaMx, scaleMx, sigma, expMx, baselines, sd_mh, priors) {
    .Call('_serrsBayes_marginalMetropolisUpdate', PACKAGE = 'serrsBayes', spectra, n, conc, wavelengths, peakWL, betaMx, scaleMx, sigma, expMx, baselines, sd_mh, priors)
}

#' Sum log-likelihoods of Gaussian.
#' 
#' This is an internal function that is only exposed on the public API for unit testing purposes.
#' 
#' The sum of the log-likelihoods (log of the product of the likelihoods)
#' for independent, identically-distributed, Gaussian random variables.
#' Note: this Rcpp function is thread-safe, unlike the equivalent alternatives. 
#' 
#' @param x Vector of i.i.d. Gaussian random varibles
#' @param mean Vector of means
#' @param sd Vector of standard deviations
#' @return log-likelihood of x
#' @seealso \code{sum(dnorm(x, mean, sd, log=TRUE))}
#' @examples
#'   x <- rnorm(100)
#'   mu <- rep(0,length(x))
#'   sd <- rep(1,length(x))
#'   sumDnorm(x,mu,sd)
sumDnorm <- function(x, mean, sd) {
    .Call('_serrsBayes_sumDnorm', PACKAGE = 'serrsBayes', x, mean, sd)
}

#' Sum log-likelihoods of i.i.d. lognormal.
#' 
#' This is an internal function that is only exposed on the public API for unit testing purposes.
#' 
#' The sum of the log-likelihoods (log of the product of the likelihoods)
#' for independent, identically-distributed, lognormal random variables. 
#' Note: this Rcpp function is thread-safe, unlike the equivalent alternatives. 
#' 
#' @param x Vector of i.i.d. lognormal random varibles
#' @param meanlog mean of the distribution on the log scale
#' @param sdlog standard deviation on the log scale
#' @return log-likelihood of x
#' @seealso \code{sum(dlnorm(x, meanlog, sdlog, log=TRUE))}
#' @examples
#' x <- rlnorm(100)
#' sumDlogNorm(x,0,1)
sumDlogNorm <- function(x, meanlog, sdlog) {
    .Call('_serrsBayes_sumDlogNorm', PACKAGE = 'serrsBayes', x, meanlog, sdlog)
}

#' Sum log-likelihoods of i.i.d. exponential.
#' 
#' This is an internal function that is only exposed on the public API for unit testing purposes.
#' 
#' The sum of the log-likelihoods (log of the product of the likelihoods)
#' for independent, identically-distributed, exponential random variables. 
#' Note: this Rcpp function is thread-safe, unlike the equivalent alternatives. 
#' @param x Vector of i.i.d. exponential random varibles
#' @param rate parameter of the exponential distribution
#' @return log-likelihood of x
#' @seealso \code{sum(dexp(x, rate, log=TRUE))}
sumDexp <- function(x, rate) {
    .Call('_serrsBayes_sumDexp', PACKAGE = 'serrsBayes', x, rate)
}

#' Compute the spectral signature using Voigt peaks.
#' 
#' Calculates the value of the pseudo-Voigt broadening function at the given wavenumbers,
#' given the parameters of the peaks. This function is thread-safe.
#' 
#' @param location Vector of location parameters of the peaks (\eqn{cm^{-1}})
#' @param scale_G Vector of standard deviations \eqn{\sigma_j} of the Gaussian components.
#' @param scale_L Vector of scale parameters \eqn{\phi_j} of the Lorentzian components.
#' @param amplitude Vector of amplitudes of the peaks (a.u.)
#' @param wavenum Vector of wavenumbers at which to compute the function.
#' @return The value of the pseudo-Voigt function at the given wavenumbers.
#' @examples
#'   Cal_V <- seq(300,400,by=5)
#'   loc <- c(320,350,375)
#'   scG <- c(10,5,1)
#'   scL <- c(3,20,7)
#'   amp <- c(100,500,200)
#'   mixedVoigt(loc,scG,scL,amp,Cal_V)
#' @references
#' Thompson, Cox & Hastings (1987) "Rietveld refinement of Debye--Scherrer synchrotron X-ray data from \eqn{Al_2 O_3},"
#' J. Appl. Crystallogr. 20(2): 79--83, DOI: \doi{10.1107/S0021889887087090}
mixedVoigt <- function(location, scale_G, scale_L, amplitude, wavenum) {
    .Call('_serrsBayes_mixedVoigt', PACKAGE = 'serrsBayes', location, scale_G, scale_L, amplitude, wavenum)
}

#' Compute the pseudo-Voigt mixing ratio for each peak.
#' 
#' Calculates the mixing parameter \eqn{\eta_j} from the scales of the Gaussian/Lorentzian
#' components.
#' 
#' First, calculate a polynomial average of the scale parameters according to
#' the approximation of Thompson et al. (1987):
#' \deqn{f_{G,L} = (\sigma_j^5 + 2.69\sigma_j^4\phi_j + 2.42\sigma_j^3\phi_j^2 + 4.47\sigma_j^2\phi_j^3 + 0.07\sigma_j\phi_j^4 + \phi_j^5)^{1/5} }
#' 
#' Then the Voigt mixing parameter \eqn{\eta_j} is defined as:
#' \deqn{\eta_j = 1.36\frac{\phi_j}{f_{G,L}} - 0.47(\frac{\phi_j}{f_{G,L}})^2 + 0.11(\frac{\phi_j}{f_{G,L}})^3}
#' 
#' @param scale_G Vector of standard deviations \eqn{\sigma_j} of the Gaussian components.
#' @param scale_L Vector of scale parameters \eqn{\phi_j} of the Lorentzian components.
#' @return The Voigt mixing weights for each peak, between 0 (Gaussian) and 1 (Lorentzian).
#' @references
#' Thompson, Cox & Hastings (1987) "Rietveld refinement of Debye--Scherrer synchrotron X-ray data from \eqn{Al_2 O_3},"
#' J. Appl. Crystallogr. 20(2): 79--83, \doi{10.1107/S0021889887087090}
getVoigtParam <- function(scale_G, scale_L) {
    .Call('_serrsBayes_getVoigtParam', PACKAGE = 'serrsBayes', scale_G, scale_L)
}

#' Initialise the vector of Metropolis-Hastings proposals.
#' 
#' This is an internal function that is only exposed on the public API for unit testing purposes.
#' 
#' @param nPK number of Raman peaks in the spectral signature
#' @param T_Prop_Theta Vector of logarithms of the MH proposals
#' @return Vector of proposals
copyLogProposals <- function(nPK, T_Prop_Theta) {
    .Call('_serrsBayes_copyLogProposals', PACKAGE = 'serrsBayes', nPK, T_Prop_Theta)
}

#' Compute the log-likelihood.
#' 
#' This is an internal function that is only exposed on the public API for unit testing purposes.
#' It computes the log-likelihood of the spline and the noise, once the spectral signature has
#' been subtracted from the observed data. Thus, it can be used with either Lorentzian, Gaussian,
#' or pseudo-Voigt broadening functions.
#' 
#' @param obsi Vector of residuals after the spectral signature has been subtracted.
#' @param lambda smoothing parameter of the penalised B-spline.
#' @param prErrNu hyperparameter of the additive noise
#' @param prErrSS hyperparameter of the additive noise
#' @param basisMx Matrix of B-spline basis functions
#' @param eigVal eigenvalues of the Demmler-Reinsch factorisation
#' @param precMx precision matrix for the spline
#' @param xTx sparse matrix cross-product
#' @param aMx orthoganal matrix A from the Demmler-Reinsch factorisation
#' @param ruMx product of Ru from the Demmler-Reinsch factorisation
#' @return The logarithm of the likelihood.
computeLogLikelihood <- function(obsi, lambda, prErrNu, prErrSS, basisMx, eigVal, precMx, xTx, aMx, ruMx) {
    .Call('_serrsBayes_computeLogLikelihood', PACKAGE = 'serrsBayes', obsi, lambda, prErrNu, prErrSS, basisMx, eigVal, precMx, xTx, aMx, ruMx)
}

#' Update the parameters of the Voigt peaks using marginal Metropolis-Hastings.
#' 
#' Updates all of the parameters (location, amplitude, std. dev., and scale) using a single Metropolis-
#' Hastings step, such that the baseline cancels out in the MH ratio, using the marginalisation identity
#' of Chib (1995).
#' Note: if \code{npart > 1}, then multiple MCMC chains will be executed independently in parallel using
#' OpenMP. This means that all functions used for the proposal distributions and to evaluate the MH ratio
#' need to be thread-safe. Specifically, no calls to \code{R::rnorm}, \code{R::dnorm}, nor their
#' Rcpp equivalents, can be made from within the parallel portion of the code.
#' 
#' @param spectra \code{n_y * nwl} Matrix of observed Raman spectra.
#' @param n number of observations to use in calculating the likelihood.
#' @param kappa likelihood tempering parameter.
#' @param conc Vector of \code{n_y} nanomolar (nM) dye concentrations
#' @param wavenum Vector of \code{nwl} wavenumbers at which the spetra are observed.
#' @param thetaMx \code{(4+npeaks*4) x npart} Matrix of parameter values for each peak.
#' @param logThetaMx \code{(4+npeaks*4) x npart} Matrix of logarithms of the parameters.
#' @param mhChol lower-triangular Cholesky factorisation of the covariance matrix for the random walk proposals.
#' @param priors List of hyperparameters for the prior distributions.
#' @return The number of RWMH proposals that were accepted.
#' @references
#' Chib (1995) "Marginal Likelihood from the Gibbs Output," JASA 90(432): 1313--1321,
#' \doi{10.1080/01621459.1995.10476635}
#' 
#' Rosenthal (2000) "Parallel computing and Monte Carlo algorithms" Far East J. Theor. Stat. 4(2): 207--236,
#' URL: \href{https://www.pphmj.com/abstract/1961.htm}{https://www.pphmj.com/abstract/1961.htm}
mhUpdateVoigt <- function(spectra, n, kappa, conc, wavenum, thetaMx, logThetaMx, mhChol, priors) {
    .Call('_serrsBayes_mhUpdateVoigt', PACKAGE = 'serrsBayes', spectra, n, kappa, conc, wavenum, thetaMx, logThetaMx, mhChol, priors)
}