Year: 2020

Abstract:
This paper proposes a stochastic variational inference (SVI) method for computing an approximate posterior path measure of a Cox process. These processes are widely used in natural and physical sciences, engineering and operations research, and represent a non-trivial model of a wide array of phenomena. In our work, we model the stochastic intensity as the solution of a diffusion stochastic differential equation (SDE), and our objective is to infer the posterior, or smoothing, measure over the paths given Poisson process realizations. We first derive a system of stochastic partial differential equations (SPDE) for the pathwise smoothing posterior density function, a non-trivial result, since the standard solution of SPDEs typically involves an Itô stochastic integral, which is not defined pathwise. Next, we propose an SVI approach to approximating the solution of the system. We parametrize the class of approximate smoothing posteriors using a neural network, derive a lower bound on the evidence of the observed point process samplepath, and optimize the lower bound using stochastic gradient descent (SGD). We demonstrate the efficacy of our method on both synthetic and real-world problems, and demonstrate the advantage of the neural network solution over standard numerical solvers.

heading: 1 INTRODUCTION
text: Cox processes (Cox, 1955; Cox & Isham, 1980), also known as doubly-stochastic Poisson processes, are a class of stochastic point processes wherein the point intensity is itself stochastic and, conditional on a realization of the intensity process, the number of points in any subset of space is Poisson distributed. These processes are widely used in the natural and physical sciences, engineering and operations research, and form useful models of a wide array of phenomena.
We model the intensity by a diffusion process that is the solution of a stochastic differential equation (SDE). This is a standard assumption across a range of applications (Susemihl et al., 2011; Kutschireiter et al., 2020). The measure induced by the solution of the SDE serves as a prior measure over sample paths, and our objective is to infer a posterior measure over the paths of the underlying intensity process, given realizations of the Poisson point process observations over a fixed time horizon. This type of inference problem has been studied in the Bayesian filtering literature (Schuppen, 1977; Bain & Crisan, 2008; Särkkä, 2013), where it is of particular interest to infer the state of the intensity process at any past time given all count observations till the present time instant (the resulting posterior is called the smoothing posterior measure).
In a seminal paper, Snyder (1972) derived a stochastic partial differential equation (SPDE) describing the dynamics of the corresponding posterior density for Cox processes. The solution of this smoothing SPDE requires the computation of an Itô stochastic integral with respect to the counting process. It has long been recognized (Clark, 1978; Davis, 1981; 1982) that for stochastic smoothing (and filtering) theory to be useful in practice, it should be possible to compute smoothing posteriors conditioned on a single observed sample path. However, Itô integrals are not defined pathwise and deriving a pathwise smoothing density is remarkably hard. 30 years after Synder’s original work Elliott & Malcolm (2005) derived a pathwise smoothing SPDE in the form of a coupled system of forward and backward pathwise SPDEs. Nonetheless, solving the system of pathwise SPDEs, or sampling from the corresponding SDE, is still challenging and intractable in general. It is well known, for example, that numerical techniques for solving these SPDEs, such as the finite element
method (FEM), suffers from the curse of dimensionality (Han et al., 2018). Therefore, it is of considerable interest to find more efficient methods to solve the smoothing SPDE.
We take a variational inference approach to computing an approximate smoothing posterior measure. Variational representations of Bayesian posteriors in stochastic filtering and smoothing theory have been developed in considerable generality; see (Mitter & Newton, 2003) for a rigorous treatment. There are a number of papers that consider the computation of an approximate posterior distribution over the paths of the underlying intensity process that is observed with additive Gaussian noise (Archambeau et al., 2007; 2008; Cseke et al., 2013; Susemihl et al., 2011; Sutter et al., 2016). Susemihl et al. (2011) studied Bayesian filtering of Gaussian processes by deriving a differential equation characterizing the evolution of the mean-square error (MSE) in estimating the underlying Gaussian process. On the other hand, Sutter et al. (2016) compute a variational approximation to the smoothing posterior density when the underlying diffusion intensity is observed with additive Brownian noise. They choose their variational family to be a class of SDEs with an analytically computable marginal density. This setting is considerably different from our setting, where the observed process is a point process. Nonetheless, Sutter et al. (2016) provides methodological motivation for our current study. In the context of the computation of approximate smoothing/filtering posteriors for point process observations, Harel et al. (2015) developed an analytically tractable approximation to the filtering posterior distribution of a diffusion modulated marked point processes under specific modeling assumptions suited for a neural encoding/decoding problem. In general, however, analytical tractability cannot be assured without restrictive assumptions.
We present a stochastic variational inference (SVI) (Hoffman et al., 2013) method for computing a variational approximation to the smoothing posterior density. Our approach fixes an approximating family of path measures to those induced by a class of parametrized SPDEs. In particular, we parametrize the drift function of the approximating SPDEs by a neural network with input and output variables matching the theoretical smoothing SPDE. Thereafter, using standard stochastic analysis tools we compute a tractable lower bound to the evidence of observing a sample path of count observations, the so-called evidence lower bound (ELBO). A sample average approximation (SAA) to the ELBO is further computed by simulating sample paths from the stochastic differential equation (SDE) corresponding to the approximating SPDE. Finally, by maximizing the ELBO, the neural network is trained using stochastic gradient descent (SGD) utilizing multiple batches of sample paths of count observations. Note that each sample path of the count observations entails the simulation of a separate SDE. We note that there are many problems in the natural and physical sciences, engineering and operations research where multiple paths of a point process (over a finite time horizon) may be obtained. For instance, we present an example in Section 5 modeling the demand for bikes rented during a 24 hour, one day time period in a bike-sharing platform, where the underlying driving intensity is subject to stochastic variations, and demand information is collected over multiple days.
In contrast to the variational algorithm developed in Sutter et al. (2016), where the variational lower bound must be re-optimized for new sample paths of the observation process, our variational method is more general and our approximation to the smoothing posterior can be used as a map for another (unobserved) sample path of count observations. Our computational approach can also be straightforwardly adapted to solve the problem of interest in Sutter et al. (2016).
In the subsequent sections, we describe our problem and method in detail and demonstrate the utility of our method with the help of numerical experiments. In particular, we show how the choice of approximating family enables us to use the trained neural network and in turn, the variational Bayesian smoothing posterior (VBSP), to compute smoothing SPDE in almost (3/4)th of the computational time required to compute the original smoothing SPDE using FEM. Moreover, we also efficiently generate Monte Carlo samples from the learned VBSP and use them for inference on the bike-sharing dataset, whereas FEM failed to compute either VBSP or the true smoothing density for the given time-space discretization.

