# Statistical Estimation of Mean Lorentzian Line Width in Spectra by Gaussian
Processes

Erik Kuitunen  Matthew T. Moores  Teemu Härkönen<EMAIL_ADDRESS>Department of Computational Engineering, School of Engineering Sciences, LUT
University, Yliopistonkatu 34, Lappeenranta, FI-53850, Finland  National
Institute for Applied Statistics Research Australia, University of Wollongong,
Wollongong, NSW 2522, Australia

###### Abstract

We propose a statistical approach for estimating the mean line width in
spectra comprising Lorentzian, Gaussian, or Voigt line shapes. Our approach
uses Gaussian processes in two stages to jointly model a spectrum and its
Fourier transform. We generate statistical samples for the mean line width by
drawing realizations for the Fourier transform and its derivative using Markov
chain Monte Carlo methods. In addition to being fully automated, our method
enables well-calibrated uncertainty quantification of the mean line width
estimate through Bayesian inference. We validate our method using a simulation
study and apply it to an experimental Raman spectrum of $\beta$-carotene.

###### keywords:

Spectral broadening , Raman spectroscopy , Fourier transform , Gaussian
process , Markov chain Monte Carlo

###### MSC:

[2020] 62P35, 78-10

††journal: Chemometrics and Intelligent Laboratory Systems

##  1 Introduction

An observed spectrum consists of a number of individual bands at various
locations, each with an intensity and line width. Together, the positions of
the bands characterize the sample, while their intensities can be used for
quantitative analysis. However, these bands often appear within close
proximity of each other and can overlap, making them difficult to distinguish
and hence complicating the analysis [1]. Several methods have been proposed to
address this problem, such as band narrowing and curve fitting. In particular,
Fourier deconvolution [2, 3, 4] uses a kernel function with a kernel width
parameter to narrow spectral bands in order to better identify the individual
spectral bands while conserving the areas of the bands.

For the Fourier self-deconvolution to perform optimally, the kernel width
parameter needs to be manually selected and to be close to the true width of
the spectral bands. As a more reasoned approach, a method for computing a mean
Lorentzian bandwidth estimate of the individual bands has been proposed [5].
In the paper, the authors studied an interesting analytical property of
spectra consisting Lorentzian, Gaussian, or Voigt line shapes that the area-
weighted average of the Lorentzian line width parameter can be estimated via
the Fourier transform properties of the spectrum. The property has very
attractive features, it does not require any a priori information on the
number of line shapes, their exact functional form, or their parameters.
Additionally, the implementation is straight-forward with fast Fourier
transforms (FFT). However, the algorithm is hampered by the requirement to
manually select points for a polynomial fit in the Fourier domain, making it
impractical for processing large numbers of spectra. Furthermore, the lack of
uncertainty quantification makes it difficult to judge the reliability of the
result.

In this paper, we extend the methodology in [5] by automating the estimation
process while simultaneously providing well-calibrated estimates of the
associated uncertainty. First we fit a Gaussian process to measured spectrum
data. We draw realizations from the Gaussian process to generate artificial
synthetic spectra which are statistically similar to the measured spectrum.
Next, we compute Fourier transforms for each of the sampled realizations,
yielding data to which we fit a second Gaussian process. The second Gaussian
process is used to sample realizations for the Fourier transforms and their
derivatives, which are required for estimating the mean Lorentzian width of
the spectral line shapes. This sampling process is repeated until a desired
number of samples have been generated. We employ Markov chain Monte Carlo
(MCMC) methods to generate these samples, thereby obtaining fully-Bayesian
estimates of the unknown parameters of the two Gaussian processes. Using MCMC,
we also estimate a Bayesian posterior probability distribution for the area-
weighted mean Lorentzian line width. By propagating uncertainty in a
consistent and rigorous manner, we are able to provide well-calibrated
estimates of posterior intervals for the parameters and the line width. In
addition to the immediate result of estimating the line width, the results can
be used as an informative prior in subsequent Bayesian spectrum analysis, see
for example [6, 7].

The remainder of this paper is structured as follows. In Section 2 we
introduce the key mathematical properties of the original approach and
introduce our statistical model for ultimately forming the mean line width
posterior distribution. Section 3 contains computational details of model,
such as specification of the prior distributions that were used and MCMC
settings. In Section 4 we present results of the proposed approach for three
synthetic spectra consisting of Lorentzian, Gaussian, and Voigt line shapes.
In Section 5 we present results for an experimental spectrum of
$\beta$-carotene. In Section 6, we conclude the study, point out possible
alternative options, and discuss possibilities for future extensions.

##  2 Theory

###  2.1 Preliminaries

An observed spectrum consisting of Lorentzian, Gaussian, or Voigt line shapes
can be represented mathematically as

$\bm{S}_{\text{N}}\equiv\bm{S}(\bm{\nu}_{\text{N}})=\bm{f}(\bm{\nu}_{\text{N}};\bm{a},\bm{\ell},\bm{\gamma},\bm{\sigma})+\bm{\epsilon}_{\text{N}},$
(1)

where
$\bm{\nu}_{\text{N}}\equiv(\nu_{1},\dots,\nu_{N})^{\top}\in\mathbb{R}^{N}$ are
the discrete measurement locations, usually in units of wavenumbers (cm-1);
$\bm{S}_{\text{N}}\equiv\left(S(\nu_{1}),\dots,S(\nu_{N})\right)^{\top}\in\mathbb{R}^{N}$
are the measurements, often in scientific arbitrary units (a.u.);
$f(\nu;\bm{a},\bm{\ell},\bm{\gamma},\bm{\sigma})$ is a continuous mathematical
model of the underlying spectral signature, with vectors of parameters
$\bm{a}$ for the areas of the line shapes, $\bm{\ell}$ for their locations,
$\bm{\gamma}$ for the Lorentzian line widths, and $\bm{\sigma}$ for the
Gaussian line widths; and
$\bm{\epsilon}_{\text{N}}\equiv(\epsilon_{1},\dots,\epsilon_{\text{N}})^{\top}\in\mathbb{R}^{N}$
are additive measurement errors, where
$\epsilon_{\text{n}}\sim\mathcal{N}(0,\sigma^{2}_{\epsilon})$ is assumed to be
white noise with unknown variance $\sigma^{2}_{\epsilon}$. We follow [6, 7] in
modelling the spectral signature as a linear combination of Lorentzian,
Gaussian, or Voigt line shapes,

