Evaluating Community Preferences for Waste-to-Energy Development in Jakarta: An Analysis Using the Choice Experiment Method

1.   Introduction

Historically, researchers and clinicians interested in tracking alcohol consumption and metabolism in the field would require data from either a drinker's self-report or from having them use a breath alcohol analyzer. Because both methods require active participation by the subject, the data they produce are often plagued by inaccuracies. Self-report often leads to misrepresentation as (1) subjects may deviate from naturalistic behaviors due to the reporting requirement seeming unnatural, and (2) alcohol directly impairs subjects' ability to be an active participant ^[1]. Using a breath alcohol analyzer correctly requires specialized training and can produce erroneous measurements due to mouth alcohol and/or a reading based on a shallow breath by the subject. Dating back to the 1930's, ethanol, the type of alcohol in alcoholic beverages, has been known to be excreted from the human body through the skin ^[2,3,4,5]. This is due to the fact that water and ethanol are highly miscible ^[6] and the ethanol finds its way into all of the water in the body. More recently, this observation paved the way for the development of a device to measure the amount of alcohol excreted transdermally through the skin ^[8,9]. The benefits derived from such a device include the availability of near continuous measurements and the ability to collect them passively (i.e., without the active participation of the subject). This gives researchers and clinicians the potential to continuously observe naturalistic drinking behavior and patterns. There is also the possibility of making these devices available on the consumer market (e.g., wearable body system monitoring technology like Fitbits, Apple watches, etc.). In addition, the ideas we discuss here may also be applicable to the monitoring of other substances once the appropriate sensor hardware has been developed.

The challenge in using transdermal alcohol sensors is that they provide transdermal alcohol concentration (TAC), whereas alcohol researchers and clinicians have always based their studies and treatments on measurements of breath alcohol concentration (BrAC) and blood alcohol concentration (BAC). Thus, a means to reliably and accurately convert TAC to BrAC or BAC would be desirable. At levels up to approximately 0.08 (see, for example, ^[10,11]) BrAC correlates well with BAC via a simple linear relationship based on an empirical relationship known as Henry's law ^[12,13]: $\text{BAC} = \rho_{B/Br} \times \text{BrAC}$, where the constant $\rho_{B/Br}$ is known as the partition coefficient of ethanol in blood and breath.

More generally, according to Henry's law, when a liquid is in contact with a gas, the concentrations, $C_L$ and $C_G$, of a compound present in both the liquid and the gas will come to equilibrium according to the linear relationship $C_L = \rho_{L/G} C_G$, where the empirical determined constant $\rho_{L/G}$ is known as the partition coefficient for the that compound in that liquid and gas. Not surprisingly, the partition coefficient, $\rho_{L/G}$, is temperature dependent and of course its actual value will vary depending on the choice of units for $C_L$ and $C_G$. It has been shown (see, for example, ^[14]) that at $34^{\circ}$C, the partition coefficient for ethanol in blood and air is $\rho_{B/A} = 2157 \pm 9.6$ for men and $\rho_{B/A} = 2195 \pm 10.9$ for women, at $37^{\circ}$C, the partition coefficient for ethanol in blood and air is $\rho_{B/A} = 1783 \pm 8.1$ for men and $\rho_{B/A} = 1830 \pm 7.8$ for women. Using a regression model, Jones ^[14] found that at $37^{\circ}$C the partition coefficient for ethanol in water and air is $\rho_{W/A} = 2133$, in blood and air is $\rho_{B/A} = 1756$, and between plasma (all of the components of blood with the exception of the oxygen carrying red blood cells) and air is $\rho_{P/A} = 2022$. All of these values are for the case when the concentration of ethanol in air is given in units of grams per liter, and in water, blood, or plasma in units of grams per deciliter. We note that it is generally accepted that a BrAC reading of 0.08 percent alcohol corresponds to.008 grams of ethanol per 210 liters of breath and a BAC of 0.08 grams of ethanol per 100 milliliters (equal to 1 deciliter (dL) or 0.1 liters (L)) of blood.

Unfortunately, however, the correlation between TAC and BrAC/BAC, on the other hand, can vary due to a number of confounding factors. These factors include, but are not limited to, stable features of the skin like its thickness, tortuosity, and porosity, particularly as they apply to the epidermal layer of the skin, which does not have an active blood supply. Environmental factors such as ambient temperature and humidity can also affect both perspiration and vasodilation, and can thus alter skin conductance, blood flow, the amount of alcohol passing below the skin in the blood, and the amount and rate of alcohol diffusing through the skin. One would also expect there to be manufacturing and operational variations among different TAC sensors.

Earlier attempts to investigate the relationship between TAC and BrAC/BAC have used deterministic models ^{[15,16,17,18,19,20,21]}. Some utilized regression-based models ^[16], whereas others utilized first principles physics-based models that on occasion included modeling the transport of alcohol all the way from ingestion to excretion through the skin ^[22,23]. In our group's initial efforts, we modeled the transport of alcohol from the blood in the dermal layer through the epidermal layer and its eventual measurement by the sensor using a one-dimensional diffusion equation ^[15,21]. The parameters in the diffusion equation model then had to be fit or tuned (i.e., calibrated) to each individual subject, the environmental conditions, and the device through the use of simultaneous BrAC/TAC training data collected in the laboratory or clinic through a procedure known as an alcohol challenge. Once the model was fit, it could then be used to deconvolve BrAC from TAC collected in the field. This two-pass approach and the related studies were relatively successful ^{[15,19,20,21,24,25]}. However, this calibration procedure is quite burdensome to researchers, clinicians, subjects and patients, and because the models were fit to a single uni-modal drinking episode, unaccounted for variation and uncertainty in the relationship between BrAC and TAC frequently arose, making it difficult to accurately convert TAC collected in the field to BrAC ^[26,27].

More recently, to eliminate the need for calibration, deconvolution of BrAC from TAC was effected using population models fit to BrAC/TAC training data from drinking episodes across a cohort of subjects, devices, and environmental conditions ^[24,25,28]. These population models took the form of the deterministic transport models but now the parameters appearing in the model equations were considered to be random. Then in fitting the models, instead of estimating the actual values of these parameters, it was their joint distributions that were estimated. Once the models were fit, they could be used to deconvolve an estimate of the BrAC input, and by making use of the distribution of the population parameters, conservative error bands could also be generated which quantified the uncertainty in the estimated BrAC ^[24,25]. The results in these studies were based on a naive pooled data statistical model and a non-linear least squares estimator.

In this paper, we seek to build on the approach described in the previous paragraph by now using a Bayesian approach to account for the underlying uncertainty and variation in the alcohol diffusion and measurement process. We obtain posterior distributions for the transport model parameters conditioned on the training BrAC/TAC data and regularized by prior distributions based on deterministic fits. Being Bayesian based, our approach yields credible sets for the estimated parameters and what we shall refer to as conservative credible or error bands for the deconvolved estimated BrAC. What is meant by the term conservative credible band will be made precise later.

An outline of the remainder of the paper is as follows. In the next section of the paper we provide a description of our method including a derivation of a new abstract parabolic hybrid PDE/ODE model for the transdermal transport of alcohol through he epidermal layer of the skin and its capture in the vapor collection bay of the sensor. Then using linear semigroup theory we obtain an input/output model in the form of a discrete time convolution. A discussion of finite dimensional approximation and convergence issues related to the use of our model to carry out the requisite computations is also included. Then in the results section of the paper we first construct our Bayesian estimator and present two theoretical results related to it: convergence of the finite dimensional approximation and consistency. We then show how our population model based on Bayesian estimates for the random parameters can be used as part of a deconvolution scheme that yields estimated BrAC curves and conservative credible or error bands from a biosensor provided TAC signal. In this section we also present and discuss a sample of our numerical findings demonstrating the efficacy of our approach. Our numerical studies were based on human subject data collected in the Luczak laboratory in the Department of Psychology at USC. A final section contains some discussion of our theoretical and numerical results along with a few concluding remarks and avenues for possible future research.

2.   Methods

2.1.   A distributed parameter model for the transdermal transport of alcohol

As in ^[21] and ^[24], and making use of an idea recently introduced in ^[28], we model the alcohol biosensor problem described in Section 1 using a one dimensional diffusion equation to describe alcohol transport through the epidermal layer of the skin coupled with an inflow/outflow compartment model to describe the perspiration vapor collection chamber of the TAC biosensor.

The epidermal layer of the skin sits atop the dermal layer. The dermal layer has an active blood supply while the epidermal layer does not. The latter consists of both dead (the stratum corneum layer which is closest to the surface) and living (the deeper layers closer to the dermal layer) cells surrounded by interstitial fluid. Not having an active blood supply, the cells in the epidermal layer obtain nourishment primarily from $\text{O}_2$ that diffuses in from the environment beyond the skin.

The SCRAM TAC biosensor (see 1 in section 3.4 below) has a perspiration vapor collection chamber on the bottom of the sensor that sits atop, and is in direct contact with, the stratum corneum layer of the skin's epidermal layer. Perspiration in vapor form collects in the chamber. A small pump extracts a sample of the vapor from the collection chamber approximately once every 30 minutes. This sample is then electro-chemically analyzed based on an oxidation-reduction (redox) reaction in much the same way that a fuel cell produces a current (and heat and water) from hydrogen and oxygen. In the TAC sensor, ethanol molecules in the sample are oxidized producing electrons in the form of an electrical current. This current is converted into the TAC measurement based on an a priori bench calibration.

To make this more precise, we let $\Lambda$ denote the thickness of the epidermal layer (units: cm) of the skin at the location of the sensor and let $\eta$ denote the depth in the epidermal layer (units: cm), $0 \leq \eta \leq \Lambda$, $\eta = 0$ denoting the skin surface and $\eta = \Lambda$ denoting the boundary between the epidermal and dermal layers. Let $t$ denote time (units: hrs) and let $x(t, \eta)$ denote the concentration of ethanol at time $t$ and depth $\eta$ in the epidermal layer (units: grams per milliliter of interstitial fluid). Let $\tilde{w}(t)$ denote the concentration of ethanol in the TAC sensor collection chamber at time $t$ (units: grams per milliliter of air), and let $u(t)$ denote the BrAC at time $t$ (units: grams per milliliter of air). Let $y(t)$ denote the TAC at time $t$ (units: grams per milliliter of air), and let $\tilde{w}_0$ (units: grams per milliliter of air) and $\varphi_0$ (units: grams per milliliter of interstitial fluid) denote the initial conditions for $\tilde{w}$ and $x$, respectively. We will assume that there is no ethanol in either the epidermal layer or the TAC biosensor collection chamber at time $t = 0$ so $\tilde{w}_0 = 0$ and $\varphi_0 = 0$. Let $T$ denote the duration of the drinking episode (units: hrs). Then, With these definitions, our model takes the form

where $\alpha > 0$ denotes the effective diffusivity of ethanol in the interstitial fluid in the epidermal layer (units: $\text{cm}^2/\text{hr}$), $\beta > 0$ denotes the effective linear flow rate at which capillary blood plasma from the dermal layer replenishes the interstitial fluid in the epidermal layer (units: cm/hr), and $\rho_{P/A}$ denotes the partition coefficient for ethanol in plasma and air with respect to the concentration units of grams per milliliter of plasma and grams per milliliter of air at $37 ^{\circ}$C (normal body temperature).