heading: 2 PROBLEM DESCRIPTION
text: Let Nt be a Cox process with unknown stochastic intensity {zt ∈ R+, t ∈ [0, T ]}. We use Nt′,t to denote a sample path realization of Nt restricted to the interval [t′, t], and useNt to denoteNt−N0; recall that N0 = 0 by definition. As noted before, a Cox process conditioned on the intensity is a
Poisson process. Therefore, given a realized sample path {zt, t ∈ [0, T ]} of the intensity, and for any 0 ≤ t′ < t ≤ T , the marginal likelihood of observing Nt −Nt′ ∈ N counts in (t′, t] is
Nt −Nt′ ∼ L(Nt −Nt′ = Nt −Nt′ |{zs}t′<s≤t) := (∫ t
t′ zsds )Nt−Nt′ e− ∫ tt′ zsds (Nt −Nt′)! , (1)
where L denotes the Poisson likelihood. Rather than directly modeling the intensity z, we will bring a little more flexibility to our setting, and assume that zt is a deterministic transformation of an another stochastic process xt through a known mapping h : Rd 7→ R+: is zt = h(xt). Note that the non-negative range of h ensures that the Poisson intensity zt = h(xt) is non-negative. Unless xt ∈ R+, the mapping h cannot be an identity function. We use the term intensity process to refer to either zt or xt.
We model the intensity process {xt ∈ Rd,∀t ∈ [0, T ]} with the following SDE,
dxt = b(xt)dt+ σ(xt)dBt,∀t ≤ T and x0 = 0, (2)
where b : Rd 7→ Rd is the drift function, σ(·) : Rd 7→ Rd×d is the diffusion coefficient, and Bt is the d−dimensional Brownian motion (or Wiener process). We assume that there exists a strong solution to the SDE above (Oksendal, 2013, Chapter 5). Moreover, we assume that b(·), h(·), and σ(·) are fixed by the modeler apriori, and we are interested in inferring the unknown intensity process with their fixed definitions. Incorporating them will obscure our main contribution, and we leave it for future work.
The model of the count observations above forms a diffusion modulated Cox process. Diffusion modulated Cox processes are widely used to model the arrival process in various service systems such as call centers, hospitals, airports etc. (Zhang et al., 2014; Wang et al., 2020). Zhang & Kou (2010) use a Gaussian process modulated Cox-process to infer proteins’ conformation, in particular, they model the arrival rates of the photons collected from a laser excited protein molecule as a Gaussian process. Schnoerr et al. (2016) model spatio-temporal stochastic systems from systems biology and epidemiology using Cox process where intensity is modelled with diffusions.
As stated in the introduction, we seek to infer the smoothing posterior measure over the unknown intensity process {xt, t ∈ [0, T ]} using the count observations upto time T . Following terminology from the Bayesian filtering theory (Särkkä, 2013), we use smoothing to refer to inferring the unobserved intensity process at any past time given the observations upto the current time. Mathematically, the smoothing posterior is defined using the conditional expectation of the form E[f(xt)|s(Nu, u ∈ [0, T ])], where s(Nu, u ∈ [0, T ]) is the smallest sigma algebra (or filtration) generated by the Cox process {Nt} from time 0 to T . For brevity we write E[f(xt)|s(Nu, u ∈ [0, T ])] as E[f(xt)|N0,T ]. Interested readers may refer to Kutschireiter et al. (2020) for more details on non-linear filtering theory.
We now provide a formal derivation of the smoothing posterior using Bayes’ theorem (Bain & Crisan (2008); Elliott & Malcolm (2005)). Observe the conditional expectation satisfies
E[f(xt)|N0,T ] = E†[Λ0,T f(xt)|N0,T ]
E†[Λ0,T |N0,T ] (3)
for any measurable function f(·) and Λs,t := L(Ns,t)L†(Ns,t) for any 0 ≤ s < t ≤ T , where L † is the unit intensity Poisson likelihood and E†[·] denotes the expectation with respect to L†. Note that L† does not depend on the stochastic intensity process x and forms a reference measure. The marginal smoothing posterior density is defined as
pt(x|N0,T ) := P(xt ∈ dx|N0,T ), (4)
which can be formally obtained from equation 3 by setting f(xt) = I{A}(xt) for any A ∈ Rd, where I{A}(y) is an indicator function that equals 1 when y ∈ A, otherwise 0. Now, define the unnormalized filtering density function q̄t(x) as the function satisfying
P(xt ∈ dx|N0,t) = q̄t(x)dx∫ Rd q̄t(ξ)dξ , (5)
and also define v̄t(x) := E†[Λt,T |N0,T ]. Then, it can be shown (Elliott & Malcolm (2005)) that for any measurable function f ,
E[f(xt)|N0,T ] = E†[Λ0,T f(xt)|N0,T ]
E†[Λ0,T |N0,T ] =
∫ Rd f(ξ)q̄t(ξ)v̄t(ξ)dξ∫
Rd q̄t(ξ)v̄t(ξ)dξ . (6)
Next, recalling that h(·) is the mapping to ensure the intensity process is positive, define the function Ψt for a given sample path of count observations (i.e., pathwise) as
Ψt := Ψ(h(x), t, Nt) = exp [(1− h(x))t+Nt log h(x)] ,∀x ∈ Rd. Following Elliott & Malcolm (2005, Theorem 4) one may use Ψt to derive a coupled system of pathwise SPDEs that characterize q̄t(x) and v̄t(x). In particular, they show that qt = Ψ−1t q̄t is a solution to the following SPDE
∂tqt(x) = Ψ −1 t L ∗[Ψtqt(x)],∀t ≤ T, q0(x) = δx0(x), (7) where L∗ is the adjoint of L[F (x)] = 12 ∑ i,j ai,j(x)∂xixjF (x) + ∑ i bi(x)∂xiF (x), which is the infinitesimal generator of the prior process for any twice-differentiable, continuous, and bounded function F : Rd 7→ R and a(x) = σ(x)σ(x)T , and δx0(x) is the Dirac delta distribution at x0. Moreover, they also show that vt(x) = Ψtv̄t(x) satisfies the following backward parabolic equation
∂tvt(x) = −ΨtL[Ψ−1t vt(x)], (8) with terminal condition vT (x) = ΨT (x).
Now it follows from equation 6 that using the solution of these two SPDEs, the marginal smoothing posterior density for any t ∈ [0, T ] satisfies
pt(x|N0,T ) = qt(ξ)vt(ξ)dξ∫ Rd qt(ξ)vt(ξ)dξ . (9)
Using the SPDEs in equation 7, and 8, together with 9, it can be shown that the marginal smoothing posterior density pt(x|N0,T ) satisfies its own SPDE: for any t ∈ [0, T ],
∂tpt(x|N0,T ) = − ∑ i ∂xi [ { (a(x)[∇ log(Ψ−1t vt(x))])i + bi(x) } pt(x|N0,T ) ] + 1
2 ∑ i,j ∂xixjai,j(x)pt(x|N0,T ) (10)
and pt(x|N0,T ) = δx0(x) with x0 = 0. We present a detailed derivation in Appendix A.1. Corresponding to this SPDE, there exists a smoothing posterior SDE, defined as
dx̄t = { a(x̄t)[∇ log(Ψ−1t vt(x̄t))] + b(x̄t) } dt+ σ(x̄t)dB̄t and x̄0 = 0, (11)
where {x̄t} is a modification of the process {xt} such that B̄t is independent of the Cox process Nt (and thus Bt).
Observe that the entire sample path of the count observations N0,T is summarized through the pathwise function Ψt and vt together in the drift term of this SDE. Also note that the diffusion coefficient of the smoothing posterior SDE is precisely the same as that of the prior SDE. The computation of the drift term in the smoothing posterior SDE requires the solving equation 8 for vt(x) which, in turn, makes the posterior computation challenging and computationally intractable in general. Consequently, the computation of the marginal posterior density (and hence the path measure). Therefore, we propose a variational inference-based method to compute an approximation to the solution of the smoothing posterior SPDE, by computing an approximate solution to the smoothing posterior SDE in equation 11.