$f(\nu;\bm{a},\bm{\ell},\bm{\gamma},\bm{\sigma})=\sum\limits_{m=1}^{M}a_{m}\mathcal{K}(\nu-\ell_{m};\gamma_{m},\sigma_{m}),$
(2)

where $a_{m}\in\bm{a}>0$, $\ell_{m}\in\bm{\ell}>0$,
$\gamma_{m}\in\bm{\gamma}\geq 0$, and $\sigma_{m}\in\bm{\sigma}\geq 0$ are the
parameters of the $M$ constituent line shapes, and where the line-shape
function $\mathcal{K}(\nu;\gamma_{m},\sigma_{m})$ can be represented as

$\mathcal{K}(\nu;\gamma_{m},\sigma_{m})=\frac{1}{\pi\gamma_{m}}\frac{\gamma_{m}^{2}}{\nu^{2}+\gamma_{m}^{2}}*\frac{1}{\sqrt{2\pi\sigma_{m}^{2}}}\exp\left(-\frac{\nu^{2}}{2\sigma_{m}^{2}}\right),$
(3)

where $\ast$ denotes convolution with respect to the wavenumbers $\nu$. This
is the Voigt line-shape function that corresponds to the Lorentzian line shape
in the limit when $\sigma_{m}=0$ and to the Gaussian line shape when
$\gamma_{m}=0$ [8]. As the Voigt line-shape function does not have a closed-
form expression, we use the pseudo-Voigt approximation, which is defined as

$\widetilde{\mathcal{K}}(\nu;\gamma_{m},\sigma_{m})=\eta\cdot\frac{1}{\pi\gamma_{m}}\frac{\gamma_{m}^{2}}{\nu^{2}+\gamma_{m}^{2}}+(1-\eta)\cdot\frac{1}{\sqrt{2\pi\sigma_{m}^{2}}}\exp\left(-\frac{\nu^{2}}{2\sigma_{m}^{2}}\right),$
(4)

where $\eta\in[0,1]$ is the mixing parameter, which is determined by the line-
width parameters $\gamma$ and $\sigma$. For further details, see [9, 10, 11].

For spectra such as above, we have the interesting property that

$\overline{\gamma}=\lim_{\xi\searrow
0}\widetilde{\gamma}(\xi)=\lim_{\xi\searrow
0}-\,\frac{1}{\pi\left|\widehat{S}(\xi)\right|}\times\frac{\partial\left|\widehat{S}(\xi)\right|}{\partial\xi}=\frac{\sum^{M}_{m=1}a_{m}\gamma_{m}}{\sum^{M}_{m=1}a_{m}},$
(5)

where $\overline{\gamma}$ denotes the area-weighted average of the Lorentzian
line widths and $\widehat{S}(\xi)$ denotes the Fourier transform of the
spectrum $S(\nu)$ [5]. It is notable that the property in Eq. (5) does not
require any a priori information on the line shape parameters, nor on the
number of them. In the following section, we formulate our Bayesian model for
estimating the mean Lorentzian line width with a two-stage, Gaussian process
approach.

###  2.2 Statistical model

We assume the spectrum in Eq. (1) to be distributed according to a Gaussian
process:

$\bm{S}_{\text{N}}\sim\mathcal{GP}\left(\bm{\mu}(\bm{\nu}_{\text{N}};\alpha),\,\Sigma(\bm{\nu}_{\text{N}},\bm{\nu}_{\text{N}};\bm{\theta})+\sigma_{\epsilon}^{2}I_{\text{N}}\right),$
(6)

where $\bm{S}_{\text{N}}$, $\bm{\nu}_{\text{N}}$ and $\sigma_{\epsilon}^{2}$
are as given above; $\bm{\mu}(\bm{\nu}_{\text{N}};\alpha)$ is the mean
function and
$\Sigma(\bm{\nu}_{\text{N}},\bm{\nu}_{\text{N}};\bm{\theta})\in\mathbb{R}^{N\times
N}$ is the symmetric, positive-definite covariance matrix of the Gaussian
process, parameterized according to $\alpha$ and $\bm{\theta}$, respectively;
and $I_{\text{N}}$ is the $N\times N$ identity matrix. Each $i,j$th element of
the covariance matrix $\Sigma(\bm{\nu},\bm{\nu};\bm{\theta})$ is defined as

$\left[\Sigma(\bm{\nu},\bm{\nu};\bm{\theta})\right]_{ij}\coloneqq
k(\nu_{i},\nu_{j};\bm{\theta}),$ (7)

where $i,j\in 1,\dots,M$. We use a squared-exponential kernel function:

$k(\nu_{i},\nu_{j};\bm{\theta})=\sigma_{\text{s}}^{2}\exp\left\\{-\frac{1}{2}\frac{(\nu_{i}-\nu_{j})^{2}}{\varphi^{2}}\right\\},$
(8)

where $\bm{\theta}\equiv(\sigma_{s}^{2},\varphi)^{\top}$ is a vector of the GP
signal variance and length scale parameters. For the GP mean function, we use
a simple constant mean, $\mu(\nu;\alpha)=\alpha$.

We formulate the task of estimating the model parameters $\alpha$,
$\bm{\theta}$, and $\sigma_{\epsilon}$ from the observed spectrum
$\bm{S}_{\text{N}}$ as a statistical inference problem in terms of Bayes’
theorem. At its simplest, a Bayesian posterior distribution comprises two
parts: the likelihood
$\mathcal{L}(\bm{S}_{\text{N}}\mid\alpha,\bm{\theta},\sigma_{\epsilon})$ of
the model fitting the data; and a prior distribution
$p(\alpha,\bm{\theta},\sigma_{\epsilon})$ encoding a priori known information
regarding the model parameters. More explicitly, a posterior distribution for
the model parameters in Eq. (6) can be given as