In modeling the TAC collection chamber, we assume that the inflow of ethanol is proportional to the flux out of (i.e., from right to left) the epidermal layer at the surface of the skin (i.e., at $\eta = 0$), $\alpha \dfrac{\partial x}{\partial \eta}(t, 0)$, with constant of proportionality $\gamma$ (units: $\text{cm}^{-1}$), and the outflow is simply proportional to the concentration of ethanol in the collection chamber (i.e., a simple linear model) with constant of proportionality $\delta$ (units: $\text{hr}^{-1}$. Finally, the output gain, $\theta$, represents the bench calibration factor for the TAC sensor that converts the concentration or ethanol in the collection chamber into a TAC (units: dimensionless).

Since the thickness of the epidermal layer, $\Lambda$, is in general difficult to measure and can be mathematically difficult to estimate computationally due to it determining the spatial domain of the diffusion equation, it is desirable to transform the system eq. $(2.1)$ to a domain of fixed length, $\Lambda = 1$. We make the standard change of variable $\eta \mapsto \frac{\eta}{\Lambda}$ thus rendering $\eta$ dimensionless. For $t \geq 0$, We also set $w(t) = \rho_{P/A}\tilde{w}(t)$. Then, recalling our assumption of zero initial conditions, the following hybrid ordinary-partial differential equation input/output system results

where $q_1 = \frac{\alpha}{\Lambda^2}$, $q_2 = \frac{\theta\beta\rho_{P/A}}{\Lambda}$, $q_3 = \frac{\gamma \alpha}{\Lambda}$, and $q_4 = \delta$. We note that since the only observable and observed quantities are BrAC, $u$, and TAC, $y$, the physiological interpretations of the variables and parameters in between that define our model in the form of an input/output map from BrAC to TAC are of little interest to us. Although we have relied on first principles modeling to derive the system of equations given in eq. $(2.2)$, our motivation was not to gain a deeper understanding of the transdermal transport of ethanol. Rather, it was to be able to keep the dimension of the space of unknown parameters as low as possible by capturing the underlying physics and physiology of the transport process, albeit in a greatly simplified form. Indeed, our primary objective here is to first fit the parameters (or, more precisely, their distributions) in the model to observed input/output BrAC/TAC training pairs and to then use the resulting population model to obtain an estimate for the BrAC and associated error bars corresponding to a given TAC signal collected in the field from a member of the cohort or population that provided the data which was used to train the model.

Let $\mathit{\boldsymbol{q}} = [q_1, q_2]^T$ denote the unknown, un-measurable, and, in general, subject-dependent physiological parameters. The parameters $q_3$ and $q_4$ are device (i.e., hardware) dependent parameters and as such, we consider them to be bench-measurable empirically in the lab. We do note however, with simple changes of variable, the theory and methods we develop below apply, and their distributions could also be estimated along with those of $q_1$ and $q_2$ with the same techniques we use here to estimate the distributions of $q_1$ and $q_2$. In addition, $q_3$ and $q_4$ could also be estimated using a deterministic scheme such as a regularized nonlinear least squares approach. For clarity and ease of exposition, we will focus our attention here on the development of a population model for a cohort of subjects by estimating the distribution of the un-measurable physiological parameters $q_1$ and $q_2$.

We consider this system on a finite-time horizon, $0 \leq T < \infty$, and we assume zero-order hold input, $u(t) = u_k$, $t \in [k\tau, (k+1)\tau)$, $k = 0, 1, 2, ...$, where $\tau$ denotes the sampling time of the biosensor. We set $x_k = x_k(\eta) = x(k\tau, \eta)$, $w_k = w(k\tau)$, and $y_k = y(k\tau)$, $k = 0, 1, \ldots, K$, where we assume $T = K\tau$. For $k = 0, 1, 2, ...$ we consider the system eq. $(2.2)$ on the interval $[k\tau, (k+1)\tau]$ and make the change of variable: $v(t, \eta) = x(t, \eta) - \xi(\eta)u_k$ where $\xi(\eta) = \frac{q_2}{q_1}\eta$. It is then easily verified that $w$ and $v$ satisfy the following hybrid system

with initial conditions $v(k\tau, \cdot) = x(k\tau, \cdot) - \xi(\cdot)u_k = x_k - \xi u_k$ and $w(k\tau) = w_k$ on $[0, 1]$.

We then use linear semigroup theory to rewrite the system eq. $(2.3)$ in state space form in an appropriately chosen Hilbert space and subsequently obtain a discrete time evolution system for $(w_k, x_k)$, $k = 0, 1, 2, ....K$ which is equivalent to eq. $(2.2)$. Let $Q$ denote a compact subset of the positive orthant of $\mathbb{R}^4$, and for $\mathit{\boldsymbol{q}} = [q_1, q_2, q_3, q_4]^T \in Q$ we define the Hilbert spaces

with respective corresponding inner products $\langle (\theta_1, \psi_1), (\theta_2, \psi_2) \rangle_{ \mathit{\boldsymbol{q}}} = (q_1/q_3)\theta_1 \theta_2 + \int_0^1 \psi_1(\eta)\psi_2(\eta)d\eta$ and $\langle (\psi_1(0), \psi_1), (\psi_2(0), \psi_2) \rangle_V = \psi_1(0) \psi_2(0) + \int_0^1 \psi_1'(\eta)\psi_2'(\eta)d\eta$, and respective norms $\vert \cdot\vert_{ \mathit{\boldsymbol{q}}}$, and $\lVert \cdot\rVert_V$. Note that the Sobolev Embedding Theorem ^[29] yields that the norm induced by the $V$ inner product is equivalent to the standard $H^1$ norm on $V$. It is not difficult to show that $V$ is densely and continuously embedded in $H_{ \mathit{\boldsymbol{q}}}$ and that we have the Gelfand triple of dense and continuous embeddings $V \hookrightarrow H_{ \mathit{\boldsymbol{q}}} \hookrightarrow V^*$.

Then based on the weak formulation of the system eq. $(2.3)$, for each $\mathit{\boldsymbol{q}} \in Q$ define the bilinear form $a(\mathit{\boldsymbol{q}}, \cdot, \cdot): V \times V \rightarrow \mathbb{R}$ by $a(\mathit{\boldsymbol{q}}, \hat{\varphi}, \hat{\psi}) = (q_1 q_4/q_3) \varphi(0)\psi(0) + q_1 \int_0^1 \varphi'(\eta)\psi'(\eta)d\eta$ for $\hat{\varphi}, \hat{\psi} \in V$, where $\hat{\varphi} = (\varphi(0), \varphi)$, $\hat{\psi} = (\psi(0), \psi)$, and $\mathit{\boldsymbol{q}} = [q_1, q_2, q_3, q_4]^T \in Q$. Standard arguments can be used to argue that the form $a(\mathit{\boldsymbol{q}}, \cdot, \cdot)$ satisfies the following three properties.

1. Boundedness There exists a constant $\alpha_0 > 0$ such that $\vert a(\mathit{\boldsymbol{q}}, \hat{\psi}_1, \hat{\psi}_2)\vert \leq \alpha_0 \lVert \hat{\psi}_1 \rVert_V \lVert \hat{\psi}_2 \rVert_V,$ $\hat{\psi}_1, \hat{\psi}_2 \in V$,

2. Coercivity There exists constants $\lambda_0 \in \mathbb{R}$ and $\mu_0 > 0$ such that $a(\mathit{\boldsymbol{q}}, \hat{\psi}, \hat{\psi}) + \lambda_0\vert\hat{\psi}\vert_{ \mathit{\boldsymbol{q}}}^2 \geq \mu_0\lVert \hat{\psi} \rVert_V^2, \quad \hat{\psi}\in V$.

3. Continuity For all $\hat{\psi}_1, \hat{\psi}_2 \in V$, we have that $\mathit{\boldsymbol{q}} \mapsto a(\mathit{\boldsymbol{q}}, \hat{\psi}_1, \hat{\psi}_2)$ is a continuous mapping from $Q$ into $\mathbb{R}$.

Note that in (1)–(3) above, $Q$ compact implies that the constants $\alpha_0$, $\lambda_0$, and $\mu_0$ may all be chosen independent of $\mathit{\boldsymbol{q}} \in Q$. Furthermore, properties (1)–(3) immediately yield that the form $a(\mathit{\boldsymbol{q}}, \cdot, \cdot)$ defines a bounded, elliptic (i.e., $\lambda_0 = 0)$ operator $\mathcal{A}(\mathit{\boldsymbol{q}}) \in \mathcal{L}(V, V^\ast)$ given by $\langle \mathcal{A}(\mathit{\boldsymbol{q}})\hat{\varphi}, \hat{\psi} \rangle = \langle \mathcal{A}(\mathit{\boldsymbol{q}})(\varphi(0), \varphi), (\psi(0), \psi) \rangle = -a(\mathit{\boldsymbol{q}}, \varphi, \psi)$ for $\hat{\varphi} = (\varphi(0), \varphi), \hat{\psi} = (\psi(0), \psi) \in V$. If we define the set $\mathcal{D} = \{\hat{\varphi} \in V: \mathcal{A}(\mathit{\boldsymbol{q}})\hat{\varphi} \in H_{ \mathit{\boldsymbol{q}}}\} = \{ (\varphi(0), \varphi) \in V : \varphi \in H^2(0, 1), \varphi'(1) = 0\}$ which is independent of $\mathit{\boldsymbol{q}}$ for $\mathit{\boldsymbol{q}} \in Q$, we obtain the closed, densely defined linear operator $A(\mathit{\boldsymbol{q}}):\mathcal{D} \subset H_{ \mathit{\boldsymbol{q}}} \rightarrow H_{ \mathit{\boldsymbol{q}}}$ given by $A(\mathit{\boldsymbol{q}})\hat{\varphi} = A(\mathit{\boldsymbol{q}})(\varphi(0), \varphi) = (q_3\varphi'(0) - q_4\varphi(0), q_1\varphi"), \quad \hat{\varphi} = (\varphi(0), \varphi) \in \mathcal{D}$. The operator $A(\mathit{\boldsymbol{q}}):\mathcal{D} \subset H_{ \mathit{\boldsymbol{q}}} \rightarrow H_{ \mathit{\boldsymbol{q}}}$ is regularly dissipative and (see, for example, ^[30]) is the infinitesimal generator of a holomorphic semigroup of bounded linear operators $\{ e^{A(\mathit{\boldsymbol{q}})t} : t \geq 0 \}$ on $H_{ \mathit{\boldsymbol{q}}}$ and $V^\ast$. Moreover, the system eq. $(2.3)$ then has a state space form where $\frac{d\hat{v}}{dt}(t) = A(\mathit{\boldsymbol{q}})\hat{v}(t) + (\frac{q_3q_2}{q_1}, 0)u_k$ and $\hat{v}(k\tau) = (w_k, x_k - \xi u_k)$ for $k\tau < t < (k+1)\tau$. Then for time step $\tau > 0$ and $k = 0, 1, \ldots, K$, letting $\hat{x}_k = (w_k, x_k)$, $\hat{A}(\mathit{\boldsymbol{q}}) = e^{A(\mathit{\boldsymbol{q}})\tau} \in \mathcal{L}(H_{ \mathit{\boldsymbol{q}}}, H_{ \mathit{\boldsymbol{q}}})$, and $\hat{B}(\mathit{\boldsymbol{q}}) = (I - \hat{A}(\mathit{\boldsymbol{q}}))\big\{(0, \xi) - A(\mathit{\boldsymbol{q}})^{-1}(\dfrac{q_3q_2}{q_1}, 0)\big\} \in \mathcal{L}(\mathbb{R}, H_{ \mathit{\boldsymbol{q}}})$, it follows that

where for $\hat{B}(\mathit{\boldsymbol{q}})$ we have used the fact that the operator $A(\mathit{\boldsymbol{q}})^{-1}$ commutes with the semigroup generated by $A(\mathit{\boldsymbol{q}})$. Note that the operator $\hat{B}(\mathit{\boldsymbol{q}})$ is in fact an element in $H_{ \mathit{\boldsymbol{q}}}$ and that $(0, \xi) \in V$, but that $(\dfrac{q_3q_2}{q_1}, 0)$ is only an element in $H_{ \mathit{\boldsymbol{q}}}$. From eq. $(2.5)$ the state space form of our discrete time model is

where $\hat{x}_k = (w_k, x_k)$, $k = 0, 1, 2, \ldots$ and the operator $\hat{C} \in \mathcal{L}(H_{ \mathit{\boldsymbol{q}}}, \mathbb{R})$ is given by $\hat{C}(\theta, \psi) = \theta$, where $(\theta, \psi) \in H_{ \mathit{\boldsymbol{q}}}$. Note that the ellipticity of $\mathcal{A}(\mathit{\boldsymbol{q}})$ guarantees the existence of $A(\mathit{\boldsymbol{q}})^{-1}$.

From eq. $(2.6)$ it is immediately clear that if we assume zero initial conditions, $\hat{x} = (w_0, \varphi_0) = (0, 0)$, the output $y$ can be written as a discrete time convolution of the input, $u$, with a filter, $h(\mathit{\boldsymbol{q}})$, as

where for $\mathit{\boldsymbol{q}} \in Q$, $h_i(\mathit{\boldsymbol{q}}) = \hat{C}\hat{A}(\mathit{\boldsymbol{q}})^i\hat{B}(\mathit{\boldsymbol{q}})$, $i = 0, 1, 2, \ldots$. Using the Trotter-Kato semigroup approximation theorem (see, for example, ^[31]), the following result can be shown (see, for proof, ^[32]).

Lemma 2.1. For $Q$ a compact subset of the positive orthant or $\mathbb{R}^4$, $K = K\tau$ for constant time step $\tau > 0$, and $h_i$ as defined in eq. $(2.7)$, we have that the mapping $\mathit{\boldsymbol{q}} \mapsto h_i(\mathit{\boldsymbol{q}})$ from $Q$ into $\mathbb{R}$ is continuous, uniformly in $\mathit{\boldsymbol{q}}$ and $i$, for $\mathit{\boldsymbol{q}} \in Q$ and $i \in \{0, 1, 2, \ldots, K\}$.

Now although the input/output model given in eq. $(2.7)$ is a standard convolution in $\mathbb{R}$, the filter, $\{h_k(q)\}$ involves the semigroup $\{ e^{A(\mathit{\boldsymbol{q}})t} : t \geq 0 \}$ which is defined on the infinite dimensional Hilbert space $H_{ \mathit{\boldsymbol{q}}}$. Consequently, finite dimensional approximation is required. For $n = 1, 2, \ldots$ let $\{ \varphi^n_i (\eta) \}_{i = 0}^n$ denote the standard linear B-splines on the interval $[0, 1]$ defined with respect to the uniform mesh $\{0, \frac{1}{n}, \frac{2}{n}, \ldots, \frac{n-1}{n}, 1\}$, $\varphi_i^n(\eta) = (n\eta - i +1)\mathbb{1}_{[\frac{i-1}{n}, \frac{i}{n}]} + (1 - n\eta +i)\mathbb{1}_{[\frac{i}{n}, \frac{i+1}{n}]}$. Set

and let $\mathcal{P}^n_{ \mathit{\boldsymbol{q}}}:H_{ \mathit{\boldsymbol{q}}} \rightarrow V^n$ denote the orthogonal projection of $H_{ \mathit{\boldsymbol{q}}}$ on to $V^n$ along $(V^n)^\perp$. Standard arguments from the theory of splines (see, for example, ^[33]) can be used to argue that $\vert\mathcal{P}^n_{ \mathit{\boldsymbol{q}}}(\theta, \psi) - (\theta, \psi)\vert_{ \mathit{\boldsymbol{q}}} \rightarrow 0$, as $n \rightarrow \infty$, for all $(\theta, \psi) \in H_{ \mathit{\boldsymbol{q}}}$, and that $\lVert\mathcal{P}^n_{ \mathit{\boldsymbol{q}}}\hat{\varphi} - \hat{\varphi}\rVert_V \rightarrow 0$, as $n \rightarrow \infty$, for all $\hat{\varphi} \in V$ with the convergence uniform in $\mathit{\boldsymbol{q}}$ for $\mathit{\boldsymbol{q}} \in Q$.

For $n = 1, 2, ...$ and $k = 0, 1, 2, ...$ we set $\hat{x}_k^n(\eta) = \sum\limits_{i = 0}^n X^{n, k}_i\hat{\varphi}_i^n(\eta),$ and we approximate the operator $A(\mathit{\boldsymbol{q}})$ using a Galerkin approach. That is, we define the operator $A^n(\mathit{\boldsymbol{q}}) \in \mathcal{L}(V^n, V^n)$ by restricting the form $a(\mathit{\boldsymbol{q}}, \cdot, \cdot)$ to $V^n\times V^n$. We then set

where $\hat{B}^n(\mathit{\boldsymbol{q}}) \in \mathcal{L}(\mathbb{R}, V^n) = V^n$. The matrix representations for these operators with respect to the basis $\{\hat{\varphi}^n_i\}_{i = 0}^n$ are then given by $[A^n(\mathit{\boldsymbol{q}})] = -[M^n(\mathit{\boldsymbol{q}})]^{-1}K^n(\mathit{\boldsymbol{q}})$, $[\hat{A}^n(\mathit{\boldsymbol{q}})] = e^{-[M^n(\mathit{\boldsymbol{q}})]^{-1}K^n(\mathit{\boldsymbol{q}})\tau}$, and $[\hat{B}^n(\mathit{\boldsymbol{q}})] = (I - [\hat{A}^n(\mathit{\boldsymbol{q}})])\big\{\Xi^n - [A^n(\mathit{\boldsymbol{q}})]^{-1}[M^n(\mathit{\boldsymbol{q}})]^{-1}[q_2, 0, 0, \ldots, 0]^T\big\} = (I - [\hat{A}^n(\mathit{\boldsymbol{q}})])\big\{\Xi^n + [K^n(\mathit{\boldsymbol{q}})]^{-1}[q_2, 0, 0, \ldots, 0]^T\big\}$, where $[M]_{i, j = 0}^n = \langle (\varphi_i^n(0), \varphi_i^n), (\varphi_j^n(0), \varphi_j^n) \rangle_{ \mathit{\boldsymbol{q}}}$, $[K]_{i, j = 0}^n = a(\mathit{\boldsymbol{q}}, \varphi_i^n(0), \varphi_i^n), (\varphi_j^n(0), \varphi_j^n))$, and $\Xi^n = \frac{q_2}{q_1}[0, \frac{1}{n}, \frac{2}{n}, \ldots, \frac{n-1}{n}, 1]^T$. Letting $[\hat{C}^n] = [1, 0, 0, \ldots, 0] \in \mathbb{R}^{1 \times(n+1)}$, we consider the discrete time dynamical system in $V^n$ given by

or equivalently in $\mathbb{R}^{n+1}$ given by the system where $X^{n, k+1} = [\hat{A}^n(\mathit{\boldsymbol{q}})] X^{n, k} + [\hat{B}^n(\mathit{\boldsymbol{q}})]u_k$, $y^n_k = [\hat{C}^n]X^{n, k}$, and $X^{n, 0} = [0, 0, \ldots, 0]^T \in \mathbb{R}^{n+1}$, we obtain that

for $k = 0, 1, 2, \ldots, K$ where $h^n_i(\mathit{\boldsymbol{q}}) = [\hat{C}^n][\hat{A}^n(\mathit{\boldsymbol{q}})]^{i}[\hat{B}^n(\mathit{\boldsymbol{q}})] \in \mathbb{R}$, $i = 0, 1, 2, \ldots, K-1$.

Using linear semigroup theory (see, for example, ^[21,34,35]) and in particular the Trotter-Kato semigroup approximation theorem (see, for example, ^[36] and ^[31]) the following results can be established (for proof, see ^[32]).

Theorem 2.1. For $Q$ a compact subset of the positive orthant of $\mathbb{R}^4$, $n = 1, 2, \ldots$, $\{ \varphi^n_i (\eta) \}_{i = 0}^n$ the standard linear B-splines on the interval $[0, 1]$ defined with respect to the uniform mesh $\{0, \frac{1}{n}, \frac{2}{n}, \ldots, \frac{n-1}{n}, 1\}$, $V^n = \mathit{\text{span}}\{(\varphi^n_i(0), \varphi^n_i) \}^n_i$, $\mathcal{P}^n_{ \mathit{\boldsymbol{q}}}:H_{ \mathit{\boldsymbol{q}}} \rightarrow V^n$ the orthogonal projection of $H_{ \mathit{\boldsymbol{q}}}$ on to $V^n$ along $(V^n)^\perp$, and $\hat{A}^n(\mathit{\boldsymbol{q}})$ and $\hat{B}^n(\mathit{\boldsymbol{q}})$ defined as in eq. $(2.9)$, we have that $\vert\hat{A}^n(\mathit{\boldsymbol{q}})\mathcal{P}^n_{ \mathit{\boldsymbol{q}}}(\theta, \psi) - \hat{A}(\mathit{\boldsymbol{q}})(\theta, \psi)\vert_{ \mathit{\boldsymbol{q}}} \rightarrow 0$, as $n \rightarrow \infty$, for all $(\theta, \psi) \in H_{ \mathit{\boldsymbol{q}}}$, that $\lVert\hat{A}^n(\mathit{\boldsymbol{q}})\mathcal{P}^n_{ \mathit{\boldsymbol{q}}}\hat{\varphi} - \hat{A}(\mathit{\boldsymbol{q}})\hat{\varphi}\rVert_V \rightarrow 0$, as $n \rightarrow \infty$, for all $\hat{\varphi} \in V$, and that $\lVert\hat{B}^n(\mathit{\boldsymbol{q}}) - \hat{B}(\mathit{\boldsymbol{q}})\rVert_V \rightarrow 0$, as $n \rightarrow \infty$, with the convergence in all cases uniform in $\mathit{\boldsymbol{q}}$ for $\mathit{\boldsymbol{q}} \in Q$.

Theorem 2.2. Under the same hypotheses as Theorem 2.1, we have that $\lVert\hat{x}^n_k(\mathit{\boldsymbol{q}}) - \hat{x}_k(\mathit{\boldsymbol{q}})\rVert_V \rightarrow 0$, that $\vert\hat{x}^n_k(\mathit{\boldsymbol{q}}) - \hat{x}_k(\mathit{\boldsymbol{q}})\vert_{ \mathit{\boldsymbol{q}}} \rightarrow 0$, that $\vert y^n_k - y_k\vert \rightarrow 0$, and that $\vert h^n_k(\mathit{\boldsymbol{q}}) - h_k(\mathit{\boldsymbol{q}})\vert \rightarrow 0$ as $n \rightarrow \infty$ uniformly in $k$ for $k \in \{0, 1, 2, ..., K\}$ and uniformly in $\mathit{\boldsymbol{q}}$ for $\mathit{\boldsymbol{q}} \in Q$.

Finally, we will assume that we have training data, $\big\{\{u_k^i\}_{k = 0}^{K}, \{y_k^i\}_{k = 0}^K\big\}_{i = 1}^R$, from $R$ participants or subjects where without loss of generality (i.e., by padding with zeros) we have assumed that all training input/output datasets have the same number, $K$, of observations. In this case, for $i = 1, \ldots, R$ we have,

where $h_j(\mathit{\boldsymbol{q}}) = \hat{C}\hat{A}(\mathit{\boldsymbol{q}})^{j}\hat{B}(\mathit{\boldsymbol{q}}) \in \mathbb{R}$ and $h^n_j(\mathit{\boldsymbol{q}}) = [\hat{C}^n][\hat{A}^n(\mathit{\boldsymbol{q}})]^{j}[\hat{B}^n(\mathit{\boldsymbol{q}})] \in \mathbb{R}$, for $j = 0, 1, \ldots, K-1$. This formulation facilitates the estimation of the population parameters $\mathit{\boldsymbol{q}}$. If one wishes to find the parameters $\mathit{\boldsymbol{q}}$ for a specific individual, the methods outlined in Section 2.2 can still be applied by letting the indices $i = 1, \ldots, R$ refer to different measured BrAC/TAC events each with $k = 0, \ldots, K$ denoting the measurement times for the desired individual subject.

2.2.   Bayesian estimation of dynamical system parameters

In this section we develop a Bayesian framework to estimate the unknown parameters $\mathit{\boldsymbol{q}} = [q_1, q_2]^T$ in the system eq. $(2.2)$. To illustrate our approach, for simplicity but without loss of generality, we have assumed that the sensor parameters $q_3$ and $q_4$ have been bench-measured and are therefore known and concentrate our effort on estimating the physiological subject-dependent parameters $q_1$ and $q_2$. All of what follows below can easily be extended to estimating all four of the parameters in the model. Our underlying statistical model incorporating noise is based on the observation of $\{y^{i}_{j}\}$ as in eq. $(2.12)$ and is given by

where $V^i_j$ are our measured TAC values, and $\varepsilon^i_{j}$ are the i.i.d. noise terms corresponding to person $i$ at time $j\tau$ with $\sigma > 0, \tau > 0$. Commonly, as we will assume in Section 3.2 and beyond, $\varepsilon^i_{j} \sim N(0, \sigma^2)$. In order to be able to carry out the requisite computations, we consider the approximating statistical model based on eq. $(2.12)$

where once again the $V^i_j$'s are assumed to be the measured TAC values. We consider $\mathit{\boldsymbol{\mathcal{q}}}$ to be a random vector on some probability space $\{\Omega, \Sigma, P\}$ with support $Q$ and assume that the prior distribution of $\mathit{\boldsymbol{\mathcal{q}}}$ is given by the push forward measure $\pi_0$. That is for $A\subset Q$, $P(\{ \mathit{\boldsymbol{\mathcal{q}}} \in A\}) = P(\mathit{\boldsymbol{\mathcal{q}}}^{-1}(A)) = \int_{ \mathit{\boldsymbol{\mathcal{q}}}^{-1}(A)}dP = \int_Ad\pi_0(\mathit{\boldsymbol{q}})$.

We assume independence across both $i$ (individuals) and $j$ (sampling times for each individual), for each $i$ and $j$ we have $V^{i}_{j} - y^{i}_{j} = \varepsilon^i_{j}$ (commonly distributed $N(0, \sigma^2)$) and similarly, $V^{i}_{j} - y^{n, i}_{j} = \varepsilon^i_{j}$ (again commonly distributed $N(0, \sigma^2)$). Letting $\varphi$ denote the density of $\varepsilon^i_j$'s, for $\mathit{\boldsymbol{q}} \in Q$ the likelihood and the approximating likelihood functions are given respectively by (see, for example, ^{[37,38,39,40]})

An application of Bayes' Theorem (see, for example, Theorem 1.31 in ^[41]) yields that the posterior distribution of $\mathit{\boldsymbol{\mathcal{q}}}$ or the conditional distribution of $\mathit{\boldsymbol{\mathcal{q}}}$ conditioned on the data, $\{V^{i}_{j}\}$, is a push forward measure $\pi = \pi(\cdot|\{V^{i}_{j}\})$ that is absolutely continuous with respect to $\pi_0$ and whose Radon-Nikodym derivative, or density, for $\mathit{\boldsymbol{q}} \in Q$ is given by

In this way, for $A \subset Q$, we have $P(\mathit{\boldsymbol{\mathcal{q}}} \in A|\{V^{i}_{j}\}) = \int_Ad\pi(\mathit{\boldsymbol{q}}) = \int_Ad\pi(\mathit{\boldsymbol{q}}|\{V^{i}_{j}\}).$ If, in addition, we have $\pi_0\ll \lambda$ with density $\dfrac{d\pi_0}{d\lambda} = f_0$ where $\lambda$ denotes Lebesgue measure on $Q$, then $\pi\ll \lambda$ with conditional density $f$ given by

Analogously, in the case of the approximating likelihood eq. $(2.15)$, eq. $(2.16)$, eq. $(2.17)$, eq. $(2.18)$, and eq. $(2.19)$ respectively become

$Z^n = \int\limits_Q L^n(\mathit{\boldsymbol{q}}|\{V^{i}_{j}\}) f_0(\mathit{\boldsymbol{q}})d\lambda(\mathit{\boldsymbol{q}}) = \int\limits_Q \prod\limits_{i = 1}^{R}\prod\limits_{j = 1}^{K}\varphi(V^{i}_{j} - y^{n, i}_{j}(\mathit{\boldsymbol{q}}))f_0(\mathit{\boldsymbol{q}})d\lambda(\mathit{\boldsymbol{q}})$, and $P(\mathit{\boldsymbol{\mathcal{q}}} \in A|\{V^{i}_{j}\}) = \int_A f^n(\mathit{\boldsymbol{q}})d\lambda(\mathit{\boldsymbol{q}}) = \int_A f^n(\mathit{\boldsymbol{q}}|\{V^{i}_{j}\})d\lambda(\mathit{\boldsymbol{q}})$.

3.   Results

3.1.   Convergence in distribution

Consider the random variable $\mathit{\boldsymbol{\mathcal{q}}}$ with posterior distribution described by the measure $\pi$ given in eq. $(2.15)$ and eq. $(2.16)$, and let $\mathit{\boldsymbol{\mathcal{q}}}^n$ denote the random variable $\mathit{\boldsymbol{\mathcal{q}}}$ but with posterior distribution $\pi^n$ given by eq. $(2.20)$ and eq. $(2.21)$. In this section we establish that $\mathit{\boldsymbol{\mathcal{q}}}^n \xrightarrow{\text{Dist}} \mathit{\boldsymbol{\mathcal{q}}}$ as $n \rightarrow \infty$; that is that $\mathit{\boldsymbol{\mathcal{q}}}^n$ converges in distribution to $\mathit{\boldsymbol{\mathcal{q}}}$. Recall that due to the physical constraints based on our model for the alcohol biosensor problem, eq. $(2.2)$, we require that the parameters $\mathit{\boldsymbol{q}}$ lie in the interior of the positive orthant of $\mathbb{R}^2$.

Theorem 3.1. For $Q$ a compact set in the interior of the positive orthant of $\mathbb{R}^2$, a prior $\pi_0$ with compact support $Q$ and a density that is continuous on $Q$, and a noise distribution with bounded density function $\varphi$ and support on $\mathbb{R}$, $\mathit{\boldsymbol{\mathcal{q}}}^n$, the random variable with posterior distribution $\pi^n$ given by eq. $(2.20)$ and eq. $(2.21)$ converges in distribution to the random variable $\mathit{\boldsymbol{\mathcal{q}}}$ with posterior distribution $\pi$ given by eq. $(2.15)$ and eq. $(2.16)$.

Proof. For $S$ a subset $Q$, the triangle inequality yields

where $\varphi$ is the normal density describing the distribution of the noise term in eq. $(2.13)$, and $Z$ and $Z^n$ are as in eq. $(2.16)$ and eq. $(2.21)$, respectively.

Focusing first on the limit of $\vert 1/Z - 1/Z^n \vert$ as $n\rightarrow \infty$, by Lemma 2.1 we have that the $y_i^j(\mathit{\boldsymbol{q}})$ are continuous in $\mathit{\boldsymbol{q}}$ for $\mathit{\boldsymbol{q}}\in Q$, $i \in \{0, 1, \ldots, R \}$, and $j \in \{0, 1, \ldots, K \}$. Since $Q$ is compact, the $\{y_j^i(\mathit{\boldsymbol{q}})\}$ are bounded and thus $0 < Z < \infty$. By Theorem 2.2, since $y_j^{n, i}(\mathit{\boldsymbol{q}}) \rightarrow y_j^i(\mathit{\boldsymbol{q}})$ uniformly in $\mathit{\boldsymbol{q}}$ for $\mathit{\boldsymbol{q}} \in Q$ as $n \rightarrow \infty$, $0 < Z^n < \infty$ for $n$ large enough. Again by Theorem 2.2 it follows from eq. $(2.16)$, eq. $(2.21)$, and the Bounded Convergence Theorem that $Z^n \rightarrow Z$ as $n \rightarrow \infty$ and therefore that $\left|1/Z - 1/Z^n \right| \rightarrow 0$ as $n \rightarrow \infty$. Then, essentially the same arguments yield that $\int_S | L(\mathit{\boldsymbol{q}}|\{V^{i}_{j}\}) - L^n(\mathit{\boldsymbol{q}}|\{V^{i}_{j}\}) | d\pi_0(\mathit{\boldsymbol{q}}) \rightarrow 0,$ from which it then immediately follows that $\left(\dfrac{1}{Z^n}\right)\int_S | L(\mathit{\boldsymbol{q}}|\{V^{i}_{j}\}) - L^n(\mathit{\boldsymbol{q}}|\{V^{i}_{j}\}) | d\pi_0(\mathit{\boldsymbol{q}}) \rightarrow 0,$ and therefore from eq. $(3.1)$ that $\mathit{\boldsymbol{\mathcal{q}}}^n \xrightarrow{\text{Dist}} \mathit{\boldsymbol{\mathcal{q}}}$ as $n \rightarrow \infty$ and the theorem has been proved.

For the push forward measures $\pi$ and $\pi^n$ from eq. $(2.15)$ and eq. $(2.20)$, respectively, we are commonly interested in their respective expected values and the convergences there within. Since $Q$ is compact, a direct invocation of the Portmanteau theorem (see, ^[41]) establishes the following corollary which guarantees the convergence of all moments described by $\pi^n$ to those of $\pi$.

Corollary 3.1. Under the hypotheses of Theorem 3.1, and for any continuous function $g: Q \rightarrow \mathbb{R}$ we have that $\mathbb{E}_{\pi^n}[g(\mathit{\boldsymbol{\mathcal{q}}})] = \int\limits_Q g(\mathit{\boldsymbol{q}}) d\pi^n(\mathit{\boldsymbol{q}}) \rightarrow \int\limits_Q g(\mathit{\boldsymbol{q}}) d\pi(\mathit{\boldsymbol{q}}) = \mathbb{E}_{\pi}[g(\mathit{\boldsymbol{\mathcal{q}}})]$ as $n\rightarrow \infty$.

3.2.   Consistency

In this section we demonstrate the strong consistency of the posterior distribution with respect to the parameters, $\mathit{\boldsymbol{q}}$, by imposing stronger assumptions on the distribution of the noise terms $\varepsilon^i_j$ in eq. $(2.13)$ and on the prior, $\pi_0$, by restricting $Q$ to a rectangle in the positive orthant of $\mathbb{R}^2$, and by applying the framework summarized in ^[42].

As in ^[42], we show consistency of the posterior distribution $\pi$ as in eq. $(2.17)$, rather than consistency of a point estimator based on the posterior distribution. As such, for prior $\pi_0$ over $Q$, posterior $\pi(\cdot\vert\{V_j^i\})$ as in eq. $(2.15)$, and i.i.d. noise $\varepsilon_j^i \sim N(0, \sigma^2)$ for $\sigma > 0$, we also consider that for $\mathit{\boldsymbol{q}} \in Q$ assumed known we have random variables $\mathcal{V}_j^i \sim N(y_j^i(\mathit{\boldsymbol{q}}), \sigma^2)$ as determined by eq. $(2.13)$ for $i = 1, 2, \ldots, R$ and $j = 1, 2, \ldots, K$. Further, we have that $\{ \mathcal{V}_j^i\}_{i, j}$ are independent in $i$ and $j$, but are non-identically distributed (i.n.i.d).

For clarity and brevity we consider the random vector $\mathit{\boldsymbol{\mathcal{V}}}^i = (\mathcal{V}^i_0, \mathcal{V}^i_1, \ldots, \mathcal{V}^i_K)^T$ with values in $\mathbb{R}^{K+1}$ and independent entries derived from the matrix equivalent to eq. $(2.13)$ and eq. $(2.14)$, namely $\mathit{\boldsymbol{V}}^i = \mathit{\boldsymbol{y}}^i +\mathit{\boldsymbol{\varepsilon}}^i = H\mathit{\boldsymbol{u}}^i + \mathit{\boldsymbol{\varepsilon}}^i$ and $\mathit{\boldsymbol{V}}^{i} = \mathit{\boldsymbol{y}}^{n, i} +\mathit{\boldsymbol{\varepsilon}}^i = H^n\mathit{\boldsymbol{u}}^i + \mathit{\boldsymbol{\varepsilon}}^i$, with noise vectors $\mathit{\boldsymbol{\varepsilon}}^i = (\varepsilon^i_0, \varepsilon^i_1, \ldots, \varepsilon^i_K)^T$, observed TAC vectors $\mathit{\boldsymbol{V}}^i = (V^i_0, V^i_1, \ldots, V^i_K)^T$, theoretical TAC vectors $\mathit{\boldsymbol{y}}^i = (y^i_0, j^i_1, \ldots, y^i_K)^T$ and analogous $\mathit{\boldsymbol{y}}^{n, i}$, BrAC data vectors $\mathit{\boldsymbol{u}}^i = (u^i_0, u^i_1, \ldots, u^i_K)^T$ and analogous $\mathit{\boldsymbol{u}}^{n, i}$, and kernel matrices $H$ and $H^n$ with entries $[H]_{i, j} = h_{i-j}\mathbb{1}_{j\leq i}$ and $[H^n]_{i, j} = h^n_{i-j}\mathbb{1}_{j\leq i}$, respectively. In this way, for every $i\in \{ 1, 2, \ldots, R \}$, by independence in $j$ we have the family of joint distributions, $\{ f_{i, \mathit{\boldsymbol{q}}}(\mathit{\boldsymbol{\cdot}}) = \prod_{j = 0}^K \varphi(\cdot_j - y^i_j(\mathit{\boldsymbol{q}})) : \mathit{\boldsymbol{q}} \in Q \}$, representing the possible densities of $\mathit{\boldsymbol{\mathcal{V}}}^i$ for $\varphi$ the noise density and $y_j^i$ as in eq. $(2.12)$. We are interested in the scenario where the number of subjects $R\rightarrow \infty$.

With our reframing, for $A \subset Q$ by independence in $i$ we may rewrite eq. $(2.15)$ as

where for all $i$ we have data vectors $\mathit{\boldsymbol{V}}^i$, and for our purposes we will be interested in the equivalent form on the right-hand side of the equation above where $\mathit{\boldsymbol{q}}_0\in Q$ is the true value of our parameters $[q_1, q_2]^T$.

We first formalize the results discussed in Section 7 of ^[42] that handle the i.n.i.d case of posterior consistency. As such, we say that our posterior distributions $\{ \pi(\cdot \vert \{ \mathit{\boldsymbol{V}}^i \}) \}$ as in eq. $(3.2)$ are strongly consistent at $\mathit{\boldsymbol{q}}_0$ if $\{ \pi(U \vert \{ \mathit{\boldsymbol{V}}^i \}) \} \rightarrow 1$ a.s $P_{q_0}^{\infty}$ for every neighborhood $U$ of $\mathit{\boldsymbol{q}}_0$ as $R\rightarrow \infty$, where $P_{ \mathit{\boldsymbol{q}}_0}^{\infty} = \prod_{i = 1}^{\infty}P_{i, \mathit{\boldsymbol{q}}_0}$ with $P_{i, \mathit{\boldsymbol{q}}_0}$ the probability distribution generated by $f_{i, \mathit{\boldsymbol{q}}_0}$ with data samples $\{ \mathit{\boldsymbol{V}}^i \}$.

For this we show that for sets $A\subset Q$ with $\mathit{\boldsymbol{q}}_0 \notin A$, $J_A(\{ V_k \})\rightarrow 0$ and $J(\{ \mathit{\boldsymbol{V}}^i \})\rightarrow \infty$ as $R\rightarrow \infty$ in some appropriate manner to be made precise below. For $J_A(\{ \mathit{\boldsymbol{V}}^i \})\rightarrow 0$ we take the same approach as expressed in ^[42] and thus state the following definition below without motivation, where we note that for any two densities $f$, $g$ on some space $\mathcal{X}$ their affinity, denoted $\text{Aff}(f, g)$, is given by $\text{Aff}(f, g) = \int_{\mathcal{X}} \sqrt{f(x)g(x)} dx$.

Definition 3.1. Let $A\subset Q$ and $\delta > 0$. The set $A$ and $\mathit{\boldsymbol{q}}_0$ are strongly $\delta$ separated if for every probability measure $\nu$ on $A$, $\text{Aff}(f_{1, \mathit{\boldsymbol{q}}_0}, \mathcal{v}_{\nu}^1) < \delta$ where $f_{1, \mathit{\boldsymbol{q}}_0}(\mathit{\boldsymbol{x}}) = \prod_{j = 0}^K \varphi(x_j - y^1_j(\mathit{\boldsymbol{q}}_0))$ for $\mathit{\boldsymbol{x}}\in \mathbb{R}^{K+1}$ as in the work surrounding eq. $(3.2)$, and $\mathcal{v}_{\nu}^R$ is the marginal density of $\{ \mathit{\boldsymbol{\mathcal{V}}}^i \}_{i = 1}^R$ given by $\mathcal{v}_{\nu}^R (\{ \mathit{\boldsymbol{V}}^i \}_{i = 1}^R) = \int_A L(\mathit{\boldsymbol{q}} | \{ \mathit{\boldsymbol{V}}^i \}_{i = 1}^R) d\nu(\mathit{\boldsymbol{q}})$ for any $R = 1, 2, \ldots$. We will say that $A$ and $\mathit{\boldsymbol{q}}_0$ are strongly separated if they are strongly $\delta$ separated for some $\delta > 0$.

From these definitions we provide the following theorem without proof. For proof see Sections 3 and 7 of ^[42] (or, for examples, ^[43])

Theorem 3.2. Let $\pi_0$ be a prior over parameter space $Q$, $\{ \mathit{\boldsymbol{V}}^i \}_{i = 1}^{\infty}$ be independent but not identically distributed data with distribution generated by $f_{i, \mathit{\boldsymbol{q}}}$ for $\mathit{\boldsymbol{q}} \in Q$, $\mathit{\boldsymbol{q}}_0 \in Q$ the true value of the parameters $[q_1, q_2]^T$, and $A \subset Q$ with $\mathit{\boldsymbol{q}}_0 \notin A$. If $A = \bigcup\limits_{i \geq 1} A_i$ such that

1. For some $\delta > 0$, all $A_i$'s are strongly $\delta$ separated from $\mathit{\boldsymbol{q}}_0$ for the model $\mathit{\boldsymbol{q}} \mapsto f_{i, \mathit{\boldsymbol{q}}}$, and

2. $\sum\limits_{i \geq 1} \sqrt{\pi_0(A_i)} < \infty$,

then for some $\beta_0 > 0$, $e^{R\beta_0}J_A(\{ \mathit{\boldsymbol{V}}^i \}_{i = 1}^R) \rightarrow 0 \quad \mathit{\text{a.s.}} P_{ \mathit{\boldsymbol{q}}_0}^{\infty}$ as $R \rightarrow \infty$ for $J_A(\{ \mathit{\boldsymbol{V}}^i \})$ as in eq. $(3.2)$.

To show $J(\{ \mathit{\boldsymbol{V}}^i \})\rightarrow \infty$ as $R\rightarrow \infty$ we utilize the approach as outlined in the proof Theorem 1 of Appendix A.2 in ^[44] (specifically the proof of (8) in Appendix A.2). For a direct proof see ^[32]. For a similar approach see ^[45].

Theorem 3.3. Let $\pi_0$ be a prior over parameter space $Q$, $\{ \mathit{\boldsymbol{V}}^i \}_{i = 1}^{\infty}$ be independent but not identically distributed data with distribution generated by $f_{i, \mathit{\boldsymbol{q}}}$ for $\mathit{\boldsymbol{q}} \in Q$, and $\mathit{\boldsymbol{q}}_0 \in Q$ the true value of the parameters $[q_1, q_2]^T$. For $\mathit{\boldsymbol{q}}\in Q$ define $\Lambda_i(\mathit{\boldsymbol{q}}_0, \mathit{\boldsymbol{q}}) = \log\dfrac{f_{i, \mathit{\boldsymbol{q}}_0}(\mathit{\boldsymbol{V}}^i)}{f_{i, \mathit{\boldsymbol{q}}}(\mathit{\boldsymbol{V}}^i)}$, $\mathcal{K}_i(\mathit{\boldsymbol{q}}_0, \mathit{\boldsymbol{q}}) = \mathbb{E}_{ \mathit{\boldsymbol{q}}_0}[\Lambda_i(\mathit{\boldsymbol{q}}_0, \mathit{\boldsymbol{q}})]$, and $S_i(\mathit{\boldsymbol{q}}_0, \mathit{\boldsymbol{q}}) = \mathrm{Var}_{ \mathit{\boldsymbol{q}}_0}[\Lambda_i(\mathit{\boldsymbol{q}}_0, \mathit{\boldsymbol{q}})]$. If there exists a set $B \subseteq Q$ with $\pi_0(B) > 0$ such that

1. $\sum\limits_{i \geq 1} \dfrac{S_i(\mathit{\boldsymbol{q}}_0, \mathit{\boldsymbol{q}})}{i^2} < \infty \quad\forall \mathit{\boldsymbol{q}} \in B$, and

2. For every $\varepsilon > 0$, $\pi_0(B \cap \{ \mathit{\boldsymbol{q}} : \mathcal{K}_i(\mathit{\boldsymbol{q}}_0, \mathit{\boldsymbol{q}}) < \varepsilon \quad \forall i \}) > 0$.

Then $\forall \beta > 0$, $e^{R\beta}J(\{ \mathit{\boldsymbol{V}}^i \}_{i = 1}^R) \rightarrow \infty \quad \mathit{\text{a.s.}} P_{ \mathit{\boldsymbol{q}}_0}^{\infty}$ as $R\rightarrow \infty$ for $J(\{ \mathit{\boldsymbol{V}}^i \})$ as in eq. $(3.2)$.

Before moving on to our main theorem we apply the following theorem and subsequent corollary to prove a lemma that will be of use to us later. For proof of the following see Theorem 5.3 of ^[46].

Theorem 3.4. For parameters $\mathit{\boldsymbol{q}} \in Q$ a subset of the positive orthant of $\mathbb{R}^n$, and $\mathit{\boldsymbol{q}}$-dependent semigroup $\{ T(\cdot; \mathit{\boldsymbol{q}}) : t > 0 \}$ with infinitesimal generator $A(\mathit{\boldsymbol{q}})$ defined in terms of a $\mathit{\boldsymbol{q}}$-dependent sesquilinear form $\sigma(\mathit{\boldsymbol{q}}): V\times V \rightarrow \mathbb{C}$ on a Hilbert space $V$ satisfying items 1 and 2 in Section 2.1, assume

1. (Affine) The map $\mathit{\boldsymbol{q}} \mapsto \sigma(\mathit{\boldsymbol{q}})$ is affine, in the sense that for any $u, v \in V$, $\sigma(\mathit{\boldsymbol{q}}, u, v) = \sigma_0(u, v) + \sigma_1(\mathit{\boldsymbol{q}}, u, v)$ where $\sigma_0$ is independent of $\mathit{\boldsymbol{q}}$ and the map $\mathit{\boldsymbol{q}} \mapsto \sigma_1(\mathit{\boldsymbol{q}}, \cdot, \cdot)$ is linear, and

2. (Continuous) For any $\mathit{\boldsymbol{q}}, \bar{ \mathit{\boldsymbol{q}}} \in Q$ with metric $d_Q(\cdot, \cdot)$ we have the bound $\vert\sigma(\mathit{\boldsymbol{q}}, u, v) - \sigma(\bar{ \mathit{\boldsymbol{q}}}, u, v)\vert \leq d_Q(\mathit{\boldsymbol{q}}, \bar{ \mathit{\boldsymbol{q}}}) \lVert u \rVert_V\lVert v \rVert_V$, for all $u, v \in V$.

Then, the semigroup $T(\cdot; \mathit{\boldsymbol{q}})$ is (Fréchet) differentiable in $\mathit{\boldsymbol{q}}$ in the interior of $Q$, where for $t > 0$, $\bar{ \mathit{\boldsymbol{q}}} \in Q$, and acting on $\delta \mathit{\boldsymbol{q}} \in Q$ the derivative is given by

for $R(\lambda, \mathit{\boldsymbol{q}}) = (A(\mathit{\boldsymbol{q}})-\lambda I)^{-1}$ the resolvent of $A(\mathit{\boldsymbol{q}})$, and the obtuse sector $\Sigma_{\gamma} = \{ \lambda \in \mathbb{C} : \arg(\lambda - \lambda_0) \leq \frac{\pi}{2} + \arctan((1+\alpha_0\mu_0)(1-\gamma)) \}$ with $\gamma \in (0, 1)$, $\alpha_0$ as in item 1, and $\lambda_0$, $\mu_0$ as in item 2.

The following corollary is an immediate consequence of the work in ^[6]. Specifically that the map $\mathit{\boldsymbol{q}} \mapsto R(\lambda, \mathit{\boldsymbol{q}})$ is analytic as a map from $Q$ to $\mathscr{L}(V^{\ast}, V)$, and from eq. $(3.3)$ the map $\bar{ \mathit{\boldsymbol{q}}} \mapsto T_{ \mathit{\boldsymbol{q}}}(t, \bar{ \mathit{\boldsymbol{q}}})$ depends continuously on $R(\lambda, \bar{ \mathit{\boldsymbol{q}}})$. Note that $\Sigma_{\gamma}$ is independent of $\mathit{\boldsymbol{q}}$ as the constants $\alpha_0$ and $\lambda_0$, $\mu_0$ from items 1 and 2, respectively, of Section 2.1 are independent of $\mathit{\boldsymbol{q}}$.

Corollary 3.2. Under the same hypotheses as Theorem 3.4, we have that the map $\bar{ \mathit{\boldsymbol{q}}} \mapsto T_{ \mathit{\boldsymbol{q}}}(t, \bar{ \mathit{\boldsymbol{q}}})$ is continuous in the operator norm on $\mathscr{L}(Q, \mathscr{L}(V^{\ast}, V))$ for $\bar{ \mathit{\boldsymbol{q}}}$ in the interior of $Q$.

Lemma 3.1. For $Q$ a rectangle in the positive orthant of $\mathbb{R}^2$, Hilbert spaces $H_q$ and $V$ as in eq. $(2.4)$, bilinear form $a(\mathit{\boldsymbol{q}}, \cdot, \cdot): V \times V \rightarrow \mathbb{R}$ as in Section 2.1 with $q_3$ and $q_4$ assumed known, and induced infinitesimal generator $A(\mathit{\boldsymbol{q}})$ as in Section 2.1, then the generated holomorphic semigroup of bounded linear operators $\{ e^{A(\mathit{\boldsymbol{q}})t} : t \geq 0 \}$ on $H_{ \mathit{\boldsymbol{q}}}$ and $V^\ast$ is (Fréchet) differentiable and Lipschitz in $\mathit{\boldsymbol{q}}$ in the interior of $Q$. Further, for $i = 1, \ldots, R$ and $j = 1, \ldots, K$, $y_j^i$ and $y_j^{n, i}$ as in eq. $(2.12)$ and eq. $(2.12)$ are (Fréchet) differentiable and Lipschitz in $\mathit{\boldsymbol{q}}$ with Lipschitz constants independent of $i$ and $j$.

Proof. First, by item 3 of Section 2.1 we have that for $\mathit{\boldsymbol{q}} \in Q$, $\mathit{\boldsymbol{q}} \mapsto a(\mathit{\boldsymbol{q}}, \cdot, \cdot)$ is continuous in $\mathit{\boldsymbol{q}}$. Second, notice that for $\hat{\psi}, \hat{\varphi} \in V$ with $\hat{\psi} = (\psi(0), \psi)$ and $\hat{\varphi} = (\varphi(0), \varphi)$, we have $a(\mathit{\boldsymbol{q}}, \hat{\psi}, \hat{\varphi}) = \dfrac{q_1 q_4}{q_3}\psi(0)\varphi(0) + q_1 \int\limits_0^1 \psi'(x)\varphi'(x)dx = a_0(\mathit{\boldsymbol{q}}, \psi, \varphi) + a_1(\mathit{\boldsymbol{q}}, \psi, \varphi)$, where $\mathit{\boldsymbol{q}} \mapsto a_0(\mathit{\boldsymbol{q}}, \psi, \varphi)$ and $\mathit{\boldsymbol{q}} \mapsto a_1(\mathit{\boldsymbol{q}}, \psi, \varphi)$ are clearly linear in $\mathit{\boldsymbol{q}}$ for $\mathit{\boldsymbol{q}} \in Q$. Hence the bilinear form $a(\mathit{\boldsymbol{q}}, \cdot, \cdot)$ is affine and continuous in $\mathit{\boldsymbol{q}}$ Thus, by Theorem 3.4, we have that the semigroup generated by $a(\mathit{\boldsymbol{q}}, \cdot, \cdot)$, $\{ e^{A(\mathit{\boldsymbol{q}})t} : t \geq 0 \} = \{ T(t, \mathit{\boldsymbol{q}}) : t \geq 0 \}$, is (Fréchet) differentiable in $\mathit{\boldsymbol{q}}$ for $\mathit{\boldsymbol{q}}\in Q$. Denote the derivative in $\mathit{\boldsymbol{q}}$ and $\bar{ \mathit{\boldsymbol{q}}}$ acting on $\delta \mathit{\boldsymbol{q}} \in Q$ by $T_{ \mathit{\boldsymbol{q}}}(t, \bar{ \mathit{\boldsymbol{q}}})\delta \mathit{\boldsymbol{q}}$.

Moreover for $t \in (0, T]$ we have for $\mathit{\boldsymbol{q}}_1, \mathit{\boldsymbol{q}}_2 \in Q$ and line segment $S = \{ s \mathit{\boldsymbol{q}}_1 + (1-s) \mathit{\boldsymbol{q}}_2 : 0 \leq s \leq 1 \}$,

with $\lVert\cdot\rVert_{\mathscr{L}(V^{\ast}, V)}$ the operator norm, where the first inequality follows from the Mean Value Theorem on Banach spaces (see Theorem 4 of Section 3.2 in ^[47]), and the second and third inequalities follow from the compactness of $Q$, Corollary 3.2, and the continuity of the map $T_{\mathit{\boldsymbol{q}}}(t, \bar{\mathit{\boldsymbol{q}}})) \mapsto \lVert T_{\mathit{\boldsymbol{q}}}(t, \bar{\mathit{\boldsymbol{q}}})) \rVert_{\mathscr{L}(Q, \mathscr{L}(V^{\ast}, V))}$.

