Log-Gaussian Gamma Processes for Training Bayesian Neural Networks in
  Raman and CARS Spectroscopies
We propose an approach utilizing gamma-distributed random variables, coupled
with log-Gaussian modeling, to generate synthetic datasets suitable for
training neural networks. This addresses the challenge of limited real
observations in various applications. We apply this methodology to both Raman
and coherent anti-Stokes Raman scattering (CARS) spectra, using experimental
spectra to estimate gamma process parameters. Parameter estimation is performed
using Markov chain Monte Carlo methods, yielding a full Bayesian posterior
distribution for the model which can be sampled for synthetic data generation.
Additionally, we model the additive and multiplicative background functions for
Raman and CARS with Gaussian processes. We train two Bayesian neural networks
to estimate parameters of the gamma process which can then be used to estimate
the underlying Raman spectrum and simultaneously provide uncertainty through
the estimation of parameters of a probability distribution. We apply the
trained Bayesian neural networks to experimental Raman spectra of
phthalocyanine blue, aniline black, naphthol red, and red 264 pigments and also
to experimental CARS spectra of adenosine phosphate, fructose, glucose, and
sucrose. The results agree with deterministic point estimates for the
underlying Raman and CARS spectral signatures.
Teemu Härkönen, Erik M. Vartiainen, Lasse Lensu, Matthew T. Moores, Lassi Roininen
# Log-Gaussian Gamma Processes for Training Bayesian Neural Networks in Raman and CARS Spectroscopies

###### Abstract

We propose an approach utilizing gamma-distributed random variables, coupled with log-Gaussian modeling, to generate synthetic datasets suitable for training neural networks. This addresses the challenge of limited real observations in various applications. We apply this methodology to both Raman and coherent anti-Stokes Raman scattering (CARS) spectra, using experimental spectra to estimate gamma process parameters. Parameter estimation is performed using Markov chain Monte Carlo methods, yielding a full Bayesian posterior distribution for the model which can be sampled for synthetic data generation. Additionally, we model the additive and multiplicative background functions for Raman and CARS with Gaussian processes. We train two Bayesian neural networks to estimate parameters of the gamma process which can then be used to estimate the underlying Raman spectrum and simultaneously provide uncertainty through the estimation of parameters of a probability distribution. We apply the trained Bayesian neural networks to experimental Raman spectra of phthalocyanine blue, aniline black, naphthalene red, and red 264 pigments and also to experimental CARS spectra of adenosine phosphate, fructose, glucose, and sucrose. The results agree with deterministic point estimates for the underlying Raman and CARS spectral signatures.

## 1 Introduction

Raman and coherent anti-Stokes Raman scattering (CARS) spectroscopies are vital tools used in chemistry, physics, and biomedical research [1, 2, 3]. The insights they offer into molecular vibrations, structural dynamics, and chemical compositions are invaluable. However, working with their data presents challenges. Measurement artifacts including noise, and, especially, background signals in Raman and CARS spectra often obscure crucial molecular information. Traditional methods for data correction are typically manual and may fall short in capturing the full complexity of the data. For instance, standard approaches used for removing the background signals include asymmetric least squares polynomial fitting, wavelet-based methods, optimization with Tikhonov regularization, and Kramers-Kronig relations [4, 5, 6, 7, 8, 9, 10, 11, 12]. While appealing, these methods suffer from practical drawbacks such as the need for manual tuning of the model or regularization parameters. The need for automated, robust, and statistically sound solutions to enhance our spectroscopic analyses is evident.

Deep neural networks offer a compelling solution for automatic spectral correction across various applications, from weather predictions [13, 14, 15] to medical imaging [16, 17, 18] and many others [19, 20, 21, 22, 23]. In the realm of Raman spectroscopy, deep neural networks have been used in chemical species identification and background removal [24, 25, 26, 27, 28]. Similarly, they have been applied to extract the underlying Raman spectra from CARS measurement [29, 30, 31, 32, 33, 34, 35]. Despite their efficacy, non-Bayesian neural networks lack a critical feature: the ability to quantify uncertainty in Raman spectrum estimation. Bayesian inference, on the other hand, provides an avenue to solve this problem.

Bayesian inference treats the parameters of a given model as random variables. These models consist of a likelihood function that is combined with prior distributions for the parameters to produce posterior estimates. The likelihood function is analogous to a utility function in an optimization context. It quantifies how well the model fits the observed data. The aforementioned prior distributions for the model parameters represent the information known beforehand, including any constraints dictated by the physical nature of the parameters, such as non-negativity. In spectroscopic analysis, the model parameters can be, for example, amplitudes, locations, and widths of Gaussian, Lorentzian, or Voigt line shape functions. The combination of the likelihood and the priors results in a posterior distribution over the model parameters. The posterior is a probabilistic representation of the uncertainty in the parameter estimates. Bayesian approaches have been considered for estimating spectrum parameters, where the authors used sequential Monte Carlo algorithms to numerically sample from the posterior distribution [36, 37]. While the uncertainty quantification provided by Bayesian modeling and Markov chain Monte Carlo (MCMC) methods is compelling, the approach
is known to be computationally expensive, see for example [38]. This becomes a major issue particularly with hyperspectral data sets, where an image can contain millions of individual spectra which are to be analyzed.

Bayesian neural networks are a synthesis of the aforementioned two ideas. Bayesian neural networks model the weights and biases of standard neural networks as random variables, which can be assigned prior distributions. When combined with a likelihood according to Bayes' theorem, the resulting utility function corresponds to the posterior for the neural network parameters. Advantages of this Bayesian neural network approach in comparison to non-Bayesian neural networks include robustness in terms of overfitting, providing uncertainty estimates instead of only point estimation, sequential learning, and better generalization [39]. In particular, uncertainty quantification has seen widespread research covering many application areas and topics, for example [40].