$p(\alpha,\bm{\theta},\sigma_{\epsilon}\mid\bm{S}_{\text{N}})\propto\mathcal{L}(\bm{S}_{\text{N}}\mid\alpha,\bm{\theta},\sigma_{\epsilon})\;p(\alpha,\bm{\theta},\sigma_{\epsilon}),$
(9)

up to an unknown normalizing constant, $p(\bm{S}_{\text{N}})$. The natural
logarithm of the GP model likelihood is given by

$\begin{split}\log\mathcal{L}(\bm{S}_{\text{N}}\mid\alpha,\bm{\theta},\sigma_{\epsilon})=&-\frac{1}{2}\left(\bm{S}_{\text{N}}-\mu(\bm{\nu};\alpha)\right)\left(\Sigma(\bm{\nu},\bm{\nu};\bm{\theta})+\sigma_{\epsilon}^{2}I_{\text{N}}\right)^{-1}\left(\bm{S}_{\text{N}}-\mu(\bm{\nu};\alpha)\right)^{\top}\\\
&-\frac{1}{2}\log\left|\Sigma(\bm{\nu},\bm{\nu};\bm{\theta})+\sigma_{\epsilon}^{2}I_{\text{N}}\right|-\frac{N}{2}\log
2\pi,\end{split}$ (10)

where
$\left|\Sigma(\bm{\nu},\bm{\nu};\bm{\theta})+\sigma_{\epsilon}^{2}I_{\text{N}}\right|$
is the matrix determinant of
$\Sigma(\bm{\nu},\bm{\nu};\bm{\theta})+\sigma_{\epsilon}^{2}I_{\text{N}}$. The
posterior distribution
$p(\alpha,\bm{\theta},\sigma_{\epsilon}\mid\bm{S}_{\text{N}})$ is analytically
intractable and thus we need to resort to numerical methods. We generate
random samples from the aforementioned posterior distribution using Markov
chain Monte Carlo, in particular with the delayed-rejection adaptive
Metropolis (DRAM) algorithm [12]. We provide details of the prior distribution
in Section 3.

Given samples from the posterior
$p(\alpha,\bm{\theta},\sigma_{\epsilon}\mid\bm{S}_{\text{N}})$, we can draw
samples, or realizations, from the predictive distribution of the GP in Eq.
(6) using the predictive mean and covariance of the associated GP fit. The
predictive mean and covariance at prediction locations $\bm{\nu}^{*}$ are
constructed with

$\bm{\mu}^{*}(\bm{\nu}^{*};\alpha,\bm{\theta},\sigma_{\epsilon})=\Sigma(\bm{\nu}^{*},\bm{\nu};\bm{\theta})\left(\Sigma(\bm{\nu},\bm{\nu};\bm{\theta})+\sigma_{\epsilon}^{2}I\right)^{-1}\left(\bm{S}_{\text{N}}-\mu(\bm{\nu};\alpha)\right)+\mu(\bm{\nu}^{*};\alpha),$
(11)

and

$\Sigma^{*}(\bm{\nu}^{*};\bm{\theta},\sigma_{\epsilon})=\Sigma(\bm{\nu}^{*},\bm{\nu}^{*};\bm{\theta})-\Sigma(\bm{\nu}^{*},\bm{\nu};\bm{\theta})\left(\Sigma(\bm{\nu},\bm{\nu};\bm{\theta})+\sigma_{\epsilon}^{2}I\right)^{-1}\Sigma(\bm{\nu}^{*},\bm{\nu};\bm{\theta})^{\top}.$
(12)

The GP is defined as a continuous random function in $\mathbb{R}$, therefore
the prediction locations $\bm{\nu}^{*}$ can be freely chosen. We use the
measurement locations as our prediction locations,
$\bm{\nu}^{*}=\bm{\nu}_{\text{N}}$. Realizations with noise
$\widetilde{\bm{S}}$ can then be drawn with

$\widetilde{\bm{S}}=\bm{\mu}^{*}(\bm{\nu}^{*};\alpha,\bm{\theta},\sigma_{\epsilon})+L\bm{w}+\sigma_{\epsilon}\bm{e},$
(13)

where $L$ is the lower-triangular Cholesky decomposition of the predictive
covariance $\Sigma^{*}(\bm{\nu}^{*};\bm{\theta},\sigma_{\epsilon})$, with
$\bm{w}=(w_{1},\dots,w_{N})^{\top}$ and $\bm{e}=(e_{1},\dots,e_{N})^{\top}$
being vectors of independent, standard normal random variables such that
$w_{n},e_{n}\sim\mathcal{N}(0,1)$. This approach can be used to generate an
arbitrary amount of artificial “measurements” which are statistically similar
to the original observed spectrum $\bm{S}_{\text{N}}$. An illustration is
shown in Figure 1. The above formulations constitute the first stage of our
statistical method.

Figure 1: On top, a spectrum consisting of 3 Lorentzian line shapes in blue
along with an example Gaussian process fit. The solid black line and the
shaded area are the predictive mean and 95% confidence intervals of the
Gaussian process. On the bottom, realizations drawn from the Gaussian process
in blue. The realizations can be used to generate additional artificial data.