Further, under zero-order hold the differentiability and Lipschitz properties of $\hat{A}(\mathit{\boldsymbol{q}}) = e^{A(\mathit{\boldsymbol{q}})\tau} = T(\tau, \mathit{\boldsymbol{q}})$ remain. Considering $\hat{B}(\mathit{\boldsymbol{q}})$ from Section 2.1, $\hat{B}(\mathit{\boldsymbol{q}}) = (I - \hat{A}(\mathit{\boldsymbol{q}}))\big\{(0, \xi) - A(\mathit{\boldsymbol{q}})^{-1}(\dfrac{q_3q_2}{q_1}, 0)\big\}$, we find that it is a sum and product of $\mathit{\boldsymbol{q}}$-differentiable and $\mathit{\boldsymbol{q}}$-Lipschitz terms and thus is differentiable and Lipschitz in $\mathit{\boldsymbol{q}}$. Since $h$ and $y$ as in eq. $(2.12)$ are a composition and sum of $\mathit{\boldsymbol{q}}$-differentiable terms they remain differentiable. Further, using eq. $(3.4)$ we have the following Lipschitz bound for all $j \in \{0, 1, \ldots, K\}$ and $\mathit{\boldsymbol{q}}, \bar{ \mathit{\boldsymbol{q}}}$ in the interior of $Q$,