One of the challenges of Bayesian neural networks is that they typically contain an enormous number of parameters. For instance, our network comprises over 11 million parameters, far beyond what is commonly considered high-dimensional for MCMC [41, 42]. Some neural networks, such as large language models (LLMs), can have billions of parameters [43]. Thus, it can be challenging to establish convergence of such a large number of parameters in a statistically rigorous manner. To combat this, partially-Bayesian neural networks have been used as a practical tool to provide uncertainty estimation with neural networks. In addition to empirical validation through practice, studies have provided compelling analytical and numerical evidence that partially-Bayesian neural networks are indeed capable of providing posterior estimates on par or even superior performance to fully-Bayesian neural networks [44]. The above points lead us to construct our neural network for this study as a partially-Bayesian neural network.

Neural networks typically require large volumes of training data. This has been noted to be a problem also in spectroscopic applications as it is difficult to acquire large sets of independent data sets [28]. Therefore, many studies mentioned above use synthetic data to train the neural networks. The synthetic data is usually generated using random linear combinations of Lorentzian line shapes, where the amplitudes, locations, and widths are sampled from predefined probability distributions, see for example [25, 29, 45]. The background data is generated similarly. The backgrounds are modeled explicitly using a parametric functional form, such as a polynomial or a sigmoidal function, and the parameters of the model are again sampled from a predefined probability distribution [25, 24, 46]. An extension to this is to use experimental Raman spectra on top of the randomly generated spectra [34].

Stochastic processes can be used to draw samples of random functions. A typical example of a stochastic process is the widely-used Gaussian process (GP). Properties of the drawn samples such as differentiability are governed through kernel functions, which are used to model dependencies between data points. For readers unfamiliar with GPs, we recommend the book by Rasmussen and Williams [47]. Instead of using explicit, parametric functions to model the spectroscopic features, we propose using stochastic processes as a more flexible tool for the purpose. In this study, we use GPs as a generative model for the additive and multiplicative backgrounds of Raman and CARS spectra, see Fig. 1.

For the purpose of generating synthetic Raman spectral signatures, we propose a specific type of doubly-stochastic Levy process which we call a log-Gaussian gamma process. Our construction of the log-Gaussian gamma process is inspired by log-Gaussian Cox process which the authors have previously used as a model for spectra [48]. While it makes sense to model spectra as a Cox process where the relaxation from higher energy levels happens at a constant rate and results in counts of photons, the data is often available in scaled floating-point numbers which prevents direct application of the log-Gaussian Cox process model. Gamma-distributed variables have direct connections to Poisson-distributed variables, which constitute the Cox process, making the extension to a log-Gaussian gamma process intuitive as a model for Raman spectroscopy. The log-Gaussian gamma process can be used to generate arbitrary amounts of synthetic spectra once parameters of the stochastic process have been estimated. We perform the estimation using MCMC methods which allow us to construct a Bayesian posterior distribution for the model parameters, thereby including uncertainty of the parameter estimates in our data generation. This also applies to our GP-based background model. We present a high-level diagram of our stochastic process method for data

Fig. 1: Structure of our generative spectrum model using GPs and log-Gaussian gamma processes. On top, an experimental CARS spectrum of adenosine phosphate in blue and an example multiplicative background in red. We model the backgrounds as a GP. At the bottom, an example underlying Raman spectral signature in blue. We assume the Raman peaks to be distributed according to our proposed log-Gaussian gamma process model. The stochastic processes are parameterized according to \(\mu_{e}\), \(\mathbf{\theta}_{e}\), \(\alpha\), and \(\beta(\nu)\). We further model \(\beta(\nu)\) using GPs which are parameterized according to \(\mu_{\beta}\) and \(\mathbf{\theta}_{\beta}\). We construct statistical samples with MCMC for the model parameters which allow us to generate synthetic spectra for training our Bayesian neural network.
generation in Fig. 1. Fig. 2 shows an example of the aim of this paper, a Raman spectral signature extracted from a CARS spectrum using a Bayesian neural network. We provide a pseudo-code description of our approach in Algorithm 1.

The key contributions of this paper are the following. We propose using log-Gaussian gamma processes for modeling Raman spectral signatures and GPs to model additive or multiplicative background signals. The aforementioned doubly-stochastic processes are sampled randomly, enabling us to generate an arbitrary number of synthetic spectra that are statistically similar to experimental spectra. Finally, we present a partially-Bayesian neural network for analyzing Raman and CARS spectral measurements, which we train using the sampled synthetic spectra. Once trained, we use these neural networks to estimate the spectral signatures for experimental Raman spectroscopy measurements of phthalocyanine blue, naphthol red, aniline black, and red 264 pigments and for experimental CARS spectra of adenosine phosphate, fructose, glucose, and sucrose in addition to synthetic test spectra.

```
Step 1: Fit a log-Gaussian gamma process to Raman spectrum data. Step 2: Fit a GP to background data. Step 3: Draw a large number of realizations from the fitted log-Gaussian gamma process. Step 4: Draw a large number of realizations from the fitted GP. Step 5: Use a forward model to combine the realizations to form a data set of synthetic spectra. Step 6. Train a Bayesian neural network using the data set of synthetic spectra.
```

**Algorithm 1** Log-Gaussian gamma process data generation for training Bayesian neural networks