Next, we randomly select $J$ samples from the MCMC output for the posterior
distribution in Eq. (9). We denote these samples with
$\left\\{(\alpha,\bm{\theta},\sigma_{\epsilon})_{1},\dots,(\alpha,\bm{\theta},\sigma_{\epsilon})_{J}\right\\}$.
For each sample $(\alpha,\bm{\theta},\sigma_{\epsilon})_{j}$, we compute the
predictive mean and covariance using Eqs. (11) and (12) and use Eq. (13) to
draw a realization $\widetilde{\bm{S}}_{j}$. This results in a set of $J$
realizations at $N$ wavenumber locations $\bm{\nu}^{*}$. Then we compute the
FFT of each sample, which we use to model the Fourier transforms in Eq. (5).
As we are only interested in the behaviour of Eq. (5) in the limit as $\xi$
approaches zero, we additionally truncate the FFTs of the realizations to an
arbitrary length, $P$. This also eases the computational burden by reducing
the size of the data set. We denote the data set consisting of the truncated
FFTs by
$\bm{Z}=\left(\widehat{\bm{S}}_{1}^{1:P},\dots,\widehat{\bm{S}}_{J}^{1:P}\right)^{\top}\in\mathbb{R}^{JP}$
where $\widehat{\bm{S}}_{j}$ is the FFT of $\widetilde{\bm{S}}_{j}$. The
corresponding frequency locations are denoted by
$\bm{\xi}=(0,\dots,\xi_{P},\dots,0,\dots,\xi_{P})^{\top}\in\mathbb{R}^{JP}$
where the locations $(0,\dots,\xi_{P})^{\top}$ are repeated $J$ times, as the
locations are identical for each $\widehat{\bm{S}}_{j}^{1:P}$ in $\bm{Z}$.
With the above, we are ready to formulate the statistical model for the
Fourier-transformed realizations, similarly to what we did for the measurement
spectrum.

We assume $\bm{Z}$ to be distributed according to a Gaussian process:

$\bm{Z}\sim\mathcal{GP}(\phi(\bm{\xi};\bm{\beta}),C(\bm{\xi},\bm{\xi};\bm{\psi})+\sigma_{z}^{2}I),$
(14)

where $\phi(\bm{\xi};\bm{\beta})\in\mathbb{R}^{JP}$ denotes the mean function
and $C(\bm{\xi},\bm{\xi};\bm{\psi})\in\mathbb{R}^{JP\times JP}$ denotes the
covariance matrix, parameterized according to $\beta$ and $\bm{\psi}$; and
with independent and identically-distributed white noise, $\sigma_{z}^{2}I$.
We again use the squared exponential kernel

$c_{0,0}(\xi_{i},\xi_{j};\bm{\psi})=\sigma_{c}^{2}\exp\left\\{-\frac{1}{2}\frac{(\xi_{i}-\xi_{j})^{2}}{\lambda^{2}}\right\\},$
(15)

where $i,j\in(1,\dots,JP)^{T}$ $\bm{\psi}=(\sigma_{c}^{2},\lambda)^{T}$ is a
vector of the GP signal and length scale parameters. Each element of the
covariance matrix $C(\bm{\xi},\bm{\xi};\bm{\psi})$ is constructed with the
covariance function $c_{0,0}(\xi_{i},\xi_{j};\bm{\psi})$ as in Eq. (7). For
the GP mean function, we now use an exponential function
$\phi(\xi;\bm{\beta})=\beta_{0}\exp(\beta_{1}\xi)$ where
$\bm{\beta}=(\beta_{0},\beta_{1})^{T}$. The posterior distribution for the
model parameters in Eq. (14) is given as

$p(\bm{\beta},\bm{\psi},\sigma_{z}\mid\bm{Z})\propto\mathcal{L}(\bm{Z}\mid\bm{\beta},\bm{\psi},\sigma_{z})p(\bm{\beta},\bm{\psi},\sigma_{z}),$
(16)

where $\mathcal{L}(\bm{Z}\mid\bm{\beta},\bm{\psi},\sigma_{z})$ is the model
likelihood and $p(\bm{\beta},\bm{\psi},\sigma_{z})$ is the prior distribution
for the model parameters. The log-likelihood of the model is now given as

$\begin{split}\log\mathcal{L}(\bm{Z}\mid\bm{\beta},\bm{\psi},\sigma_{z})=&-\frac{1}{2}\left(\bm{Z}-\phi(\bm{\xi};\bm{\beta})\right)\left(C(\bm{\xi},\bm{\xi};\bm{\psi})+\sigma_{z}^{2}I\right)^{-1}\left(\bm{Z}-\phi(\bm{\xi};\bm{\beta})\right)^{T}\\\
&-\frac{1}{2}\log\left|\Sigma(\bm{\xi},\bm{\xi};\bm{\psi})+\sigma_{z}^{2}I\right|-\frac{JP}{2}\log
2\pi.\end{split}$ (17)

We document the prior distributions in Section 3. To obtain samples from the
posterior distribution in Eq. (16), we again employ the DRAM algorithm.

Next, we draw samples in order to model the Fourier transform and its
derivative in Eq. (5). However, this requires some additional mathematical
machinery, see for example [13, 14], in comparison to Eq. (13) which we will
define below. The predictive mean for the Fourier transform and its derivative
at prediction locations $\bm{\xi}^{*}$ are now given as

$\begin{bmatrix}\phi^{*}(\bm{\xi}^{*};\bm{\beta},\bm{\psi},\sigma_{z})\\\
\phi^{\prime*}(\bm{\xi}^{*};\bm{\beta},\bm{\psi},\sigma_{z})\end{bmatrix}=C(\bm{\xi}^{*},\bm{\xi};\bm{\psi})\left(C(\bm{\xi},\bm{\xi};\bm{\psi})+\sigma_{z}^{2}I\right)^{-1}\left(\bm{Z}-\phi(\bm{\xi};\bm{\beta})\right)+\begin{bmatrix}\phi(\bm{\xi}^{*};\bm{\beta})\\\
\phi^{\prime}(\bm{\xi}^{*};\bm{\beta})\end{bmatrix}$ (18)

where $\phi^{\prime}(\bm{\xi}^{*};\bm{\beta})$ is the derivative of
$\phi(\bm{\xi}^{*};\bm{\beta})$ with respect to $\xi$ and the corresponding
predictive covariance is given as

$C^{*}(\bm{\xi}^{*};\bm{\psi},\sigma_{z})=C(\bm{\xi}^{*},\bm{\xi}^{*};\bm{\psi})-C(\bm{\xi}^{*},\bm{\xi};\bm{\psi})\left(C(\bm{\xi},\bm{\xi};\bm{\psi})+\sigma^{2}_{z}I\right)^{-1}C(\bm{\xi}^{*},\bm{\xi};\bm{\psi})^{T},$
(19)