for $C_1$ the max of the operator norms for $\hat{A}(\mathit{\boldsymbol{q}})$, $\hat{B}(\mathit{\boldsymbol{q}})$, $\hat{C}$ over $Q$, and $\tilde{C}_{\hat{A}}$, $\tilde{B}_{\hat{A}}$ the max Lipschitz constants of $\hat{A}(\mathit{\boldsymbol{q}})$ and $\hat{B}(\mathit{\boldsymbol{q}})$ over all $k$ and $Q$. The final inequality above follows from eq. $(3.4)$ by noticing that for all $\hat{\varphi} \in H_{ \mathit{\boldsymbol{q}}}$, $\lVert \hat{A}(\mathit{\boldsymbol{q}})\hat{\varphi} - \hat{A}(\bar{ \mathit{\boldsymbol{q}}})\hat{\varphi} \rVert_{H_{ \mathit{\boldsymbol{q}}}} \leq C_V \lVert \hat{A}(\mathit{\boldsymbol{q}})\hat{\varphi} - \hat{A}(\bar{ \mathit{\boldsymbol{q}}})\hat{\varphi} \rVert_{V}$ and that by identification the supremum over $V^{\ast}$ is larger than that over $H_{ \mathit{\boldsymbol{q}}}$. For $\hat{A}^n(\mathit{\boldsymbol{q}})$ and $\hat{B}^n(\mathit{\boldsymbol{q}})$ as in eq. $(2.9)$ by a repetition of the above arguments we maintain differentiability in $\mathit{\boldsymbol{q}}$, and thus $h^n$ and $y^n$ as in eq. $(2.12)$ are differentiable and Lipschitz in $\mathit{\boldsymbol{q}}$.

