A Probabilistic State Space Model for Joint Inference from Differential Equations and Data
Authors: Jonathan Schmidt
Nicholas Krämer
Philipp Hennig
Venue: NeurIPS 2021 Poster (2021)
Keywords: state space models
differential equations
gaussian processes

Abstract: Mechanistic models with differential equations are a key component of scientific applications of machine learning. Inference in such models is usually computationally demanding because it involves repeatedly solving the differential equation. The main problem here is that the numerical solver is hard to combine with standard inference techniques. Recent work in probabilistic numerics has developed a new class of solvers for ordinary differential equations (ODEs) that phrase the solution process directly in terms of Bayesian filtering. We here show that this allows such methods to be combined very directly, with conceptual and numerical ease, with latent force models in the ODE itself. It then becomes possible to perform approximate Bayesian inference on the latent force as well as the ODE solution in a single, linear complexity pass of an extended Kalman filter / smoother — that is, at the cost of computing a single ODE solution. We demonstrate the expressiveness and performance of the algorithm by training, among others, a non-parametric SIRD model on data from the COVID-19 outbreak.

A Probabilistic State Space Model for Joint Inference from Differential Equations and Data
Jonathan Schmidt
Nicholas Krämer
University of Tübingen
University of Tübingen
Tübingen, Germany
Tübingen, Germany
jonathan. schmidt@un huwettttttttttttteenneeee -tuebingen.dde
nicholas kww.kkraaaeear Camneerrunitttubbii -tt.bttgebnnnnnde de
Philipp Hennig
University of Tübingen
Max Planck Institute for Intelligent Systems
Tübingen, Germany
philipp
de
Abstract
Mechanistic models with differential equations are key component of scientific applications machine learning. Inference in such models usually computation- ally demanding, because involves repeatedly solving the differential equation. The main problem here is that the numerical solver is hard combine with standard inference techniques. Recent work in probabilistic numerics has developed new class solvers for ordinary differential equations (ODEs) that phrase the solution process directly in terms of Bayesian filtering. We here show that this allows such methods to combined very directly, with conceptual and numerical ease, with latent force models in the ODE itself. It then becomes possible to perform approximate Bayesian inference on the latent force as well as the ODE solution single, linear complexity pass of an extended Kalman filter / smoother that is, at the cost of computing a single ODE solution. We demonstrate the expressiveness and performance of the algorithm by training, among others, non-parametric SIRD model on data from the COVID-19 outbreak.
Introduction
Mcchanistic models based on ordinary differential equa- tions (ODEs) are popular across wide range of scientific disciplines. To increase the descriptive power such models, it is common consider parametrized versions of ODEs and find a set of parameters such that the dynamics reproduce empirical observations accurately possi- ble. Algorithms for this purpose typically involve repeated forward simulations in the context of, e.g., Markov-chain Monte Carlo or optimization. The need for iterated com- putation of ODE solutions may demand simplifications in the model to meet limits in the computational budget.
Only data Only ODE
Both (ours)
Data Truth
Figure 1: Inferring unknown func- tion with Gaussian Process and differ- ent sources of information.
This work describes algorithm that merges mechanistic knowledge in the form of an ODE with a non- parammertttic model over the parameters controlling the ODE latent force that represents quantities of interest. The algorithm then infers a trajectory that is informed
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