where the covariance matrices $C(\bm{\xi}^{*},\bm{\xi};\bm{\psi})$ and
$C(\bm{\xi}^{*},\bm{\xi}^{*};\bm{\psi})$ are constructed as

$C(\bm{\xi}^{*},\bm{\xi};\bm{\theta})=\begin{bmatrix}c_{0,0}(\bm{\xi}^{*},\bm{\xi};\bm{\psi})\\\
c_{1,0}(\bm{\xi}^{*},\bm{\xi};\bm{\psi})\end{bmatrix},$ (20)

and

$C(\bm{\xi}^{*},\bm{\xi}^{*};\bm{\theta})=\begin{bmatrix}c_{0,0}(\bm{\xi}^{*},\bm{\xi}^{*};\bm{\psi})&c_{0,1}(\bm{\xi}^{*},\bm{\xi}^{*};\bm{\psi})\\\
c_{1,0}(\bm{\xi}^{*},\bm{\xi}^{*};\bm{\psi})&c_{1,1}(\bm{\xi}^{*},\bm{\xi}^{*};\bm{\psi})\end{bmatrix},$
(21)

where the $ij$th elements of each of the individual covariance matrices are
defined as

$\left[C_{\bm{\cdot},\bm{\cdot}}(\bm{\xi},\bm{\xi};\bm{\psi})\right]_{ij}\coloneqq
c_{\bm{\cdot},\bm{\cdot}}(\xi_{i},\xi_{j};\bm{\psi}).$ (22)

The three additional covariance functions can be obtained by differentiating
the covariance function $c_{0,0}(\xi_{i},\xi_{j};\bm{\psi})$ defined in Eq.
(15):

$\displaystyle c_{0,1}(\xi_{i},\xi_{j};\bm{\psi})$
$\displaystyle=\frac{\partial
c_{0,0}(\xi_{i},\xi_{j};\bm{\psi})}{\partial\xi_{j}}$
$\displaystyle=\frac{\sigma^{2}_{c}}{\lambda^{2}}(\xi_{i}-\xi_{j})\exp\left\\{-\frac{1}{2}\frac{(\xi_{i}-\xi_{j})^{2}}{\lambda^{2}}\right\\},$
(23) $\displaystyle c_{1,0}(\xi_{i},\xi_{j};\bm{\psi})$
$\displaystyle=\frac{\partial
c_{0,0}(\xi_{i},\xi_{j};\bm{\psi})}{\partial\xi_{i}}$
$\displaystyle=\frac{\sigma^{2}_{c}}{\lambda^{2}}(\xi_{i}-\xi_{j})\exp\left\\{-\frac{1}{2}\frac{(\xi_{i}-\xi_{j})^{2}}{\lambda^{2}}\right\\},$
$\displaystyle c_{1,1}(\xi_{i},\xi_{j};\bm{\psi})$
$\displaystyle=\frac{\partial^{2}c_{0,0}(\xi_{i},\xi_{j};\bm{\psi})}{\partial\xi_{j}\partial\xi_{i}}$
$\displaystyle=\frac{\sigma_{c}^{2}}{\lambda^{2}}\left(\lambda^{2}-(\xi_{i}-\xi_{j})^{2}\right)\exp\left\\{-\frac{1}{2}\frac{(\xi_{i}-\xi_{j})^{2}}{\lambda^{2}}\right\\},$

where $c_{0,1}(\xi_{i},\xi_{j};\bm{\psi})$ along with
$c_{1,0}(\xi_{i},\xi_{j};\bm{\psi})$ model the covariance between a
realization and its derivative and $c_{1,1}(\xi_{i},\xi_{j};\bm{\psi})$ models
the covariance of a derivate realization similarly to how
$c_{0,0}(\xi_{i},\xi_{j};\bm{\psi})$ models the covariance of a realization.

Finally, we are able to draw realizations to model the Fourier transforms in
Eq. (5) with

$\begin{bmatrix}\widetilde{\bm{Z}}\\\
\widetilde{\bm{Z}}^{\prime}\end{bmatrix}=\begin{bmatrix}\phi^{*}(\bm{\xi}^{*};\bm{\beta},\bm{\psi},\sigma_{z})\\\
\phi^{\prime*}(\bm{\xi}^{*};\bm{\beta},\bm{\psi},\sigma_{z})\end{bmatrix}+Q\bm{r},$
(24)

where $\widetilde{\bm{Z}}$ and $\widetilde{\bm{Z}}^{\prime}$ denote
realizations for the Fourier transform and its derivative, respectively, $Q$
is the lower triangular Cholesky decomposition matrix of the predictive
covariance $C^{*}(\bm{\xi}^{*};\bm{\psi},\sigma_{z})$ and
$\bm{r}=(r_{1},\dots,r_{N})^{T}$ is a vector of standard normal variables such
that $r_{n}\sim\mathcal{N}(0,1)$. We can then compute an estimate for the mean
Lorentzian line width with

$\overline{\gamma}\sim-\frac{\widetilde{\bm{Z}}^{\prime}_{0}}{\pi\widetilde{\bm{Z}}_{0}},$
(25)

where $\widetilde{\bm{Z}}_{0}$ and $\widetilde{\bm{Z}}^{\prime}_{0}$ are the
realizations evaluated at $\xi=0$. We sample estimates for the mean gamma line
width $\overline{\gamma}$ by repeatedly sampling
$(\bm{\beta},\bm{\psi},\sigma_{z})$ from the MCMC chain for the posterior
distribution in Eq. (16) and by computing Eq. (25). This yields us our
ultimate goal, a posterior distribution for $\overline{\gamma}$ which we
denote by $p(\overline{\gamma}\mid\bm{Z})$. We show an illustration of
$\bm{Z}$ with a respective GP fit, the respective GP predictive distribution
for the derivative, mean line width function $\widetilde{\gamma}(\xi)$, and
the posterior distribution $p(\overline{\gamma}\mid\bm{Z})$ in Figure 2. We
summarize the above in Algorithm 1.