Lastly, for all $i = 1, \ldots, R$, $j = 1, \ldots, K$ we have that for $\mathit{\boldsymbol{q}}, \bar{ \mathit{\boldsymbol{q}}} \in Q$

where $\{u_j^i\}$ are BrAC values bounded by definition to be in $[0, 1]$, $\tilde{M}$ is the Lipschitz constant from eq. $(3.5)$, and $K+1$ is the fixed upper bound on the number of temporal observations, $j$. Hence, the Lipschitz constant for $y_j^i$ is independent of $(i, j)$. By noticing that the previous statement holds for $y^{i, n}_j$ with a repetition of the work leading to eq. $(3.6)$, our lemma has been proved.

A direct consequence of Lemma 3.1 is that for all $i = 1, 2, \ldots, R$ with $\mathcal{K}_i(\bar{ \mathit{\boldsymbol{q}}}, \mathit{\boldsymbol{q}})$ and $\Lambda_i(\bar{ \mathit{\boldsymbol{q}}}, \mathit{\boldsymbol{q}})$ as in the statement of Theorem 3.3, we have

where $\sigma > 0$ is the standard deviation of the $N(0, \sigma^2)$ noise density, and $\bar{M}$ is the Lipschitz constant from eq. $(3.6)$ that is independent of $i$ (and $j$). Thus, for any $\delta^{\ast} > 0$, $i \in \{1, \ldots, R\}$, and $\bar{ \mathit{\boldsymbol{q}}}, \mathit{\boldsymbol{q}} \in \{ \mathit{\boldsymbol{q}} \in Q : \lVert \mathit{\boldsymbol{q}}_0 - \mathit{\boldsymbol{q}} \rVert_1 < \delta^{\ast} \}$, we have that $\lVert f_{i, \bar{ \mathit{\boldsymbol{q}}}} - f_{i, \mathit{\boldsymbol{q}}} \rVert_{L_1} \leq (2 \vert \mathcal{K}_i(\bar{ \mathit{\boldsymbol{q}}}, \mathit{\boldsymbol{q}}) \vert)^{\frac{1}{2}} \leq (2 \tilde{\ell} \lVert \bar{ \mathit{\boldsymbol{q}}} - \mathit{\boldsymbol{q}} \rVert_1)^{\frac{1}{2}} < 2(\tilde{\ell}\delta^{\ast})^{\frac{1}{2}}$ by the relationship between total variation and Kullback-Leibler distances, for $\tilde{\ell}$ as in eq. $(3.7)$. If we let $\mathit{\boldsymbol{q}}^{\ast} \in Q$ be such that $\lVert f_{i, \mathit{\boldsymbol{q}}_0} - f_{i, \mathit{\boldsymbol{q}}^{\ast}} \rVert_{L_1} > \delta^{\ast}$ and consider the set $\mathcal{G} = \{ \mathit{\boldsymbol{q}} \in Q : \lVert \mathit{\boldsymbol{q}}^{\ast} - \mathit{\boldsymbol{q}} \rVert_1 < (\delta^{\ast})^2/(16\tilde{\ell}) \}$, then $\mathcal{G}$ is strongly separated from $\mathit{\boldsymbol{q}}_0$ (see Definition 3.1). This follows from the relationship between Affinity and total variation distance (via the Hellinger distance) as well as by noticing that for any density $\nu$ on $\mathcal{G}$, the marginal density of $\mathit{\boldsymbol{V}}^1$ satisfies $\lVert f_{i, \mathit{\boldsymbol{q}}^{\ast}} - \mathcal{v}_{\nu}^i(\mathit{\boldsymbol{V}}^i) \rVert_{L_1} \leq \delta^{\ast}/2$ (for full example, see ^[32] or Example 3.5 in ^[42]). If this holds for all $i$ then $\mathcal{G}$ and $\mathit{\boldsymbol{q}}_0$ are strongly $\bar{\delta}$ separated with $\bar{\delta}$ independent of $i$.

With this example in mind we note that items 1 and 2 of Theorem 3.2 are satisfied if the following special condition is met:

1. For every $\delta^{\ast} > 0$, there exist sets $A_1, A_2, \ldots$ with $L_1$ diameter less than $\delta^{\ast}$, $\text{diam}(A_i) < \delta^{\ast}$, $\bigcup\limits_{i \geq 1}A_i = Q$, and $\sum\limits_{i \geq 1} \sqrt{\pi_0(A_i)} < \infty$ for the mappings $\mathit{\boldsymbol{q}} \mapsto f_{i, \mathit{\boldsymbol{q}}}$

where $\pi_0$ is the prior over $Q$. This follows from the fact that if special item 1 holds then we may take an $\varepsilon^{\ast}$-neighborhood of $\mathit{\boldsymbol{q}}_0$, $U = \{ \mathit{\boldsymbol{q}} \in Q : \lVert f_{i, \mathit{\boldsymbol{q}}_0} - f_{i, \mathit{\boldsymbol{q}}} \rVert_{L_1} < \varepsilon^{\ast} \quad \forall i\}$. As discussed above, since $\lVert f_{i, \bar{ \mathit{\boldsymbol{q}}}} - f_{i, \mathit{\boldsymbol{q}}} \rVert_{L_1}$ is independent of $i$, $U$ is non-empty and contains the set $\{ \mathit{\boldsymbol{q}} \in Q : \lVert \mathit{\boldsymbol{q}}_0 - \mathit{\boldsymbol{q}} \rVert_1 < (\varepsilon^{\ast})^2/(4\tilde{\ell}) \}$. Now set $\delta^{\ast} = (\varepsilon^{\ast})^2 / (16\tilde{\ell})$, and by compactness cover $Q$ with a finite number of disjoint sets $A_i$ determined by the balls $\{ \mathit{\boldsymbol{q}} \in Q : \lVert \bar{ \mathit{\boldsymbol{q}}}_i - \mathit{\boldsymbol{q}} \rVert_1 < \delta^{\ast} \}$ with model $\mathit{\boldsymbol{q}} \mapsto f_{i, \mathit{\boldsymbol{q}}}$, where $\{\bar{ \mathit{\boldsymbol{q}}}_i\}_{i = 1}^{\gamma}$ represents a finite set of points in $Q$ chosen so that $\bigcup\limits_{i \geq 1}A_i = Q$. From these $A_i$'s we have that the finite subset that intersect with $U^{ \mathit{c}}$ must cover $U^{ \mathit{c}}$. This finite subset of $A_i$'s subsequently satisfies the assumptions of Theorem 3.2. Specifically, the strong separation condition is satisfied as per the discussion leading up to special item 1 by noticing that on each $A_i$ we have $\lVert f_{i, \mathit{\boldsymbol{q}}_0} - f_{i, \bar{ \mathit{\boldsymbol{q}}}_i} \rVert_1 > \varepsilon^{\ast}$, and the convergent sum condition is satisfied by the fact that the $A_i$'s can be considered (made) mutually exclusive with union contained in $Q$. We now state and prove our main theorem.

Theorem 3.5. For $Q$ a rectangle in the interior of the positive orthant of $\mathbb{R}^2$, a prior $\pi_0$ with compact support $Q$ and a density that is continuous on $Q$, i.i.d noise distributed as $N(0, \sigma^2)$ for $\sigma > 0$, data $\{ \mathit{\boldsymbol{V}}^i \}_{i = 1}^{R}$ drawn from independent but not identically distributed distributions generated by $f_{i, \mathit{\boldsymbol{q}}}$ as in eq. $(3.2)$, Hilbert spaces $H_q$ and $V$ as in eq. $(2.4)$, bilinear form $a(\mathit{\boldsymbol{q}}, \cdot, \cdot): V \times V \rightarrow \mathbb{R}$ as in Section 2.1 with $q_3$ and $q_4$ assumed known, induced infinitesimal generator $A(\mathit{\boldsymbol{q}})$ as in Section 2.1, and true parameter $\mathit{\boldsymbol{q}}_0 \in Q$, we have that our posterior $\pi(\cdot \vert \{ \mathit{\boldsymbol{V}}^i\})$ as in eq. $(3.2)$ is consistent for $\mathit{\boldsymbol{q}}_0$ as $R \rightarrow \infty$.

Proof. For any set $A \subset Q$ with $\mathit{\boldsymbol{q}}_0 \notin A$ we will use the form of $\pi(A \vert \{ \mathit{\boldsymbol{V}}^i\})$ as in eq. $(3.2)$ and handle the numerator, $J_A$ and denominator, $J$ separately.

First, as $Q$ is compact, for any $\delta > 0$ we may cover $Q$ by a finite number of sets $A_i$, $i = 1, 2, \ldots, \gamma$ where each $A_i$ is a subset of an $L_1$ ball in $Q$. That is, for every $i$ and $\mathit{\boldsymbol{q}}, \bar{ \mathit{\boldsymbol{q}}} \in A_i$ we have that $\lVert \mathit{\boldsymbol{q}} - \bar{ \mathit{\boldsymbol{q}}} \rVert_1 < \delta$. For $R$ large enough, if on each $A_i$ we consider the model $\mathit{\boldsymbol{q}} \mapsto f_{i, \mathit{\boldsymbol{q}}}$ for $i \in \{ 1, \ldots, R \}$ and $f_{i, \mathit{\boldsymbol{q}}}$ the density of the random variable $\mathit{\boldsymbol{\mathcal{V}}}^i$ with $\mathit{\boldsymbol{q}}$ assumed known, then special item 1 is satisfied for prior $\pi_0$. Hence, by Theorem 3.2 we have that for some $\beta_0 > 0$, $e^{R\beta_0} J_A(\{ \mathit{\boldsymbol{V}}^i \}) \rightarrow 0$ a.s. $P_{ \mathit{\boldsymbol{q}}_0}^{\infty}$ as $R\rightarrow \infty$.