The rest of the paper is structured as follows. We detail the steps used to generate the synthetic training data in three stages in the following three sections. We first present the log-Gaussian gamma process as a model for Raman spectral signatures and explain how to draw realizations of this doubly-stochastic process. This is followed by a description of our GP-based additive and multiplicative background models. We finalize the explanation of our synthetic data generation method with definitions of the forward models used to simulate synthetic training data for Raman and CARS measurements with additive and multiplicative backgrounds, respectively. Next, we present our partially-Bayesian neural network architecture, which we train against the synthetic data sets that we have generated. We document computational details and prior distributions in the next section, followed by a presentation of our results for both artificial and real experimental data. Finally, we conclude with a discussion of the significance and other potential applications for our method.

## 2 Log-Gaussian gamma process spectrum model

We model a Raman spectral signature as a collection of conditionally-independent, gamma-distributed random variables

\[r_{k}:=r(\nu_{k})\sim\text{Gamma}\left(\alpha,\beta(\nu_{k})\right), \tag{1}\]

where \(r_{k}\) denotes a Raman measurement at wavenumber location \(\nu_{k}\) with \(\alpha\) and \(\beta(\nu_{k})\) being the shape and scale parameters of the gamma distribution, respectively. The above construction is motivated by log-Gaussian Cox processes [49] but without the restriction of modeling of only integer-valued data and with an additional parameter in the stochastic process allowing for more flexible modeling of uncertainty. Poisson-distributed random variables, which constitute the Cox process, have a single parameter to control both the mean and variance of the distribution. Very often in real data, this assumption is found to be too restrictive, leading to a model that is either under- or over-dispersed [50]. In contrast, the gamma distribution has two parameters which together allow for a range of different variances for a given mean.

We extend Eq. (1) by modeling the log-scale as a GP, resulting in a hierarchical model

\[\log\beta(\mathbf{v})\sim\text{GP}(\mu_{\beta},\Sigma_{\beta}(\mathbf{v},\mathbf{v},\mathbf{ \theta}_{\beta})), \tag{2}\]

where \(\mathbf{v}:=(\nu_{1},\dots,\nu_{K})^{T}\) is a vector of the wavenumber locations with \(\mu_{\beta}\in\mathbb{R}\) and \(\Sigma_{\beta}(\mathbf{v},\mathbf{v},\mathbf{\theta}_{\beta})\in\mathbb{R}^{K\times K}\) being a constant mean and a covariance matrix parameterized according to hyperparameters \(\mathbf{\theta}_{\beta}\). This doubly-stochastic model introduces dependence between values \(r_{i}\) and \(r_{j}\) at different wavenumbers \(\nu_{i}\) and \(\nu_{j}\). For the covariance function of the log-scale GP, we use the squared exponential kernel

\[\left[\Sigma_{\beta}(\mathbf{v},\mathbf{v},\mathbf{\theta}_{\beta})\right]_{ij}=\sigma_{ \beta,f}^{2}\exp\left(-\frac{1}{2}\frac{\left(\nu_{i}-\nu_{j}\right)^{2}}{l_{ \beta}^{2}}\right)+\sigma_{\beta}^{2}\delta(\nu_{i}-\nu_{j}), \tag{3}\]

where \(\left[\Sigma_{\beta}(\mathbf{v},\mathbf{v},\mathbf{\theta}_{\beta})\right]_{ij}\) denotes the \(ij\)th element of the covariance matrix \(\Sigma_{\beta}(\mathbf{v},\mathbf{v},\mathbf{\theta}_{\beta})\), \(\sigma_{\beta,f}^{2}\) is the signal variance, \(l_{\beta}\) is the length scale, \(\sigma_{\beta}^{2}\) denotes the noise variance, and \(\delta(\nu_{i}-\nu_{j})\) is the Dirac delta function with \(\mathbf{\theta}_{\beta}=\left(\sigma_{\beta,f}^{2},l_{\beta},\sigma_{\beta}^{2} \right)^{T}\). The GP construction yields an analytical form for the log-scale \(\log\beta(\mathbf{v})\) which we will

Fig. 2: On top, an experimental CARS spectrum of adenosine phosphate in blue. With a trained Bayesian neural network, we can extract the underlying Raman spectral signature from the data along with an uncertainty estimate for the spectrum. At the bottom, the corresponding Bayesian neural network median Raman spectrum estimate and 90% confidence interval of the estimate for the adenosine phosphate data.
detail below as we construct the posterior distribution according to Bayes' theorem. This log-GP parameterization is identical to the log-intensity model for Poisson variables that features in log-Gaussian Cox processes. For more details on the log-Gaussian Cox process, see [49] and for example [51].

The posterior distribution involves the likelihood function \(\mathcal{L}(\mathbf{r}\mid\alpha,\beta(\mathbf{\nu}))\), a log-GP prior for the scale \(\pi_{0}(\beta(\mathbf{\nu})\mid\mu_{\beta},\mathbf{\theta}_{\beta})\), and a joint prior distribution \(\pi_{0}(\alpha,\mu_{\beta},\mathbf{\theta}_{\beta})\) for rest of the model parameters. Given a measured Raman spectrum \(\mathbf{r}:=(r(\nu_{1}),\ldots,r(\nu_{K}))^{\intercal}\), we can formulate the likelihood as a product of conditionally-independent, gamma-distributed random variables

\[\mathcal{L}(\mathbf{r}\mid\alpha,\beta(\mathbf{\nu}))\propto\prod_{k=1}^{K}\frac{r_{k} ^{\alpha-1}\exp(-r_{k}/\beta_{k})}{\Gamma(\alpha)\beta_{k}^{\alpha}}, \tag{4}\]

where \(\beta_{k}:=\beta(\nu_{k})\), and \(\Gamma(\alpha)\) is the gamma function. The hierarchical prior for \(\beta(\mathbf{\nu})\) can be evaluated as