Figure 2: On the top left, truncated FFTs computed for Gaussian process
realizations, in blue, shown in Figure 1 along with the associated Gaussian
process fit. The solid black line corresponds to the predictive mean
$\phi^{*}(\bm{\xi}^{*};\bm{\beta},\bm{\psi},\sigma_{z})$ and the shaded area
is the 95% confidence interval. On top right, the same GP is shown in black
and gray along with the the predictive mean for the derivative
$\phi^{\prime*}(\bm{\xi}^{*};\bm{\beta},\bm{\psi},\sigma_{z})$ and its 95%
confidence interval. On the bottom left, the expectation and 95% confidence
interval for $\widetilde{\gamma}(\nu)$ obtained by sampling Eq. (25). The
behaviour of $\widetilde{\gamma}(\nu)$ at $\xi=0$ corresponds to the posterior
distribution for the mean Lorentzian line width,
$p(\overline{\gamma}\mid\bm{Z})$. The true value used to generate the data in
Figure 1 is illustrated with a solid vertical black line. Algorithm 1 Two-
stage Gaussian process estimation of mean Lorentzian lime width

Step 1: Fit a GP to the measurement spectrum $\bm{S}$ with MCMC using Eq. (9).

Step 2: Draw $J$ realizations $\widetilde{\bm{S}}$ via Eq. (13).

Step 3: Form the truncated FFT data set $\bm{Z}$ and fit a GP using Eq. (16)
and MCMC.

Step 4: Draw $J_{z}$ Fourier-domain realizations with Eq. (24).

Step 5: Compute the mean Lorentzian line width for the $J_{z}$ realizations
with Eq. (25).

##  3 Computational details

The posterior distributions defined in Eqs. (9) and (16) have prior
distributions $p(\alpha,\bm{\theta},\sigma_{\epsilon})$ and
$p(\bm{\beta},\bm{\psi},\sigma_{z})$. We use simple uniform priors for all our
parameters to enforce positivity of the parameters which are either physically
or mathematically positive such as the mean Lorentzian line width
$\overline{\gamma}$ or standard deviation parameters such as
$\sigma_{\epsilon}$ and $\sigma_{z}$. Additionally, all our prior
distributions are modelled as jointly independent distributions such that
$p(\alpha,\bm{\theta},\sigma_{\epsilon})=p(\alpha)p(\sigma_{s})p(\ell)p(\sigma_{\epsilon})$
and
$p(\bm{\beta},\bm{\psi},\sigma_{z})=p(\beta_{1})p(\beta_{2})p(\sigma_{c})p(\lambda)p(\sigma_{z})$.
We detail these individual distributions in Table 1. We also enforce
positivity on the posterior $p(\overline{\gamma}\mid\bm{Z})$.

Table 1: Prior distributions for the GP parameters. Parameter | Distribution | Parameter | Distribution
---|---|---|---
$p(\sigma_{s})$ | $\mathcal{U}(0,\infty)$ | $p(\sigma_{c})$ | $\mathcal{U}(0,\infty)$
$p(\ell)$ | $\mathcal{U}(0,2\left(\nu_{N}-\nu_{1}\right))$ | $p(\lambda)$ | $\mathcal{U}(0,3\xi_{P})$
$p(\sigma_{\epsilon})$ | $\mathcal{U}(0,\infty)$ | $p(\sigma_{z})$ | $\mathcal{U}(0,\infty)$
$p(\alpha)$ | $\mathcal{U}(0,\infty)$ | $p(\beta_{1})$ | $\mathcal{U}(0,10\times\max\bm{Z})$
|  | $p(\beta_{2})$ | $\mathcal{U}(-\infty,\infty)$

We run the DRAM algorithm to sample both posterior distributions
$p(\alpha,\bm{\theta},\sigma_{\epsilon}\mid\bm{S})$ and
$p(\bm{\beta},\bm{\psi},\sigma_{z}\mid\bm{Z})$ in Eqs. (9) and (16) with 3
proposal steps per each MCMC iteration and with a chain length of $50000$. We
use a burn-in of $25000$.

##  4 Numerical examples

We apply the method proposed in Section 2 to three synthetic spectra with
additive noise. The three synthetic data sets consist of a spectrum simulated
with $M=8$ Lorentzian line shapes, one simulated with $M=10$ Gaussian line
shapes, and lastly of a spectrum with $M=6$ Voigt line shapes. All parameters,
$(\bm{a},\bm{l},\bm{\gamma},\bm{\sigma})$, and noise level of the synthetic
Lorentzian and Gaussian spectra are randomly sampled from uniform
distributions which detail in Table 2. The synthetic Voigt spectrum is
calculated using pseudo-Voigt approximation. A log-normal prior
$\delta_{m}\sim\mathrm{Lognormal}(\mu_{V},\sigma_{V}^{2})$ with
$\sigma_{V}^{2}=0.16$ and $\mu_{V}=\ln(25)-\sigma_{V}^{2}/2$ is used for Voigt
line width, and uniform prior $\mathcal{U}(0,\delta_{m})$ is used for the
Lorentzian width needed for the calculation of pseudo-Voigt approximation. The
same priors for parameters $\bm{a}$ and $\bm{l}$ are used as with Lorentzian
and Gaussian spectra.

The synthetic Lorentzian spectrum along with the obtained area-weighted mean
Lorentzian line width posterior distribution are illustrated in Figure 3. The
Gaussian spectrum and the corresponding posterior distribution is shown in
Figure 6. In Figure 7, the Voigt spectrum and the posterior of the mean
Lorentzian line width can be seen. The posterior distribution for the
Lorentzian spectrum in Figure 3 contains the true value for the area-weighted
mean Lorentzian width. The same is true for the Gaussian spectrum in Figure 6
which exhibits also a far higher degree of uncertainty for the estimate. The
Voigt spectrum posterior in Figure 7 likewise contains the true value of the
mean width, with similar uncertainty to the Lorentzian spectrum. The results
for the Gaussian and Voigt spectrum are presented in Appendix A.