heading: 3 VARIATIONAL BAYES FOR APPROXIMATING THE SMOOTHING DENSITY
text: Observe that the posterior path measure is the minimizer of the following variational optimization problem (Mitter & Newton (2003, Proposition 2.1)),
max Π(x0,T )∈P(C)
{ KL(Π‖Π0) + ∫ dΠ(x0,T ) logL(N0,T |h(x0,T )) } , (12)
where P(C) is the space of all absolutely continuous measures with respect to Π0, the measure induced by a solution of the intensity SDE (equation 2) on the space C[0, T ] of continuous functions with support [0, T ], and KL denotes the Kullback-Leibler divergence between two absolutely continuous measures.Note that P(C) also contains Π0. Solving this optimization problem over all measures in P(C) is intractable. Therefore, we choose the subset of absolutely continuous measures Qb̄ ⊂ P(C) induced by solutions of the following SDE:
dxt =b̄(xt, Nt, t)dt+ σ(xt)dBt, for t ≤ T and x0 = 0, (13)
where b̄(·, ·, ·) : Rd × N × [0, T ] 7→ Rd is the drift function. We term this space of measures Qb̄ as the variational family. The measures in Qb̄ are absolutely continuous with respect to Π0 as they have the same diffusion coefficient, therefore Qb̄ ⊂ P(C). This choice of the variational family is not arbitrary, but rather motivated by the smoothing posterior SDE derived in equation 11, where the diffusion coefficient is σ(·) and the drift coefficient has an intractable form (that depends on the prior drift b(·) and diffusion coefficient σ(·), and Nt through Ψt and vt(·)). Notice that the choice of drift function spans the space of measures in the variational family Qb̄. Since, Qb̄ ⊂ P(C), it follows from equation 12 that
max Π(x0,T )∈P(C)
{ −KL(Π‖Π0) + ∫ dΠ(x0,T ) logL(N0,T |h(x0,T )) } ,
≥ max Q∈Qb̄
{ −KL(Q‖Π0) + ∫ dQ(x0,T ) logL(N0,T |h(x0,T )) } (14)
The right hand side above is known as the evidence lower bound (ELBO). The corresponding ELBO maximization problem to compute the optimal Q ∈ Qb̄ (for a given sample path N0,T ) is simply
Q∗(·|N0,T ) = arg max Q∈Qb̄ EQ {logL(N0,T |x0,T )} − KL(Q‖Π0). (15)
Note that absolutely continuous measures on path space correspond to changes in the drift function, for a fixed diffusion coefficient (else, the measures are singular). As a consequence of Girsanov’s theorem (Oksendal, 2013, Theorem 8.6.8) (see Appendix A.2 for the proof) we have,
KL(Q‖Π0) = 1
2 EQ [∫ T 0 ‖σ−1(xt)(b(xt)− b̄(xt, Nt, t))‖2dt ] , (16)
where recall b(·, ·) is the drift of the prior SDE defined in equation 2 and b̄(·, ·, ·) is the drift of the variational SDE. Substituting this into equation 15 yields
Q∗(·|N0,T ) = arg max Q∈Qb̄ EQ
{ logL(N0,T |x0,T )− 1
2 ∫ T 0 ‖σ−1(xt)[b(xt, t)− b̄(xt, Nt, t)]‖2dt } .
(17)
We denote Q∗(·|N0,T ) as the variational Bayesian smoothing posterior (VBSP) path measure . Next, we lay down the details of the SVI algorithm to solve the above optimization problem to compute the VBSP.