\[\pi_{0}(\beta(\mathbf{\nu})\mid\mu_{\beta},\mathbf{\theta}_{\beta})=\frac {1}{\sqrt{(2\pi)^{k}}}\left|\Sigma_{\beta}(\mathbf{\nu},\mathbf{\nu};\mathbf{\theta}_{\beta })\right|^{-1/2} \tag{5}\] \[\times\exp\left(-\frac{1}{2}\left(\beta(\mathbf{\nu})-\mu_{\beta} \right)^{\intercal}\Sigma_{\beta}(\mathbf{\nu},\mathbf{\nu};\mathbf{\theta}_{\beta})^{-1} \left(\beta(\mathbf{\nu})-\mu_{\beta}\right)\right),\]

where \(\left|\Sigma_{\beta}(\mathbf{\nu},\mathbf{\nu};\mathbf{\theta}_{\beta})\right|\) denotes the determinant of the covariance matrix. With the above and a joint prior \(\pi_{0}(\alpha,\mu_{\beta},\mathbf{\theta}_{\beta})\), we can construct the posterior distribution for the model parameters conditioned on the measured spectrum data \(\mathbf{r}\) as

\[\pi(\alpha,\beta(\mathbf{\nu}),\mu_{\beta},\mathbf{\theta}_{\beta}\mid\bm {r})\propto\mathcal{L}(\mathbf{r}\mid\alpha,\beta(\mathbf{\nu}))\pi_{0}(\beta(\mathbf{ \nu})\mid\mu_{\beta},\mathbf{\theta}_{\beta}) \tag{6}\] \[\times\pi_{0}(\alpha,\mu_{\beta},\mathbf{\theta}_{\beta}).\]

In the posterior in Eq.(6), the dimension of \(\beta(\mathbf{\nu})\) is \(K\).

The scale is a vector of the same dimension as the data, \(\beta(\mathbf{\nu})\in\mathbb{R}_{+}^{K\times 1}\). MCMC methdos are known to struggle estimating high-dimensional parameters. At a minimum, the high-dimensional parameters incur a computational cost for inference with MCMC. To amend these issues and to simplify the inference, we perform dimension reduction for the scale \(\beta(\mathbf{\nu})\). To achieve this, we observe that our data \(\mathbf{r}\) should be a reasonable estimate for the expectation of the gamma process in Eq.(1), \(\mathbf{r}\approx\mathbb{E}[\text{Gamma}(\alpha,\beta(\mathbf{\nu}))]=\alpha\beta( \mathbf{\nu})\). This implies that the shape of the data \(\mathbf{r}\) is close to the shape of the scale function, \(\beta(\mathbf{\nu})\). Thus, we approximate the scale \(\beta(\mathbf{\nu})\) as a convolution between a Gaussian kernel and the data

\[\beta(\mathbf{\nu})\approx c_{\beta}G(\mathbf{\nu};\sigma_{G})*\mathbf{r}, \tag{7}\]

where \(*\) denotes convolution, \(c_{\beta}\) is a scaling constant, and \(G(\mathbf{\nu};\sigma_{G})\) is Gaussian smoothing kernel with width \(\sigma_{G}\). By this, we reduce the inference of the scale \(\beta(\mathbf{\nu})\in\mathbb{R}_{+}^{K\times 1}\) to inference of two parameters, \(c_{\beta}\) and \(\sigma_{G}\). With this smoothing approximation, we formulate an approximate posterior for Eq. (6) as

\[\pi(\alpha,c_{\beta},\sigma_{G},\mu_{\beta},\mathbf{\theta}_{\beta} \mid\mathbf{r})=\widetilde{\mathcal{L}}(\mathbf{r}\mid\alpha,c_{\beta},\sigma_{G})\pi_ {0}(\beta(\mathbf{\nu})\mid\mu_{\beta},\mathbf{\theta}_{\beta}) \tag{8}\] \[\times\pi_{0}(\alpha,c_{\beta},\sigma_{G},\mu_{\beta},\mathbf{\theta} _{\beta}),\]

where \(\widetilde{\mathcal{L}}(\mathbf{r}\mid\alpha,c_{\beta},\sigma_{G})=\mathcal{L}( \mathbf{r}\mid\alpha,c_{\beta}G(\mathbf{\nu};\sigma_{G})*\mathbf{r})\) and \(\pi_{0}(\alpha,c_{\beta},\sigma_{G},\mu_{\beta},\mathbf{\theta}_{\beta})\) is the prior distribution augmented with \((c_{\beta},\sigma_{G})^{\intercal}\). We detail the prior distribution \(\pi_{0}(\alpha,c_{\beta},\sigma_{G},\mu_{\beta},\mathbf{\theta}_{\beta})\) in the section on computational details and prior distributions. We perform inference of the posterior in Eq. (8) by sampling all the model parameters simultaneously using the DRAM algorithm [52].

Given samples from the posterior distribution \(\pi(\alpha,c_{\beta},\sigma_{G},\mu_{\beta},\mathbf{\theta}_{\beta}\mid\mathbf{r})\) obtained with MCMC, we can sample realizations for the synthetic spectra to generate an arbitrary amount of synthetic data in the following way. First, we sample the GP parameters \((\widetilde{\mu}_{\beta},\widetilde{\mathbf{\theta}}_{\beta})^{\intercal}\) from the MCMC chain. Next, we use \((\widetilde{\mu}_{\beta},\widetilde{\mathbf{\theta}}_{\beta})^{\intercal}\) sample a GP realization \(\widetilde{\beta}(\mathbf{\nu}^{*}\mid\widetilde{\mu}_{\beta},\widetilde{\mathbf{ \theta}}_{\beta})\) at prediction locations \(\mathbf{\nu}^{*}:=(\nu_{1}^{*},\ldots,\nu_{\widetilde{K}})^{\intercal}\) modeling the scale \(\beta(\mathbf{\nu}^{*})\) with