Table 2: Distributions of the spectrum parameters $(\bm{a},\bm{l},\bm{\gamma},\bm{\sigma})$ and measurement standard deviation $\sigma_{\epsilon}$ used to generate $M$ line shapes for generating the synthetic Lorentzian and Gaussian spectra. The top row corresponds to the prior distributions for the Lorentzian spectrum and the bottom row to the prior distributions for the Gaussian spectrum. Area | Location | Width | Noise level
---|---|---|---
$a_{m}\sim\mathcal{U}(1,30)$ | $l_{m}\sim\mathcal{U}(1625,1675)$ | $\gamma_{m}\sim\mathcal{U}(2.5,20)$ | $\sigma_{\epsilon}=0.05\times\max\bm{S}$
$a_{m}\sim\mathcal{U}(1,30)$ | $l_{m}\sim\mathcal{U}(1625,1675)$ | $\sigma_{m}\sim\mathcal{U}(10,30)$ | $\sigma_{\epsilon}=0.05\times\max\bm{S}$

Figure 3: On top, a spectrum consisting of $M=8$ Lorentzian line shapes in
blue along with the individual overlapping line shapes illustrated with black
dashed lines. On bottom, the obtained posterior distribution
$p(\overline{\gamma}\mid\bm{Z})$ for the area-weighted mean Lorentzian line
width. The true value for $\overline{\gamma}=17.12$ used to generate the data
is shown with a vertical solid line.

##  5 Experimental results

In addition to the synthetic validation spectra, we apply the method to an
experimental Raman spectrum of $\beta$-carotene [15]. We found that the
results were sensitive to the selected region of data used for the analysis.
Wide regions of data where the spectrum has multiple individual areas of peaks
exhibit blowup for the posterior confidence intervals. Smaller, more
constrained intervals with less features present are better handles by the
algorithm, producing worthwhile estimates for the mean Lorentzian line width
and its confidence intervals. We present results for different selections of
data in Table 3. As a showcase, we illustrate the results for a single
selection of data in Figure 5. Illustrations for the other selections of data
are shown in A.

Table 3: Results obtained for different choices for the experimental $\beta$-carotene Raman spectrum for different regions of data. Interval [cm-1] | Interval width [cm-1] | Posterior mean | 95% Confidence interval
---|---|---|---
$[1240,1060]$ | $180$ | $23.2$ | $[19.5,27.0]$
$[1300,1060]$ | $240$ | $36.2$ | $[26.3,46.1]$
$[1360,900]$ | $460$ | $61.5$ | $[44.2,78.8]$
$[1615,900]$ | $715$ | $2.9\cdot 10^{6}$ | $[0,1.4\cdot 10^{7}]$
$[1615,1060]$ | $555$ | $260.4$ | $[0,3.8\cdot 10^{3}]$
$[1615,1360]$ | $255$ | $24.1$ | $[22.3,26.0]$
Figure 4: An experimental spectrum of $\beta$-carotene in blue. Black
horizontal lines show the selected region ($[1240,1060]$) of the spectrum.

Figure 5: On top, the region between $[1240,1060]$ of the experimental
spectrum of $\beta$-carotene shown in Figure 4. On bottom, the obtained
posterior distribution $p(\overline{\gamma}\mid\bm{Z})$ for the area-weighted
mean Lorentzian line width.

##  6 Conclusions

We have proposed and implemented a statistical algorithm for estimating the
area-weighted mean Lorentzian line width in spectra consisting of Lorentzian,
Gaussian, or Voigt line shapes. This extends previously introduced methodology
by automating the estimation process and by enabling uncertainty
quantification with Gaussian processes and Markov chain Monte Carlo methods in
a two-stage implementation.

The method is validated with three synthetic spectra. In all three cases, the
proposed method yields posterior distributions which contain the true
parameter value used to generate the synthetic spectra. Additionally, the
method is applied to an experimental spectrum of $\beta$-carotene.

The method allows for some straight-forward modifications. The squared
exponential covariance function in Eq. (8) can be replaced with the more
general Matérn covariance function with a fixed regularity parameter $\nu$ or
so that the regularity parameter is also estimated with MCMC. The squared
exponential covariance function is a special case of the Matérn covariance
function when $\nu\rightarrow\infty$. Discussion on different covariance
functions can be found for example in [13]. For computational speed-up, MCMC
sampling can be replaced with optimization. This would yield a single point
estimate for the GP parameters. This was however found to have an effect on
the resulting mean Lorentzian line width posterior distribution,
$p(\overline{\gamma}\mid\bm{Z})$. An example of a more involved extension
would be using mixtures of Gaussian process experts where multiple Gaussian
processes are used to model different parts of data with individual Gaussian
processes [16, 17, 18]. This could come into play with spectra with wide tails
or where the spectrum is not adequately modelled by a single GP.

The statistical approach also enables alternative approaches to further
validate the approach with simulation-based calibration [19, 20]. This has
been previously considered in the context of spectrum analysis by the authors
[21].

## Acknowledgements

This work was supported by the Research Council of Finland through the
Flagship of Advanced Mathematics for Sensing, Imaging and Modelling (decision
number 359183).