by the observations but also follows sensible dynamics, as defined by the ODE, in the absence of observations (Figure IT. The main insight enabling this approach is that probabilistic ODE solvers use the language of (extended) Kalman filters, conditioning on observations and solving the ODE itself is possible in one and the same process of Bayesian filtering and smoothing. Instead of iterated computation of ODE solutions, posterior distribution arises from single forward simulation, which has complexity equivalent to numerically computing an ODE solution, once, with a filtering-bssed probabilistic ODE solver [3)]. Intuitively, one can think of this opening up the black box ODE solver and acknowledging that each task - solving the ODE and discovering a latent force is probabilistic inference in state-space model.
The main contribution of this work is formalizing this intuition. Several experiments empirically prove the efficiency and the expressivity of the resulting algorithm. In particular, a practical model for the dynamics the COVID-19 pandemic is considered, in which non- non-parametric latent force captures the effect of policy measures that continuously change the contact rate among the population.
2 Problem setting
Let ar /max] Rdd be a process that is observed a discrete set of points T := (10 OBS
through a sequence of measurements Yo:N := (yo, .........
....... WN) E [R(++)xX Assume that these measure- ments are subject to additive i...d. Gaussian noise, according to the observation model
Un Hw((kn) + Ens en ~ NN,, R),
for 0, ..: N and matrices H E &Roxd and R Rexk Further suppose that sr() solves the ODE
a(t) = f(z());
and satisfies the initial condition a=(o) Ty E R.. The vector field f Rd x B'' + Rd is assumed to be autonomous, which is no loss of generality (c.g. (99) but simplifies the notation. The latent force [oo mmxx] + Re parametrizes f and shall be unknown.
SIR-type models (e.g. 7 are a common choice to describe the evo- lution of the COVID-19 pandemic. In SIR-type models, a population partitions into a discrete set of compartments. The differential equation then describes the transition of counts of individuals between these com- partments. For example, the SIRD model [T0] formulates the transitions between susceptible, infectious, recovered, and deceased people as
Susceptible
Infectious
(t) = (()) f(t) = (())((()))))))/))))))))))))))))) (-)(((()))/))())))) (T((), D(1) = TT((t),
R(()= = (t),
(3)
Recovered
Deceased
governed by contact rate B(t) : [tto mmaxl - (0. 1), recovery rate y E [0, 1]. and mortality rate 7) E [0, 1] (Figure 2}. S, I, Rt and D evolve over time, but the total population P (as the sum of the compartments) assumed remain constant. In this context, the contact rate ((1) is the latent force and varies over time (in the notation from Eq. (2). 3 is u). A time-varying contact rate provides a model for the impact of governmental measures on the dynamics of the pandemic. The experiments in Section 5 Iisolate the impact of the contact rate on the course of the infection counts, by assuming that and are fixed and known. The method by no means restricted to inference over a single latent force, as will also be shown in Section 5.1 In this SIRD setting, the goal is to infer an (approximate) joint posterior over 3(t) and the dynamics of S(t), I(t). R(t), and D(t) as well as to use the reconstructed dynamics to extrapolate into the future. Section Bexplains the conceptual details, Section 4 distinguishes the method from related work, and Section eevallates the performance.
Figure 2: SIRD dynamics.
3 Method
This section explains how to infer the unknown process uut) and the ODE solution .r(t) in a single forward solve. Section 3.1 defines the prior model, Section 3.2 describes the probabilistic numerical ODE inference setup, and Section 3.3 describes approximate Gaussian filtering and smoothing in this context. Section 3.4 summarizes the resulting algorithm. The exposition of classic concepts here necessarily compact. In-depth introductions can be found, c.g., in the book by Särkkä and Solin [28).

---

3.1 Gauss-Markoy prior
Let V N. Define two independent
processes U: (f0 fmax] + R and
X: [to, {max] + 1d(++) that solve the linear, time-invariant stochastic differential equations 225 dU(t) = (UU(tt dt + Lu dWu(t),
dx = (x) dt Lx (W)(t),
(4)
with drift matrices Fu Rext and Fx E [ed((+))xd(+))) well as dispersion matrices Lu Rxx and Lx E [d(x+1)xd Wu [to, maax] + IK and Wx [to, fmax] + R are Wiener processes, U and X satisfy the Gaussian initial conditions,
U(to) (mm. X((t) Nm Px).
defined by my R², Pu E R my E IRd(+t1), and Pu Rdl++t)x(e+t)) U(t) models the unknown function u(t) and can be any Mauss-Markov process that admits representation as the solution of a linear SDE with Gaussian initial conditions. X(t) = (x (t), d'
(t) E Rd((+11 dr (() Rd, i O, ........ In other words, the first element in (() is an estimate for x()), the second element is an estimate for (t), et cetera. Encoding that the state X consists of a model for sxt)) as well as its first v derivatives imposes structure on Fx and Lx (see c.g. [88). Examples include the Matérn, integrated Uaiieinn-UUlebbeck, and integrated Wiener processes; the canonical choice for probabilistic ODE solvers would be integrated Wiener processes 29
models the ODE dynamics, in light of which require (1)
19].
The class of Maruko priors inherits its wide generalizability from Gaussian process models; processes like U and X are Gaussian processes with the Markov property.
recall that While not every Gaussian process with one-dimensional input space is Markovian, a large number of descriptions of Gauss- processes emerge by translating covariance function into an (approximate) SDE representation (28) Chapter 12]. For example, this applies to (quasi-)periodic,
equared-epponentil, or rational quadratic kernels; in particular, sums and products of Gauss-Markov processes admit state-space representation [30 28). Recent research has considered approximate SDE representations of general Gaussian processes in one dimension (20). With these tools, prior knowledge over U or X can be encoded straightforwardly into the model.
3.2 Two likelihoods: for observations and for the ordinary differential equation
A functional relationship between the processes U(t), X(t) and the data yo:N emerges by combining two likelihood functions: one for the observations yo:N (recall Equation [.) and one for the ordinary differential equation. The present section formalizes both. Let T T U ToPE be the union of the observation-ges -rridonn-rrd Th which has been introduced in Section 2 and an ODE-grid TM (1g0e .... ) The name "ODE-grid" expresses that this grid contains the locations on which the ODE information will enter the inference problem, as described below.
T8as contains the locations of No in light of which the first of two observation models is
Y, ((1 ~ (HX( (1TTh)) R
(6)
for n = O.....,., N. This is a reformulation of the relationship between process x and observations Yo:N in Eq. ① in terms of X (instead of x, which is modeled by x ())). Including this first measurement model ensures that the inferred solution remains close to the data points. Tipe contains the locations on which U(t) connects to (() through the ODE. Specifically, the set of random variables Zo.M E (M++))x defined as
Zm X(1 U((oRR)) ~ 8 (x()) - -5(x (x(0) (Lon (x(0)(10))) U((22R)
where 8 is the Dirac delta, describes the discrepancy between the current estimate of the derivative of the ODE solution (i.c. (()) and its desired value (i.c. f(x(); U)), as prescribed by the vector field f. If the random variables Zo.M realize small values everywhere, (u) solves the ODE as parametrized by U. This motivates introducing artificial data points z0.M E (M)++)xd that are equal to zero, 7m 0 E Rd m 0, .... M. Conditioning the stochastic processes X and U on attaining this (artificial) zero data ensures that the inferred solution follows ODE dynamics throughout the domain. Figure 3 shows the discretized state-space model.
3

---

Data
YA
Xa
ODE
Un
U.
Up
Figure 3: Instance of the described state-space model, visualized as directed graphical model. Shaded variables are observed. Either only data, only mechanistic knowledge, or both sources of information can be conditioned on during inference (recall Figure
3.3 Approximate inference with extended Kalman filter
Both X and enter the likelihood in Eq. 7 through possibly non-linear vector field f. Therefore, the posterior distribution (recall 20:M
0)
((((t), X(t) Zo.M = 20:M: Yo. = Yo.N)
(8)
intractable, but can be approximated efficiently. Even though the problem is discretized, the posterior distribution is continuous [28] Chapter 10). There are mainly two approaches to computing tractable approximation of the intractable posterior distribution Eq. (8): approximate Gaussian filtering and smoothing (27). which computes a cheap, Gaussian approximation of this posterior, and sequential Monte Carlo methods [24], whose approximate posterior may be more descriptive, but also more expensive to compute. Like the literature on probabilistic ODE solvers [32 4 this work uses approximate Gaussian filtering and smoothing techniques for their low computational complexity.
The continuous-discred state-space model inherits its non-linearity from the ODE vector field f. Linearizing f with a first-order Taylor series expansion creates a tractable inference problem; more specifically, it gives rise to the extended Kalman filter (EKF) 13 22). Loosely speaking, if the random variable Z is large in magnitude, then and U are poor estimates for the ODE and its parameter. An EKF update, based on the first-order lincarization of f. approximately corrects this misalignment. If sufficiently many ODE measurements 20:M are available, sequence of such updates preserves sensible ODE dynamics over time. An alternative to a Taylor-series linearization the unscented transform, which yields the unscented Kalman filter [4 [i]. The computational complexity of both algorithms is linear in the number of grid points and cubic in the dimension of the
(4
state-space. Detailed implementation schemes can found, for instance, the book by Särkkä [27].
The EKF approximates the filtering distribution
p (U(t), Zoom 20m. Yo. 30 such that 4DDE 40BS 5t)..
It describes the current state of the system given all the previous measurements and allows updates in an online fashion as soon as new observations emerge. If desired, the Rauch-Tung- Tuntt-t- Striebel smoother turns the filtering distribution into an approximation of the full (smoothing) posterior (in Eq. (8)) In doing so, all observations - that is, measurements according to both Eq. (6) and Eq. 7 - are taken into account for the posterior distribution at each location l. As special cases, this setup recovers: (i) a Kalman filter/Rauch- Tung-Striebe smoother [16] if the ODE likelihood (Eq. (7)) is omitted; (ii) an ODE solver (32), if the data likelihood (Eq. (6) is omitted. In the present setting, however, both likelihoods play an important role.
3.4 Algorithm and implementation
The procedure is summarized in Algorithm The prediction step is determined by the prior and available closed-form (Appendix A.2 times which data is observed according to linear Gaussian measurement model in Eq. 6). the update step follows the rules of the standard Kalman filter. Before updating on pseudo-observation the non-linear measurement model is linearized at the predicted mean. More details are provided
according to the ODE likelihood (Eq. (7)),

---

Algorithm 1 Compute the filtering distribution by conditioning on both yooN and To:M-
Input: data yo:N, time grid T -Ns U U
vector field J, mxx Px, my, Pu
Output:
Filtering distribution
Eq.
99
Initialize Xo = - [m)x,Px) Px) and Up- = (m,, Po) for łj E T do Predict X, from Xj-1 and predict U from U if t, E Thes then update X, on y/j end if
EEq.
[Eq. (6 [Eq.
if t E TMP then linearize measurement model and update X, and U on 2j end if
end for
in Appendix A The filtering distribution can be turned into a smoothing posterior by running backwards-pass with a SruuchTuun-Strienn Striebel smoother (e.g. (27).
The computational cost of obtaining either, the filtering or the smoothing posterior, are both linear in the number of grid points and cubic in the dimension of the state-space, i.e. O((N (dd + E*)) Only single borward-backwart pass is required. If desired, the approximate Gaussian posterior can be refined iteratively by means of posterior linearization and iterated Gaussian filtering and smoothing. which yields the maximum-a cerin-a-posste a-posteriori (MAP) estimate [1]. The experiments presented in Section 5 show how single howard-backward pass already approximates the MAP estimate accurately.
4 Related work
The explained method closely relates to probabilistic ODE solvers and latent force models (37), especially the kind of latent force model that exploits the state-space formulation of the prior (9. The difference is that, in the spirit of probabilistic numerical algorithms, the mechanistic knowledge in the form of ODE is injected through the likelihood function instead of the prior. A similar approach of linking observations to mechanistic constraints has previously been used in the literature on constrained Gaussian processes [1] and gradient matching (5] (3). Probabilistic ODE solvers have been used by Kersting et al. [17] for efficient ODE inverse problem algorithms, but their approach is different the present algorithm, in which the need for iterated optimization or sampling is avoided altogether.
Latent forces and ODE solvers:
Monte Carlo methods: (Markov-chail Monte Carlo methods are also able to infer time-dependent ODE latent force from set of state observations. Options that are compatible with a setup similar the present work would include sequential Monte Carlo techniques (24), elliptical slice sampling 23). or Hamiltonian Monte Carlo [3] (for instance realized as the No-U-Turn sampler [12).. The shared disadvantage of Monte Carlo methods applied the resulting ODE inverse problem is that the complexity of obtaining a single Monte Carlo sample is of the same order of magnitude as computing the full Gaussian approximation of the posterior distribution. In Appendix we show results from parametric version of the SIRD-latent force model (using the No-U-Turn sampler as provided by NumPyro (26). This sampler requires thousands of numerical ODE solutions, compared to the single solve of our method. This fact is also reflected in the wall-clock time needed for both types of inference. While the MCMC experiment in Appendix takes in the order of hours, each experiment with our approach takes under one minute to complete. In other words, the algorithm in the present work poses an efficient yet expressive alternative to Monte Carlo methods for approximate inference with dynamical systems.
Experiments
This section describes three blocks of experiments. The implementation based on ProbNum and all experiments use a conventional, consumer-level CPU. First, a range of artificial datasets is generated by sampling ODE parameters from prior state-space model and simulating a solution of the corresponding ODE. Inference in such a controlled environment allows comparing to the ground truth, thereby assessing the quality of the approximate inference. We consider three ODE models to this end. Second, COVID-19 dataset will probe the predictive performance of the probabilistic model and the resulting approximate posterior distribution. Third, some changes to the model from the COVID-19 experiments, for instance, ensuring that the number of case counts must be positive,
5

---

will improve the interpretability (for example, of the credible intervals). Controlling the range of values that the prior state-space can realize introduces additional non-linearity into the model which can also be locally approximated by the EKF - and makes the solution more physically meaningful.
5.1 Simulated environments
As a first test for the capabilities of the proposed method, we consider three simulated environments. To this end, the training data is generated as follows. The starting point is always an initial value problem with dynamics defined by a vector field f and a Gauss-Markov prior over the dynamics x and the unknown parameters of the vector field. Then, (i) sample the vire-ngrying parameter trajectories from the Moruss prior; (ii) we solve the ODE, as parametrized by the sampled trajectories from (i), using LSODA [T] with adaptive step sizes using SciPy [33: (iii) we subsample
the ground-truth solution on uniform grid (which will become " to generate artificial state observations No:N: (iv) we add Gaussian i.i.d. noise the observations.
The procedure described above generates both a ground truth to compare to and a noisy, artificially observed data set. Given such a set of observations, Algorithm C computes posterior distribution over the true trajectories under appropriate model assumptions. In this posterior, we look for the proximity of the mean estimate to the underlying ground truth; the closer, the better. We measure this proximity in the root-mean-square error. Furthermore, the width of the posterior (expressed by the posterior covariance) should deliver an appropriate quantification of the mismatch. We report the -ttstttittcc [.. which suggests that the posterior distribution is well-calibrated if the x -tasttiit is close to the dimension d of the ground truth. Three mechanistic models serve as examples.
The first of three test problems the van-der-Pol oscillator g It has one parameter μ (sometimes referred to stiffness constant, because for large Jr, the van-der-Pol system is stiff).
Van-der-Pol:
As a prior over the dynamics we choose a xx 300. The stiffness parameter j is modeled as a Matérn-% process with lengthscale ty 10 and diffusion intensity of = 0.3. The posterior is computed on a grid from to 0 to mmax 25 units of time with step size AL = 0.025.
Wiener process with diffusion intensity
Volteraa The Lotka- Volterra equations [2] describe the change in the size of two populations, predators and prey. There are four parameters, which we call a, 6, c, and d, which describe the interaction and deatt/reproduucion rates of the populations. As a prior over the dynamics we choose a twice-integrated Wiener process with diffusion intensity ox 10. The four parameters
modeled as Matérn-V2 processes with lengthscales Éu,, = Lu, = fu. = fu. = 40. The diffusion intensities ot. o = 0.01 and ot. = ot. 0.001. The posterior is computed on grid from to 0 to max = 60 units of time with step size At- = 0.1.
SIRD: As detailed in Section 2. the SIRD model is governed by a contact rate B().. Recall that we assume a time-dependent /3 to account for governmental measures in reaction to the spread of CO The recovery rate y and fatality rate 17 are fixed 7 0.06 and 17 0.002, like they will be in the experiments with real data in Sections 5.2 and 5.3 below. As a prior over the dynamics we choose twice-integrated Wiener process with diffusion intensity ox 50. The contact rate 8 is modeled as a Matérn-3 process with lengthscale L = 14 and diffusion intensity of = 0.1. The posterior is computed on a grid from to 0 to fmax 100 units of time with step size At 0.1.
The model allows for straightforward restriction of parameter values by using link functions. The natural support for the SIRD-contact rate is the interval [0, 1]. but U()). as Gauss-Markov process, takes values on the real line. A change in the basis of (() with a logistic sigmoid function before it enters the likelihood fixes this misspecification. Similarly, the inferred in log-space to ensure positivity. is an appealing aspect of the EKF that these non- aincar transformations do not require significant adaptation of the algorithm. Instead, the EKF treats it as merely another level of lincarization of Eq. 7 Section 5.3 extends this to the state dynamics.
parameters are
The results are shown in Figure 4 On all test problems, the algorithm recovers the true states and the true latent force accurately. The recovery is not exact, which shows how the Gaussian posterior is only an approximation of the true posterior. The -ttsttiiit for the van-der- Per- stiffness parameter ju is 1.11, which lies in (0.0039, 3. 3.8415), the 90% confidence interval of the x" distribution with 1 degree of freedom. The moot-mean- square error (RMSE) to the truth is 0. 14. The ²"statistic for the Lotka- Volterra parameters 8.06, which lies (0.7107, 9.4877), the 90% confidence interval of the x² distribution with 4 degrees of freedom. The RMSE to the truth is 0.04 in log space and 0.018 in

---

State trajectories
Inferred parameters
VDP
12
10
20
25
Time
T Time
LV Population
0.06
0.06
0.03
0.03
30
30
60
30
60
Time
Time
Time
SIRD Counts
1000
0.50
0.00
Time
Time
Data
Posterior states
Posterior parameters
MAP estimate
Truth
Figure 4: State recovery simulated environments. The stiffness parameter of the van-der-Pod oscillator (top row) and the VolterVolerrr parameters (middle row) inferred accurately with appropriately high certainty. For the SIRD experiment (bottom row), the uncertainty is high, where low case counts provide little information about the latent contact rate. With more fluctuations in the observed counts, the approximated contact rate displays more certainty.
linear space. The -tssttiittic for the contact rate 3 is 0.91, which lies in (0.0039, 3.8415). the 90% confidence interval of the x" distribution with I degree of freedom. The RMSE to the truth is 0.2 in logit space and 0.033 in linear space.
5.2 COVID-19 data
We continue with the SIRD model introduced in Eq. 33, now using data collected in Germany over the period from January 22, 2020, May 27, 2021. Throughout the pandemic, the German government has imposed mitigation measures of varying severity. Together with seasonal effects, summer vacations, etc., they caused continual change in the contact rate. The next experiments aim to recover said contact rate (and the SIRD counts) from the German dataset.
The Center for Systems Science and Engineering at the Johns Hopkins University publishes daily counts of confirmed (yfaammee) recovered (ecooodred) and deceased (gfcocasee) individuals 6.. One transform this data to suit the SIRD model
Jn = Vn confirmed
R, D.. Rr == Un
recovered
D, = yn
deccased
(10)
The counts In Rn, and Dn are available for each day, starting with January 22, 2020. Assuming a constant population over time, the numbers of susceptible individuals S, are always evident from the other quantities, thus left out of the visualizations. We fix the population at / 83 783 945, based on public record. We rescale the data to cases per thousand people (CPT).
As a prior over (t)., due to its popularity in constructing probabilistic ODE solvers 32 we assume twice-integrated Wiener process. ((()) is modelled
Matémn3/2 process with length scale 75 and diffusion intensity of 0.05. The state-space model is straightforwardly extendable to sums and products of (more) processes [30] (2). Inferring parameters that are constant over time, however, is not straightforward due to potentially singular transition models (27] Section 12.3.1].
As described in Section 5.1 the contact rate is inferred in logit space. We shift the logistic sigmoid function such that it fulfills 8(0) 0.1 in which case the stationary mean U 0 translates to
stationary mean ((U) 3 0.1 of the Matém process that models the contact rate. The recovery
7

---

Data (train//val)
Posterior
Posterior contact rate
MAP estimate
Counts [cot]
0.0
Jon
May
Aug
Nov
Feb
Jan
Jann
May
Aug
Nov
Feb
Jun
Figure 5: Estimated counts of infectious cases and contact rate based on real COVID-19 data. The case counts of infectious people are scaled cases per thousand (cpt). The uncertainty over the contact rate mereases when the case counts are low. After single forward solve, the inferred mean is already close the MAP estimate. The shaded areas show the 95 credible interval and the dotted black lines are samples from the posterior.
Table I: List of selected governmental measures imposed in Germany with the aim to contain the spread of COVID-19. These events are depicted in Figures and (see column 'Mark'.. Links to the sources are provided in Appendix C
Mark
Governmental Measures
Stringent contact restrictions, partial shutdown of public life
- Continual relaxations of measures
Partial shutdown of public life ("lockdown light')
Hard lockdown, stringent contact restrictions
First nationwide decree of restrictions, increased intensification of measures
rate and mortality rate are considered known and fixed effect of the inference procedure on recovering the evolution of the contact rate C(t) = 3()). We set the mean of the Gaussian initial conditions the first data point that is available. The diffusion intensity of the prior process X(t) is set to ox = 10. The latent process and all derivatives are initialized at zero. Note that due to the logistic sigmoid transform, an initial value U, 0 amounts to an initial contact rate 9n = 0.1.
0.06 and 77 0.002 to isolate the
In the present scenario, we cannot take the SIRD model as accurate description of the underlying data but merely as a tool that aids the inference engine in recovering physically meaningful states and forces. In order to account for this model mismatch, the Dirac likelihood from Eq. is relaxed towards a Gaussian likelihood with measurement noise X2 0.01. This equals the data observation noise and thus balances the respective impact of either (misspecified) source of information. Intu- itively, adding ODE measurement noise reduces how strictly the vector field dynamics are enforced during inference and therefore avoids overconfident estimates of (t)).
The mesh-size of the ODE At = 1/24 days, i.e. ODE updates are computed on an hourly basis. The final 14 observations are excluded from the training set to serve as validation data for evaluating extrapolation behavior of the proposed method. Figure the results. The mean of the state X estimates the case counts accurately in both interpolation and extrapolation tasks. The estimated contact rate rapidly decreases around late March, remains low until fall, increases momentarily, and dampened again (compare Figure to Table Tp. The uncertainty over the estimated contact rate is high in the early beginning when the case counts are still low. It then increases again in summer and with the beginning of the extrapolation phase.
after. This aligns with set of political measures imposed by the govermment
If the experiment is taken as-is, the credibility intervals of the posterior over X(t) include negative numbers (mostly where the case counts are low and the uncertainty high, and when extrapolating). Of course, in system that models counts of people in different stages of a disease, negative numbers should be excluded altogether. The proposed method provides straightforward means address this issue. Section 5.3 explains the details.

---

scH00o Data (train/val)
Posterior 1, R.
Posterior contact rate
MAP estimate
Counts [cpt]
Counts [spt]
25
Jan
May
Aug
Nov
Feb
Ju
J60
Aug
Feb
1
2
3
5
Counts D [cpt]
0.0
0.0
Jan
May
Aug
Nov Feb
Jun
Jan
Aug Feb
Figure 6: Estimated case counts and contact rate, inferred in the logarithmic basis on real COVID- 19 and vaccination data. The case counts of infectious people are scaled to cases per thousand (cpt). Again, the uncertainty of the contact rate increases where the case counts are low. Now, the posterior credible interval is restricted to the positive reals. The shaded areas show the 95 % credible interval and the dotted black lines are samples from the posterior.
5.3 Non-negative state estimates
The following experiment evaluates how the proposed method performs in combination with state- space model that constrains the support of the dynamics. Concretely, let X(() model the logarithm of SIRD dynamics and the respective derivatives. With a slight abuse of notation, we will continue writing "X" even though it lives in a different space than in the previous sections. The structure of dynamic model is the same. The diffusion intensity the prior process X(I) is o = 0.05. The diffusion is not comparable the value in the previous section because the state dynamics moved to d log-space. Using de exp(x(t)) exp(x(t)) the ODE likelihood becomes
7mm X UDDE
N(i f(2: Sa), A2Ta),
(11)
with auxiliary quantities (recall the logistic sigmoid 9)
Ci = exp (0) (1
(( C2 = exp (1
Ca = 30==(UUmm ((EE)))))
(12)
The exponential function introduces an additional non-linearity into the state-space model, which necessitates smaller step-sizes for the ODE measurements (see below).
The observed case count data Yo.N transformed into the log-space, too, in which we assume additive, i.i.d. Gaussian noise. On the one hand, transforming the measurements into log-space implies that the measurement model for the counts remains linear; on the other hand, imposes log-normal noise model (if viewed back in "linear space"). Log-normal noise underlines how the estimated states cannot negative. Again, we scale the counts to cases per thousand.
As depicted in Figure 6, the reconstruction of the driving processes in this setting yields results that first glance, look similar the previous experiment. The states match the data points well. However, the extrapolation is more realistic in that the credible intervals encode that negative values are impossible (which is due to the log-trai trgntrraffoom.. The mean of the recovered contact rate closely resembles the estimate of the previous experiment, Again, upon implementation of strict governmental measures, the uncertainty decreases, whereas the context of relaxations, the uncertainty high.
Statement on Societal Impact
This work performs methods research to develop an efficient numerical algorithm to infer latent forces governing ordinary differential equations. As a testbed, we use data from the COVID-19 pandemic. We do so motivate and visualize the practical value of our methods. The results of this

---

algorithm, however, should not be taken as policy advice. The model used in the paper is deliberately simplistic. The presented work therefore should be misunderstood as epidemiological research. The machine learning community has, over time, frequently used data of contemporary societal concern to motivate and test new algorithmic concepts (well-knowe examples from the UCT collection include the Wisconsin Breast Cancer Dataset, the mushroom classification dataset, and the German credit data set), Our work follows in this line. Of course, if this algorithm. or any competitor, is used to derive policy advice, the underlying differential equation and latent states must be carefully considered by domain experts, which we are not.
Conclusion
By coupling mechanistic and data-driven inference so directly, the algorithm builds on the premise of probabilistic numerics that computation itself data source that does not differ, formally, from observational data. Information from observations and mechanistic knowledge (in the form of an ODE) can thus be described in the same language Bayesian filtering and smoothing. This removes the need for outer loop over multiple forward solves and thus drastically reduces the computational cost. Our experimental evaluation corroborates that the resulting approximate posterior close to the ground truth and drastically reduces computational cost over Monte Carlo alternatives. It faithfully captures multiple sources of uncertainty from the data, numerical (discretization) error, and epistemic uncertainty about the mechanism. We hope this framework helps empower practitioners, not just by reducing computational burden but also by providing a more flexible modelling platform.
Acknowledgements
The authors gratefully acknowledge financial support by the European Research Council through ERC StG Action 757275 / PANAMA; the DFG Cluster of Excellence "Machine Learning New Perspectives for Science", EXC 2064/1, project number 390727645; the German Federal Ministry of Education and Research (BMBF) through the Tübingen AI Center (FKZ: 011S1 011818039A); and funds from the Ministry of Science, Research and Arts of the State of Baden-Württente authors thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for supporting Krämer. Moreover, the authors thank Nathanael Bosch and Marvin Pförtner for valuable discussions.
The
References
[1] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan. Estimation With Applications to Tracking and Navigation: Theory Algorithms and Software, John Wiley & Sons, 2004.
[2] B. M. Bell. The iterated Kalman smoother as Gauss-Newton method. SIAM Journal on Optimization, 4(3):626-634) 1994.
[3] M. Betancourt. A conceptual introduction to Hamiltonian Monte Carlo. arXiv:1701.0243. 2017.
[4] N. Bosch, P. Hennig, and Tronarp. Calibrated adaptive probabilistic ODE solvers. In International Conference on Artificial Intelligence and Statistics, 2021.
[5] B. Calderhead, M. Girolami, and N. Lawrence. Accelerating Bayesian inference over nonlinear differential equations with Gaussian processes. In Advances Neural Information Processing Systems, 2009.
[6] E. Dong, H. Du, and L. Gardner. An interactive web-based dashboard to track COVID-19 in real time. The Lancet Infectious Diseases, 20(5):533-534, 2020.
[7] G. Giordano, F. Blanchini, R. Bruno, P. Colaneri, A. Di Filippo, Matteo, and M. Colaneri. Modelling the COVID-19 epidemic and implementation of population-we wide interventions in Italy. Nature Medicine, pages 1-6, 2020.
[8] Guckenheimer. Dynamics of the van der Pol equation. IEEE Transactions on Circuits and Systems, 27(11):983-988 1980.
10

---

[9] J. Hartikainen, M. Seppänen, and S. Särkkä. State-space inference for non-linear latent force models with application satellite orbit prediction. In International Conference on Machine Learning, 2012.
[10] H. W. Hethcote. The mathematics of infectious diseases. SIAM Review, 42(4))599-653 2000.
[11] A. Hindmarsh and L.. Petzold. LSODA, ordinary differential equation solver for stiff or non-stiff system. 2005.
[12] M. Hoffman and A. Gelman. The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(1):1993-122 2014.
[13] A. H. Jazwinski. Stochastic Processes and Filtering Theory. Academic Press, 1970.
[14) C. Jidling, N. Wahlström, A. Wills, and T. B. Schön. Linearly constrained gaussian processes. In Advances Neural Information Processing Systems, 2017.
[15] Julier and K. Uhlmann. Unscented filtering and nonlinear estimation. Proceedings of the IEEE, 92(3) 4-222, 2004.
[16] R. E. Kalman. A new approach to linear filtering and prediction problems. Transactions of the of Basic Engineering, 82(Series D):35-45, 1960.
[17] H. Kersting, N. Krämer, M. Schiegg, C. Daniel, M. Tiemann, and Hennig. Differentiable like- lihoods for fast inversion of 'likeiihood-free" dynamical systems. In International Conference on Machine Learning, 2020.
[18] H. Kersting, T. Sullivan, and P. Hennig. Convergence rates of Gaussian ODE filters. Statistics and Computing, 30(6):1791-1811 2020.
[19] N. Krämer and P. Hennig. arXiv:2012. 10106, 2020.
Stable implementation of probabilistic ODE solvers.
[20] Loper, D. M. Blei, J. Cunningham, and L. Paninski. General linear-time inference for Gaussian processes one dimension. 1xXi::2003.05554 2020.
[21] A. J. Lotka. The growth of mixed populations: two species competing for common food supply. In The Golden Age of Theoretical Ecology: 1923-1940, pages 274-286. Springer, 1978.
[22] S. Maybeck. Stochastic Models, Estimation, and Control. Academic Press, 1982.
[23] I. Murray, R. Adams, and D. MacKay. Elliptical slice sampling. International Conference on Artificial Intelligence and Statistics, 2010.
[24] C. Naesseth, F. Lindsten, and T. B. Schön. Elements of sequential Monte Carlo. Foundations and Trends@ in Machine Learning, 12(3):307-392, 2019.
[25] B. Øksendal. Stochastic Differential Equations. Springer, 2003.
[26] D. Phan, N. Pradhan, and M. Jankowiak. Composable effects for flexible and accelerated probabilistic programming in NumPyro. arXiv:1912.11554 2019.
[27] S. Särkkä. Bayesian Filtering and Smoothing. Cambridge University Press, 2013.
[28] S. Särkkä and A. Solin. Applied Stochastic Differential Equations. Cambridge University Press, 2019.
[29] M. Schober, S. Särkkä, and P. Hennig. probabilistic model for the numerical solution of initial value problems. Statistics and Computing, 29:99-122, 2019.
[30] A. Solin and S. Särkkä. Explicit link between periodic covariance functions and state space models, In International Conference on Artifcial Intelligence and Statistics, 2014.
[31] F. Tronarp, Á. F. García-Fermnndez, and S. Särkkä. Iterative filtering and smoothing in nonlinear and non-Gaussian systems using conditional moments. IEEE Signal Processing Letters, 25(3): 408-412, 2018.

---

[32] F. Tronarp. H. Kersting, S. Särkkä, and P. Hennig, Probabilistic solutions to ordinary differential equations as nonlinear Bayesian filterig: new perspective. Statistics and Computing, 29(6): 1297-1315, 2019.
[33] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, Peterson, W. Weckesser, Bright, et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261-274 2020.
[34] E. A. Wan and R. Van Der Merwe. The unscented Kalman filter for nonlinear estimation. In IEEE 2000 Adaptive Systems for Signal Processing, Communications, and Control Symposium, pages 153-158, 2000.
[35] Wenger, N. Krämer, M. Pförtner, Schmidt, N. Bosch, Effenberger, Zenn, A. Gessner, T. Karvonen, F.X. Briol, et al. ProbNum: Probabilistic numerics in Python. arXiv preprint arXiv:211..02100, 2021.
[36] Wenk, G. Abbati, M. A. Osborne, Schölkopf, A. Krause, and Bauer. ODIN: ODE- informed regression for parameter and state inference in time-continuous dynamical systems. AAAI Conference on Artificial Intelligence, 34(04):6344-6371, 2020.
[37] M. Álvarez, D. Luengo, and N. D. Lawrence. Latent force models. In International Conference on Artificial Intelligence and Statistics, 2009.
12


Reviews:
[{"review_id": "fdXhrZncKMH", "reviewer": "Reviewer_SfYp", "summary": "This paper presents a technique for computing approximate ordinary differential equation (ODE) solution posteriors. As shown in previous work [31], the state solutions of mechanistic ODE systems can be approximated by an extended Kalman filter. The authors further extend this methodology by conditioning the filter on observed sequences. Unlike its Monte Carlo alternatives, the resulting inference scheme computes approximate posteriors in a single forward-backward pass. The synthetic data experiments demonstrate that both the state and parameter trajectories are accurately inferred. The vanilla method and a variant that allows non-negative state estimation are tested on a COVID-19 dataset, implying interpretable findings that agree with government measures.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "Overall, the paper is written clearly and very easy to follow. The background is explained in sufficient detail. Figures are highly illustrative and communicate the ideas well (figure-1 might be slightly misleading as the curves seem like point estimates while the inference is probabilistic). Both the simulated and COVID experiments nicely show how high level of accuracy and calibration of uncertainty estimates. \n\nTo my understanding, the model largely builds upon [31]. Since [31] has already defined the prior process and introduced artificial data points $z_{0:M}$, the main theoretical contributions of the manuscript are (i) the additional prior on $u(t)$ and likelihood, and (ii) the inference in log space. In my opinion, these extensions are straightforward and thus NeurIPS might not be a suitable venue. That being said, I believe the COVID experiments pertain to the challenges we globally face today, which boosts the impact of this manuscript. Therefore, if the NeurIPS organizing committee particularly welcomes COVID-related submissions (about which I couldn't find any information) or I miss any major aspect of the paper, I recommend a reject and re-submission to a more suitable venue. Below are more detailed comments:\n\n- The observation model in eq. (1) is linear, which might be too restrictive for more complicated systems. A discussion on how to extend the method into non-linear observation mappings would be nice.\n- The performance of the proposed framework should be analyzed thoroughly. Although the comparison against the MCMC baseline is impressive, I wonder how difficult the addressed problems are. For instance, how does the estimation error change with (intermediate) gaps in the observed sequences? Since the inference is already based on discrete filter update equations, it would be interesting to see how simple discrete-time approximations (to the underlying true ODE system) would behave (this is not to show the merits of continuous-time modeling but to show when/if discrete approximations fail).\n- How are the prior parameters selected? An ablation study showing how, e.g., length-scale and M impact the overall performance would be nice. Similarly, does $\\lambda$ provide a trade-off between data likelihood and accuracy of ODE solutions?\n\nMinor comments:\n- Assuming fixed $\\gamma$ and $\\eta$ parameters may be too simplistic. In the context of the global pandemic, for instance, different COVID variants and the ratio of people who have developed immunity directly affect these parameters. This should at least be addressed verbally.\n- In Figures 5-6, what do the thin gray dashed curves represent (look like the \"truth\" curve in Figure 4)?\n- (minor) The state notation $X(t) = (X^{(0)}(t), ..., X^{(\\nu)}(t))$ might be updated to reflect that $X(t)$ is concatenated, e.g., $\\mathbf{X}(t) = (X^{(0)}(t), ..., X^{(\\nu)}(t))$.\n\n**Post-rebuttal update:** Thanks the authors for the rebuttal. After reading the rebuttal, other reviews/responses, and a little bit of extra thinking, I started to think the paper should be accepted. Utilizing Kalman filters to infer unknown latent forces is indeed one idea that could impact further research and developments. That being said, I have three additional comments that should provide the ground for my final rating:\n- The approach assumes a dense observed trajectory. The gradient matching-based \"ODE likelihood\" (7) would start failing with more and more sparse data. I believe the approach should be stress tested.\n- The methodology is tested only on very low-dimensional systems with very simple data likelihood models. So, I wonder how generalizable this framework is.\n- The inference heavily depends on [31]. This is not a problem but makes me question the originality of the contribution.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The limitations and potential negative societal impact are addressed. Maybe an additional sentence reminding that the findings are highly related to the assumptions made before inference would be nice.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "b93aoK4icc6", "reviewer": "Reviewer_K7YV", "summary": "Paper proposes a methodology to combine ODE with latent force models to efficiently solve ODEs and ground them with observed data. ", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "Authors utilize the recent developments in probabilistic numerics to solve  ODE with a latent force in single forward pass. Authors show the usefulness of their methodology on a non-parametric SIRD model applied to COVID19 spread in Germany. \n\nIn general experiments on simulated and real data are done properly and show usefulness of the approach.\nThe ability to solve ODE model in a single pass is admittedly very useful and something that can make the technique adoptable. However, I would like to see comparisons being done on more sophisticated models and being compared against a MCMC approach in much more detail. Right now authors do a comparison as something of an after thought. Also, all novelty is useful only if first the runtimes for the proposed approach is really short than say NUTS. Secondly, the estimates from the approach are comparable to NUTS or in case of simulated data near to real data. \n\nFinally, it seems everything authors have proposed here is already done in previous works they have referenced apart from specifically applying it to SIRD models. It might be case of misunderstanding on my side as most of the technical details have been taken away from paper and gives me impression that they have been done before. If not, I would at least like to see in SI a full blown description that can justify or outline the exact methodological contribution by authors.  \n\nAlso, I would like authors to claim less novelty about modeling time-varying contact rate. It is something standard in literature. In case of COVID  I would encourage authors to see work done by Walker, Patrick G. T., et al. \"The impact of COVID-19 and strategies for mitigation and suppression in low- and middle-income countries.\" Science, vol. 369, no. 6502, 24 July 2020, pp. 413-22, doi:10.1126/science.abc0035. The supplement has the model description Which was first published as a pre-print in March 2020. The work along this line is used as a standard model for COVID for quite long. This is essentially a problem with our field that we don't look at work done outside the echo chamber. Even the framework byWalker, Patrick G. T., et al. is opensourced here https://github.com/mrc-ide/squire. Even a more sophisticated version than the author's and Walker, Patrick G. T., et al. is used here https://www.medrxiv.org/content/10.1101/2021.01.11.21249564v1. Which is now a Science Translational Medicine article here DOI: 10.1126/scitranslmed.abg4262. All this models have a time varying contact rate in them. \n\nI have updated my score based on replies from authors and reading comments from other reviewers. My concerns have been mostly meet. So I am updating the score. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Authors have not even included them. I would first like them to even attempt it before giving any suggestion. They say we do have a model but we don't use it for policy making. This is weird the whole point about this section is authors thinking a bit more constructively about their work rather than passing the buck that it is not our concern. I would encourage authors to at least acknowledge the limitations and assumptions and then say output derived from their work might only be valid under certain specific conditions and any work deriving policy from these estimates should keep them in mind. ", "review_text": "", "extra_scores": {"ethical_concerns": "no", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "RSp1soY_22r", "reviewer": "Reviewer_RRdP", "summary": "The paper \"A Probabilistic State Space Model for Joint Inference from Differential Equations and Data\" describes a method of how to estimate latent functions in ordinary differential equations using measuremt data on observed parts of the model. This is achieved by combining a (single!) run through an ODE solver with a suitable probabilistic framework (Bayes filter).", "questions": "", "limitations": "", "rating": 8, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "The research question is highly relevant, not only due to COVID-19, the methods are suitable for this problem, and the results are good. The paper is written very clearly, both from the language aspect and from a formal point of view. Suitable experiments are conducted and clearly described. Previous work is cited. This paper is ready to be published at NeurIPS. \n\nSome minor criticism:\n - The title of the paper could be more what kind of differential equations are being used and what is inferred.\n - It might have been helpfull to mention further possible areas of application.\n - Is it possible to work with a non-autonomous differential equation? (lines 49-50)\n - For the wider machine learning audience: What is a latent force? What are latent force models? \n - Due to potential ethical concerns: is the \"SIR-type Modell\" the canonical/only choice to model COVID-19? What are the alternatives?\n - Line 64: how to estimate gamma and eta might be worth a footnote. \n - How is the ODE grid chosen? equi-distant? adaptive?\n - I would have liked the related work section earlier in the paper.\n - Don't you think that \"\n    (a) Did you state the full set of assumptions of all theoretical results? [Yes]\n    (b) Did you include complete proofs of all theoretical results? [N/A]\n   \" is a little bit contradictory? I mean one could have said more about function spaces which allow for solutions of the ODEs, existence of solution of ODEs and more aspects of the usual problems when dealing with differential equations. I do not think that this would have helped the paper.\n\nIt is a minefield to adequately cite all papers relating differential equations and machine learning. However, the following paper (and references therein) might be suitable to also be cited:\n - Chen et al., Neural ordinary differential equations\n - Dupont et al., Augmented neural odes\n - Jidling et al., Linearly constrained Gaussian processes\n - Lange-Hegermann, Algorithmic linearly constrained Gaussian processes\n - Rubanova et al., Latent odes for irregularly-sampled time series\n - Yang et al., Inference of dynamic systems from noisy and sparse data via manifold-constrained Gaussian processes", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "See ethical concerns.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "12"}, "extra_text": {"ethical_concerns": "Regarding \n\n\"Did you discuss any potential negative societal impacts of your work? [No] While the work includes an experimental setting of societal interest (the COVID pandemic), it only serves as an extended example for the utility of the algorithm. It is not used to derive policy recommendations.\"\n\nOne fundamental application of this paper seems to be the usage in modeling spread of disease. This is cause for potential ethical concerns, as political decisions and lives of many people might depend on such models. In particular, I disagree with the authors that in heated online debates such an extended example can not influece policy. However, problems in such a model could lead to real world problems: e.g. detecting changes in viral infectivity too slowly can lull politicians into safety and refusing to take timely countermeasures or e.g. predicting an unreasonaly high rate of infectivity could lead to civil liberties being taken away unnecessary. I have no direct ethical concerns with these models. I even find it helpfull that variance bands are part of the methodology. However, the impact of this paper on society is much higher than for other ML papers, and this should be clearly acknowledged.\n\nI should add that I personally feel very comfortable assessing the technical content of the paper (Gaussian processes, differential equations, solvers, ...), but have no prior knowledge of disease modelling. This justifies my downgrade in reviewer confidence: I cannot bear the responsibility of high confidentiality for such a paper, when lives are at stake and I had no prior point of contact to its application domain."}}, {"review_id": "J3KrLce1go_", "reviewer": "Reviewer_qgjY", "summary": "The paper discusses how to combine data with known differential equations of an observed system, to infer the parameters of these differential equations. To this end, the authors combine the state-space framework with ODE solvers and provide an approximate inference procedure based on extended Kalman filters.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "Though the paper doesn’t provide new techniques itself, the combination of recent developments in probabilistic ODE solvers, and combining them with a state-space framework is interesting for inference of dynamical systems. The paper is mostly well written. Integrating knowledge of the underlying dynamic system into state spaced models, and combining it with observations, is important in many applications. The paper lines out nicely the methodology, that achieves this.\nThe main critique I have, is that the paper is written very high level, and from the main text it was quite difficult to follow what e.g. is happening in Algorithm 1. I think a bit more detail on the Filter equations couldn’t harm here in Sec. 3.4. Also in the experiments some details, e.g. explicitly writing what is $X(t)$, $U(t)$ and $Y$ for the different cases would help readability in my opinion. Then another important point is that there is no baseline in the experiments. I guess it would be simple to at least show the improvement over a normal state space model, which just uses observations. (Much like it is done schematically in Fig. 1, but providing numbers). Otherwise, I think it is a good paper, and with some changes, it could be even better.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "+ It would be nice to provide some details on the Gauss-Markov priors in the appendix.\n+ I miss a discussion, how parameters could be estimated in this framework. For example, if some parameters on the system are unknown (but constant), how would you estimate them. I guess, one would think about some EM procedure, like done in standard state-space models. Also, it would help to say what emission parameters H and R (from Eq. 1) are in the experiments. I couldn’t find any details on them in the text, neither, what are observations. \n+ I noticed, that the values for the diffusion parameter $\\sigma_X^2$ were chosen quite high. How were these parameters selected? Is there a heuristic that can be followed?\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4h"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors present a method for inferring latent functions in ODE model from observations. Reviewers widely praised the manuscript for clarity and motivation. 

The discussion focussed on a technical aspect: that the method only requires a single pass of an ODE solver. This is achieved through a reworking of the problem of ODE solving as a state-space inference problem. this is an exciting insight that could change how the community thinks about solving these sorts of inverse problems. 

One reviewer raised a serious concern about the lack of a societal impact statement: the authors have provided their statement in the rebuttal with which I am satisfiedL this must make it into the camera-ready edition. 
Decision: Accept (Poster)