\[\widetilde{\beta}(\mathbf{\nu}^{*}\mid\widetilde{\mu}_{\beta},\widetilde{\mathbf{ \theta}}_{\beta})=\exp\left(\widetilde{\mu}_{\beta}+L(\mathbf{\nu}^{*}\mid \widetilde{\mathbf{\theta}}_{\beta})\mathbf{n}\right), \tag{9}\]

where \(L(\mathbf{\nu}^{*}\mid\widetilde{\mathbf{\theta}}_{\beta})\) is the lower triangular Cholesky decomposition matrix of \(\Sigma_{\beta}(\mathbf{\nu}^{*},\mathbf{\nu}^{*};\widetilde{\mathbf{\theta}}_{\beta})\) and \(\mathbf{u}:=(u_{1},\ldots,u_{\widetilde{K}})^{\intercal}\) is Gaussian white noise such that \(u_{\widetilde{K}}\sim\mathcal{N}(0,1)\). Finally, by sampling \(\widetilde{\alpha}\), we can draw a spectrum realization \(\widetilde{r}(\mathbf{\nu})\) from the gamma process, \(\text{Gamma}(\widetilde{\alpha},\widetilde{\beta}(\mathbf{\nu}))\).

We normalize the realizations \(\widetilde{r}(\mathbf{\nu})\) such that \(\max\left\{\widetilde{r}(\mathbf{\nu})\right\}=1\) and introduce an additional parameter to control amplitudes of the realizations. With an amplitude parameter \(A\), we sample a normalized shape \(\widetilde{r}_{N}(\mathbf{\nu}\mid\mathbf{\psi})\) of the spectrum and multiply this by a sampled amplitude \(\widetilde{A}\). This procedure results in the following statistical model

\[r(\mathbf{\nu}\mid A,\mathbf{\psi}) \sim A\,\frac{r_{N}(\mathbf{\nu}\mid\mathbf{\psi})}{\max r_{N}(\mathbf{\nu}\mid\mathbf{ \psi})}, \tag{10}\] \[r_{N}(\mathbf{\nu}\mid\mathbf{\psi}) \sim\text{Gamma}(\alpha,\beta(\mathbf{\nu})),\] \[A \sim\pi_{0}(A),\] \[\mathbf{\psi} \sim\pi_{0}(\mathbf{\psi}),\]

where \(\mathbf{\psi}:=(\alpha,c_{\beta},\sigma_{G},\mu_{\beta},\mathbf{\theta}_{\beta})^{\intercal}\) is a shorthand for the gamma process parameters and \(\pi_{0}(A)\) is a prior distribution for the amplitude \(A\). Example realizations from the above statistical model are shown in Fig. 3. In the following section, we detail how we model additive and multiplicative backgrounds for Raman and CARS spectra using GPs.

## 3 Additive and multiplicative background models

We propose GPs as a flexible way to randomly draw additive and multiplicative background functions for Raman and CARS spectrum modeling. This is in contrast to more standard polynomial models such as the ones used in [32].

As noted above, we model additive or multiplicative spectral backgrounds as a GP

\[e(\mathbf{\nu})\sim\text{GP}(\mu_{e},\Sigma_{e}(\mathbf{\nu},\mathbf{\nu},\mathbf{\theta}_{e})), \tag{11}\]

with \(\mu_{e}\in\mathbb{R}\) and \(\Sigma_{e}(\mathbf{\nu},\mathbf{\nu},\mathbf{\theta}_{e})\in\mathbb{R}^{K\times K}\) being a constant mean and the covariance matrix of the GP parameterized according to hyperparameters \(\mathbf{\theta}_{e}\). For the background GP covariance function,
we use again the squared exponential kernel

\[[\Sigma_{e}(\mathbf{v},\mathbf{\psi},\mathbf{\theta}_{e})]_{i,j}=\mathbf{\sigma}_{e}^{2}\exp\left( -\frac{1}{2}\frac{\left(\nu_{i}-\nu_{j}\right)^{2}}{l_{e}^{2}}\right)+\sigma_{e} ^{2}\delta(\nu_{i}-\nu_{j}) \tag{12}\]

where \([\Sigma_{e}(\mathbf{v},\mathbf{\psi},\mathbf{\theta}_{e})]_{i,j}\) denotes the \(ij\)th element of the covariance matrix, \(\mathbf{\sigma}_{e}^{2}\), is the signal variance, \(l_{e}\) is the length scale, and \(\mathbf{\sigma}_{e}^{2}\) denotes the noise variance with \(\mathbf{\theta}_{e}:=(\mathbf{\sigma}_{e,f},l_{e},\mathbf{\sigma}_{e})^{\intercal}\).

Given a measurement of the background process, \(\mathbf{e}:=(e(\nu_{1}),\ldots,e(\nu_{K}))^{\intercal}\), we can formulate a posterior distribution for the background GP parameters \((\mu_{e},\mathbf{\theta}_{e})^{\intercal}\) as

\[\pi(\mu_{e},\mathbf{\theta}_{e}\mid\mathbf{e})\propto\mathcal{L}(\mathbf{e}\mid\mu_{e}, \mathbf{\theta}_{e})\pi_{0}(\mu_{e},\mathbf{\theta}_{e}), \tag{13}\]

where \(\mathcal{L}(\mathbf{e}\mid\mu_{e},\mathbf{\theta}_{e})\) is the GP likelihood and \(\pi_{0}(\mu_{e},\mathbf{\theta}_{e})\) denotes the prior distribution for the GP parameters. The log-likelihood is given as