heading: 4 STOCHASTIC VARIATIONAL INFERENCE OF THE VBSP
text: It is evident from the ELBO in equation 17 and the definition of the variational familyQb̄ that computing the VBSP measure entails the computation of the unknown drift function b̄(·, ·, ·) in equation 13. We further restrict the family of measures Qb̄ by assuming the drift functions belong to a class of parametrized, smooth functions. A feasible way to model this class of drift functions is through a neural network. However, it is possible to use simpler approximation function classes as done in, for example, Sutter et al. (2016), who fixed b̄ to ensure that the marginal distributions of the variational smoothing SDE belong to a specific exponential family of distributions. We note
that in choosing a parametrized class, we must still ensure that the resulting drift functions are Lipschitz continuous and satisfy sufficient regularity so that a solution to the SDE equation 13 exists. Furthermore, restricting the drift functions in this way entails a further restriction of the class of (approximating) path measures. An open question is here is how much of a loss this entails (in terms of the Kullback-Leibler divergence from the ‘true’ posterior path measure in this instance).
To fix the idea, we assume that b̄(·, ·, ·) lies in a general class of functions parametrized by θ. Henceforth, we write the drift coefficient as b̄(·, ·, ·, θ) to make its dependence on θ explicit. We use stochastic gradient descent (SGD) to maximize the ELBO, which requires the computation of stochastic gradients of the ELBO with respect to θ. To compute the gradients, we generate sample paths of xθ, the solution of the variational SDE equation 13 for a given θ, using a first-order EulerMaruyama integration of the SDE. We do this for convenience, though higher order approximations could be used. Specifically, we partition the time interval [0, T ] in K equal sub-intervals of length ∆ = T/K, {t0, t1, . . . tK}, where t0 = 0 and tK = T and then generate the sequence of {xθti} K i=1 using the following recursive equation and initial condition xθt0 = 0:
xθti = x θ ti−1 + b̄(x θ ti−1 , Nti−1 , ti−1, θ)∆ + ∆σ(x θ i−1)Zi,∀ i ∈ {1, 2, . . .K} (18)
where {Zi}Ki=1 is the sequence of K independent and identically distributed (i.i.d.) d−dimensional standard Gaussian random vectors.
We generate M independent sample paths of the discrete-time process in equation 18 denoted as {xθ,mti }, for m ∈ {1, 2, . . .M}, to compute a sample average approximation (SAA) of the ELBO over the partition {t0, t1, . . . tK} as
ÊLBO = 1
M M∑ m=1 [ K−1∑ i=0 [ logL(Nti,ti+1 |h(x θ,m ti )∆)− 1 2 ‖σ(xθ,mti ) −1[b(xθ,mti , ti)− b̄(x θ,m ti , Nti , ti, θ)]‖ 2∆ ]]
= 1
M M∑ m=1 [ K−1∑ i=0 [ Nti,ti+1 log(h(x θ,m ti )∆)− h(xθ,mti )∆− 1 2 ‖σ(xθ,mti ) −1[b(xθ,mti , ti)− b̄(x θ,m ti , Nti , ti, θ)]‖ 2∆ ]] + C,
(19)
where we used the definition of the Poisson likelihood L from equation 1 and C is a constant independent of θ. Now, to compute gradients of ELBO with respect to θ, observe that the gradient operator can be exchanged with the expectation, since the only source of randomness in each sample path of xθ,mt are K i.i.d. Gaussian random vectors {Zmi }Ki=1, which are independent of θ. Notice too that this is a pathwise analog of the reparametrization trick, and has been used recently to learn deep latent models (Tzen & Raginsky, 2019; Li et al., 2020). Also, note that thus far the ÊLBO is defined for a single sample path of the count observations N0,T . However, as noted before, in our method we will also take a sample average of ÊLBO over multiple batches of sample paths of count observations at each epoch of the training algorithm.