## References

  * [1] L. Antonov, D. Nedeltcheva, Resolution of overlapping UV–Vis absorption bands and quantitative analysis, Chem. Soc. Rev. 29 (2000) 217–227.  doi:10.1039/A900007K.
  * [2] J. K. Kauppinen, D. J. Moffatt, M. R. Hollberg, H. H. Mantsch, A new line-narrowing procedure based on Fourier self-deconvolution, maximum entropy, and linear prediction, Applied Spectroscopy 45 (3) (1991) 411–416.  doi:10.1366/0003702914337155.
  * [3] J. K. Kauppinen, D. J. Moffatt, H. H. Mantsch, D. G. Cameron, Fourier self-deconvolution: A method for resolving intrinsically overlapped bands, Applied Spectroscopy 35 (3) (1981) 271–276.  doi:10.1366/0003702814732634.
  * [4] J. K. Kauppinen, J. Partanen, Fourier Transforms in Spectroscopy, Wiley, Berlin, 2001.
  * [5] V. A. Lóenz-Fonfría, E. Padrós, Method for the estimation of the mean Lorentzian bandwidth in spectra composed of an unknown number of highly overlapped bands, Appl. Spectrosc. 62 (6) (2008) 689–700.  doi:10.1366/000370208784658129.
  * [6] M. T. Moores, K. Gracie, J. Carson, K. Faulds, D. Graham, M. Girolami, Bayesian modelling and quantification of Raman spectroscopy, arXiv preprint 1604.07299 (2016).  arXiv:1604.07299.
  * [7] T. Härkönen, L. Roininen, M. T. Moores, E. M. Vartiainen, Bayesian quantification for coherent anti-Stokes Raman scattering spectroscopy, The Journal of Physical Chemistry B 124 (32) (2020) 7005–7012.  doi:10.1021/acs.jpcb.0c04378.
  * [8] M. Diem, Modern Vibrational Spectroscopy and Micro-Spectroscopy: Theory, Instrumentation and Biomedical Applications, John Wiley & Sons, 2015.
  * [9] J. F. Kielkopf, New approximation to the Voigt function with applications to spectral-line profile analysis, J. Opt. Soc. Am. 63 (8) (1973) 987–995.  doi:10.1364/JOSA.63.000987.
  * [10] J. Olivero, R. Longbothum, Empirical fits to the Voigt line width: A brief review, Journal of Quantitative Spectroscopy and Radiative Transfer 17 (2) (1977) 233–236.  doi:https://doi.org/10.1016/0022-4073(77)90161-3.
  * [11] T. Ida, M. Ando, H. Toraya, Extended pseudo-Voigt function for approximating the Voigt profile, Journal of Applied Crystallography 33 (6) (2000) 1311–1316.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1107/S0021889800010219, doi:https://doi.org/10.1107/S0021889800010219.
  * [12] H. Haario, M. Laine, A. Mira, E. Saksman, DRAM: Efficient adaptive MCMC, Statistics and Computing 16 (4) (2006) 339–354.  doi:10.1007/s11222-006-9438-0.
  * [13] C. E. Rasmussen, C. K. I. Williams, Gaussian Processes for Machine Learning, The MIT Press, 2005.  doi:10.7551/mitpress/3206.001.0001.
  * [14] H. Susmann, Derivatives of a Gaussian process, http://herbsusmann.com/2020/07/06/gaussian-process-derivatives/, [Accessed: 2023-07-04] (2020).
  * [15] B.-K. Hsiung, T. A. Blackledge, M. D. Shawkey, Spiders do have melanin after all, Journal of Experimental Biology 218 (22) (2015) 3632–3635.  doi:10.1242/jeb.128801.
  * [16] V. Tresp, Mixtures of Gaussian processes, in: T. Leen, T. Dietterich, V. Tresp (Eds.), Advances in Neural Information Processing Systems, Vol. 13, MIT Press, 2000, pp. 654–660.
  * [17] M. M. Zhang, S. A. Williamson, Embarrassingly parallel inference for Gaussian processes, Journal of Machine Learning Research 20 (169) (2019) 1–26.
  * [18] T. Härkönen, S. Wade, K. Law, L. Roininen, Mixtures of Gaussian process experts with SMC2 (2022).  arXiv:2208.12830.
  * [19] S. Talts, M. Betancourt, D. Simpson, A. Vehtari, A. Gelman, Validating Bayesian inference algorithms with simulation-based calibration, arXiv preprint 1804.06788 (2020).  arXiv:1804.06788.
  * [20] J. McLeod, F. Simpson, Validating Gaussian process models with simulation-based calibration, in: 2021 IEEE International Conference on Artificial Intelligence Testing (AITest), 2021, pp. 101–102.  doi:10.1109/AITEST52744.2021.00028.
  * [21] T. Härkönen, E. Hannula, M. T. Moores, E. M. Vartiainen, L. Roininen, A log-Gaussian Cox process with sequential Monte Carlo for line narrowing in spectroscopy, Foundations of Data Science (2023) 0–0doi:10.3934/fods.2023008.

##  Appendix A Additional results

Figure 6: On top, a spectrum consisting of $M=10$ Gaussian line shapes in blue
along with the individual overlapping line shapes illustrated with black
dashed lines. On bottom, the obtained posterior distribution
$p(\overline{\gamma}\mid\bm{Z})$ for the area-weighted mean Lorentzian line
width. The true value for $\overline{\gamma}=0$ used to generate the data is
shown with a vertical solid line.

Figure 7: On top, a spectrum consisting of $M=6$ Voigt line shapes in blue
along with the individual overlapping line shapes illustrated with black
dashed lines. On bottom, the obtained posterior distribution
$p(\overline{\gamma}\mid\bm{Z})$ for the area-weighted mean Lorentzian line
width. The true value for $\overline{\gamma}=33.50$ used to generate the data
is shown with a vertical solid line.

Figure 8: On top, an experimental spectrum of $\beta$-carotene with selection
$[1300-1060]$ in the middle. On bottom, the obtained posterior distribution
$p(\overline{\gamma}\mid\bm{Z})$ for the area-weighted mean Lorentzian line
width.

Figure 9: On top, an experimental spectrum of $\beta$-carotene with selection
$[1360-900]$ in the middle. On bottom, the obtained posterior distribution
$p(\overline{\gamma}\mid\bm{Z})$ for the area-weighted mean Lorentzian line
width.

Figure 10: On top, an experimental spectrum of $\beta$-carotene with selection
$[1615-900]$ in the middle. On bottom, the obtained posterior distribution
$p(\overline{\gamma}\mid\bm{Z})$ for the area-weighted mean Lorentzian line
width.

Figure 11: On top, an experimental spectrum of $\beta$-carotene with selection
$[1615-900]$ in the middle. On bottom, the obtained posterior distribution
$p(\overline{\gamma}\mid\bm{Z})$ for the area-weighted mean Lorentzian line
width.

Figure 12: On top, an experimental spectrum of $\beta$-carotene with selection
$[1615-900]$ in the middle. On bottom, the obtained posterior distribution
$p(\overline{\gamma}\mid\bm{Z})$ for the area-weighted mean Lorentzian line
width.