Now for $\Lambda_i$, $\mathcal{K}_i$, and $S_i$ as in the statement of Theorem 3.3, for $i = 1, 2, \ldots, R$ and $\mathit{\boldsymbol{q}} \in Q$, we have $\vert \mathcal{K}_i(\mathit{\boldsymbol{q}}_0, \mathit{\boldsymbol{q}}) \vert \leq \tilde{\ell} \lVert \mathit{\boldsymbol{q}}_0 - \mathit{\boldsymbol{q}} \rVert_1$ and

for $\tilde{\ell}$ as determined by eq. $(3.7)$.

Thus, for $\mathit{\boldsymbol{q}}_0, \mathit{\boldsymbol{q}} \in Q$ we find that $\sum\limits_{i \geq 1} \dfrac{S_i(\mathit{\boldsymbol{q}}_0, \mathit{\boldsymbol{q}})}{i^2} \leq \left(\dfrac{\tilde{\ell}^2}{4} \lVert \mathit{\boldsymbol{q}}_0 - \mathit{\boldsymbol{q}} \rVert^2_1 \right) \sum\limits_{i \geq 1} \dfrac{1}{i^2} < \infty$. Further, by the bounds above we have that for every $\varepsilon > 0$ and $i$, $\{ \mathit{\boldsymbol{q}} : \vert \mathcal{K}_i(\mathit{\boldsymbol{q}}_0, \mathit{\boldsymbol{q}})\vert < \varepsilon \}$ is non-empty and our choice in such $\mathit{\boldsymbol{q}}$ does not depend on $i$. Hence the set $\{ \mathit{\boldsymbol{q}} : \vert \mathcal{K}_i(\mathit{\boldsymbol{q}}_0, \mathit{\boldsymbol{q}})\vert < \varepsilon \quad \forall i\}$ is non-empty. Thus, for $B = Q$ we satisfy the assumptions of Theorem 3.3 and therefore find that $\forall \beta > 0$, $e^{R\beta}J(\{ \mathit{\boldsymbol{V}}^i \}) \rightarrow \infty$ a.s. $P_{ \mathit{\boldsymbol{q}}_0}^{\infty}$ as $k\rightarrow \infty$.

So for any set $A \subset Q$ with $\mathit{\boldsymbol{q}}_0 \notin A$, from eq. $(3.2)$ we have that $\pi(A \vert \{ \mathit{\boldsymbol{V}}^i\}) \rightarrow 0$ a.s. $P_{ \mathit{\boldsymbol{q}}_0}^{\infty}$ as $R\rightarrow \infty$ and thus the theorem has been proved.

From Lemma 3.1 we find that we maintain the differentiability and Lipschitz properties of the finite-dimensional semigroup as in eq. $(2.9)$ and respective kernel as in eq. $(2.12)$. Thus, with a straightforward rewriting of eq. $(3.2)$ and eq. $(3.2)$ in terms of the finite-dimensional posterior eq. $(2.20)$, and repetition of the work following Lemma 3.1 through the proof of Theorem 3.5 we have the following corollary.

Corollary 3.3. Under the same hypotheses as Theorem 3.5, for fixed positive integer $n$ we have that our finite-dimensional posterior $\pi^n(\cdot \vert \{ \mathit{\boldsymbol{V}}^i\})$ as in eq. $(2.20)$ is consistent at $\mathit{\boldsymbol{q}}_0$ as $R \rightarrow \infty$.

3.3.   Deconvolution of BrAC from TAC

In this section we consider the problem of using the biosensor measured TAC signal to estimate BrAC. We do this by deconvolving it; to wit we invert the convolution given in eq. $(2.12)$ subject to a positivity constraint and regularization to mitigate the inherent ill-posedness of the inversion. Recall that the convolution given in eq. $(2.12)$ was found by solving the finite-dimensional discrete time system eq. $(2.10)$ derived from eq. $(2.2)$. We employ the method originally described in ^[25], wherein the problem is formulated as a constrained, regularized, optimization problem (see, for example, ^[48]).

We first briefly summarize the treatment in ^[25] and then follow by showing how our work is able to make direct use of this theory. Let $\tilde{V}$ and $\tilde{H}$ be Hilbert spaces forming a Gelfand Triple, $\tilde{V} \hookrightarrow \tilde{H} \hookrightarrow \tilde{V}^{\ast}$. For an admissible set $Q$, a compact subset of the positive orthant of $\mathbb{R}^{2}$, with $\mathit{\boldsymbol{q}} \in Q$, let $A(\mathit{\boldsymbol{q}})$ be an abstract parabolic operator defined by a sesquilinear form $a(\mathit{\boldsymbol{q}}, \cdot, \cdot): V\times V \rightarrow \mathbb{R}$ (i.e., one that satisfies items 1 to 3 in Section 2.1) that when restricted to $\{ \varphi \in \tilde{V} : A(\mathit{\boldsymbol{q}})\varphi \in \tilde{H} \}$ generates a holomorphic semigroup on $\tilde{H}$, $\{ e^{A(\mathit{\boldsymbol{q}})t} : t\geq 0\}$. For bounded operators $B(\mathit{\boldsymbol{q}}) \in \mathscr{L}(\mathbb{R}, \tilde{H})$ and $C(\mathit{\boldsymbol{q}}) \in \mathscr{L}(\tilde{H}, \mathbb{R})$ consider the input/output system where $\dot{x}(t) = A(\mathit{\boldsymbol{q}})x(t) + B(\mathit{\boldsymbol{q}})u(t)$, $x(0) = x_0 \in \tilde{H}$, and $y(t) = C(\mathit{\boldsymbol{q}})x(t)$ where on the interval $[0, T]$, $u\in L_2(0, T)$ is the input, $y$ the output, and $x$ is the state variable.

For sampling interval $\tau > 0$ and zero-order hold input $u(t) = u_j$, $t \in [j\tau, (j+1)\tau)$, $j = 0, 1, 2, ...$ the corresponding sampled-time system is given by

where $\hat{A}(\mathit{\boldsymbol{q}}) = e^{A(\mathit{\boldsymbol{q}})\tau} \in \mathscr{L}(\tilde{H}, \tilde{H})$, $\hat{B}(\mathit{\boldsymbol{q}}) = \int_0^{\tau} e^{A(\mathit{\boldsymbol{q}})s}B(\mathit{\boldsymbol{q}}) ds \in \mathscr{L}(\mathbb{R}, \tilde{H})$, and $\hat{C}(\mathit{\boldsymbol{q}}) = C(\mathit{\boldsymbol{q}}) \in \mathscr{L}(\tilde{H}, \mathbb{R})$, $x_j = x(j\tau) \in \tilde{H}$, $y_j = y(j\tau) \in \mathbb{R}$ $j = 0, 1, 2, \ldots$, and for all $j$, $\{u_j\} \subset \mathbb{R}$ are zero-order hold input values.

Now let $\mathit{\boldsymbol{\mathcal{q}}}$ be a random variable with support the parameter space $Q$. For $\tilde{\pi}$ the probability measure of $\mathit{\boldsymbol{\mathcal{q}}}$, define the Bochner spaces $\mathcal{V} = L^2_{\tilde{\pi}}(Q; \tilde{V})$, $\mathcal{H} = L^2_{\tilde{\pi}}(Q; \tilde{H})$, and $\mathcal{U} = L^2_{\tilde{\pi}}(Q; \mathbb{R})$. It is easily shown that the spaces $\mathcal{V}$ and $\mathcal{H}$ form a Gelfand triple $\mathcal{V} \hookrightarrow \mathcal{H} \hookrightarrow \mathcal{V}^{\ast}$. Define $\mathcal{a}(\cdot, \cdot): \mathcal{V} \times \mathcal{V}\rightarrow \mathbb{R}$ by $\mathcal{a}(\varphi, \psi) = \mathbb{E}_{\tilde{\pi}}[a(\mathit{\boldsymbol{\mathcal{q}}}, \varphi(\mathit{\boldsymbol{\mathcal{q}}}), \psi(\mathit{\boldsymbol{\mathcal{q}}}))] = \int_{Q} a(\mathit{\boldsymbol{q}}, \varphi(\mathit{\boldsymbol{q}}), \psi(\mathit{\boldsymbol{q}})) d\tilde{\pi}(\mathit{\boldsymbol{q}})$ for $\varphi, \psi \in \mathcal{V}$. Then, as in Section 2.1, the form $-\mathcal{a}(\cdot, \cdot)$ satisfies items 1 to 3 and therefore defines a linear map $\mathcal{A}$ that when restricted to $\{ \varphi \in \mathcal{V} : \mathcal{A}\varphi \in \mathcal{H} \}$, generates an analytic semigroup on $\mathcal{H}$, $\{ e^{\mathcal{A}t} : t\geq 0\}$.

Assume further that the map $\mathit{\boldsymbol{q}} \mapsto B(\mathit{\boldsymbol{q}})$ is in $L_\infty(Q, \mathscr{L}(\mathbb{R}, \tilde{H}))$ and that the map $\mathit{\boldsymbol{q}} \mapsto C(\mathit{\boldsymbol{q}})$ is in $L_2(Q, \mathscr{L}(\tilde{H}, \mathbb{R}))$ with respect to the measure $\tilde{\pi}$. (Note that since the domain space in $\mathscr{L}(\mathbb{R}, \tilde{H})$ and the co-domain space in $\mathscr{L}(\tilde{H}, \mathbb{R})$ are both $\mathbb{R}$, it follows that in fact the mapping $\mathit{\boldsymbol{q}} \mapsto B(\mathit{\boldsymbol{q}}) \in L_\infty(Q, \tilde{H})$, and by the Riesz Representation Theorem, that effectively the mapping $\mathit{\boldsymbol{q}} \mapsto C(\mathit{\boldsymbol{q}}) \in L_2(Q, \tilde{H}) = \mathcal{H}$). Then define bounded linear operators $\mathcal{B} \in \mathscr{L}(\mathcal{U}, \mathcal{H})$ and $\mathcal{C} \in \mathscr{L}(\mathcal{H}, \mathbb{R})$ by $\langle \mathcal{B}\mathcal{u}, \psi \rangle_{\mathcal{H}} = \mathbb{E}_{\tilde{\pi}}[\langle B(\mathit{\boldsymbol{\mathcal{q}}})\mathcal{u}(\mathit{\boldsymbol{\mathcal{q}}}), \psi(\mathit{\boldsymbol{\mathcal{q}}}) \rangle_{\tilde{H}}] = \int_{Q} \langle B(\mathit{\boldsymbol{q}})\mathcal{u}(\mathit{\boldsymbol{q}}), \psi(\mathit{\boldsymbol{q}}) \rangle_{\tilde{H}} d\tilde{\pi}(\mathit{\boldsymbol{q}})$ and $\mathcal{C}\psi = \mathbb{E}_{\tilde{\pi}}[C(\mathit{\boldsymbol{\mathcal{q}}})\psi(\mathit{\boldsymbol{\mathcal{q}}})] = \int_{Q} C(\mathit{\boldsymbol{q}})\psi(\mathit{\boldsymbol{q}}) d\tilde{\pi}(\mathit{\boldsymbol{q}})$, respectively, for $\mathcal{u} \in \mathcal{U}$ and $\psi \in \mathcal{H}$. It can then be shown ^[49,50] that for $\mathcal{u} \in L_2([0, T], \mathcal{U})$ the solution to the input/output system from above with $u(t) = u(t, \mathit{\boldsymbol{q}}) = \mathcal{u}(t)$ agrees with the solution to the system where $\dot{ \mathcal{x}}(t) = \mathcal{A} \mathcal{x}(t) + \mathcal{B}\mathcal{u}(t)$, $\mathcal{x}(0) = \mathcal{x}_0 \in \mathcal{H}$, and $\mathcal{y}(t) = \mathcal{C} \mathcal{x}(t)$ for $\tilde{\pi}-$almost every $\mathit{\boldsymbol{q}} \in Q$. Then with sampling interval $\tau > 0$ as in eq. $(3.8)$, and zero-order hold input $\mathcal{u}(t) = \mathcal{u}_j$, $t \in [j\tau, (j+1)\tau)$, $\mathcal{u}_j \in \mathcal{U}$, $j = 0, 1, 2, ...$, this system becomes

for $j = 0, 1, 2, \ldots$ where $\hat{\mathcal{A}} = e^{\mathcal{A}\tau} \in \mathscr{L}(\mathcal{H}, \mathcal{H})$, $\hat{\mathcal{B}} = \int_0^{\tau} e^{\mathcal{A}s} \mathcal{B} ds \in \mathscr{L}(\mathcal{U}, \mathcal{H})$, and $\hat{\mathcal{C}} \in \mathscr{L}(\mathcal{H}, \mathbb{R})$.

Now, with eq. $(3.9)$, note that $\{ \mathcal{u}_j\} \subset \mathcal{U}$ is obtained by zero-order hold sampling a continuous time signal. That is, the input to eq. $(3.9)$ is $\mathcal{u}(j\tau) = \mathcal{u}_j \in \mathcal{U}$ with $\mathcal{u}$ at least continuous on $[0, T]$. We seek an estimate for the input based on this model, wherein the input estimate $\mathcal{u}$ is a function of both time and the random parameters $\mathit{\boldsymbol{\mathcal{q}}}$. For optimization purposes (more precisely, to be able to include regularization) we require additional regularity. Given the time interval $[0, T]$, let $\mathcal{u} \in S(0, T) = H^1(0, T; \mathcal{U})$ and let $\mathcal{U}$ be a compact subset of $S(0, T)$.

The input estimation or deconvolution problem is then given by

where $\lVert \cdot \rVert^2_{S(0, T)}$ is a norm on $S(0, T)$ that will be defined below, $\{\hat{y}_k\}$ are measured output values, the term $\lVert \mathcal{u} \rVert^2_{S(0, T)}$ serves as regularization, and $\mathcal{y}_k(\mathcal{u}) = \sum\limits_{j = 0}^{k-1}\langle h_{k-j-1}, \mathcal{u}_j\rangle_{\mathcal{U}},$ for $k = 1, 2, \ldots, K$ with $\mathcal{u}_j = \mathcal{u}(j\tau)$ for $j = 1, 2, \ldots, K$ the zero-order hold input to the discrete time system eq. $(3.9)$, and convolution filter $h_{\ell} = \hat{\mathcal{C}}\hat{\mathcal{A}}^{\ell -1} \hat{\mathcal{B}} \in \mathscr{L}(\mathcal{U}, \mathbb{R}) = \mathcal{U}^\ast$ (which is equal to $\mathcal{U}$, by the Riesz Representation Theorem) where $\hat{\mathcal{C}} \in \mathscr{L}(\mathcal{H}, \mathbb{R})$, $\hat{\mathcal{A}} \in \mathscr{L}(\mathcal{H}, \mathcal{H})$, and $\hat{\mathcal{B}} = \int_{Q} e^{\mathcal{A}s} \mathcal{B} ds \in \mathscr{L}(\mathcal{U}, \mathcal{H})$.