\[\begin{split}\log\mathcal{L}(\mathbf{e}\mid\mu_{e},\mathbf{\theta}_{e}) =&-\frac{1}{2}(\mathbf{e}-\mu_{e})^{\intercal}\Sigma_{e}(\mathbf{v}, \mathbf{\psi},\mathbf{\theta}_{e})^{-1}(\mathbf{e}-\mu_{e})\\ &-\frac{1}{2}\log|\Sigma_{e}(\mathbf{v},\mathbf{\psi},\mathbf{\theta}_{e})|- \frac{K}{2}\log 2\pi,\end{split} \tag{14}\]

where \(|\Sigma_{e}(\mathbf{v},\mathbf{\psi},\mathbf{\theta}_{e})|\) is the determinant of the covariance matrix. Again, we perform the posterior estimation for Eq. (13) by sampling all the model parameters simultaneously using DRAM.

Given a posterior \(\pi(\mu_{e},\mathbf{\theta}_{e}\mid\mathbf{e})\), we construct realizations for the spectrum by drawing realizations from the GP predictive distribution. We sample starting and ending points for the background function from priors \(\pi_{0}(e_{\text{start}})\) and \(\pi_{0}(e_{\text{step}})\), \(\pi_{0}(e_{\text{start}},e_{\text{stop}})=\pi_{0}(e_{\text{start}})\pi_{0}(e_{ \text{stop}})\). Next, we compute the predictive mean

\[e^{*}(\mathbf{v}\mid\widetilde{\mu}_{e},\widetilde{\mathbf{\theta}}_{e},\widetilde{ \mathbf{\varepsilon}}_{\text{ss}})=\Sigma_{e}(\mathbf{v},\mathbf{v}_{\text{ss}};\widetilde {\mathbf{\theta}}_{e})\Sigma_{e}(\mathbf{v}_{\text{ss}},\mathbf{v}_{\text{ss}};\widetilde {\mathbf{\theta}}_{e})^{-1}(\widetilde{\mathbf{\varepsilon}}_{\text{ss}}-\widetilde{ \mu}_{e})+\widetilde{\mu}_{e}, \tag{15}\]

and the predictive covariance

\[\begin{split}\Sigma_{e}^{*}(\mathbf{v},\mathbf{\psi};\widetilde{\mathbf{ \theta}}_{e})=\Sigma_{e}(\mathbf{v},\mathbf{\psi};\widetilde{\mathbf{\theta}}_{e})-\Sigma _{e}(\mathbf{v},\mathbf{\psi}_{\text{ss}};\widetilde{\mathbf{\theta}}_{e})\Sigma_{e}(\bm {v}_{\text{ss}},\mathbf{v}_{\text{ss}};\widetilde{\mathbf{\theta}}_{e})^{-1}\\ \times\Sigma_{e}(\mathbf{v},\mathbf{v}_{\text{ss}};\widetilde{\mathbf{\theta} }_{e})^{\intercal},\end{split} \tag{16}\]

where \((\widetilde{\mu}_{e},\widetilde{\mathbf{\theta}}_{e})\) are samples from the posterior distribution \(\pi(\mu_{e},\mathbf{\theta}_{e}\mid\mathbf{e})\) obtained via MCMC, and \(\mathbf{\psi}_{\text{ss}}=(\nu_{\text{start}},\nu_{\text{stop}})^{\intercal}\) are the wavenumber locations corresponding to the sampled starting and ending points \(\widetilde{\mathbf{\varepsilon}}_{\text{ss}}=(\widetilde{\varepsilon}_{\text{ start}},\widetilde{\varepsilon}_{\text{stop}})^{\intercal}\). Elements of the covariance matrix \(\Sigma_{e}(\mathbf{v},\mathbf{\psi};\widetilde{\mathbf{\theta}}_{e})\) are given as defined in Eq. (12) and elements of the covariance matrices \(\Sigma(\mathbf{v},\mathbf{\psi}_{\text{ss}};\widetilde{\mathbf{\theta}}_{e})\) and \(\Sigma(\mathbf{v}_{\text{ss}},\mathbf{v}_{\text{ss}};\widetilde{\mathbf{\theta}}_{e})\) are given by otherwise the same covariance function but without the diagonal elements produced by the Dirac delta function.

With the above mathematical machinations, we can sample realizations for the background function by

\[\widetilde{e}(\mathbf{v}\mid\widetilde{\mu}_{e},\widetilde{\mathbf{\theta}}_{e}, \widetilde{\mathbf{\varepsilon}}_{\text{ss}})\sim e^{*}(\mathbf{v}\mid\widetilde{\mu}_{e },\widetilde{\mathbf{\theta}}_{e},\widetilde{\mathbf{\varepsilon}}_{\text{ss}})+L(\mathbf{v }\mid\widetilde{\mu}_{e},\widetilde{\mathbf{\theta}}_{e})\mathbf{w}, \tag{17}\]

where \(L(\mathbf{v}\mid\widetilde{\mu}_{e},\widetilde{\mathbf{\theta}}_{e})\) is the lower triangular Cholesky decomposition matrix of \(\Sigma^{*}(\mathbf{v},\mathbf{\psi};\widetilde{\mathbf{\theta}}_{e})\) and \(\mathbf{w}\in\mathbb{R}^{K\times 1}\) is a Gaussian white noise vector. This is compiled into the following statistical model for the background function sampling:

\[\begin{split}\widehat{e}(\mathbf{v}\mid\mu_{e},\mathbf{\theta}_{e}, \mathbf{\varepsilon}_{\text{ss}})&\sim\mathcal{N}(\mathbf{e}^{*},\Sigma^{*}), \\ (\mu_{e},\mathbf{\theta}_{e},\mathbf{\varepsilon}_{\text{ss}})&\sim \pi(\mu_{e},\mathbf{\theta}_{e}\mid\mathbf{e})\pi_{0}(\mathbf{\varepsilon}_{\text{ss}}), \end{split} \tag{18}\]

