Estimating Differential Equations from Temporal Point Processes
Authors: Shuichi Miyazawa
Daichi Mochihashi
Venue: Accepted by TMLR ()

Abstract: Ordinary differential equations (ODEs) allow interpretation of phenomena in various scientific fields. They have mostly been applied to numerical data observed at regular intervals, but not to irregularly observed discrete events, also known as point processes. In this study, we introduce an ODE modeling of such events by combining ODEs with log-Gaussian Cox processes (Møller et al., 1998). In the experiments with different types of ODEs regarding infectious disease, predator-prey interaction, and competition among participants, our method outperformed existing baseline methods assuming regularly observed continuous data with respect to the accuracy of recovering the latent parameters of ODEs. Through both synthetic and actual examples, we also showed the ability of our method to extrapolate, model latent events that cannot be observed, and offer interpretability of phenomena from the viewpoint of the estimated parameters of ODE.

Published in Transactions on Machine Learning Research (09/2023)
Estimating Differential Equations from Temporal Point Pro- cesses
Shuichi Miyazawa
miyazawa
Deportment of Statistical Science, The Groduate Universily for Advanced Studies
Daichi Mochihashi
daichi Waechii@eem.co dam.hi@iem.ac.pp
The Instibute of Statistical Mathematics
Reviewed on OvenReview: https://openrevieet http:p/reeeeeeeeeeooootttaee---e
Abstract
Ordinary differential equations (ODEs) allow interpretation of phenomena in various scien- tific fields. They have mostly been applied to numerical data observed at regular intervals, but not to irregularly observed discrete events, also known point processes. In this study, we introduce an ODE modeling of such events by combining ODEs with log- Gaussian Cox processes (Moller et al. 1998). In the experiments with different types of ODEs regard- ing infectious discase, predator-prey interaction, and competition among participants, our method oulperformed existing baseline methods assuming regularly observed continuous data with respect to the accuracy of recovering the latent parameters of ODEs. Through both synthetic and actual examples, we also showed the ability of our method to extrapo- late, model latent events that cannot be observed, and offer interpretability of phenomena from the viewpoint of the estimated parameters of ODF.
1 Introduction
Ordinary differential equations (ODEs) have been used in various scientific fields to model naturally occurring phenomena (Jones el al.] 2009 Simmons 2016). ODEs can be derived from first principles mathematical models developed by domain expert to explain phenomena (Murray & Murray 2003 Brauer 2017). Using domain knowledge of the phenomena, one can specify the form of the ODE. However, the parameters of the ODE are often unknown. Since these parameters can help understand phenomena, the inverse problem, i.e., estimation of the parameters from observed data, has been studied as one of the main problems in the scientific literature (Kryazhimskiy & Osipov] 1995). Primary data for this purpose have been numerical observations that regularly collected; however, in practice, the data often not collected under such experimental control, limiting the application of ODE modeling. To broaden the application scope of ODEs, here we focus on event data, which are discrete phenomena that occur irregularly in continuous time, also known as discrete events. This type of data has previously not been the subject of ODE modeling. Thus, we address the problem estimating the parameters of ODEs from such discrete data.
Event data are ubiquitous in various domains, from science to industry, e.g., carthquakes, blog posts, tweet posts, and business transactions. The properties of event data are significantly different from those of numerical data commonly used for ODE modeling: event data only consist of the time of their occurrences, and while they may include some metadata, they do not come with numerical information by themselves. More importantly, event data are often observed irregularly, leading to the application of existing ODE modeling methods that assume regularly observed data to be a non-trivial problem.
To model such discrete events, temporal point processes (TPPs) have been studied extensively (Daley Vere-Jones 2007 Yan 2019). TPPs model the probability of such event occurrences using a latent intensity
function; in particular, model where the intensity varies time is called Cox process (Cox 1955). To

---

Published in Transactions on Machine Learning Research (09/2023)
Observations
Observations
Estimated ODE
Estimated ODE
das/ dt = -1.992521,
dzs/dt - 1.882s21,
dzy/ dt = 19992g21-244442
d2y/ dt 1.882221-22214
dz dt 2.40z1.
d dt 2.2121.
Latent intensity
Latent intensity
z z= = exeppx
s
10
R
10
R
5
0.0
0
0.0
0.2
0.4
0.6 0.8 1.0
0.0
0.2
0.4 0.6 0.8 1.0
time
time
Figure 1: Illustaatonn of the proposed method using SIR equations in epidemiology as an ODE. The Param- eters of ODE and latent intensities are estimated from discrete events (point processes) for two cases: when all types of observations are available (left column), and when only type observations are available (right column). The shaded areas represent the 75% quantiles of the estimated intensities.
estimate the parameters from discrete events, this study defines an ODE-guided Poisson process, where the latent intensity is assumed to follow the dynamies of an ODE. An illustrative result of ODE-guided Poisson process is shown in Fig- [ To infer this model from discrete events, we propose inference method that combines the log- Gaussian Cox process (LGCP) (Moller et al. 1998 with the fast Gaussian process-based gradient matching (Wenk et al. 2019), which is an efficient ODE parameter estimation method. Because the posterior probability density distribution of the proposed model is not given in a closed form, an approximate inference in Bayesian framework is necessary. We used the Markov Chain Monte Carlo (MCMC) method, an algorithm widely used approximate posterior inference, to obtain samples of latent variables from complex, high-dimensional probability distributions (Gelman
1995 MacKay 2003).
Recently, Neural ODEs CChen et al. 2018) have been proposed for flexible time series modeling using deep neural networks with ODEs. Subsequently, extensions of Neural ODE incorporating event modeling have been proposed (Jia & Benson 2019 Rubanova et al. 2019 Chen et al. 2020). Although these studies aim to model discrete events as well as this study, they employ implicit ODEs that do not have an explicit form, thus prechuding the interprctability of the data from their parameters.
In contrast, the parameters of explicit ODEs allow the interpretation of phenomena related event data, which is the aim of this study, Most existing ODE parameter estimation methods assume a discrete time model with regularly observed continuous data; however, TPPs continnous time models that handle irregularly observed discrete events. The main issue is to fill such gap between the parameter estimation of ODEs and TPPs, providing new flexibility and applicability in ODE modeling, In this study, we
define an ODE'-gided Poisson process, which is a generative model of discrete event sequences governed by an explicit ODE,
develop an approximate inference method for ODE-guided Poisson process by incorporating LGCP with ODE parameter ostimation mothod,
firstly demonstrate on synthetic data simulated by using ODEs regarding infectious disease, predator- prey interaction, and competition, that our method outperforms bascline ODE parameter estimation methods developed for regularly observed numerical data with respect to the accuracy of estimating parameters of ODE, showing the effectiveness of combining TPPs with ODE parameter estimation,
2

---

Published in Transactions on Machine Learning Research (09/2023)
show that our method improves extrapolation performance by introducing mechanistic knowledge through ODFs and can be used even when some types of events are unavailable as shown in the right column of Fig. D
apply the proposed method to real-world data to present examples of modeling the generation process of event data with ODEs, and interpreting phenomena from the viewpoint of estimated ODE parameters.
2 Backgrounds
2.1 Modeling with Differential Equations
Consider a dynamical system of ODE consisting of K component states. In general, the state of the system
at time t E R can be represented by a vector ztt) = {+)(t) where 2x(t) denotes the state quantity of the k-th component. Assuming that the observations y() = {yx (1)] are available for N time points t1....,ty tu (regularly spaced), the relationship between y(t) and z(t) can be described with the additive i.idd. Gaussian noise c(t) NO, o² I), as follows:
y(4) 2(4)) <(ti), 1=
(1)
Given z(t) and ODE parameters 0, the derivatives of the states are computed over known function f(-) as
0z(1)/0l = ((2)))
(2)
In the inverse problem for ODEs, we need to estimate 0 given observations y; however, this is challenging for general nonlinear function f(-) due to the lack of a closed form solution.
In this study, focus ODEs with non-negative states, z. Generally, in ODEs, it is not guaranteed that the state variables are non-negative, and whether they are non-negative depends the phenomenon being modeled and the problem setting. For example, in physics, the states of harmonic oscillators can have negative values. However, in mathematical models of population dynamics or infectious diseases, it is natural to define states as non-negative. In the following, we introduce examples of ODEs addressed this study.
SIR equations. The SIR (susceptible, infected, and recovered) equations in epidemiology describe the spread of infectious discases (Kermack & McKendrick 1927). This system handles three types of states 2s: 21, and zR, which denote the sizes of susceptible, infected, and recovered populations, respectively. Using {a, b) as the parameters of ODE, the equations are expressed as follows:
das -02521, dzj
azgzi bāg, dzn bzg.
In mathematical biology, predator-prey equations, also known the Lotka- Volterra predator-prey model, describe the population dynamics of two species, one as prey and the other The equations are given as follows:
Predator-prey equations.
as predator Lotka 1925
d 021 b2122, dt
dz2
-222 d 22,
dt
where and denote the population of prey and predator, and {a, b, c, d) are the parameters of ODE.
Competition equations. Competition equations, also known as the Lotka-Voltern competition model, describe the population dynamics among multiple species in competition (Smale, 1976). Under competition among K species, the time derivative of 24, the population of species i, is described as follows:
dz; = r; 1 " 03.525/Th))
K
dt
where r = (rx), == {n3, and a } for i= I,...,K, j 1,...,K K are the parameters of ODE. Ta denotes the internal growth rate of species i. n denotes the carrying capacity of species i. 0g denotes the degree to which the population of species j decreases the derivative of the population of species i. In particular, 0g is set to 1.
3