Solving eq. $(3.10)$ requires finite dimensional approximations. For index $M$, let $\mathcal{U}^M$ define an approximating family of closed subsets of $\mathcal{U}$, where each subset is contained within a corresponding finite dimensional subspace, $\mathcal{S}^M$ of $S(0, T)$. Further we require that for each $\mathcal{u} \in \mathcal{U}$ there exists a sequence $\{ \mathcal{u}^M \}$ with $\mathcal{u}^M \in \mathcal{U}^M$ such that $\mathcal{u}^M \rightarrow \mathcal{u}$ in $S(0, T)$ as $M\rightarrow \infty$. For index $N$, let $\mathcal{V}^N$ be an element of an approximating family of finite-dimensional subspaces of $\mathcal{V}$, and let $\mathscr{P}^N_{\mathcal{H}}: \mathcal{H} \rightarrow \mathcal{V}^N$ be the orthogonal projection of $\mathcal{H}$ onto $\mathcal{V}^N$. We also require of the spaces $\mathcal{V}^N$ that for each $\mathcal{v}\in\mathcal{V}$, $\mathscr{P}^N_{\mathcal{H}}\mathcal{v} \rightarrow \mathcal{v}$ in $\mathcal{V}$ as $N\rightarrow \infty$.

We next specify finite-dimensional operators $\hat{\mathcal{A}}^N \in \mathscr{L}(\mathcal{V}^N, \mathcal{V}^N)$, $\hat{\mathcal{B}}^N \in \mathscr{L}(\mathcal{U}, \mathcal{V}^N)$, and $\hat{\mathcal{C}}^N \in \mathscr{L}(\mathcal{V}^N, \mathbb{R})$ that define the finite-dimensional system analogous to eq. $(3.9)$. That is, let $\mathcal{A}^N: \mathcal{V}^N \rightarrow \mathcal{V}^N$ be given by $\langle \mathcal{A}^N\varphi^N, \psi^N \rangle_{\mathcal{H}} = -\mathcal{a}(\varphi^N, \psi^N)$ for $\varphi^N, \psi^N \in \mathcal{V}^N$, $\hat{\mathcal{A}}^N = e^{\mathcal{A}^N\tau}$, $\mathcal{B}^N = \mathscr{P}^N_{\mathcal{H}} \mathcal{B}$, $\hat{\mathcal{B}}^N = \int_0^{\tau} e^{\mathcal{A}^Ns}\mathcal{B}^N ds$, and $\hat{\mathcal{C}}^N = \mathcal{C}$. In this way we obtain a doubly-indexed sequence of approximating finite-dimensional optimization or deconvolution problems given by

where $L = (M, N)$, $h^N_{\ell} = \hat{\mathcal{C}}^N(\hat{\mathcal{A}}^N)^{\ell}\hat{\mathcal{B}}^N \in \mathcal{U}$, and

Using the approximation properties of the subspaces $\mathcal{V}^N$ and $\mathcal{U}^M$ (that is, that for each $\mathcal{u} \in \mathcal{U}^M$ there exists a sequence $\{ \mathcal{u}^M \}$ with $\mathcal{u}^M \in \mathcal{U}^M$ and $\lVert \mathcal{u}^M - \mathcal{u} \rVert_{S(0, T)} \rightarrow 0$ as $M\rightarrow \infty$, and that for each $v \in \mathcal{V}$, $\lVert \mathscr{P}^N_{\mathcal{H}}v - v \rVert_{\mathcal{V}} \rightarrow 0$ as $N \rightarrow \infty$), and the corresponding operators $\hat{\mathcal{A}}^N \in \mathscr{L}(\mathcal{V}^N, \mathcal{V}^N)$, $\hat{\mathcal{B}}^N \in \mathscr{L}(\mathcal{U}, \mathcal{V}^N)$, and $\hat{\mathcal{C}}^N \in \mathscr{L}(\mathcal{V}^N, \mathbb{R})$, it can be shown that 1) for each multi-index $L$, eq. $(3.11)$ admits a solution $\mathcal{u}_L^{\ast}$, and 2) there exists a subsequence of $\{ \mathcal{u}_L^{\ast} \}$, $\{ \mathcal{u}_{L_k}^{\ast} \} \subset \{ \mathcal{u}_L^{\ast} \}$ with $\mathcal{u}_{L_k}^{\ast} \rightarrow \mathcal{u}^{\ast}$ strongly as $k\rightarrow \infty$ with $\mathcal{u}^{\ast}$ a solution of eq. $(3.10)$. Further, if in addition $\mathcal{U}$ is assumed to be a closed and convex subset of $S(0, T)$, for each $M$, $\mathcal{U}^M$ is a closed and convex subset of $\mathcal{U}$, and the optimization problem given in eq. $(3.10)$ admits a unique (with respect to sampling) solution, then the sequence of solutions to eq. $(3.11)$, $\{ \mathcal{u}_L^{\ast} \}$ converges strongly, or in $S(0, T)$ to the unique solution of eq. $(3.10)$, $\mathcal{u}^{\ast}$. For the proofs of these results see Section 5 of ^[25]

To numerically carry out the requisite computations to actually determine $\mathcal{u}_L^{\ast}$ for given values of $M$, $N$ and $L = (M, N)$, we continue to apply the results in ^[25] while also connecting them to our treatment in Sections 2.1 and 2.2 above. We assume that the feasible parameter set $Q$ is a compact rectangle in the positive orthant of $\mathbb{R}^2$, we set $\tilde{H} = H_{ \mathit{\boldsymbol{q}}}$ and $\tilde{V} = V$ as in eq. $(2.4)$, and we identify the operators in eq. $(3.8)$ with those in eq. $(2.6)$. Our distribution over $\mathit{\boldsymbol{\mathcal{q}}}$, $\tilde{\pi}$, is the finite-dimensional posterior $\pi^n(\cdot \vert \{V_j^i\})$ for fixed $n$ as in eq. $(2.20)$ and we proceed with the Bochner spaces $\mathcal{V} = L^2_{\pi^n(\cdot \vert \{V_j^i\})}(Q; V)$ and $\mathcal{H} = L^2_{\pi^n(\cdot \vert \{V_j^i\})}(Q; H_{ \mathit{\boldsymbol{q}}})$ to achieve eq. $(3.9)$.

For the state variables $\mathcal{x}_j(\eta, \mathit{\boldsymbol{q}})$ we have that $\eta \in [0, 1]$ and $\mathit{\boldsymbol{q}} \in Q = [a_1, b_1] \times [a_2, b_2]$ for $0 < a_i < b_i$, $i = 1, 2$. Further, for the inputs $\mathcal{u}(t, \mathit{\boldsymbol{q}})$ we have that $t \in [0, T]$ and $\mathit{\boldsymbol{q}} \in Q$. Let $n$ be as in eq. $(2.8)$ and $m$ a positive integer, and we discretize $[0, 1]$ and $[0, T]$ using the sets of linear B-splines, $\{ \psi_j^{n} \}_{j = 0}^{n}$ and $\{ \zeta_j^{m} \}_{j = 0}^{m}$, respectively, on the uniform meshes, $\{ \frac{j}{n} \}_{j = 0}^{n}$ and $\{ \frac{i T}{m} \}_{j = 0}^{m}$, respectively. Further, for positive integers $m_1$ and $m_2$, we discretize $Q$ with the $0^{\text{th}}$-order B-splines $\{ \chi_{i, j}^{m_i} \}_{j = 1}^{m_i}$, $i = 1, 2$ defined with respect to the uniform grids $\{ a_i - \frac{(b_i - a_i)j}{m_i} \}_{j = 0}^{m_i}$, on $[a_i, b_i]$, $i = 1, 2$.