where \(\pi_{0}(\mathbf{\varepsilon}_{\text{ss}})=\pi_{0}(e_{\text{start}},e_{\text{stop}})\). Example realizations for a multiplicative background relevant for CARS are shown in Fig. 4.

## 4 Raman and CARS spectrum models

In the preceding two sections we formulated mathematical procedures to sample synthetic spectrum and background realizations which are statistically similar to measurement data. Below, we combine these two approaches for generating arbitrary amounts of statistically realistic spectral data which are ultimately used for training our Bayesian neural networks. We present two forward models which are used to generate data for Raman measurements with an additive background and CARS measurements with a multiplicative background.

Raman spectra \(y(\mathbf{\psi})\) with an additive background \(B(\mathbf{\psi})\) are constructed using

\[y(\mathbf{\psi})\sim r(\mathbf{v}\mid A,\mathbf{\psi})+B(\mathbf{\psi}), \tag{19}\]

where \(r(\mathbf{v}\mid A,\mathbf{\psi})\) is distributed according to the model defined in

Fig. 4: Example realizations drawn from the background model defined in Eq. (18) for a multiplicative background. The starting and end points are sampled from a prior distribution and the GP predictive mean and covariance are used to sample the background shape.

Fig. 3: Example realizations drawn from the log-Gaussian gamma process model defined in Eq. (10). On the left, realizations for the scale process \(\beta(\mathbf{\psi})\), drawn from a log-Gaussian process. On the right, corresponding realizations from the gamma process. All realizations are normalized and multiplied by a sampled amplitude.
Eq. (10). The background \(B(\mathbf{\nu})\) is sampled with Eq. (18).

CARS spectra \(z(\mathbf{\nu})\) are generated similarly to the additive Raman realizations. The CARS model consists of a multiplicative background function \(\varepsilon_{\text{m}}(\mathbf{\nu}\mid\mu_{e},\mathbf{\theta}_{e})\) distorting a CARS spectrum \(S(\mathbf{\nu}\mid B_{\text{NR}},\mathbf{\psi})\) given as

\[z(\mathbf{\nu})\sim\varepsilon_{\text{m}}(\mathbf{\nu}\mid\mu_{e},\mathbf{\theta}_{e})S( \mathbf{\nu}\mid B_{\text{NR}},\mathbf{\psi}), \tag{20}\]

where the CARS spectrum \(S(\mathbf{\nu}\mid B_{\text{NR}},\mathbf{\psi})\) can be given as

\[S(\mathbf{\nu}\mid B_{\text{NR}},\mathbf{\psi})\sim|B_{\text{NR}}+(ir(\mathbf{\nu}\mid A, \mathbf{\psi})-\mathcal{H}\left\{r(\mathbf{\nu}\mid A,\mathbf{\psi})\right\})|^{2}, \tag{21}\]

and \(B_{\text{NR}}\sim\tau_{0}(B_{\text{NR}})\) is a non-resonant background inherent to the CARS phenomenon distributed according to a prior distribution \(\pi_{0}(B_{\text{NR}})\) and \(\mathcal{H}\) denotes the Hilbert transform. The model for the CARS spectrum has been previously used for example in [9, 37]. We show example realizations for the Raman model in Fig. 5 and the CARS model in Fig. 6. We use the two models defined in Eqs. (19) and (20) to generate two synthetic data sets which are used to train two separate Bayesian neural networks. In the following section, we discuss the Bayesian neural network architecture.

## 5 Bayesian neural network architecture

Our neural network architecture used in the experiments is based on the SpecNet architecture [29]. The SpecNet architecture is composed of convolutional layers encoding the input, the measurement spectrum. The encoded information is then decoded using fully-connected hidden layers, resulting in estimates for the underlying true Raman spectrum. We present our changes to the SpecNet architecture below.

To achieve a partially Bayesian neural network [44], we use a Bayesian layer for the first convolutional layer. Additionally, we augment the architecture with a probabilistic output layer. This transforms the neural network estimate into estimates of a stochastic process instead of directly estimating the Raman spectrum. We use a gamma distribution as our output layer, following our formulation of Raman spectra as a log-Gaussian gamma processes. We also found that \(L_{1}\) or \(L_{2}\) regularization was not necessary for the deterministic parts of the network and therefore only employ Dropout [53] regularization with the last dense layer of the network. This in agreement with the documented robustness of Bayesian neural networks with respect to overfitting [39]. The above results in the following partial posterior probability distribution, or cost function, used for training the neural network

\[\pi(\Psi_{\text{D}},\Psi_{\text{S}}\mid R)\propto\mathcal{L}(R\mid\Psi_{\text {D}},\Psi_{\text{S}})\pi_{0}(\Psi_{\text{S}}), \tag{22}\]

where \(R\in\mathbb{R}^{I\times J}\) is a data matrix of \(I\) synthetic spectra of length \(J\) generated using either the Raman or CARS forward models in Eqs. (19) and (20) and \(\mathcal{L}(R\mid\Psi_{\text{D}},\Psi_{\text{S}})\) denotes the likelihood of the neural network estimate and \(\pi_{0}(\Psi_{\text{S}})\) is the prior distribution for the stochastic parameters of the network. As our outputs are modeled as gamma-distributed random variables, the likelihood \(\mathcal{L}(R\mid\Psi_{\text{D}},\Psi_{\text{S}})\) is given as