heading: 5 NUMERICAL EXPERIMENTS
text: We present three experiments demonstrating the efficacy and utility of our proposed SVI method. First, we consider a setting when the underlying stochastic intensity process is 1-dimensional. We compare the SVI approximation with the ‘true’ smoothing posterior density computed using the solution of the forward and backward SPDEs, defined in equation 7 and equation 8. We solve the SPDEs using a standard finite element method (FEM) solver (Skeel & Berzins, 1990, Matlab solvers for 1-D PDEs). In a second experiment, we demonstrate the performance of our algorithm on a subset of a Bike-sharing dataset obtained from the UCI machine learning repository (Fanaee-T & Gama, 2013). In this experiment we estimate a smoothing posterior density for the observed counts of the demand for bikes in a 24 hour period, assuming that the demand process is well-modeled by a Cox process. In our third and final experiment, we apply our method to compute an approximation to a 4- dimensional smoothing posterior density. We note that despite being low dimensional, the standard FEM solver does not scale to this setting, while our method can be straightforwardly adapted.

heading: 5.1 VARIATIONAL APPROXIMATION OF UNIVARIATE SMOOTHING POSTERIOR DENSITY
text: As defined in Section 2, we are interested in learning the posterior measure over an unknown process {xt ∈ R} where the intensity process {zt} satisfies zt = h(xt). We set b(x) = −x and σ(x) = 1 in the prior SDE as defined in equation 2. Furthermore, motivated from the mathematical structure of the true smoothing SDE equation 11, we fix our variational family Q to be the class of measures induced by solutions to the class of SDEs in equation 13 with drift and diffusion coefficient set as
b̄(x, t,Nt) = − Ψ′t Ψt + V (x, t,NT −Nt; θ)− x, σ(x) = 1,
where Ψ′t is the derivative of Ψt with respect to x. Here V (x, t,NT − Nt; θ) is modeled using a neural network with 2 hidden layers whose parameters we call θ (see Appendix A.3 for more details on the architecture).