---

Published in Transactions on Machine Learning Research (09/2023)
2.2 Temporal Point Processes
Temporal point processes (TPPs) are stochastic processes of the event occurrence times. Using TPPs, we can analyze the probabilistic structure of a sot of randomly occurring time points, whose nature is determined by a non-negalive function, which expresses the rale of events to occur and is called the intensity function.
The Poisson process (Kingman 1992 is typical model in which it is assumed that the number of events in the interval A follows a Poisson distribution Poisson (f&s X(t)át). where X(t) is the intensity function. A Pois- son process with a time-varying intensity is referred to as an inhomogeneous Poisson process. Furthermore, when the time- varying intensity follows a stochastic process, it is particularly known as a Cox process (Cox, 1955). As a flexible model for fitting the Cox process, Moller et al. developed the log-Gaussian Cox process, in which the time-vare varying intensity function modeled as the exponent of a Gaussian process (GP) (Ras- mussen & Williams 2005), which is distribution over nonlinear functions. To represent the non-negative intensity function by temporal function x E Il generated from the GP, exponential transformation is used and the intensity at time t is described as follows:
X(1) Xn - exp (z(()) (>0), 0),
where A denotes the base intensity, meaning the expected number of events, and exp (z(t)) is the modulation.
An approach termed multivariate LGCP introduces correlations between GPs of multiple LGCPs to model interactions among different point processes (Waagepetersen al.) 2016 Hessellund et al. 2022a b). This is closely related to our study. However, our method differs in that it models the interactions between different point processes using ODFs.
Another major model in the literature on TPPs is a Hawkes process (Hawess 1971 self-excitatory process to model phenomena in which multiple events are concentrated over a relatively short time span. While the proposed model and Hawkes process are similar that they both have mechanisms that change intensity over time, Hawkes process limits its modeling to self-exciting events. However, the proposed model focuses on a diverse range of discrete events with properties other than self-excitation, through various ODE assumptions.
Recently, Neural TPPs, incorporating the expressive power of neural networks into TPPs, have attracted much attention (Shchur et al.) 2021). However, much of the litcrature on ODE parameter estimation ad- dressed in our study relies on MCMC inference (Gelman et. al.]
1995
MacKay 2003). Therefore, we focus
on extensions of the classical latent variable type of TPPs (Cox
1955
Moller et al. 1998) rather than the
neural approach, which often has difficulty in MCMC inference [Wenzel et al. 2020
Izmailov et al.) 2021
3 Proposed Model
3.1 An ODE-guided Poisson Process
We define an ODE-guided Poisson process as TPP with a timo-varying intensity, which follows the dynamics of an ODE. Suppose the ODE consists of K equations with parameters 0, and consider that the events of the k-th type are generated from a k-th inhomogeneous Poisson process with intensity function (t). Based
on Eq. 2 and Eq. (3), considering (t) exp(zx(t)) (x) as the k-th state quantity in the ODE, the derivatives of the modulation exp(x(()) = (exp(=x(1))) are written as
i @exp(x(t)) = f (exp (x(1). 0).
(4)
Øt
Here, the state of the ODE corresponds to the modulation of TPP. Because the intensity of the TPP is proportional to the state of the ODE, the higher the state value, the more frequent the event. This relationship between the state and event frequency must hold the ODE. In practico, ODEs, such as those in the population dynamics literature presented in Section 2.1 need to be chosen or developed based on domain knowledge and assumptions about the phenomena associated with discrete events.
4

---

Published in Transactions on Machine Learning Research (09/2023)
3.2 -CCP-based Gradient Matching
Since the inference of ODE-guided Poisson processes requires joint inference of the latent intensities and the parameters of ODE, we propose new approximate inference framework that combines ODE parameter ostimation with LGCP (Moller ct al., 1998). The basic approach to ODE parameter estimation, f.c., solving the inverse problem, is numerical integration (Tarantola, 2005). However, this approach is practically inap- plicable to complex problems due to computational complexity. Since the scalability of the ODE parameter estimation method with respect to problem complexity essential for handling a variety of ODEs, we incor- porale the gradient matching method, a Bayesian approach originally devised for efficient ODE parameter estimation (Varah 1982
In gradient matching, regression function approximates the trajectory of an ODE by aligning the gradients of the function with the derivatives of the ODE at series of discrete time points. Inspired by the success of studies on gradient matching using Gaussian processes Rasmusssn & Williams 2005 as regressors (Calder- head et al. 2008 Dondelinger et al. 2013 Gorbach et al. 2017] Wenk et al. 2019) we modify such GP-based gradient matching (GPGM) share GPs with LGCPs and develop gradient matching method applicable to discrete events. The resulting inference framework
extension of fast Gaussian procoss-based gradicnt matching (FGPGM) (Wenk ct al. 2019), which is statc-of-the-art gradient match- ing method. We call this framework log- Gaussian Cox process-basel gradient matching (LGCP GM). In the following, we begin with a description of FGPGM and then derive LGCP-CL CCMGMMMM by making three modifications to FGPGM.
be regarded
3.2.1 FGPGM (Wenk et al., 2019)
In the context of modeling with ODEs, Eq. [ which describes the observational model, can be described in probabilistic form, as follows:
plylz,o) NIyolaz,,
(5)
Then, we build a standard zero-mean GP prior on each latent state zs, as follows:
p(z II N[ Cou):
(6)
where o- { and Co. denotes the kernel matrix of GP with parameters di. Note that the latent variables z representing the states of the ODE are treated in GPGM as real variables, sctting in which z can take negative values even in ODEs defined with non-negative states.
Suppose the kernel function is differentiable, then the distribution over the derivatives conditioned on the state is given as follows:
p()|,, II p(zz Zkk Φbk)
II N22 (DZz, Ax). k
(7)
Dj TrT C-1
Φc &b.
A C" -CTT C- C"
tb
where C% is the autocovariance of each state derivative and C¹.. is the cross-covariance between the state and its derivative.
FGPGM introduces additional latent variables, fade and fyp: which have different origins, derivatives of ODE and gradients of GPs, respectively, and considers the equivalence between them First, given the state and the ODE parameters 0, fode can be directly calculated by ODE, i.e., fode f(z, 0) and particularly edd = (z, 0) for the k-th component. This deterministic relationship between fode and f(z, 0) is represented
In Ween et 20191, fode Eop are notated as F1 and F2, respectively.

---

Published in Transactions on Machine Learning Research (09/2023)
using Dirac's delta functions, follows:
(fode z, 0) = =I k
unde 2, 8)
((()
(8)
II
ode 5e(2..0)
p((k)
Second, given the state 2, the distribution over fgp = { as follows:
p((k)
K
is given with additional noise scaled by
gp
k=1
p(ceple, 22 7)
II W
p(k) EP
(9)
Finally, considering fade and Ep match, this equivalence is represented by Dirac's delta function 8 (fnde fep) Thus, joint density of latent variables and observation y can be written as
(z, 2, fode, fgp. 010,0,11 = p(y/z, ptllo//ppp//pppppzzl/// ()(())
(10)
x (fode z, p) (2, (fode - -pp
where p(0) denotes some prior on the parameters of ODE.
Following Theorem of Wenk and Eqs.
10
the joint posterior density of z and 0 in GP-based
gradient matching is written
0i,,0ly, x
N (yulze, T
II N (zx10, Cau)
x p()) II O( As 721)..
(11)
Note that é, fnde; and Eg are now marginalized out and do not appear in Eq. 11 In the last Gaussian factors on the right-hand side of Eq. (IIp, the variances small, thus inducing the gradients of GP, 2µ, to match the derivatives of ODE, (z, 0). From a different perspective, for each component k, the covariance of the gradient of GP, AE, along with additional noise scaled by 7, probabilistically models the expected value of the mismatch between the gradients of both GP and ODE, fg (z, 0) DD as a multivariate normal distribution. In other words, these variances control the degree to which violations of the deterministic (i.e., no misalignment allowed) constraints of the ODE are permitted.
3.2.2 Derivation of LGCP-GM
In this section, we derive LGCP-GM by making three modifications to FGPGM.
Modification Switching the likelihood. continuous Gaussian likelihood is typically used to model the observation noise in ODE modeling methods, including gradient matching. However, discrete likelihood, such as Poisson distribution, is necessary for TPPs. To bridge this gap, use the LGCP's hierarchical structure (Mller et al. 1998), which combines Poisson likelihoods and Gaussian prior. For approximate inference on LGCP, time is finely discretized and event data are converted to count data by counting the event occurrences over each interval. Let m&,e and Ik, be the event count and the state of the k-th type at the t-th discretized time with interval A, respectively. The likelihood of the state x given event counts m is evaluated using Poisson distributions as follows:
p(m|x)
p( (m..)x.)))
(12)
Poisson (mz A Ao exp(*x,)))

---

Published in Transactions on Machine Learning Research (09/2023)
Modification 2: Function completion with sparse GPs. This LGCP requires fine discretization of time, whereas gradient matching requires that the number of discrete points at which gradient evaluation is performed be kept smaller to reduce the computational complexity. To resolve this discrepancy, we use the inducing variable method for GPs (Quinonero-Canun
& Rasmussen, 2005), which also known as sparse GPs (Snelson Ghahramani 2005 While the inducing variable method was originally devised to reduce the computational complexity of the inverse covariance matrix, we use this method to handle many discrctized time points. While we set equally spaced and finely discretized T time points
{ referred to observation points in this study, for evaluating Poisson likelihoods of LGCP, also set coarsely discretized U T) time points v = {Nu)==1
inducing points Let. be the state at the inducing points and & be the approximated value of the state at the observation points. The distribution of & is given as the posterior predictive distribution of the Gaussian process conditioned on with the covariance approximated diagonally, as follows:
p Φ) = III (St& Mse Ezc),
k
(13)
Us.. C.oCCrt
Σx, diag(C Cro
where C.. denotes the kernel matrix with parameters Φv and diag(v) denotes a diagonal matrix with vector u as its diagonal elements. Subsequently, by restricting the gradient evaluation to the inducing points, the fine discretization for LGCP and the coarse discretization for gradient matching become compatible. In this modification, interpolating the states at the observation points is required for the above compatibility, whereas sparse approximation, i.e., the diagonal approximation of the covariance of the conditional GP, is not essential. However, sparse approxination reduces the computational complexity of evaluating large multivariate normal distributions for the states of the observation points, which is suitable for the requirement of a larger number of observation points by &, as follows:
3 Note that as a result of this modification, x in Eq. (12) is replaced
p(m[ pl p( (
(14)
Poisson (ma,l A An exp($xx))
Modification 3: Transformation of the ODE. According to Eq. ( the derivatives of the ODE at the inducing points are evaluated in a non-negative real space for exp(x), whereas the gradient of the GP is evaluated over the entire real space rather than the nonnegative real space. To match both gradients in the same space in the gradient matching scheme, we redefine the derivatives of the ODE function of x instead of exp(x), as in Eq. as follows:
ox
σi g(x,0) f (cxxp(x) 0)
(15)
exp(x)
3.2.3 Joint Posterior of LGCP-GM
Considering the above, we modify Eq. that of LGCP-GM as follows:
11 the unnormalized joint posterior density of GPGM, to obtain
p(x, x 0 () 7) &
p Am (x..e) x II N (Stp Mxxe Ex. ) x II N (x+l0,
k
(16)
x ( x II (gx O) A 2
If the two probability factors in the second line on the right-hand side of Eq. 16 are ignored (i.e. removal of constraints regarding ODEs), LGCP GM is reduced to LGCP with sparse GPs, as follows:
"Note that the number of induction points U must not be so small that the GP conditioned the induction points cannot capture the global pattern of the function.
3 AAthhough the computational complexity T-dimensional Gaussian distribution is O(T3), when the covariance matrix is diagonal matrix, the complexity is reduced O(T).
7

---

Published in Transactions on Machine Learning Research (09/202))
f SP
7
Egp
X
0
fode
Φ
0
Sode
X
Φ
yk.t
σk
m
T
T
K
K
(a) FGPGM (Wenk et al.
2019
(b) LGCP GMM (ours)
Figure 2: Graphical representations of gradient matching schemes for ODE parameter estimation.
x Φ) x IIII p (mmaalaaa) x II M2L Exn. x Nxx Cor).
(17)
k
For prior distributions of 0, logit normal distribution is set the prior for the scaled variable Logit of), where (0- mm - mmii) to set the range (0min maxx) on 0 to reflect prior domain knowledge. We treat Φ and 7 fixed parameters, in FGPGM (Wenk et al. 2019). In FGPGM, both Φ and are set by optimizing GPs on the observed data whereas in LGCP- GM, optimizing kernel parameters & by applying LGCP to the observed variable m can easily fall into local optima due to the complexity caused by the non-Gaussian Poisson distribution. Thus, in LGCP GM, we empirically set the kernel parameters & by grid search using the log posterior density evaluated in prior experiments as metric.
3.2.4 Comparison of Graphical Representations
For comparison, Fig. 2 shows graphical models of both the GP-based and LGCP-based gradient matching schemes. There are differences in these graphical models due to the three modifications described above: (1) the shaded node of ma.t and the edge from &R,t to mK, in Fig. 2(b) reflect the first. modification regarding Poisson likelihood, (2) the edge from Φ to tt Fig. 2(b) reflects the second modification regarding function completion with sparse GPs, and (3) the changes in variables from fode to gode and from fgp to ggp reflect the third modification regarding the Lransformation of ODE. Here, fode and gode denote the derivatives of the ODE, in Eq. (2) and Eq. [15 respectively. In addition, fgp and ggp denote the gradients of GPs with the addition of small noise & N (0, 21). The shaded edges denote deterministie relationships between variables.
3.3 Extrapolation using LGCP-GM
The proposed model allows the estimation of the parameters of ODE within normalized time 0 S 1, where data are available, referred to the interpolation time. A simple extension makes it possible to simultancously predict the dynamics at the time when no data exist, in the extrapolation time 1 < t. Let us introduce inducing points v' and observation points 7' in extrapolation time and let x' and J' be the states at v' and T' respectively. z' denotes the latent dynamics in the extrapolation time range. Additionally, let m' be the null count of event occurrence at T'. By setting p(m'[x') = 1, we can skip the likelihood evaluation for the null counts to avoid erroncous evaluation, time. As result, we estimate x' and & via the mechanistic knowledge of ODEs learned from discrete event observations at the interpolation time.
if the event count is not null but zero in the extrapolation
4 Inference
Because of the analyiilllaal form of Eq. [6]. we require an approximate inference using the Markov chain Monte Carlo method. While previous studies used the random-walk Metropolis-Hastiine algorithm
8

---

Published in Transactions on Machine Learning Research (09/2023)
used Hamiltonian Monte Carlo (HMC) method, which is suitable for high-dimensional problems (Duane et al. 1987 Neal 2011). To use HMC, we need to derive the first derivatives of the log posterior density with respect to cach latent variable to be sampled. These can be derived analytically or computed using a library of automatic differentiation, but are not detailed here.
to extract samples from the posterior distribution (Dondelinger et al. 2013 Wenk al. 2019
When HMC is applied hierarchical Bayesian model, the combination of HMC and blocked Gibbs sampling can be effective [Neal] 2011). In the block sampling strategy (Roberts & Sahu, 1997), the latent variables are partitioned into groups, and subset of the latent variables in the groups under update is sampled alternately. Therefore, we take the latent variables &, x, and 0 as different groups, probabilistically select
one or more groups for each iteration of the inference, and repeatedly update the variables in the selected groups.
The high corrclation of latent functions, x, that follow GPs makes MCMC inference difficult (Titsias et al. 2011). We use variable transformation to reduce correlations (Kuss al. 2005). Consider linear transfor- mation Lxxx, where C. LL is the Cholesky decomposition. Because is white with respect to C.. can avoid the high-correlation issue by sampling rather than Moreover, because the covariance of the GP gradient A also results in a high correlation of the posterior density, we mitigate the high correlation by approximating A by its diagonal matrix, A = diag(A) A).
Recall that the last Gaussian factors on the right-hand side of Eq. 16 are multivariate Gaussian distributions with small variances that evaluate the deviation between the gradients of GPs and ODEs. Since the latent variables randomly initialized at the beginning of the inference, the deviation between the gradients will not be small, and the likelihood of these Gaussian factors will be low. Thercfore, in the carly iterations of inference, exploration of the state space may occur, such that the likelihood of these Gaussian factors is improved, even as the Poisson likelihoods involving directly observed data are reduced. This exploration may cause divergence of the state from the region of high probability density in the target distribution, which we originally intend to explore. To address this issue, we use the idea of annealing (Kirkpatrick et al. 1983) to apply the inverse temperature parameter & these Gaussian factors during the burn-in phase of the inference, as follows:
p x, 0l,, 2.
IIII
p p
II (8& (&x M8se Exe ) II
[ (xx10, Coc)
(18)
((0) (IE 0) 2
Initially, 8 0 is used for sampling, and the latent variables repeatedly updated while gradually ap- proaching 3 = in the burn-in iterations. Here, when B = 0, the Gaussian factors for evaluating the gradient error are completely ignored; when 3 = 1, Eq. (8) is consistent with Eq. (16). This allows the inference to begin without diverging the state from the high probability density region of the target distribution.
5 Experiments
5.1
Settings
This section describes the detailed setup of both the model and the MCMC algorithm. Times were normalized to values between O and 1. For observation points, T 100 and r = |0.005, 0.015, 0.995). For inducing points, U = 21 and v = [0.00, o.o.....,....0.00 1.00). additive noise Lerm,
We used the squared exponential kernel function with an
k(t, = of exp (- [ + 030(t= -
203
where σ1, d2, and d3 denote the amplitude, length, and white noise scales, respectively. As mentioned in Scction 3.2.3 through preliminary experiments using grid search, we empirically set Φ1 5.0, d2 0.15 for SIR and competition equations and d2 0.1 for prodator-prey equations, d3 = 0.1, and γ = 0.1, hyperparameters.
9

---

Published in Transactions on Machine Learning Research (09/2023)
Table 1: Root mean squared deviation between the posterior mean of 0 and ground truth
Type of ODEs
SIR
Predator- Prey
Competition
base intensity X0
50
100
1000
50
1000
50
100
1000
GPGM (T=20)
0.657
0.735
0.537
7.685
6.091
13.811
3.183
3.350
2.646
GPGM (T=100)
0.643
0.588
0.442
4.884
1.863
2.872
3.345
3.049
2.636
LGCP GM (ours)
0.379
0.070
0.119
4.643
2.235
1.786
2.305
2.208
2.047
Table 2: Mean of the negative log-likelihood of posterior samples for extrapolation time (1 <t<<>> 1.5)
Type of ODEs
SIR
Predator Prey
Competition
base intensity A
50
100
1000
50
100
1000
50
100
1000
LGCP LGCP-GM (ours)
619.4
464.8
5356.7
118.1
174.7
818.5
242.6
338.2
1741.8
166.5
213.2
433.2
112.8
144.1
361.7
195.4
268.6
602.3
For the HMC setup, the number of leapfrog steps was fixed at 10, and the step size was adjusted within the burn-in MCMC simulations. After 10, 000 burn-in iterations, 1, 1,000 posterior samples were collected from the subsequent 20, 000 simulations thinned to every 20 samples. Unless an exception mentioned, these settings were used for the proposed model and baselines in all experiments.
5.2
Simulated Study
5.2.1 Procedure to Synthesize Experimental Event Data
First, we fixed the initial states and ground-cruth parameters of the ODE and simulated the ground-truth dynamics by numerically solving the ODE. Next, we multiplied the simulated states by the same Å for all k = ,,..., K to obtain the intensity function (t) Åo exp( (t)) for the k-th component of the ODEs. Then, we discretized the time for t in [0, 1] with interval of 0.001 and synthesized the data by drawing event counts from Poisson distributions with the parameter 0.001 x (t) at each time point t.
5.2.2 ODE Parameter Estimation
We compared LGCP GM with GPGM with respect to the accuracy of ODE parameter estimation. As GPGM, we used slightly modified versions of FGPGM (Wenk et al. 2019) incorporating the inducing variable method and HMC sampling to improve scalability and conduct a fair comparison against LGCP GM. Since GPGM aggregates event data into bins as count data and models them in discrete time, it necessary specify the number of bins for discretization. We used two versions of GPGM as basclines: observation points for the coarse discretization case, referred to as GPGM (T = 20), and the other with 100 observation points, same as in LGCP-GM, for the fine discretization case, referred to GPGM (T 100). GPGM (T' = 20) can be viewed as an alternative baseline for comparison to the original FGPGM, since it has almost the same y and z lengths, and its setup similar to that of the FGPGM with which they coincide. particular, since GPGM (T 100) has the same number of observation points as T.GCP-GM. GM, the choice of likelihood, i.c., Gaussian or Poisson likelihood, the only difference between it and LGCP-GM GM. Therefore, GPGM (T = 100) considered purely comparative benchmark for differences in likelihood. As shown in Fig. 2 GPGM is applied to the observed value of y, which is computed event counts divided by the average number of events per observation point T.
with 20
Table D shows the rool mean squared deviation between the posterior mean of 0 and the ground-truth parameters of ODE as the accuracy metrics ODE parameter recovery for each method, type of ODEs, and bascline intensity λo which indicates the rate of data size. The bold font denotes the best result for each case. In most cases, the proposed method outperformed the baseline methods for both the trajectories and parameters of the ODE. In Fig. draw the estimated ODE trajectories, recovered by numerical
10

---

Published in Transactions on Machine Learning Research (09/2023)
syothesized events
2
A -
Ap
Ap 1000
50-
GPGM (T=20)
monusatio
modutation
00
0.0
GPGM (T=100)
50
LGCP- GM (ours)
on
1.
0.0
0.5
1.0
80
0.5
00
05
10
0.0
0L5
000
0L5
time
time
(a) SIR
(b) Predator-Prey
Ag - 50
Ag 100
Ag 1001
5.0
modulatinn
00
50
50
30
50
5.0
13
2.5
2.5
00
(10)
(10
time
(c) Competition
Figure 3: Numerically simulated modulations for different ODEs and base intensities, using the ground- truth imitial state and posterior samples of parameters of ODE. Ground-truth (dashed lines), median (solid lines), and 75% quantiles (shaded arcas) of the simulated modulations indicate that LGCP-GM (bottom row) outperforms GPGM baselines with resepect to modulation recovery.
synthesized
events
(xxxxs
LGCP
(xXixx
LGCP- GM
(ours)
0.0
0.5
100
1.5
time
Figure 4: Extrapolation of modulation for predator- predator-prey prey equations and Ào 1000. The estimated modula- tions indicate that LGCP-GM can recover the true modulation (red dashed line) in the extrapolation time (the right side of the dashec-dotted line), but LGCP cannot.
integration using the postcrior sample of ODE parameters and the ground-truth initial state, with the synthesized discrete events, and compare them with the ground truth trajectory. For each type of ODE, comparing the results from GPGM with those from LGCP-GM, the recovered ODE trajectories also showed the advantage of the proposed method, especially in settings with little data.
5.2.3 Extrapolation
We evaluated extrapolation accuracy for the time < < 1.5. As a baseline, we used LGCP (Moller el al. 1998) with sparse GPs using Eq. (17) as the posterior density. The kernel parameters & for LGCP were determined by a grid search using the log postcrior density evaluated in the preliminary MCMC inference as a metric. We used the negative log- ligeliid log pmm'|>'), as a metric for extrapolation accuracy. To evaluate the negative liellhod, additionally simulated discrete events for the extrapolation time from
11

---

Published in Transactions on Machine Learning Research (09/2023)
3.5
3.0
Observations: S, I, R Only I observed Ground- truth
2.5
2.0
1.5
1.0
a
b
Figure 5:
Estimated parameters of SIR model for two cases: (1) all types of events are observed, and (2)
only type I events are available.
the ground truth intensities at several different Ao settings and counted the events to calculate m'. This was repeated 100 times to simulate different samples of m'. We evaluated the negative log-likelihood for each combination of simulated m' and the posterior sample of x from LGCP and LGCP-GM, respectively.
Table 2 shows the mean of the evaluated negative log-likelihoods for each method, the base intensity, and the type of ODE. The bold font denotes the best result for each case. In all cases, LGCP GM consistently outperformed LGCP with respect to extrapolation. The performance especially differs in the case with a large base intensity, where the estimated parameters of ODE tend to be more accurate. This suggests that the accurately estimated ODE parameters are useful for extrapolation. Fig. | illlustrates the estimated modulation for both LGCP and LGCP GM in the case of predator-prey equations and Åo = 1000. Looking at the ground-truth modulation (red dashed line), median (solid line), and 75% quantile (shaded area) of the estimated modulations, the proposed model succeeded in recovering the dynamics in the extrapolation time, on the right side of the dotted lines, whereas LGCP failed.
5.2.4 Modeling on Non-available Events
Consider a scenario where we cannot obtain discrete event observations of certain components in an ODE system consisting of K states. In such setting, several studies have developed ODE modeling methods for numerical data (Huang al. 2020 Linial et al., 2021 Yang et al. 2021 and we performed ODE modeling for event data using the proposed method. Let U be a set of unavailable component indices in the ODE. Based on the same idea of skipping the likelihood evaluation in extrapolation, we can naturally incorporate the unseen k-th components of the ODEs by sctting the likelihood at observation point as p(m,,,l = and p( ok. Φb)) = as long as & e U. We estimated the parameters of SIR model using only events for I, assuming that events for S and are unavailable. In addition, we assumed the event count for R at the first observation point to be zero, which softly bound the initial value of the modulation of R to be near zero. This is ad hoc but is natural assumption in SIR model.
In Fig. 1 we compare two cases: one when events, synthesized with Ao = 1000 for S, I, and R are available, and one when only those for I are available. Fig. 5 shows the estimated parameters for the two cases. Although the posterior variances of both the modulations and parameters of ODE were higher in the second case, the estimation was as successful as in the first case. These results suggest that there are scenarios where the proposed model can be applied, even when some types of events are unavailable.
5.3 Applications
Competition in Technological Industry. To show an application example of the proposed method to real data, we estimate the competitive relationships among companies related a specific technological theme using actual patent application data published by the United States Patent and Trademark Office T We extracted patents filed between 2012 and 2020 that had an organization name in the first applicant, had "Handling natural language daia" as the International Patent Classificatiion (IPC) label, and included the word "speech" as a keyword in their abstracts. The first applicant (organization) of patents was counted, and the patents of the top five companies in terms of number of patents were extracted. Using the proposed method and competition equations, we analyzed the patent application event series for the 210 patents extracted, as described above. Fig. shows the patent application event data and an adjacency matrix consisting of the estimated competition coefficients. Coefficient di.j denotes the competitive influence from
4Data source: USPTO patent application data
12

---

Published in Transactions on Machine Learning Research (09/202))
Google
1:114
Google
1BM
Microsoft
Microsoit
SAMSUNG
SAMSUNG
TOSHIBA
TOSHIBA
Google
IDM
Minosoft
SAMBUNG
TOSHIBA
2014
2015
2016
2017
2018
2019 2020
Patent application date
(a) "speech" related patent application events
(b) competition coefficienss
Figure 6: (a) Observed events and (b) adjacency matrix consisting of the posterior means of competition coefficients a.
Darker colors denote stronger impact of the column component on the row component.
100
Querry
jQuery
80
Interest
60
Angular
React
Angular
40
Vue.js
React
Vue js
2010
2015
2020
jQuery
Angular
React
Vue.js
Year
(a) Search trends of web frameworks
(b) competition coefficients
Figure (a) Observed event counts and (b) adjacency matrix the same style as Fig. 6(b)
the j-th component to the i-th component. The row for SAMSUNG has small competition coefficients (in- dicated by light colors). This corresponds to steady increase in the number of SAMSUNG's applications, independent of the application trends of other companies. In contrast, the for TOSHIBA has large com- petition coefficients (indicated by darker colors). This corresponds to the fact that TOSHIBA's 's applications decreased during the period when other companies' applications increased; and that it has not filed any applications since 2018. The column for SAMSUNG also has large compelition coefficients, corresponding to the fact that the trend of others' applications declined during the period when SAMSUNG's applications increased. Although these results do not directly represent the actual rivalry among these companies, they allow estimate the competition coefficients corresponding to changes patent application trends.
Competition of Web frameworks. In the industry of information technology, web application frame- works have several major tools and their popularity has changed over time. Using the proposed model and competition equations, we estimated the competition dynamics between the following representative web frameworks: JQuery, Angular, React, and Vue.js. We retrieved data from Google Trends, tool that evalu- ates search volume over time, for worldwide web searches the Internet and telecommunications category from January 1, 2005, to December 31, 2021. Note that the data obtained with the keyword " "Angular" be a mixture of results for AngularJS and Angular, which are strictly different software. Data consist. of positive integers (up to 100) proportional to the search volume, or the string, "less than 1", for each search target and month. We then converted all string values to zeros in order to make all values numeric. We considered these integers as count data aggregated for each month and applied the proposed model. For each component k, (k) was set to the total number of counts for events of the k-th typc. Fig. 7(a) shows the line plots of the data, where the legends are arranged in order of their appearance. Fig. 7(b) shows an adjacency matrix consisting of the estimated competition coefficients, in the same manner as in Fig, 6(b) Looking the larger competition coefficients represented by darker colors, we can see the impact of Angular on JQuery, React on Angular, and Vue on React. This result illustrates the impact of newer tools on the older ones.
13

---

Published in Transactions on Machine Learning Research (09/2023)
#of twaets/day
25
Minions (2015)
175
A Aengerss: End Game (2019)
20
15
100
10
A. Mad
0
30
60
90
120
30
90
120
day
day
Figure 8: Number of tweets about the two movies over time.
Table 3:
Estimated parameters (mean 22 standard deviation) of SIR model for two movies
Movie Production
b
Minions (2015)
3.03 0.49
7.13 = 1.10
Avengers: End Game (2019)
4.46 0.69
12.98
Trend Spreading on Twitter. Recently, information sharing among individuals on Twitter has become regular. To understand the trend spreading on Twitter, many attempts have been made to analyze tweet data using epidemiological mathematical models (Abdullah & Wu, 2011 Jin et al. 2013). In these studies, a user who tweeted about topic was considered as an infectious agent the mathematical model of infectious diseases with respect to the given topic. Then, using SIR model with S and R as components that cannot be observed, we attempled lo analyze movie trends from MovieTweetinng
(Dooms et al. 2013), dataset of movie review tweets [ Fig. S shows the number of tweets per day for 120 days, starting 30 days before the time of the first review tweet for the two films. For all component A including S and R, x was set the total number of tweets. The estimated parameters of SIR model for each movie presented in Table
The magnitudes of and b can be interpreted as the rates at which the trend rises and then falls. The estimated results show that the first movic generally has slower rising trend and more sustained topicality than the second movie. Although this can be observed in Fig. g because these movies have significantly different trends, the estimated parameters of ODE allow for numerical comparisons between movies with similar trends and between multiple movies, which is visually difficult.
6 Conclusion
Combining temporal point processes with gradient matching methods could increase the accuracy of ODE parameter estimation when the data are irregularly observed discrete events compared with existing gradient matching methods for numerical data, and allow for a suggestive interprctation of the event data via param- eters of ODE. Furthermore, exploiting the mechanistic knowledge from the ODEs learned by the proposed method could compensate for the absence of data, making better extrapolation modeling with unavailable components in the ODE system.
References
Saeed Abdullah and Xindong Wu. An epidemic model for news spreading on twitter. In 2011 IEEE 23rd international conference tools with artificial intelligence, pp. 163-169. IEEE, 2011.
Fred Brauer. Mathematical epidemiology: Past, present, and future. Infectious Disease Modelling, 2(2): 113-127, 2017.
Ben Calderhead, Mark Girolami, and Neil Lawrence. Accelerating bayesian inference nonlinear differ- ential equations with gaussian processes. Advances in neural information processing systems, 21, 2008.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Daaa source: MovieTwwetting. (fng:://gghhaab..mm///moo/////h
14

---

Published in Transactions on Machine Learning Research (09/2023)
Ricky TQ Chen, Brandon Amos, and Maximilian Nickel. Learning neural event functions for ordinary differential equations. arXiv preprinl arXix:201100300 2020.
D. R. Cox. Some statistical methods connected with series events. Journal of the Royal Statistical Society. Series B (Methodological), 17(2):129-164, 1955. ISSN 00359246. URL http://www.jstor.com
Saale/2983955)
DJ Daley and David Vere-Jones. An Introduction to the Theory of Point Processes: Volume II: General Theory and Structure. Springer Science & Business Media, 2007.
Frank Dondelinger, Dirk Husmeier, Simon Rogers, and Maurizio Filippone. Ode parameter inference using adaptive gradient matching with gaussian processes. AISTATS, pp. 216-228. PMLR, 2013.
Simon Dooms, Toon De Pessemier, and Lue Martens. Movietweetings: a movie rating dataset collected from twitter. In Workshop on Crowdsourcing and human computation for recommender systems, CrowdRec at RecSys, volume 2013, pp. 43, 2013.
Simon Duane, Anthony D Kennedy, Brian Pendleton, and Duncan Roweth. Hybrid monte carlo. Physics letters B, (95(2): 222, 1987.
Andrew Gelman, John B Carlin, Hal S Stern, and Donald B Rubin. Bayesion data onalysis. Chapman and Hal//CRC, 1995.
Nico Gorbach, Stefan Bauer, and Joachim M Buhmann. Scalable variational inference for dynamical systems. Advances in neural information processing systems, 30, 2017.
Alan G Hawkes. Spectra of some self-exciting and mutually exciting point processes. Biometrika, 58(1): 83-90) 1971.
Kristian Bjorn Hessellund, Ganggang Xu, Yongtao Guan, and Rasmus Waagepetersen. Second-order semi- parametric inference for multivariate log gaussian cox processes. Journal of the Royal Statistical Society Series C: Applied Statistics, 7141):244-268 2022a.
Kristian Bjørn Hessellund, Ganggang Xu, Yongtao Guan, and Rasmuu Waagepetersen. Semiparametric multinomial logistic regression for multivariate point pattern data. Journal of the American Statistical Association, 117(639)::5000-6111 2022b.
Hanwen Huang, Andreas Handel, and Xiao Song. bayesian approach to estimate parameters of ordinary differential equation. Computational statistics, 35(3):1481 1499, 2020.
Pavel Izmailov, Sharad Vikram, Matthew D Hoffman, and Andrew Gordon Gordon Wilson. What are bayesian neural network posteriors really like? In International conference on machine learning, pp. 4629-4640. PMLR, 2021.
Junteng Jia and Austin R Benson. Neural jump stochastic differential equations. Advances in Neural Information Processing Systems, 32, 2019.
Fang Jin, Edward Dougherty, Parang Saraf, Yang Cao, and Naren Ramakrishnan. Epidemiological modeling of news and rumors on twitter. In Proceedings of the 7hh workshop on social network mining and onalysis, pp. 9, 2013.
Douglas Samuel Jones, Michael Plank, and Brian Sleeman. Differential equations and mathematical biology. Chapman and Hal//CRC, 2009.
William O Kermack and Anderson G McKendrick. A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London A: mathematical, physical and engineering sciences, 115(772): 700-721, 1927.
John Frank Charles Kingman. Poisson processes, volume 3. Clarendon Press, 1992.
15

---

Published in Transactions on Machine Learning Research (09/2023)
Scott Kirkpatrick, C Daniel Gelatt Jr, and Mario P Vecchi. Optimization by simulated annealing. science, 220(4598):671-66)) 1983.
AV Kryazhimskiy and Yu Osipov. Inverse problems for ordinary differential equalions: dynammicul solutions. (Gordon and Breach) Taylor and Francis, 1995.
Malte Kuss, Carl Edward Rasmussen, and Ralf Herbrich. Assessing approximate inference for binary gaussian process classification. Journal of machine learning research, 6(10), 2005.
Ori Linial, Neta Ravid, Danny Eytan, and Uri Shalit. Generative ode modeling with known unknowns. In Proceedings of the Conference Health, Inference, and Learning, CHIL '21, pp. 79 94, New York, NY, USA, 2021. Association for Computing Machinery. ISBN 9781450383592. doi: 10.1145/3450439.0 1866.
URL
Alfred James Lotka. Elements of physical biology. Williams & Wilkins, 1925.
David J. C. MacKay, Information theory, inference, and learning algorithms. 2003. URL http: www.
inference. www.:phhyyyyccc.mmmmm www.c/aa.aa.....aa... mackay/iiiial..
Jesper Møller, Anne Randi Syversveen, and Rasmus Plenge Waagepetersen. Log gaussian processes.
Scandinavian Journal of Statistics, 25(3):451 482, 1998. doi: 1tt.//.or///..111////4/444.
URL httppe//mmmimmiiirrrryeelle..o/.
James Dickson Murray and James Dickson Murray. Matherattiaal Biology: II: Spatial Models and Biomedical Applications, volume Springer, 2003.
Radford M Neal. Meme using hamiltonian dynamics. Handbook of Markov Chain Monte Corlo, 2(11):2, 2011.
Joaquin Quinonero-Candela and Carl Edward Rasmussen. A unifying view of sparse approximate gaussian process regression. The Journal of Machine Learning Research, 6:1939 1959, 2005.
Carl Edward Rasmussen and Christopher I. Williams. Gaussian Processes for Machine Learning (Adaptive Compulation and Machine Learning). The MIT Press, 2005. ISBN 026218253X.
Gareth O Roberts and Sujit K Sahu. Updating schemes, correlation structure, blocking and parameterization for the gibbs sampler. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 59(2): 291-317, 1997.
Yulia. Rubanova, Ricky TQ Chen, and David K Duvenaud. Latent ordinary differential equations for irregularly-sammlled time series. Advances in newral information processing systems, 32, 2019.
Oleksandr Shchur, Ali Caner Türkmen, Tim Jamuschowski, and Stephan Günnemann. Neural temporal point processes: A review. In Zhi-Hua Zhou (ed.). Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21, pp. 4585
4585-4593. 4593. International Joint Conferences on Artificial Intelligence
Organization, 8 2021. doi: 10.24963/ijcaa.222 2021 2021/633. /623. Uttphtts: /ww.:////ooo////// httpd//.org....9933/////// 2ttpg//000229433jjaa2222202 /623
Survey
Track.
George F Simmons. Differential equations with applications and historical notes. CRC Press, 2016.
Stephen Smale. On the differential equations of species in competition. Journal of Mathematical Biology, 3: 5-7, 1976.
Edward Snclson and Zoubin Ghahramani. Sparse gaussian processes using pscudo-inputs. Advances in neural information processing systems, 18, 2005.
Albert Tarantola. Inverse problem theory and methods for model parameter estimation. SIAM, 2005.
Michalis K. Titsias, Magnus Rattray, and Neil D. Lawrence. Markov chain Monte Carlo algorithms for Gaussian processes, 2011. URL http://wwwsssallllloooooolllllls
16

---

Published in Transactions on Machine Learning Research (09/2023)
James M Varah. A spline least squares method for mumerical parameter estimation in differential equations. 46, 1982.
SIAM Journal on Scientific and Stalistical Compuling,
Rasmus Waagepetersen, Yongtao Guan, Abdollah Jalilian, and Jorge Mateu. Analysis of multispecies point patterns by using multivariate log-gaussian processes. Journal of the Royal Statistical Sociely Series C: Applied Statistics, 65(1):77-96, 2016.
Philippe Wenk, Alkis Gotovos, Stefan Bauer, Nico S Gorbach, Andreas Krause, and Joachim M Buhmann. Fast gaussian process based gradient matching for parameter identification in systems of nonlinear odes. In The 22nd International Conference Arlificial Intelligence and Slatistics, pp. 1351-1360 PMLR, 2019.
Florian Wenzel, Kevin Roth, Bastiaan S. Veeling, Jakub Swigtkowski, Linh Tran, Stephan Mandt, Jasper Snoek, Tim Salimans, Rodolphe Jenatton, and Sebastian Nowozin. How good the bayes posterior doep neural networks really? In Proceedings of the 37th International Conference on Machine Learning. ICML'20. JMLR.org, 2020.
Junchi Yan. Recent advance in temporal point process: from machine learning perspective. SJTU Technical Report, 2019.
Shihao Yang, Samuel WK Wong, and SC Kou. Inference of dynamic systems from noisy and sparse data via maanfold-constrainge gaussian processes. Proceedings of the National Academy of Sciences, 118(15): e2020397118, 2021.
17


Reviews:
[{"review_id": "H2KHO745E6", "reviewer": "Reviewer_db6R", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths.**\n\nThe main strengths of the work are perhaps at the beginning and at the end. I particularly liked the presentation/motivation and clarity around the importance of ODE modelling for discrete events or irregularly spaced data for dynamical systems. In that sense, I remark that the paper was clear and concise -- which makes everything more easy to be understood. On the other hand, the experimental results included are meaningful for me, and the authors did an effort to provide sufficient empirical values to prove the performance of their model. In this way, I think that the paper has a complete \"tail\" of experiments that matches what would be usually required for a submission of this sort to TMLR. However, this does not mean that the empirical results could not be improved, and I will add some comments later, particularly on the way that the method is compared to other baselines.\n\n**Weaknesses.**\n\nIn my opinion, the work might be understood as an extension of the ODE model proposed in [Wenk et al. 2019] for accepting discrete events. My main concerns and fears around the novelty/quality of the proposed method come mainly from the section 3. Using [Wenk et al. 2019] as a reference, one can easily see the similarities between both works, and how the current manuscript inherits both the notation and the structure from the other. However, in the current case, the reasons behind the use of a GP regressor and the elements playing some role in the joint distribution of Eq. 5 are not well explained. I mainly say this bc one can easily follow which distribution plays which role in [Wenk et al. 2019] (for instance, Eq. 3, Eq. 4, Eq. 5, Eq. 6, ...), but this is not happening in the current manuscript and I only understood how the model was working once I went to the mentioned reference. In that sense, I would say that the work should be self-contained, not requiring this sort of support from the main SOTA.\n\nAdditionally, what is described at the beginning (pp. 2) as a list of contributions, where the ODE-guided Cox process shines as the main one, later seems to be reduced to the three modifications included in section 3. Moreover, in this direction I am concerned about the modification #2, and how and why the inducing variables fit in the ODE method. This is not clear to me. Some questions around this point might be around how large is the computational cost, where does it come from, how the number of inducing points affect to the method, and if is it worth to make such approximation to the ODE model. I can see that [Wenk et al. 2019] chose different approximations and heuristics, that are not considered or mentioned here. I just want to remark that if the previous part was somehow incomplete in the number of details, the one on sparse GPs is even more limited and difficult to understand. Additionally, just for reference, in the GP literature, a change of likelihood plus a sparse approximation and a linking transformation is generally not considered as a novel advance. For this reason, I think is very important to understand why the current approach is novel and important in the ODE context with discrete events.\n\nAlso, I am fine with the experiments, and I see that the performance looks relatively good and improved wrt. to the chosen baseline [Wenk et al. 2019]. However, the lingering question in my mind is that if the LGCP model, (which as far as I understood considers a Gaussian likelihood) is trained on discrete events. In that direction, I would like to know just if it makes sense to run a ODE-GP model which has a continuous likelihood function on discrete data, and if this might cause the performance to be significantly worse than the proposed approach. \n\n**References.**\n[Wenk et al. 2019] -- Fast Gaussian process based gradient matching for parameter identification in systems of nonlinear ODEs. *Philippe Wenk, Alkis Gotovos, Stefan Bauer, Nico S. Gorbach, Andreas Krause, Joachim M. Buhmann Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, PMLR 89:1351-1360, 2019.*", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "ODE modelling and particularly the estimation of their parameters is usually considered from observed data that is regularly sampled and of continuous nature. Importantly, this work considers cases where the event data are discrete/irregularly spaced. This point leads to use temporal point processes (i.e. Cox processes) to model such discrete events. As a consequence, a new ODE-guided Cox process is presented, where the latent intensity function is assumed to follow the dynamics of an ODE. In practice, the work builds on the previous model of [Wenk et al. 2019], where a system of parametric ODEs is combined with GP regression in order to match the time points to the corresponding state values. This one is extended in three specific ways: 1) A switch of likelihood that now is Poisson instead of Gaussian, 2) the introduction of inducing points and 3) a transformation of the ODE to match to the real space (instead of the non-negative).", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": ">**C1.** Section 3 should be -- *at least* -- as clear as [Wenk et al. 2019] and self-contained. Including an Algorithm would be also helpful and the main differences should be remarked. In some way, the point around if the current work is an extension of Wenk et al. 2019 or not should be extremely clear for the reader since the beginning.\n\n>**C2.** Motivation, details and consequences of including sparse approximations should be discussed and clarified also in section 3 or 4.\n\n>**C3.** The comparison of the proposed model wrt baseline should be as equal as possible, in the sense of having a similar likelihood model or using discrete/continuous events for each one if needed. Just to perceive if the difference of performance is only produced by this difference in the likelihood function or not.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "aPkSyLhVw1", "reviewer": "Reviewer_A6CU", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The idea to use ODE to model interactions between different point processes is novel and highly practical.\nThe paper is well written and the presentation is clear.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces an innovative ODE-based method for concurrently estimating the intensities of multiple point processes. By employing ordinary differential equations controlled by a set of parameters, the model captures the interplay between various point processes, leading to easily interpretable estimates. To estimate intensities in a nonparametric fashion, the approach utilizes Gaussian process regression with MCMC sampling. This method proves beneficial for extrapolation and handling non-available events. The efficacy of the proposed approach is validated through simulation studies and real-world data examples.\n\n\n\n\n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. Some statements in the paper is inaccurate. For example, one page 1, it is stated that \"TPPs model the probability of such event occurrences using a latent intensity function; in particular, a model where the intensity varies over time is called a Cox process (Cox, 1955).\" The intensity of a Cox process is not only time-varying but also stochastic (as opposed to a deterministic function). A TPP with a time-varying but deterministic intensity is usually called inhomogeneous/nonhomogeneous Poisson process. Based on my understanding of the paper, the intensities of the point processes under consideration (i.e., exp(x_k(t))) are still deterministic. Therefore, the following term \"ODE-guided Cox process\" is a bit misleading. I would suggest changing this to \"ODE-guided Poisson Process\" to be more accurate.\n\n2. The statement in the abstract \"our method outperformed existing baseline methods assuming regularly observed continuous\ndata with respect to the accuracy of recovering the latent parameters of ODEs\", is also a bit mis-leading. Based on the likelihood given in (7), the proposed approach also discretize the point process and obtain counts data on a regular grid, which is the same as the baseline method. I think the difference is that the baseline approach assumes a Gaussian likelihood for the counts data and therefore is mis specified. Please clarify this issue.\n\n3. Using ODE to model interactions between different point processes is an interesting idea, but there are other existing approaches for the same purpose. For example, one can assume that each point process is a LGCP and introduce correlations between these latent Gaussian processes, which is commonly referred to as multi-variate LGCP. See, e.g., [1], [2], [3]. The paper seems to completely disregard this line of closely related research. Some comments on the differences between these two approaches are warranted here.\n\n4. For the Gaussian process regression, the choice of the parameters of the covariance kernel function of the latent Gaussian process is important. How do you choose these parameters in your simulation study and real data examples? How sensitive are the performances to the choice of these parameters?\n\n\nReferences:\n\n[1] Waagepetersen, R., Guan, Y., Jalilian, A., & Mateu, J. (2016). Analysis of multispecies point patterns by using multivariate log-Gaussian Cox processes. Journal of the Royal Statistical Society Series C: Applied Statistics, 65(1), 77-96.\n\n[2] Hessellund, K. B., Xu, G., Guan, Y., & Waagepetersen, R. (2022). Semiparametric multinomial logistic regression for multivariate point pattern data. Journal of the American Statistical Association, 117(539), 1500-1515.\n\n[3] Hessellund, K. B., Xu, G., Guan, Y., & Waagepetersen, R. (2022). Second-order semi-parametric inference for multivariate log Gaussian Cox processes. Journal of the Royal Statistical Society Series C: Applied Statistics, 71(1), 244-268.\n\n\n", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "W4ZClwC9io", "reviewer": "Reviewer_5J9q", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper is well-written and is easy to follow. The method is well derived. The results are sufficient and convincing, and show improvements, which are expected since the underlying mechanistic ODE matches the usecases. The paper has application-focused impact.\n\nThe method is incremental, and there is little machine learning novelty. The method derivation also has some issues. The model is not stated in full, but built gradually which makes the final model somewhat implicit. For instance, eq 4 doesn't even define a LGCP, but just a deterministic Cox process. The paper relies heavily on the FGPGM, which also decreases its novelty. \n\nThe bigger question I have towards this work is why do we need the gradient matching at all. The paper is mostly about the SIR model, which only has 2 parameters to learn. To fit the SIR to event data, one only then needs to fit 2 parameters and 3 initial states of z's. Surely this is easily doable with off-the-shelf optimisers even with large datasets or long trajectories. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents a log-Gaussian Cox process where the intensity is driven by a mechanistic ODE, whose parameters are estimated. The paper adapts a gradient matching based inference towards the model.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "No issues", "broader_impact_concerns": "No issues", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Reviewers commended the paper for its clarity as well as the interesting application studies. They commented the method could be useful for many different applications in event data analysis, and they found the results on non-conventional ML applications are worth publishing at TMLR.

The main concern is about novelty, where reviewers pointed out a missing related work (Wenk et al. 2019). In revision the authors included substantial text regarding discussions of Wenk et al. (2019) which has largely cleared reviewers' concern.

In the revision for camera ready, I would suggest the authors to also include some experiments on the fast GPGM method (Wenk et al. 2019) or at least clarify why the original GPGM method serves as a sufficiently good baseline.
Decision: Accept (poster)