\[\mathcal{L}(R\mid\Psi_{\text{D}},\Psi_{\text{S}})=\prod_{i=1}^{I}\prod_{j=1}^{ J}\frac{R_{i,j}^{\alpha_{\text{NN},j}-1}\exp(-R_{i,j}/\beta_{\text{NN},j})}{ \Gamma(\alpha_{\text{NN},j})R_{\text{NN},j}^{\alpha_{\text{NN},j}}}, \tag{23}\]

where \(R_{i,j}\) denotes the \(j\)th data point of the \(j\)th spectrum, \(\alpha_{\text{NN},j}\) and \(\beta_{\text{NN},j}\) are neural network outputs for the gamma distribution parameters. For the prior distribution \(\pi_{0}(\Psi_{\text{S}})\), we use an independent normal distribution \(\mathcal{N}(0,1)\) for all the weights and biases of the first layer, \(\pi_{0}(\Psi_{\text{S}})\propto\prod_{p=1}^{P}\mathcal{N}(\Psi_{\text{S},p};0,1)\) where \(P\) is the total number of parameters in the first layer and \(\mathcal{N}(\Psi_{\text{S},p};0,1)\) denotes the evaluation of the probability density at the parameter value \(\Psi_{\text{S},p}\). We illustrate the neural network architecture in Fig. 7.

In the log-Gaussian gamma process section, we estimate parameters of a doubly-stochastic process via MCMC. The Bayesian neural network architecture proposed here can be seen as an estimate of a _triply_-stochastic process where the neural network outputs are two stochastic process realizations \(\alpha_{\text{NN}}(\mathbf{\nu})\) and \(\beta_{\text{NN}}(\mathbf{\nu})\), an extension to the analytical log-Gaussian gamma process in Section 2 where the log-Gaussian parameterization of the scale process \(\beta_{\text{NN}}(\mathbf{\nu})\) is used for mathematical convenience due to the closed form of the probability density in Eq. (5).

Fig. 5: Example realizations for the Raman spectrum model defined in Eq. (19). The realizations correspond to the log-Gaussian gamma process realizations in Fig. 3.

Fig. 6: Example realizations for the CARS spectrum model defined in Eq. (20). The realizations correspond to the log-Gaussian gamma process realizations in Fig. 3.
## 6 Computational details and prior distributions

We use 4 experimental Raman spectra and 4 CARS spectra to generate the synthetic training data sets. We use a wavelet-based approach [11] to obtain point estimates for the underlying Raman spectra in all 8 cases. Additionally, the method provides point estimates for the additive and multiplicative background signal which we use to estimate the parameters of the background GP model defined in Eq. (18). We show the obtained Raman data point estimates for the Raman spectra and additive backgrounds in Fig. 8 and CARS point estimates for the Raman spectra and multiplicative backgrounds in Fig. 9. The four cases of measurement data are used to train their respective Bayesian neural network architectures. It should be noted that for cases with significantly different Raman spectral signatures, such as where the spectra consists of either significantly sharper or wider line shapes, the training should be done using experimental data which contain such features.

We run the DRAM algorithm with 5 proposal steps and with a length of \(100\,000\) samples for both the log-Gaussian gamma process parameters and the GP parameters. We use a burn-in of \(50\,000\) samples. The prior distributions for the log-Gaussian gamma process likelihood and the GP background likelihood are documented in Table 1. We use TensorFlow and TensorFlow Probability together with Keras to implement the neural network architecture [54, 55, 56]. We use the Adam optimizer for estimating the network parameters \(\Psi_{\text{D}}\) and \(\Psi_{\text{S}}\).
in Fig. 11. The experimental details of the CARS samples have been described in detail elsewhere, see for example [37]. The Raman spectra are from an online database of Raman spectra of pigments used in modern and contemporary art (The standard Pigments Checker v.5) [57].

Results for synthetic CARS spectra are shown in Fig. 12 and results for experimental CARS spectra of adenosine phosphate, fructose, glucose, and sucrose are presented in Fig. 13. The spectra themselves were not part of the training data set. The results show the median estimate of the Raman spectrum obtained from the trained Bayesian neural network along with the 90% confidence intervals of the Raman spectrum estimate. We overlay the Raman spectrum estimate with a scaled versions of the point estimates in Fig. 9. The point estimates are scaled such that the minima and maxima of the point estimate are equal to the minima and maxima of the median estimate of the Raman spectrum. The results coincide with the overall shape of the point estimates, supporting the validity of the data generation approach and the Bayesian neural network design.
Conclusions

We propose a novel approach utilizing log-Gaussian gamma processes and Gaussian processes to generate synthetic spectra and additive or multiplicative backgrounds that are statistically similar to experimental measurements, even when using a limited number of experimental spectra. The parameters of these stochastic processes are learned through Markov chain Monte Carlo methods, enabling the generation of extensive training data for neural networks by sampling from Bayesian posterior distributions of the parameters.

This data generation method is applied to train two Bayesian neural networks, specifically designed for correcting spectral measurements. One network is tailored for Raman spectra with additive backgrounds, while the other is optimized for coherent anti-Stokes Raman scattering (CARS) spectra with multiplicative backgrounds. Bayesian neural networks expand upon prior research involving neural networks for spectral corrections, offering not only point estimates but also the critical capability of uncertainty quantification.

Our approach is validated using synthetic test data generated from the stochastic processes and experimental Raman spectra of phthalocyanine blue, aniline black, naphthol red, and red 264 pigments, along with experimental CARS spectra of adenosine phosphate, fructose, glucose, and sucrose. The results demonstrate excellent agreement with deterministically obtained point estimates of the Raman spectra, while simultaneously providing valuable uncertainty estimates for the Raman spectrum estimates.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

The authors were supported by Academy of Finland (grant number 353095).