heading: 5.1.1 SIMULATED DATASET
text: We generate sample paths of the count observation N0,T from a non-homogeneous Poisson process, where the intensity process {z0t } is the solution of the following ordinary differential equation, dz0t dt = 20(2− t) exp(−0.85(2− t)
2), and z00 = 0. We fix the map h(a) = 5 ∗ exp(−.08 ∗ (a− 5)2) in this experiment. To train the neural network V , we use 150 samples paths of the count observation between time interval [0, 2] and optimize ÊLBO in equation 19 using Adam (Kingma & Ba, 2014).
To demonstrate the efficacy of our approach, we first generate 20 test sample paths of count observations and compute the true smoothing posterior density (defined in equation 9) using the solution of the forward and backward SPDEs (defined in equation 7 and equation 8 respectively). We do this using the FEM method. Then for the same test observations we compute the VBSP density by FEM using the trained drift coefficient of VBSP SDE (see equation 10 and 11).
Comparative results are presented in Figure 1. We clearly see from the the top row plot that these are very similar to each other, with our variational approximation capturing the sharp rises in density with high fidelity. Note from the first two plots in the bottom row that as the ELBO decreases the gap between the the true and VB smoothing posterior reduces too. Moreover, the time required to compute the smoothing density using the forward and backward SPDEs on the test data is about 2.2 seconds which is approximately fifty percent higher than the trained VBSP density, which required 1.6 seconds (on an 3.1 GHz Intel i5 CPU). This is due to the fact that we are required to solve one SPDE in the latter case instead of two in the former case.
Notice that the learnt drift of the smoothing SDE equation 13 is a map which can be used with any sample path of count observations to generate Monte Carlo samples from the an approximate smoothing posterior density. In contrast, it is challenging to sample from the true smoothing SDE as it involves computing the solution of the system of SPDEs in equation 7 and 8. Furthermore, this solution must be recomputed for each new sample path of count observations.

heading: 5.1.2 BIKE SHARING DATASET
text: In this experiment, we compute the VBSP density for the hourly counts of demand in a bike-sharing system. The experimental setting remains unaltered, albeit the diffusion coefficient is set to σ(x) = 1.1, to capture the increased variability in counts of bike demand than the variability in simulated counts in the previous experiment. Notice that fixing the map h is a modeling question, and we consider mappings of the form h(x) = a exp(−b(x − c)2) parametrized by a, b, and c. We take a simple empirical Bayes heuristic to fix their values after observing the count observations, based on the fact that h(x)∆ is the mean of the Poisson counts in the interval ∆. We set a in such a way that h(x)∆ equals the maximum of the median observed count. For the Bike-sharing data we re-scaled the problem to interval [0, 2] and fixed ∆ = 0.083 and thus choose 90/0.083 ' 1050 as a. The choice of b and c depends on the diffusion coefficient σ(x) and x0, as the SDE should be able to explore the relevant domain of h to appropriately model the actual count observations. Thus, after looking at the count data, we chose h(x) = 1050 ∗ exp(−0.001 ∗ (y − 50)2). The empirical results are summarized in Figure 2. We note here that the FEM approach (our implementation) to compute the VBSP density was numerically unstable and failed; this may be attributed
to the concentrated nature of the VBSP, as visible in the plots. It is evident from the last plot in Figure 2 that for the current h, VBSP density captures the trend in the real count observations well, however we anticipate that a different choice of the map h can better model the count observations.
Note that in a smoothing problem, the prior intensity process (specifically the drift coefficient b(·), diffusion coefficient σ(·) and the map h(·)) are sourced from an expert, and the objective is to update the modeler’s beliefs with count observations to compute the smoothing posterior density. In many settings, the functions b(·), σ(·) and h(·) are known only up to some unknown parameters. It is fairly straightforward to combine the parameters of h and b with the neural network parameters θ, and learn them all in a data-driven manner. We choose not to do this to keep the discussion simple. Learning the parameters of σ presents a slightly greater challenge since the path measures for different settings of σ are singular. This is a crucial difference from the finite dimensional setting, where the Lebesgue measure is a common reference. We leave solving this problem for future work.

heading: 5.2 VARIATIONAL APPROXIMATION OF MULTIVARIATE SMOOTHING POSTERIOR DENSITY
text: We demonstrate our method on a 4-dimensional smoothing problem. In this case, we fix the map h(a) = 25 ∗ exp(−.08 ∗ ‖a − 5‖2), where ‖ · ‖ is the L-2 norm, and a ∈ R4. We also choose the prior density to be induced by an SDE defined in equation 2 with b(x) = −x and σ(x) = I , where x ∈ Rd and I is a d × d identity matrix. Furthermore, we choose our variational family to be a family of SDE as defined in equation 13, with drift and diffusion coefficients, b̄(x, t, Nt) = −∇ΨtΨt + Vd(x, t, NT −Nt; θ)− x and σ(x) = I .
To train the neural network V , we use 150 samples paths of the count observation between time interval [0, 2] and optimize the ELBO defined in equation 19 using Adam (Kingma & Ba, 2014). We plot the empirical results in Figure 3.