Then for multi-indices $N = (n, m_1, m_2)$ and $M = (m, m_1, m_2)$ we define the approximating subspaces $\mathcal{V}^N$ and $\mathcal{S}^M$ as follows using tensor products. That is, let $\mathcal{V}^N = \text{span}\{ \hat{\psi}_i^N \}_{i = 1}^{^\mathcal{N}} = \text{span}\{(\psi^{n}_j(0) \chi^{m_1}_{1, j_1} \chi^{m_{2}}_{2, j_2}, \psi^{n}_j \chi^{m_1}_{1, j_1} \chi^{m_{2}}_{2, j_{2}}) \}_{j = 0, j_1 = 1, j_{2} = 1}^{n, m_1, m_{2}}$ and $\mathcal{S}^M = \text{span}\{(\zeta^{m}_j \chi^{m_1}_{1, j_1} \chi^{m_{2}}_{2, j_2}\}_{j = 0, j_1 = 1, j_{2} = 1}^{m, m_1, m_{2}} = \text{span}\{ \hat{\zeta}_i^M \}_{i = 1}^\mathcal{M}$, where $\mathcal{N} = (n+1)m_1m_2$ and $\mathcal{M} = (m+1)m_1m_2$. Standard approximation theoretic arguments (see, for example, ^[33]) can be used to argue that the subspaces defined in above satisfy the required approximation assumptions on $\mathcal{V^N}$ and $\mathcal{S}^M$. Then $\mathcal{x}^N \in \mathcal{V^N}$ and $\mathcal{u}^M \in \mathcal{S}^M$, can be written as $\mathcal{x}^N(\eta; \mathit{\boldsymbol{q}}) = \sum\limits_{i = 0, i_1 = 1, i_{2} = 1}^{n, m_1, m_{2}} \mathcal{x}^N_{i, i_1, i_{2}} \psi_i^{n}(\eta) \chi^{m_1}_{1, i_1}(q_1) \chi^{m_2}_{2, i_2}(q_2)$ and $\mathcal{u}^M(t; \mathit{\boldsymbol{q}}) = \sum\limits_{i = 0, i_1 = 1, i_{2} = 1}^{m, m_1, m_{2}} \mathcal{u}^M_{i, i_1, i_{2}} \zeta_i^{m}(t) \chi^{m_1}_{1, i_1}(q_1) \chi^{m_2}_{2, i_2}(q_2)$, respectively.

Then with the bases for $\mathcal{V}^N$ and $\mathcal{S}^M$ as chosen above, it is an elementary exercise to determine the matrix representations for the operators $\mathcal{A}^N$, $\mathcal{B}^N$, $\mathcal{C}^N$, $\hat{\mathcal{A}}^N$, $\hat{\mathcal{B}}^N$, and $\hat{\mathcal{C}}^N$. It then follows that eq. $(3.9)$ takes a matrix system where for $k = 1, 2, \ldots, K$, $\mathbb{M}^N\mathbb{X}^N_{k+1} = \mathbb{K}^N\mathbb{X}^N_{k} + \mathbb{B}^L\mathbb{U}^N_{k}$ and $\mathcal{y}^L_k = \mathbb{C}^N\mathbb{X}^N_k$ with $\mathbb{X}^N_k \in \mathbb{R}^\mathcal{N}$ the coefficients of the basis elements $\{ \hat{\psi}^N_i \}$, $\mathbb{U}^N_k \in \mathbb{R}^\mathcal{M}$ the coefficients of the basis elements $\{ \hat{\zeta}^M_i\}$ as in the previously mentioned approximating subspaces, $\mathbb{M}^N \in \mathbb{R}^{\mathcal{N}\times \mathcal{N}}$ a matrix with entries $[\mathbb{M}^N]_{i, j} = \langle \hat{\psi}_i, \hat{\psi}_j\rangle_{\mathcal{H}}$, $\mathbb{K}^N \in \mathbb{R}^{\mathcal{N}\times \mathcal{N}}$ a matrix with entries $[\mathbb{K}^N]_{i, j} = \langle \mathcal{A}^N\hat{\psi}^N_i, \hat{\psi}^N_j\rangle_{\mathcal{H}}$, $\mathbb{B}^L \in \mathbb{R}^{\mathcal{N}\times \mathcal{M}}$ a matrix with entries $\langle \mathcal{B}^N \hat{\zeta}^{M}_i, \hat{\psi}^N_j \rangle_{\mathcal{H}}$, and $\mathbb{C}^N \in \mathbb{R}^{1\times \mathcal{M}}$ given by $[1, 0, \ldots, 0]$. From here the matrix representation of $h_{k}^L$ (with $L = (M, N)$ in place of $N$ due to the joint dependence on the multi-indices $M$ and $N$) can be found using this matrix system.

We note that the optimization problem eq. $(3.11)$ is a constrained problem, in that $\mathbb{U}^N_k$ of the previously stated matrix system are to be non-negative. With a proper placement of $\{ h_k^L \}$ into the block matrix $\mathbb{H}^L$, the approximating deconvolution problem eq. $(3.11)$ is now given by $u_{L}^{*} = \underset{\mathscr{U}^{M}}{\arg \min } J^{L}(u)$ where we have that

where $\mathbb{U}^M$ is the $K\mathcal{M}$ dimensional column vector of the coefficients of $\mathcal{u} \in \mathscr{U}^{M}$, and $\mathbb{Y}^M$ is the $K + K\mathcal{M}$ column vector of measured output values $\{ \hat{y}_k \}$ followed by $K\mathcal{M}$ zeros. Further, $\mathbb{Q}^M_i$ for $i = 1, 2$ are matrices with entries given by the $\mathcal{U}$ inner products of the basis elements for the subspaces $\mathcal{S}^M$ as determined by the regularization term $\| \mathcal{u}\|_{S(0, T)}^{2}$ below. Note that the regularization term $\| \mathcal{u}\|_{S(0, T)}^{2}$ is derived from a weighted inner product on $S[0, T]$ and thus corresponds to a squared norm on $S(0, T)$ given by $\| \mathcal{u}\|_{S(0, T)}^{2} = r_{1} \int_{0}^{T}\| \mathcal{u}(t)\|_{\mathcal{U}}^{2} d t+r_{2} \int_{0}^{T}\|\dot{ \mathcal{u}}(t)\|_{\mathcal{U}}^{2} d t$.

The values of the regularization weights used in eq. $(3.13)$, $r_i = r^{\ast}_i > 0$ for $i = 1, 2$ are chosen optimally. Indeed, in order to find $(r_1^{\ast}, r_2^{\ast})$, BrAC-TAC input-output training data pairs, $\{(u^i_{j}, V^i_{j})\}_{i = 0, j = 0}^{P, K}$ are used to optimize $(r_1, r_2)$ via the following scheme:

where $\bar{u}_{L; j}^{i, *} = \mathbb{E}_{\tilde{\pi}(\mathit{\boldsymbol{\theta}}^{\ast})}[\tilde{u}_{L; j}^{i, *}]$, $\tilde{u}_{L; j}^{i, *}$ are the predicted BrAC values found by finding the minimum of $J^L(\ \cdot \; r_1, r_2)$ given by eq. $(3.13)$ with $r_1$ and $r_2$ candidate values for the regularization weights from a specified feasible set in the positive orthant of $\mathbb{R}^2$, $\mathbb{R}^{+}\times \mathbb{R}^+$, and $\tilde{y}_{L; j}^{i, *}$ are the TAC values found by using $\tilde{u}_{L; j}^{i, *}$ as input to eq. $(3.12)$.

3.4.   Numerical results

All of the data used in the studies detailed below, unless otherwise specifically stated (e.g., as in Section 3.4.2), were collected in USC IRB approved human subject experiments designed and run by researchers in the laboratory of one of the authors (S. E. L.) as part of a National Institutes of Health (NIH) funded investigation (see, ^[51]). These experiments were carried out in controlled environments wherein 40 participants completed one to four drinking episodes, with viable data recorded in 146 drinking episodes. BrAC was obtained using Alco-sensor IV breath analyzer devices from Intoximeters, Inc, St. Louis, MO, and participants each wore two SCRAM (Secure Continuous Remote Alcohol Monitoring) devices manufactured by Alcohol Monitoring Systems (AMS) in Littleton, Colorado (see Figure 1) simultaneously placed on the participants' left and right arms for TAC. For each separate SCRAM device, participants started their readings with a TAC and BrAC of 0.000, consumed alcohol (equivalent across all sessions per participant) in one of three different drinking patterns (single: over 15 minutes; dual: over two 15-min periods spaced 30-minutes apart; or steady: over 60 minutes), and then ended their session when their TAC and BrAC had returned to 0.000. We note that the placement of the two sensors challenges the independence assumption from Section 2.2, but for experimental purposes we will include all of the data as independently measured drinking episodes with this caveat in mind. In addition, we did not focus on any specific drinking pattern as including all possible patterns is in line with real-world, varying drinking patterns and may improve the generalizability of our model. In the calculations of Sections 3.4.1 and 3.4.2, as in eq. $(2.5)$, time is discretized by a constant sampling time $\tau$ of 5 minutes and is subject to our zero-order hold assumption. While this challenges the implications of our zero order hold assumption, namely that $\tau = .0833$ hours implies that subjects' BAC is constant for 5 minutes, this restriction is needed as computational complexity becomes unstable as $\tau$ decreases. In order to achieve this sampling time, we first linearly interpolate all of the data (both BrAC and TAC), and then re-sample at our desired rate of $\tau = 5$. For Section 3.4.3, a $\tau$ will be discussed. Further, in all sections we assume a truncated multivariate normal (tMVN) prior $\pi_0$ (as in eq. $(2.20)$) on $\mathit{\boldsymbol{q}}$ with mean $\mu$ and covariance matrix $\Sigma$ which varies from example to example.

Unfortunately, the USC IRB approved experiments for collecting human subject data were not designed around the problem of estimating the sensor collection chamber inflow and outflow parameters, $q_3$ and $q_4$ as in eq. $(2.3)$, nor do the authors have the laboratory facilities or expertise to determine them experimentally. Moreover, since the approach developed in sections 2.1 and 2.2, and in particular our underlying hybrid PDE/ODE model given in 2.3, are relatively novel, no values for $q_3$ and $q_4$ are available from either the manufacturers of the sensors or the current literature. Consequently, for the purposes of this study we have chosen values for $q_3$ and $q_4$ arbitrarily as $q_3 = q_4 = 1$. However, we note that the precise values chosen for $q_3$ and $q_4$ had no perceptible qualitative effect on the results to be presented below. Finally we note that all computational work was done in Python 3.7.2 and includes ported MATLAB code from the work of ^{[15,21,24,25]}, in particular with respect to the creation of the finite-dimensional, discrete-time kernel as in eq. $(2.12)$. Ported code was verified against the original code through the use of unit tests.

3.4.1.   Convergence in distribution

We used surface plots as well as Metropolis Hastings (MH) Markov Chain Monte Carlo (MCMC) methods to validate our convergence in distribution results. Throughout the results described here we have that from eq. $(2.14)$ for all sample times the i.i.d. noise $\varepsilon$ is distributed as $N(0, 0.005^2)$, and prior $\pi_0$ as in eq. $(2.20)$ is distributed as the optimal distribution found in Section 6 of ^[25]. Specifically, the prior is a tMVN random variable with mean, $\mu = \left(\begin{array}{l}0.6318\\ 1.0295\end{array} \right)$ and covariance matrix, $\Sigma = \left(\begin{array}{l} 0.0259 & 0.0077\\ 0.0077 & 0.1232 \end{array} \right)$ with the feasible parameter set, $Q$, taken to be $Q = [0.01, 2.2877] \times [0.01, 2.1410]$ for $\mathit{\boldsymbol{q}} = [q_1, q_2]^T$. The choice of $0.005$ for the standard deviation of $\varepsilon$ was made to limit the role of noise in our subsequent sampling algorithms so that we may focus on the role of the dimension of our approximating system in the resulting posterior distribution. In addition, when comparing this choice in standard deviation to the peak TAC values of our training dataset, we had a typical peak TAC to noise ratio of 20. For computational reasons, we limit ourselves to measurements from a random subgroup of $R = 3$ subject drinking episode measurements. Figure 2 contains the resulting surface plots for $n$ values of $1, 3$, and $25$. Further, Table 1 contains the means and credible regions for $n$ values of $1, 2, 3,$ and $25$ as determined by respective 1000 sample (1100 draws with a 100 draw burn-in period) MH MCMC sampling runs. The MCMC sample size chosen here was due to computational complexities and runtimes. Figure 3 displays deconvolution results for a randomly chosen, non-training drinking episode for different values for the dimension of the approximating system, $n$. This figure used the method from Section 3.3 along with the resulting posteriors as shown in Figure 2 and Table 1.

3.4.2.   Consistency

We again used surface plots as well as MH MCMC sampling methods to verify our consistency results. For these studies we have assumed that the noise $\varepsilon$ is now distributed as $N(0, 0.025^2)$ while our prior $\pi_0$ from eq. $(2.20)$ is still the optimal distribution found in Section 6 of ^[25]. That is, $\pi_0$ is a tMVN with $\mu = \left(\begin{array}{l} 0.6318\\ 1.0295 \end{array} \right)$ and covariance matrix, $\Sigma = \left(\begin{array}{l} 0.0259 & 0.0077\\ 0.0077 & 0.1232 \end{array} \right)$ with bounds $[0.01, 2.2877]$ and $[0.01, 2.1410]$ for $q_1$ and $q_2$, respectively. The choice for a distribution for the random noise is meant to simulate a more realistic situation where little is assumed known about any external effects that play a role in perturbing the sensor measurements. Consequently the data is assumed noisy. When comparing this choice for the standard deviation of the noise process to the peak TAC values of our training dataset, we had a typical peak TAC to noise ratio of 8.

To test Theorem 3.5, we generated 276 TAC values using subject-measured BrAC values via eq. $(2.14)$ with a predetermined $\mathit{\boldsymbol{q}}_0$ value of $[1, 1]$, $n = 24$, and noise variance of $0.025^2$. Table 2 displays the calculated means and $90\%$ credible circle radii for the posterior distribution eq. $(2.22)$ for increasing amounts of idealized (BrAC, TAC) data pairs ($R$) all generated using the "true" $\mathit{\boldsymbol{q}}_0$ value previously stated. To calculate these values, MH MCMC samples were drawn with a sample of size 1400 (1500 data points with a 100 sample burn-in phase) where the MCMC sample size was increased from that of Section 3.4.1 due to the increase in noise variance.

We now investigate the results of Section 3.2 with respect to the field-measured (BrAC, TAC) data pairs. Note that we no longer are able to know the true value of the parameters, $\mathit{\boldsymbol{q}}_0$. Surface plots for increasing amounts of subject drinking episode measurements, $R = 1, 26, 76$, and $101$ are contained within Figure 4. Table 3 displays the calculated means and $90\%$ credible circle radii for increasing numbers of subjects, and thus data (corresponding to $R$ as in Section 3.2) included in determination of the prior. To calculate these values, for each $R$, we again used 1400 MH MCMC samples (1500 draws with a 100 sample burn-in phase) generated according to our chosen prior.

3.4.3.   Deconvolution

As in Section 3.3, we rely on the treatment in ^[25] for deconvolving BrAC from TAC using a distribution for $\mathit{\boldsymbol{\mathcal{q}}}$ over $Q$. The chosen distribution was the posterior eq. $(2.22)$ with $n = 3$. To determine the posterior we elected to investigate the case where a non-informative, or what is more aptly described as an uneducated, prior was used. Thus we chose a prior of a tMVN random variable with bounds $[0.01, 10] \times [0.01, 10]$, and parameters $\mu = \left(\begin{array}{l} 5\\ 5 \end{array} \right)$ and $\Sigma = \left(\begin{array}{l} 0.7 & 0.1\\ 0.1 & 0.55 \end{array} \right)$. The noise used was distributed as $\mathcal{N}(0, 0.025^2)$. Note that this choice in prior also highlights the effects of data on the posterior by not providing any initial information to the posterior. When comparing this choice in noise standard deviation to the peak TAC values of our training dataset, we had a typical peak TAC to noise ratio of 8. Further, for the subspaces from Section 3.3 we set our discretization to be $n = 3$, time discretization as $m = 1300$, and discretized $Q$ with $m_1 = m_2 = 20$. As with Section 3.4.2, the noise distribution is meant to simulate a situation where little is known about external effects that play a role in determining noise, and so the data is assumed noisy.

In all of the numerical results presented and discussed in this section, the test dataset used consisted of five drinking episodes from four different participants. These drinking episodes were chosen heuristically so that the test dataset had two drinking episodes with peak BrAC greater than peak TAC, two drinking episodes with peak BrAC less than peak TAC, and one drinking episode with peak BrAC within 0.015 of peak TAC (deemed, "close"). The remaining drinking episodes were used as training data with the added restriction that whenever the desired number of training sets to be used was not too large, BrAC/TAC pairs from any participant who had a dataset included in the selected test data, would be excluded from being among the data used for determining the posterior. The primary exception to this restriction being Figure 6c, wherein we allowed all data that wasn't the current test data point to be included in the training set.

By linearly interpolating the BrAC and TAC data for each subject in all test and training datasets, we are able to re-sample our data with sampling interval $\tau = 45$ seconds, and the time discretization $m = 1300$ previously mentioned. The associated participant IDs, TAC device placement (left vs. right arm), type of drinking pattern used (single, dual, or steady), and number of subjects used in posterior distribution determination ($R$) are labeled in Figures 5 and 6. As in eq. $(3.14)$, we utilized all available non-test subject drinking episode measurements ($R = 136$) to determine population parameters $(r_1^{\ast}, r_2^{\ast})$ to be $(4.7733, 1.7020)$.

Figure 5, shows varying deconvolution attempts for three test data participants, whereas Figure 6 shows deconvolution attempts for the same test data participant reading (BT333), with varying amounts of training subject data within the posterior eq. $(2.22)$, $R = 25, 75,145$. In both Figures 5 and 6, gray bands represent 90% error regions that are determined by sampling respective parameter posterior distributions and utilizing these samples with (3.13) to determine estimated BrAC values. It follows that these error regions contain the 90% credible regions for the pointwise BrAC values as functions of the population parameters appearing in the model. This is the basis for our referring to them in what follows as conservative credible bands.

4.   Discussion

4.1.   Bayesian estimation of model parameters

Figure 2 illustrates rapid convergence in dimensionality of our spatial dimensions as $n$ grows, thus bolstering the results of Theorem 3.1. Within two steps ($n = 3$), we have a graph that visually differs from that of $n = 25$ in ways barely perceptible. Paired with the credible circles in Table 1, these provide evidence that after $n = 3$ the mean and radius of the $\mathit{\boldsymbol{q}}$ credible circles stay consistent. Thus one can choose a computationally efficient $n$ value that minimizes data lost when projecting eq. $(2.6)$ into finite dimensions, eq. $(2.11)$.

For the consistency results, Table 2 exemplifies the theoretical prediction in Theorem 3.5 that as the amount of subject data $R$ grows, the posterior distribution better predicts the true $\mathit{\boldsymbol{q}}$ value by localizing the true parameter $\mathit{\boldsymbol{q}}_0$ in mean with higher confidence (smaller credible circles). This increasing confidence is backed by the decreasing variance results shown in Figure 4. Notice that although the variance decreases, the mean is allowed to shift as more data are incorporated, as evident from comparing Figure 4c to Figure 4d. This shifting mean is permitted by the theoretical results and is likely due to the incorporation of 70 extra data points. Table 3 displays the shifting of the mean as more data are incorporated while quantitatively displaying a decreasing 90% credible circle radius, as expected.

As a final note, recall that TAC data were collected simultaneously from both the right and left arms of participants. For an investigation into this see ^[32].

4.2.   Deconvolution of BrAC from TAC

In Figure 5a, the deconvolved mean BrAC curve more closely resembles the overall curve of the measured TAC values rather than the desired BrAC, with its increased values towards the latter part of the curve. This is to be expected as the measured TAC plays a role in the Bayesian step, but notice that the severity of the increase in the mean value curve is attenuated when compared to that of the TAC curve (red vs. yellow curves at the five hour mark). A similar phenomenon also appears in Figure 6a. For Figures 6a to 6c, as the number of subject drinking episodes $R$ increases, we find that the mean curve grows towards the actual BrAC curve, an expected convergence phenomenon given the theoretical consistency results from Section 3.2.

Lastly, the 90% conservative credible bands about the deconvolved BrAC curves appear to always have a lower bound of zero. For the upper bound, the extreme case is shown in Figure 5c. These wide ranges in BrAC values allow us to capture the true BrAC value with high probability, but also leave us capturing far more area under the curve than needed. Thus, there are times when our two-step method would falsely signal that the TAC device wearer is far more inebriated than they actually are. This incorrect signaling might be due in part to the quantitative inaccurate readings in Figure 5c, wherein the TAC curve is greater than the BrAC curve. If our (training) data are mainly composed of the other cases (TAC following BrAC at an attenuated rate), then the algorithm will learn to "guess up" when turning the TAC back into BrAC. Lastly, this phenomenon may be due to the use of an uninformed prior as the credible regions in Table 3 do not approach zero. Hence, in the future use of an informed prior may be preferable.

4.3.   Conclusions

We believe that the i.n.i.d. assumption from Section 2.2 (specifically Section 3.2) may not reflect the realities of the data collection method wherein two sensors are worn simultaneously on participants' left and right arms. We are currently investigating the elimination of this i.n.i.d assumption. However, the results from Section 3.4 are quite reasonable and are extremely useful when seeking to use this approach computationally in practice. Further investigation is needed regarding the traveling mean exhibited in the numerical results and how it is related to the non-inclusion of other covariate data (age, height, weight, etc.). This investigation may also be aided by attempting to combine the results of Sections 3.1 and 3.2 and let both the approximating dimension of the kernel, as well as the amount of training data, go to infinity simultaneously.

We also believe that the packaging of all error sources into a single random variable in Section 2.2 may yield larger uncertainties than formulations where many additive errors are considered. Namely, mixed-effects formulations may be utilized in order to separate errors and might lower overall uncertainty. However, the results from Section 3.4 are again quite reasonable, and the usage of mixed-effects formulations can be left as a design choice when considering the main goals and implementations of the PDE model from Section 2.1.

When our approach and results are optimized for use in actual practice, some care will have to be taken in regard to the sampling methods used in Sections 3.4.1 and 3.4.2. If Markov Chain Monte Carlo methods are still the method of choice, then issues such as sample size, convergence of the chains, and randomized chain starting points will need to be taken into account. In addition, a laboratory protocol will be need to be developed to estimate the sensor-dependent values of $q_3$ and $q_4$ that appear in eq. $(2.3)$. As far as the numerical results presented in Section 3.4 are concerned in regard to the values chosen for $q_3$ and $q_4$, they primarily serve to reinforce the theoretical results in Sections 3.1 and 3.2.

Finally, Of primary interest is the direct inversion of BrAC, $u$, given TAC as in eq. $(2.12)$ without the need for a two-step process like that of the method used in this paper. We believe that a hierarchical model paired with a Gaussian Process framework may reduce the problem down to a single step (see, ^[52]). In such a framework, we place a prior on $\mathit{\boldsymbol{q}}$, as well as a function space prior over $u$. In this way, we obtain a method that statistically deconvolves BrAC from TAC while providing a distribution from which we may derive error bars on the estimated BrAC values. We are also currently examining the inclusion of another hierarchical Bayesian model that incorporates covariates in both priors placed over $\mathit{\boldsymbol{q}}$ and $u$. We believe that this will improve the accuracy of our predictions by allowing the use of all available subject and environment data.

Acknowledgments

The authors wish to acknowledge that this research was supported in part by the National Institute for Alcohol Abuse and Alcoholism (NIAAA) under grant no. R01AA026368.

Conflict of interest

The authors declare no conflicts of interest.