heading: A APPENDIX
text: A.1 DERIVATION OF SMOOTHING SPDE According to Theorem * [sic] in Elliott & Malcolm (2005), Kt := ( ∫ Rd qt(ξ)vt(ξ)dξ)
−1 is almost surely constant for t ≤ T . Using this result it follows that
∂tpS,t =Ktqt(x)∂tvt(x) +Ktvt(x)∂tqt(x)
=Ktqt(x)∂tvt(x) +Ktvt(x) [ Ψ−1t L ∗[Ψt pS,t
Ktvt(x) ] ] =− pS,t
vt(x) ΨtL[Ψ
−1 t vt(x)] + vt(x) [ Ψ−1t L
∗[Ψt pS,t vt(x) ] ] =− pS,t
Vt(x) L[Vt(x)] + Vt(x)
[ L∗[
pS,t Vt(x) ]
] . (20)
where Vt(x) = Ψ−1t vt(x) and PS,t = pt(x|N0,T ) are introduced for brevity. Now observe that
Vt(x)L ∗ [ pS,t Vt(x) ] = Vt(x) 1 2 ∑ i,j ∂xixj ( ai,j(x) pS,t Vt(x) ) − ∑ i ∂xi ( bi(x) pS,t Vt(x) ) = Vt(x) 1
2 ∑ i,j ∂xixj ( ai,j(x) pS,t Vt(x) ) + ∑ i [ ∂xiVt(x)bi(x)pS,t Vt(x) ] − ∑ i [(∂xibi(x)pS,t + bi(x)∂xipS,t)] . (21)
Consider the summand in the first term of equation 21 and observe that
∂xixj [ ai,j(x)pS,t Vt(x) ] = ∂xi [ Vt(x)∂xj (ai,j(x)pS,t)− ∂xjVt(x)ai,j(x)pS,t V 2t (x) ] = Vt(x)∂xixj (ai,j(x)pS,t)− ∂xiVt(x)∂xj (ai,j(x)pS,t)
V 2t (x)
+ 2∂xiVt(x)
V 3t (x)
[ ∂xjVt(x) (ai,j(x)pS,t) ] − 1 V 2t (x) ∂xi ( ∂xjVt(x)ai,j(x)pS,t ) Now the second term in equation 21 can be expressed as
Vt(x)∂xixj [ ai,j(x)pS,t Vt(x) ] = ∂xixj (ai,j(x)pS,t)− ∂xiVt(x)∂xj (ai,j(x)pS,t) Vt(x)
+ 2∂xiVt(x)
V 2t (x)
[ ∂xjVt(x) (ai,j(x)pS,t) ] − 1 Vt(x) ∂xi [ ∂xjVt(x)ai,j(x)pS,t ] (22)
Substituting the above expression in equation 21, we have Vt(x)L ∗ [ pS,t Vt(x) ] = 1 2 ∑ i,j { ∂xixj (ai,j(x)pS,t)− ∂xiVt(x)∂xj (ai,j(x)pS,t) Vt(x)
+ 2∂xiVt(x)
V 2t (x)
[ ∂xjVt(x) (ai,j(x)pS,t) ] − 1 Vt(x) ∂xi [ ∂xjVt(x)ai,j(x)pS,t ] } + ∑ i [ ∂xiVt(x)bi(x)pS,t Vt(x) ] − ∑ i [(∂xibi(x)pS,t + bi(x)∂xipS,t)] . (23)
Now consider the first term in the RHS of equation 20
− pS,t Vt(x) L[Vt(x)] =− 1 2 ∑ i,j ai,j(x)pS,t Vt(x) ∂xixjVt(x)− ∑ i bi(x)pS,t Vt(x) ∂xiVt(x). (24)
Now substituting equation 23 and equation 24 into equation 20, we obtain
∂ ∂t pS,t = − pS,t Vt(x) L[Vt(x)] + Vt(x)
[ L∗[
pS,t Vt(x) ] ] =− 1
2 ∑ i,j ai,j(x)pS,t Vt(x) ∂xixjVt(x)− ∑ i bi(x)pS,t Vt(x) ∂xiVt(x) + 1 2 ∑ i,j { ∂xixj (ai,j(x)pS,t)
− ∂xiVt(x)∂xj (ai,j(x)pS,t)
Vt(x) +
2∂xiVt(x)
V 2t (x)
[ ∂xjVt(x) (ai,j(x)pS,t) ] − 1 Vt(x) ∂xi [ ∂xjVt(x)ai,j(x)pS,t ] } + ∑ i [ ∂xiVt(x)bi(x)pS,t Vt(x) ] − ∑ i [(∂xibi(x)pS,t + bi(x)∂xipS,t)]
= 1
2 ∑ i,j ∂xixj (ai,j(x)pS,t)− ∑ i ∂xi [bi(x)pS,t]− 1 2 ∑ i,j {∂xixjVt(x) Vt(x) [ai,j(x)pS,t]
+ ∂xiVt(x)∂xj (ai,j(x)pS,t)
Vt(x) −
2∂xiVt(x)∂xjVt(x)
V 2t (x) [(ai,j(x)pS,t)]
+ 1
Vt(x)
[ ∂xjVt(x)∂xi(ai,j(x)pS,t) + ∂xixjVt(x)ai,j(x)pS,t ] } = 1
2 ∑ i,j ∂xixj (ai,j(x)pS,t)− ∑ i ∂xi [bi(x)pS,t]
− 1 2 ∑ i,j { ∂xi [ ∂xj log Vt(x)[ai,j(x)pS,t] ] + ∂xiVt(x)∂xj (ai,j(x)pS,t) Vt(x)
− ∂xiVt(x)∂xjVt(x)
V 2t (x) [(ai,j(x)pS,t)] +
1
Vt(x)
[ ∂xixjVt(x)ai,j(x)pS,t ] } Since, ∂xixjVt(x) = ∂xjxiVt(x) therefore
∂ ∂t pS,t = 1 2 ∑ i,j ∂xixj (ai,j(x)pS,t)− ∑ i ∂xi [bi(x)pS,t]
− 1 2 ∑ i,j { ∂xi [ ∂xj log Vt(x)[ai,j(x)pS,t] ] + ∂xj [∂xi log Vt(x)[ai,j(x)pS,t]] } = 1
2 ∑ i,j ∂xixj (ai,j(x)pS,t)− ∑ i ∂xi [bi(x)pS,t]− ∑ i,j { ∂xi [ ∂xj log Vt(x)[ai,j(x)pS,t] ] } = 1
2 ∑ i,j ∂xixj (ai,j(x)pS,t)− ∑ i ∂xi {[(a(x)[∇ log Vt(x)])i + bi(x)]pS,t} . (25)
A.2 KL-DIVERGENCE BETWEEN A MEMBER OF VARIATIONAL FAMILY AND PRIOR SDE
We derive a pathwise expression for KL(Q‖Π0) for a given count observation path N0,T . Theorem A.1. Define u(xt, t, Nt; θ) := σ(xt) −1 ( b(xt, t)− b̄(xt, t, Nt; θ) )
and suppose that u satisfies a strong Novikov’s condition:
E [ exp ( 1
2 ∫ T 0 ‖u(xt, t, Nt; θ)‖2dt )] < +∞ ∀θ, φ.
Then,
KL(Q‖Π0) = EQ
[ 1
2 ∫ T 0 ‖u(xt, t, Nt; θ)‖2ds ] . (26)
Proof. Given samples path of count observationN0,T , using the definition of u and under Novikov’s condition, using Girsanov’s theorem (Oksendal, 2013, Theorem 8.6.8), we have
dQ
dΠ0 = exp
( − ∫ t
0
u(xt, t, Nt; θ)dBs − 1
2 ∫ t 0 ‖u(xt, t, Nt; θ)‖2ds ) ,
and
B̂t := ∫ t 0 u(xt, t, Nt; θ)ds+ B(t) (27)
is a Brownian motion w.r.t. Q. Furthermore, we also have
dxt = b(xt, t)dt+ σ(xt)dB̂t. (28)
The following expression is obtained by substituting for dQdΠ0 : EQ [ log ( dQ(x0:T )
dΠ0(x0:T ) )] = −EQ [∫ T 0 ( u(xt, t, Nt; θ)dBs + 1 2 ‖u(xt, t, Nt; θ)‖2 ) ds ] .
(29)
Now, applying B̂(t) := ∫ t
0 u(xt, t, Nt; θ)ds+ Bt in Equation equation 29, we have
EQ [∫ T 0 u(xt, t, Nt; θ)dBs ] = EQ [∫ T 0 u(xt, t, Nt; θ)[dB̂s − u(xt, t, Nt; θ)ds] ]
= EQ [∫ T 0 u(xt, t, Nt; θ)dB̂s − ∫ T 0 ‖u(xt, t, Nt; θ)‖2ds ]
= EQ [ − ∫ T
0
‖u(xt, t, Nt; θ)‖2ds ] . (30)
Substituting equation 30 into equation 29 yields KL(Q‖Π0) = EQ [ log ( dQ(x0:T )
dΠ0(x0:T ) )] = EQ [ 1
2 ∫ T 0 ‖u(xt, t, Nt; θ)‖2ds ] .
(31)
and thus concludes the proof.
A.3 NEURAL NETWORK ARCHITECTURE
For all the experiments, we use the neural network architecture depicted in Figure A.3, with ReLU activation functions between fully-connected hidden layers.
A.4 COMPARING COMPUTATIONAL TIME REQUIRED TO NUMERICALLY COMPUTE VBSP AND TRUE SMOOTHING DENSITY USING FEM

Review summaries:
The paper under review proposes a variational inference procedure for a specific class of Cox processes whose intensity is derived from a stochastic differential equation. The methodology relies on a restriction of candidate solutions the the subset for which the drift depends on $x_t$, $N_t$ and $t$; the drift is then modelled with a neural network. By simulating from the candidate model, a sample average approximation of the ELBO is used to compute a stochastic gradient descent algorithm, optimize the bound and thus estimate non-parametrically the drift.