Inference in conditioned dynamics through causality restoration

Estimating observables from conditioned dynamics is typically computationally hard. While obtaining independent samples efficiently from unconditioned dynamics is usually feasible, most of them do not satisfy the imposed conditions and must be discarded. On the other hand, conditioning breaks the causal properties of the dynamics, which ultimately renders the sampling of the conditioned dynamics non-trivial and inefficient. In this work, a Causal Variational Approach is proposed, as an approximate method to generate independent samples from a conditioned distribution. The procedure relies on learning the parameters of a generalized dynamical model that optimally describes the conditioned distribution in a variational sense. The outcome is an effective and unconditioned dynamical model from which one can trivially obtain independent samples, effectively restoring the causality of the conditioned dynamics. The consequences are twofold: the method allows one to efficiently compute observables from the conditioned dynamics by averaging over independent samples; moreover, it provides an effective unconditioned distribution that is easy to interpret. This approximation can be applied virtually to any dynamics. The application of the method to epidemic inference is discussed in detail. The results of direct comparison with state-of-the-art inference methods, including the soft-margin approach and mean-field methods, are promising.


I. INTRODUCTION
The method we will present is rather general and applies to a wide family of stochastic processes.We will thus first describe it below in complete generality, and delay its description for a specific important application (namely the risk assessment problem in epidemic spreading processes) to the following section.Let us denote by P [x] = P [x (0) , ..., x (k∆t)] the probability distribution of trajectories x of a (known) dynamical model.Given a (hidden) realization x * , consider a set of observations O = (O 1 , . . ., O M ) sampled from a (known) conditional distribution P [O|x * ].The scope of this work is to devise an efficient method to infer information about x * given O, in particular, to be able to estimate averages over the posterior distribution Although it might be generally feasible to sample efficiently from the prior P [x], sampling from P [x|O] is normally difficult.A naive approach is given by importance sampling [1,2], that consists in evaluating the average of a function f by first generating M independent samples x 1 , . . ., x M from P [x] and then computing Unfortunately, this method is impractical when observations deviate significantly from the typical case, as for the case in which P [O|x µ ] becomes very small (or even zero), rendering the convergence of f to the true average value inefficient.
One reason for which sampling from P [x] is usually feasible is that the causal structure induced by the dynamical nature of the stochastic process can be exploited to efficiently generate trajectories.The causal property of the stochastic dynamics lies in the fact that the state of the system at a given time depends naturally (in a stochastic way) on states at previous times.When considering discrete time-steps or epochs 0, ∆t, 2∆t, . . .(in the following discussion, for simplicity of notation, we will assume ∆t = 1), this property implies that the distribution of trajectories of the stochastic dynamics assumes the following factorized form: P [x(0), . . ., x(T )] = T t=0 P [x(t)|x(t − 1), . . ., x(0)] , where in the t = 0 term the conditioning part is empty and thus the probability is unconditioned.In most models, it is computationally simple (or at least feasible) to sample x(t) from P [x(t)|x(t − 1), . . ., x(0)], implying that (3) can be exploited to generate trajectories by sequentially sampling x(0), then x(1), etc.The intrinsic difficulty associated with sampling from the conditioned distribution P [x|O] is a consequence of causality breaking [3] induced by the addition of the extra information in O.In general, P [x(t)|x(t − 1), . . ., x(0)] and P [x(t)|x(t − 1), . . ., x(0), O] are very different objects.For example, even if the former is time and space invariant, the latter is generally not, because this symmetry is typically broken by the observations.This difference ultimately implies that we cannot sample from the posterior distribution sequentially as in the unconditioned case (3).Although we can write an exact expression similar to (3), sampling from ( 4) is unfortunately still problematic.Indeed, the expression for P [x(t)|x(t − 1), . . ., x(0), O] is in general extremely difficult to compute and involves a marginalization over times t > t (with an exponential number of terms): P [x(t)|x(t − 1), . . ., x(0), O] ∝ x(t+1),...,x(T ) T t =0 P [x(t )|x(t − 1), . . ., x(0)] P [O|x(T ), . . ., x(0)] . ( This dependence on future times is in our opinion the real source of the causality breaking phenomenon.When dynamics are unconditioned, i.e. causality applies, information is intuitively flowing from past to future.Although it is a very intuitive concept, the study of information flow is actually rather involved and it opens to interesting insights into the collective interactions among agents in agent-based systems.We refer the interested reader to [4,5].The approach proposed here, called Causal Variational Approach (CVA), aims at providing a variational approximation of the posterior distribution P [x|O], for which causality features are restored and, therefore, independent samples can be efficiently generated from it.In particular, we propose to approximate P [x|O] ≈ Q(x), where Q (x(0), . . ., x(T )) = processes satisfy a spatial conditional independence property [7], namely that: P[x(t)|x(t − 1), . . ., x(0)] = N i=1 P[x i (t)|x(t − 1), ..., x(0)], (6) As this property is often crucial for efficient sampling from P [x], CVA maintains it on the approximating q t i.e. q t (x(t)|x(t − 1), . . ., x(0)) = i q i t (x i (t)|x(t − 1), . . ., x(0)).
2. Local interactions.Each variable of the prior process, moreover, might depend only on a restricted (local) set of variables on a given contact network; we choose to preserve this dependence in q t .Note that this property is in general not present in the posterior.
3. Markovianity.If the prior distribution defines a memory-less stochastic process [8], namely CVA extends this property to the approximating factors as well, q i t x i (t)|x(t − 1), . . ., x(0) = q i t x i (t)|x(t − 1) .Note that if additionally the observations are time-factorized, namely P[O|x] = CVA can be used to tackle some difficult problems emerging in the field of epidemic inference, such as epidemic risk assessment from partial and time-scattered observations of cases, or the detection of the sources of infection.These problems have been recently addressed within a Bayesian probabilistic framework using computational methods inspired by statistical physics [9][10][11], and generative neural networks [12].With respect to existing similar approaches based on variational autoencoders (e.g.[12,13]), the CVA ansatz for the posterior distribution does not employ neural networks, has comparatively a much smaller set of parameters and allows for much simpler physical interpretations.In particular, risk assessment from contact tracing data is of major importance for epidemic containment, because having access to an accurate measure of the individual risk can pave the way to effective targeted quarantine plans based on contact tracing devices [14][15][16].Moreover, in epidemic problems, there are quantities of interest that are not known a priori.An example is the infection rate of the disease.Our method can be used to compute such quantities, treated as hyperparameters of the CVA distribution.Being able to find the prior parameters of a distribution gives also the possibility to simplify the inference problem by adopting a simpler model.For example, in the context of epidemic inference, CVA allows one to study inference problems related to the SEIR model (introduced later) with an effective SI model, which is simpler than SEIR.This part, which we call model reduction is illustrated in detail in the Results section.After presenting the CVA in a general setting, its main features will be discussed by exploiting a conditioned random walk [17][18][19] as a toy model.Then, an application to the important problem of epidemic inference and risk assessment on dynamic contact networks is developed and analyzed in detail.We stress that the two reference cases, i.e. the epidemic inference and the conditioned random walk, represent two very different dynamical processes, the former continuous in time while the latter advances in discrete time-steps.

II. METHODS
The method is based on approximating the original constrained process by introducing an effective unconstrained causal process that is naturally consistent with the observations.Let Q θ (x) be the probability distribution of a generalized dynamics, parametrized by the vector θ of parameters.The best approximation to P [x|O] (in a precise variational sense) can be obtained by observing that Eq.(1) can The crucial point in this estimation is that both Q θ and P [x, O] = P [x] P [O|x] have explicit expressions that, due to their causal structure, can be efficiently computed using rejection-free sampling, at difference with P [x|O] and P [O] which do not benefit from this property.The fact that the samples are independent allows for trivial parallelism in implementation.For a detailed description of the gradient descent optimization adopted in this work, we refer the reader to Appendix E. When a fixed point is reached, the corresponding distribution Q θ (x) is the argument that (locally) minimizes the free energy in (9) therefore providing an approximation of the posterior distribution P [x|O].Finally, the result can be used to generate samples satisfying the constraints given by O and to compute interesting observables from them by efficiently computing sample averages.random walk, starting at site x (0) = 0.If the generating process is spatially homogeneous, the probability of every feasible trajectory x of length T is P [x] = 2 −T .Note that every possible trajectory can be directly sampled by means of a causal generative process, namely a discrete-time Markov chain in which the conditional probability of a jump is P [x (t + 1) = x (t) ± 1 |x (t)] = 1/2.Fig. 1(a) displays a space-time representation for a set of realizations of such an unbiased random walk (black paths).For this process, let us imagine a procedure that, given a time instant t µ and position x µ in space, can test if the trajectory was at time t µ to the left or right of position x µ , and denote the corresponding half-line as W µ ⊆ Z. Assume that for a given unknown trajectory, we have M observations of this kind O = (t µ , W µ ) µ with µ = 1, . . ., M .The posterior probability of a trajectory x can be written as where the numerator is 1 only if the trajectory x satisfies all observations, and zero otherwise.The denominator is the sum over all trajectories (so the variable y runs over the space of all the possible trajectories) of the numerator and plays the role of a normalization term for the posterior.The effect of O is to select (or constrain to) a subset of the trajectories of a free random walk, i.e. those compatible with the observations.One could naively sample trajectories from the free dynamics and then select only those compatible with O.However, as depicted in Figure 1, the fraction of trajectories compatible with the constraints might be very small to allow for a feasible computation: in the example of Fig. 1(b), where it is assumed that three regions at specific time steps (black horizontal barriers) cannot be crossed, several realizations of the unconstrained dynamics are discarded (red paths), while only a small fraction is kept (black paths).In this regard, the CVA provides an efficient way of generating trajectories compatible with the constraints by building up an effective probability distribution that is -by construction -compatible with the former.Within the framework provided by CVA, the following causal ansatz can be introduced for the conditioned random walk problem: with θ = {r t x } t=1,...,T x=−T,...,T being the set of site-dependent and time-dependent rates to jump to the right, and l t x(t) = 1 − r t x(t) the associated probabilities to jump to the left.We remark that Eq. ( 12) has the same functional form as the unconstrained distribution, i.e. it still represents the probability distribution of a random walk, but with heterogeneous (in general, both in space and time) jump rates.The distribution Q θ requires T (2T + 1) parameters, where 2T + 1 is the total number of sites that can be visited by the realizations of the random walk.These parameters are sought by minimizing the KL distance between Q θ and the posterior distribution Eq.(11).The resulting probability Q θ (x) obtained using the CVA is characterized by heterogeneous rates r t x whose dependence in time and space perfectly mirrors the constraints introduced by the barriers.The marginal distributions of the trajectories sampled from Q θ (x) are represented in Fig. 1(c), where the color gradient is associated with the marginal probability of occurrence of each step.

Epidemic Models and observations
From now on, we will consider a class of individual-based epidemic models describing a spreading process in a community of N individuals, interacting through a (possibly dynamic) contact network.The overall state of the system at time t (consisting of the state of each individual) is described by a vector x (t) ∈ X N , where X is a finite set of possible health conditions (called compartments).The simplest, but already non-trivial, model of epidemic spreading is the discrete-time Susceptible-Infected (SI) model [22], in which X = {S, I} (corresponding to an individual being "susceptible" and "infected", respectively) where the only allowed transition occurs from state S to state I.More precisely, each time t, if an infected individual j is in contact with a susceptible individual i, the former can infect the latter (which moves into state I) with a transmission probability λji (t), sometimes called transmissivity.Since transmissions are independent, the individual transition probabilities are and . A simple assumption is that time dependence only enters to describe the dynamic nature of the contact network (with λji (t) = 0 if there is no contact between j and i at epoch t).More realistically, the transmission probability λji (t) should also depend on the current stage of infection of the infector j (e.g.presence/absence of symptoms) and thus mainly on the time elapsed since her own infection, making the epidemic dynamics non-Markovian.
Assuming transmission probabilities proportional to ∆t in Eq.( 13) and defining transmission rates as λ ij (t) = lim ∆t→0 + λij (t) /∆t, a continuous-time version of the SI model is obtained.Both for the discrete-time and continuous-time SI models, the history of the epidemic process can be fully specified by the "infection times" {t i } N i=1 of all individuals; we conventionally set t i = 0 if individual i is already infected at the initial time and t i = +∞ if i is never infected during the whole epidemic process.In terms of the trajectory vector t := (t 1 , . . ., t N ), the probability pseudo-density of an epidemic history can be generally written as where γ denotes the probability of each individual to be a patient zero, and is the first-success distribution density of an event with rate f (t).Hence, the quantity Λ j =i I [t j ≤ t] λ ji (t) , t i is the probability density associated with the infection, at time t i , of individual i by one of its infectious contacts at previous times.Notice that when ´b −∞ Λ (f, t) dt < 1, it means that there is a non-zero probability of the individual remaining susceptible.In that case, we will formally assign the defect mass In addition to the epidemic model, a set of observations has to be defined.In real epidemics, observations mirror the outcomes of medical tests, namely the state of an individual i at time t.For the sake of simplicity, an auxiliary variable r ∈ {+.−} representing positive or negative tests, respectively, is defined, such that each observation can be encoded as a triplet (i, t, r).Given the stochastic nature of clinical tests, it is assumed that the outcome r of a test performed on individual i at time t obeys a known conditional distribution law P [r|t i ] where t i represents the infection time of individual i.When medical tests are affected by uncertainty, i.e. there exist non-zero false positive and false negative rates of the diagnostic tests, the conditional probability states For a population undergoing M individual test events, the set O of observations is then identified with the set of triplets (i µ , t µ , r µ ) for µ = 1, . . ., M .As in the random walk example, each observation constrains the dynamics: in the noise-less case (i.e.p FNR = p FPR = 0) the posterior distribution gives zero measure to all the epidemic trajectories that violate the observations.Given a realization of an epidemic model defined on a contact network, and the (possibly noisy) observation of the states of a subset of the individuals (at possibly different times), the epidemic risk assessment problem consists of estimating the epidemic risk, i.e. the risk of being infected, of the unobserved individuals at some specific time.In practice, it amounts to computing marginal probabilities from the posterior distribution where ´dt denotes the integral over all infections times t 1 , . . ., t N .A richer epidemic model, that is often used as a testing ground for more realistic scenarios (see e.g.Refs [23,24]), is the SEIR model [25], which includes also the Exposed (E) and Recovered (R) states, i.e.X = {S, E, I, R}.The only allowed transitions are S → E (representing the contagion event), E → I, and I → R; the latter ones occur for each individual independently of the others, with latency and recovery rates ν i and µ i , respectively.The previous representation in terms of transmission times can be straightforwardly generalized to the SEIR model (by introducing individual-wise infective and recovery times) as well as the definition of observations from clinical tests and the measure of the individual risk (see Appendix B for additional details).
The choice of the parameters θ of the CVA ansatz reflects and somehow generalizes the features of the generative model P [t]: in the SI case for instance, for each individual i, heterogeneous infection rates λ i (t), zero-patient probabilities γ i , and self-infection rates ω i (t) are defined.Since λ i (t) and ω i (t) are time-dependent quantities, an additional parametrization is introduced for computational purposes.Then, the trial distribution Q θ (t) is optimized with respect to the full set of parameters.The total number of parameters for the inference in the SI model in a population of N individuals is 7N , while for the SEIR model is 13N .We refer to Appendix C for a more detailed discussion of the parameter choice and of the implementation of the gradient descent.

Performances on synthetic networks
The performance of the CVA in reconstructing epidemic trajectories can be tested by measuring its ability in classifying the state of the unobserved individuals based on their predicted risk.The results are directly compared with those obtained using other inference techniques previously proposed in the literature, such as Belief Propagation [26] as implemented in [9] (sib), a Monte Carlo method (MC), a Soft Margin method (soft) adapted from [27], and simple heuristic methods based on Mean Field equations (MF) [9] or on sampling (heu).A description of the implementation of the several methods used for comparison is provided in the Appendix G.
For the sake of simplicity, we considered SI epidemic processes on proximity graphs [28], i.e. random graphs generated by proximity relationships between N = 50 individuals randomly drawn from a uniform distribution on a two-dimensional square region (a definition of proximity graphs is given in Appendix H).Each instance corresponds to a different realization of both the dynamical network and the forward epidemic propagation.Observations O are noiseless, i.e. p FNR = 0, and performed on a randomly chosen fraction of the population at a fixed time t obs = T .The comparison among the different inference methods is performed by ranking the individual marginal probabilities of being in state I at a chosen time t * and building the corresponding receiving operating characteristic (ROC) curve [29].The AUC (area under the ROC curve) at a time t * is an indicator of the accuracy of the method in reconstructing the state of the individuals at that time.The AUC at initial time (t * = 0, patient-zero problem) and at the final time (t * = T , risk assessment) are shown in Figure 2 as function of the number of observations available.
As one may expect, in both cases the average performances of all methods improve when the number of observations increases.In particular, the Soft Margin method is expected to converge to the exact results for this type of experiment when the number N of individuals is small.The results obtained with CVA are very close (and closer than any other technique) to those obtained by means of Soft Margin (soft), even in the interesting and challenging regime with only a few observations.
To further investigate the performances of CVA against other state-of-the-art techniques, the AUC associated with the prediction of individual risk is quantified as a function of time, with two different observation protocols: (i) the states of a fraction of individuals are observed at observation times scattered over the duration of the simulation or (ii) observations are performed at the last time t obs = T .More precisely, in the first protocol observation times are randomly drawn a priori with uniform distribution in the interval [1, T ], and observations are biased towards tested-positive outcomes to mimic a realistic scenario where symptomatic, i.e. infected individuals, are more likely tested than susceptible ones.For these experiments, two realistic dynamic contact network instances are considered, one generated using the Spatio-temporal Epidemic Model (StEM) in continuous-time in Ref. [30] and the other using the discrete-time OpenABM model in Ref. [14] (see Appendix H for a brief description of StEM and OpenABM models).For sake of simplicity, instead of adopting the complex epidemic dynamics described in Ref. [14] and Ref. [30], epidemic realizations are generated using a continuous-time SI model on these contact graphs.
A measure of the individual risk is computed according to all different methods (CVA, sib, soft, and MC), and the corresponding AUCs are shown as functions of time (in days), in Figure 3(a,c) and 3(b,d) for OpenABM and the StEM respectively.For the latter only, we also consider different MC parameterizations, in particular when using δ ∈ {12, 24, 48, 96} hours; a further increase of δ does not carry any improvement of the results.The quantity δ is associated with the MC move's proposal (see Appendix G for a detailed description).Panels (a) and (b) are associated with the observations scattered in time, while panels (c) and (d) use observations at the last time only.In panel (a) simulations are run for N = 2000, while in panel (c) we set N = 1000.It is easy to see that, in panels (a) and (c), CVA (blue dots) is the best-performing method in terms of AUC; only MC (pink triangles) reaches comparable AUC for t ∼ T .The results achieved by Belief Propagation are similar to those produced by CVA when the size of the graph is N = 1000 (panel (c)), and significantly deteriorate for N = 2000 (panel (a)).For the instances generated according to StEM in panels (b) and (d), the comparison reveals that CVA achieves the largest values of the AUC at all times and only Belief Propagation (sib, orange squares) performs comparably to CVA for the risk assessment problem, i.e. the inference at the last time of the dynamics.MC for δ ∈ {48, 96}, approaches CVA performances in the last days while is not able to predict the zero patient.Indeed, the AUC associated with MC predictions for all parametrizations is slightly larger than 0.5 for t < 5 when observations are performed at the last time of the dynamics.

Hyperparameters inference
In the previous numerical experiments, the parameters of the generative SI model, i.e. the (homogeneous) infection rate λ and the probability of being the zero patient γ, are assumed to be known.These quantities enter the CVA formalism as hyperparameters of the prior distribution which are often inaccessible in realistic applications, but can be estimated as those realizing the minimum of the free-energy F = − log P[O].This can be achieved by gradient descent if the number n obs of available observations is sufficiently large (see Appendix F for details).An example of the quality of the parameter inference is provided by the following experiment.For

Model reduction
For viral diseases with sufficiently known transmission mechanisms, agent-based modeling using discretestate stochastic processes has proven useful to build large-scale simulators of epidemic outbreaks and design containment strategies [24,30].Such mathematical models are much more complex than the SI model analyzed previously, as they need to include additional specific features of real-world diseases.In particular, models may assume different infected states, characterized by a different capability of transmitting the virus and diverse sensitivity to diagnostic tests.Another important feature that can emerge from realistic transmissions is that individuals can stop being infectious, even before recovering from the infected state, because of the decay of their viral load.These ingredients may be effectively included in the SI and SEIR models by introducing timedependent infection rates, which is a natural assumption in the framework of the Causal Variational Approach (see App. B-C).This property makes the latter a very suitable inference method to approximate unknown and possibly complex generative epidemic processes using classes of simpler probabilistic models.A simple test of such a potentiality is provided by the following example.Several epidemic realizations are generated with an SEIR prior model and the quality of the inference obtained by the Causal Variational Approach is evaluated when (i) the SEIR model is also used as an ansatz for the posterior distribution and (ii) when the posterior distribution is approximated with a simpler probabilistic model, such as the SI model.If the parameters of the generative SEIR model are known, the hyperparameters of the SEIR posterior are also known.The corresponding results (green diamonds) for the AUC as a function of time on a proximity graph are displayed in Fig. 5(a).Otherwise, the hyperparameters of the SEIR posterior can be inferred by means of the CVA (blue circles).Finally, the ansatz for the posterior distribution can be simplified to a SI model, and the corresponding hyperparameters can be inferred as well within the CVA (red squares).The overall quality of the inference depends on the possibly different regimes of information contained in the observations.Strikingly, when the generative model is not known, the results from SEIR-based and SI-based inference are always very close to each other.For a sufficiently large number of observations, such results are also close to those obtained with the SEIR posterior and known hyperparameters.
From a generative perspective, the inferred hyperparameters can be also interpreted as the epidemic parameters of some prior model, from which epidemic realizations can be sampled.It is natural to ask what are the statistical properties of such generative processes compared to the original one, from which the observations were sampled.Figure 5 also shows, in two different regimes, as a function of time the average number of infected individuals estimated from the original SEIR prior model (green diamond), the SEIR prior model with inferred hyperparameters (blue circles) and the SI prior model with inferred hyperparameters (red squares).The average number of infected individuals computed over the realizations from which the observations are sampled is also displayed (black line).The regimes shown in Fig. 5 correspond to unbiased observations (panel (b), for λ = 0.3), and to observations preferentially sampled from large outbreaks (panel (c), for λ = 0.15).Although the discrepancy between the different curves is significant, the moderate difference between predictions obtained using SEIR and SI prior models with inferred hyperparameters suggests that model reduction is only a minor source of information bias.

IV. CONCLUSIONS
Sampling from the posterior distribution of a conditioned dynamical process can be computationally hard.In this work, a novel computational method to accomplish this task, called the Causal Variational Approach, was put forward.The Causal Variational Approach is based on the idea of inferring the posterior with an effective, unconditioned, dynamical model, whose parameters can be learned by minimizing a corresponding free energy functional.An insight into the potential of the method is obtained by analyzing a one-dimensional conditioned random walk in which some regions of space are forbidden.The CVA produces a generalized random walk process, with space-dependent and time-dependent jump rates, whose unconditioned realizations satisfy the imposed constraints.An application of greater practical interest concerns epidemic inference, in particular the risk assessment from partial and time-scattered observations.For simple stochastic epidemic models, such as SI and SEIR, taking place on contact networks of moderately small size, the CVA performs better or as well as the best methods currently available.Moreover, the variational nature of the method allows one to estimate the parameters of the original epidemic model that generated the observations, which enter the CVA in the form of hyperparameters.Since the CVA approximates the posterior distribution of the epidemic process by learning a set of generalized individual-based, time-dependent parameters, even with a rather simple ansatz for the epidemic model, such as an SI model, inference from observations coming from more complex epidemic processes can be performed.In fact, a generalized SI model with time-dependent infection rates and self-infection rates allows one to accommodate many features of real-world epidemic diseases, such as time-varying viral load and transmissivity, incubation, and recovery.The performances of the Causal Variational Approach do not seem to suffer from model reduction from SEIR to SI, suggesting that simplified epidemic models could be effective for inference also in real-world cases.The Causal Variational Approach is very flexible and, employing sampling to perform estimates, it can be applied virtually to any dynamics for which the latter can be carried out efficiently.In particular, the method can thus be applied to inference problems involving recurrent epidemic processes, such as the SIS model [31] or other models (e.g.[32]) where immunity decays over time.There are, however, some limitations.The CVA relies on the fact that the functional form of the posterior should be similar to the one of the prior.This is not true in general.For example, let us take an epidemic SI model in which the zero-patient probability γ is infinitesimally small.If one individual is tested positive at a certain time, then the posterior distribution is substantially different from the prior.In particular, in the prior process each individual is the patient zero independently with probability γ, while in the posterior there is a strong (anti)correlation: indeed the measure will concentrate on trajectories with exactly one infected individual (and this is impossible to reproduce with independent patient zero probabilities).Of course, this example is extremely contrived, as the probability that infection occurs at all in this system, and thus such a test result can be obtained, is infinitesimally small as well.Moreover, in this case the problem can be simply solved by adopting a more natural distribution for the initial state (either by using a non-infinitesimal initial infection probability in the prior, or by adopting a single initial infection in the test distribution q, see also Appendix D ).Nevertheless, it is a simple example in which the prior functional form is substantially different from the one of the posterior.
The posterior distribution P[t|O] can be computed using Bayes theorem as where the denominator is given by It is convenient to compute the posterior probability that individuals are infected at the initial time given the observation O, For both individuals, non-causal effects arise as the infection probabilities at time 0 also depend on the infection rate λ and on the observation time T .Moreover, the expressions of Γ A , Γ B are different because the observation on A has broken the symmetry between the two individuals.
It is also possible to compute the posterior probability that individual B is infected at time t, which implies that A is initially infected, namely and the posterior probability that individual A is infected at time t by B, i.e. (A10) The two quantities differ because here the observation forces individual B to infect A before the observation time t obs = T , causing the infection rate Λ BA (t) to diverge when t → T .The exact posterior distribution (A5) can be expressed in the form of a generalized SI model with asymmetric probabilities Γ A and Γ B that the individuals are already infected at the initial time and time-dependent infection rates Λ AB (t) and Λ BA (t).A schematic representation of the role played by these parameters as compared with the ones of the original SI model is provided in Fig. 6 (right).Such a generalized SI model can be used as an ansatz distribution Q θ in the CVA.The minimum of the corresponding free energy is obtained when the generalized parameters Γ A , Γ B , Λ AB (t) and Λ BA (t) assume exactly the form derived above in Eqs.(A7),(A8),(A9) and (A10).It follows that the CVA provides a formally exact solution to the inference problem under study.This is shown in Figure 7, where the exact solution of the model is compared with the solution found using the CVA.FIG. 7. .Comparison between the exact calculations (dashed black lines) and results obtained using the CVA (red lines) for a SI model (λ = 0.05 and γ = 1/2) with two individuals A and B, where A is observed infected at the final time T = 10.Left: Marginal probabilities of being infected for individuals A and B as a function of time.Due to the observation on A at time T , the two marginals are different: in particular, individual A's marginal reaches 1 at time T in such a way as to be consistent with the observation.Right: Effective infection rate ΛBA(t) as function of time.The divergence is due to the observation on individual A at time T .

Appendix B: Probabilistic description of the SEIR model with observations
The Susceptible-Exposed-Infected-Recovered (SEIR) model is a generalization of the SI model in which incubation (E) and recovery (R) are included.The only allowed transitions are S → E, E → I, I → R. A transmission event by an infected individual j brings a susceptible individual i into the E state with a rate λ ji , then the transition to I occurs independently of the rest of the system with rate ν i .Finally, an infected individual i recovers, again independently of the others, with rate µ i .In this model, the trajectory of the epidemic process is fully specified by three times: , representing the times in which individual i enters the states E, I and R, respectively.For an initially infected individual t E i = t I i = 0.In this notation, the probability weight of the trajectory t = (t E i , t I i , t R i ) N i=1 can be written as ) Test-based observations can be defined by generalizing the argument discussed in Section III in the main text for the SI model.Once again, the simplest realization of test r consists of a noisy evaluation of the individual's infection state, whose outcome obeys the following stochastic expression, conditioned w.r.t. the time trajectory of individual i: where p FPR (resp.p FNR ) is the false positive (resp.negative) ratio of the test.Here it was implicitly assumed that the test outcome does not depend on any acquired immunity.Finally, the definition of the risk measure in Eq. ( 16) in the main text can be straightforwardly generalized to the SEIR model where ´dt denotes the integral over all transition time triplets Appendix C: Parametrization in the Causal Variational Approach A strong advantage of CVA is its flexibility in the parametrization of the solution ansatz.It is sufficient, indeed, to weakly generalize the prior model in order to obtain the form of the Q θ (x) to optimize.This section describes the parametrizations used for the three models whose results are presented in the main text, namely the Random Walk and the epidemic SI and SEIR models.
Random Walk -The Random Walk generative model is a discrete-time stochastic process, where at each time the walker chooses to jump right or left with a uniform probability 1/2.As discussed in Section II of the main text, the goal is to characterize the probability distribution of a constrained random walk.The CVA simply consists in assuming a more comprehensive parametrization of the generative model, in a way that, the conditioned distribution thus obtained still describes a Random Walk.In particular, instead of fixing the jump probability to 1/2, it is allowed to take different, heterogeneous values, depending on time and space.The parameters of the Q θ (x) are then simply identified with probabilities r t i to jump to the right for a walker that at time t is in position x(t) = i.The ansatz has therefore the following form, which describes an inhomogeneous time-dependent random walk.CVA reduces to optimize Q θ over the probabilities {r t i } t=1,...,T i=1,...,N .Homogeneous Markovian SI model -To characterize a homogeneous and Markovian SI one needs to specify the constant infection rate λ, together with the probability γ of being the sources of the infection at the initial time.Thus, the probability of the infection times in Eq. ( 14) of the main text reads The Causality ansatz for this model is obtained by introducing, for each individual i, effective infection rates {λ i (t)} i=1,...,N and zero-patient probabilities {γ i } i=1,...,N respectively.In other words, the conditioned, and constrained, SI dynamical model is approximated with an unconditioned but inhomogeneous and time-dependent SI model.Each λ i (t) refers to the 'incoming' infection rate at site i at time t, namely, we impose that λ ji (t) = λ i (t) for every j ∈ ∂i.The benefit carried by this parametrization is twofold: (i) this individual infection allows us to simplify the calculations and (ii) it guarantees 'outgoing' heterogeneous infection rates.A simple example can be used to illustrate this parametrization.Suppose to consider an individual i that has only two contacts with j 1 and j 2 , observed to be respectively I and S at time τ .Since j 2 is susceptible at time τ , none of its previous contacts has infected it at previous times; in our parametrization, this case can be encoded by setting λ j2 (t) = 0 for t < τ .To cope with the observation of the state of j 1 one can tune λ j1 (t) (for t < τ ) to sufficiently high values to guarantee that j 1 is in state I at the observation time.In the present work, self-infection rates {ω i } i=1,...,N are also introduced, that is the rates with which individuals can get infected without any contact with an infectious individual.Although this transition is not contemplated in the generative model, it is conveniently included to satisfy some constraints associated with the observation of infected individuals, that are hard to justify through the infection rates alone.In this case, the CVA ansatz reads Notice that Eq. (C3) is almost identical to (14) of the main text, except for the presence of the self-infection rates.The SI model considered here is continuous in time, so that the variational parameters λ i and ω i should in principle be treated as continuous functions over time as well.Since it is not possible to optimize over a function defined on a finite real domain, it is more convenient to express such rates using a family of functions, defined by a set of parameters, and optimize over the latter.A simple choice -adopted in all the results presented in the main text -is a Gaussian-like (not normalized) function for each site-dependent rate λ i and ω i .Gaussian rate will thus depend on three parameters, namely the value of the peak (labeled using µ), the standard deviation (σ), and a scale factor (p).For instance, the infection rate λ i (t) reads and analogous expressions hold for the self-infection rates ω i (t) with parameters (ω p i , ω µ i , ω σ i ).With this choice, the total set of parameters to optimize over is given by: Homogeneous Non-Markovian SI model -A more realistic scenario is that depicted by the homogeneous Non-Markovian SI model, where the infection rates are not individual-dependent (e.g.spatially homogeneous) but they can vary in time.In this case, Eq. ( 14) of the main text slightly modifies as This is quite a realistic hypothesis.In fact, when an individual gets infected, its infectivity is not constant in time in real situations but it depends on the viral load of the infectious individuals.Typically, infectivity is very low in a first time-window after the contagion, then increases according to time scales that depend on the type of pathogens, and finally, it decreases to zero.In particular, the generative rate λ depends on the time elapsed since infection, i.e. λ(t − t j ).This introduces a Non-Markovian character to the dynamic because in order to know the configuration at time t + dt it is not sufficient to know the state at time t, but a memory of all the infection times of each individual must be kept.Even though the generative model becomes more complex with respect to the Markovian case described above, the inference with CVA has almost no modifications.For this reason, the same inferential parameters introduced in the Markovian case are used, with the same interpretation.The only difference is that the effective interaction rate is now defined as (where λ 0 is a rate needed to preserve the correct dimension, but it can be numerically set to 1).It is worth noting that the effective infection rate encompasses both the time-dependent infectivity of individual j (which is given by the generative model) and the susceptibility of i to be infected (which instead is learned by CVA).
The formula for the CVA ansatz is therefore .
This manipulation is called variance reduction [13] and is known to facilitate the gradient descent.Averages are evaluated by sampling from Q θ using the implementation details given in Appendix D. Since sampling can be performed in parallel, also the gradient computation is parallelizable, due to the form of equation (E11).The learning process is therefore the following: 1.The parameters are initialized to θ 0 (we initialize them with the values of the generative model, where known.For example, in the random walk example the jump probabilities are all initialized to 1/2).
2. The derivatives of the KL divergence are performed by sampling from Q θ0 and through Eq. (E11).

The process stops when
) does not significantly change after two consecutive iterations.
In the evaluation of the KL divergence (or its derivatives), it is necessary to evaluate quantities that include log P (O | x).When O imposes a hard constraint on the trajectory -as it happens in the epidemic models with noiseless observations -the violation of one of the constraints (which is likely to occur especially in the first stages of the gradient descent) would result in log P (O | x) → −∞.To avoid this issue, it is convenient to relax (soften) all the constraints: a small acceptance p ∼ 10 −10 is introduced, in such a way that log P (O | x) has always non-diverging values.As a consequence, the distribution Q θ (x) resulting from the KL minimization gives a non-zero (yet very small) probability to events x which do not satisfy the constraint.

Appendix F: Inference of Hyperparameters
As explained in Section III in the main text, the term hyperparameters refers to parameters of a generative model.For example, the Random Walk presented in Section II in the main text has a unique hyperparameter, namely the probability of a right jump, which is set to 1/2 for all times and sites.In the SI model, the hyperparameters are the zero-patient probability γ and the infection rate λ while in the SEIR model, we need to add to the hyperparameters set the latency rate ν and the recovery rate µ.Usually, in inference problems, they are not known as the only source of information, in fact, are the observations.However, inferring the hyperparameters from observations is crucial to ensure an effective solution to the underlying inference problem.The CVA allows one to estimate them.Let us call the set of hyperparameters θ p .First, notice that in general, not only the generating distribution P [x], but also the ansatz Q θ depends on the hyperparameters θ p .The dependency on Q θ might be present when using part of the generative model in the ansatz distribution.The example implemented in this paper is the case of the Non-Markovian SI model (see Appendix C).In that case, the infectivity function depends on both the inferred parameters ({λ i } ∈ θ) and the relative infectivity of the generative model (λ ∈ θ p ), as shown in Eq. (C7).Thus, the KL divergence depends on θ p through both Q θ and P [x], indicated in the following using the notations Q θ,θp = Q θp and P = P θp .As it happens for the parameters of the variational ansatz, the hyperparameters can be inferred through a gradient descent.The derivative of the 2. it extends [9] by allowing the patient zero inference; 3. it allows one to infer the hyperparameters of the generative model and to perform epidemic reconstruction in the case of non-Markovian dynamics.

Monte Carlo
The results obtained for the SI model in Figs.2-3 in the main text are also compared to a standard Markov-Chain Monte-Carlo (MCMC) sampling for the posterior distribution.Since epidemic trajectories can be fully described in terms of the infection time vector t = (t 1 , . . ., t N ), the Markov chain defines dynamics on these continuous variables that eventually converge to a stationary distribution (i.e. the posterior).At each step of the MCMC, a node i is randomly selected and a new value of its infection time -denoted with t i -is proposed, by drawing it from a suitable kernel K (t i | t i ).We adopted a Gaussian kernel centered at t i , with fixed standard deviation δ.The proposed value t i sampled from K is then clamped within the interval [0, T ]: this procedure allows to sample with non-zero probability the two extreme values, associated to node i being a patient-zero (when t i = 0), or node i never being infected (in which case its infection time is formally set to infinity, as previously discussed).This procedure is equivalent to use as proposal kernel a 2-sided rectified Gaussian distribution in the window [0, T ].The acceptance probability of such a proposal is computed as where t i represents an epidemic trajectory where only t i is modified to the corresponding proposed value, namely t i = (t 1 , . . ., t i , . . ., t N ), and K is the rectified kernel given by: K (x | t i , δ) = Φ (0; t i , δ) δ (x) + I [x ∈ (0, T )] K (x; t i , δ) + [1 − Φ (T ; t i , δ)] δ (x − T ) (G2) with K (x; µ, σ) being the probability density of a Gaussian random variable and Φ (x; µ, σ) its cumulative function.The initial condition for the Markov Chain is sampled from the prior distribution of the SI model.In order to have a fair comparison, the number of samples collected during the MCMC is equal to the one used by CVA.To diminish the effect of initial equilibration time an initial number of steps is typically required to let the MC forget the initial condition and sample efficiently the posterior distribution; notice that a "step" consists in proposing a move for each node in a random permutation.
Soft-Margin -The Soft-Margin estimator is described in [27].We adapted this method by sampling from the prior probability distribution P[x] and weighting each sample with the observation likelihood P[O|x], in which we introduced a small artificial noise in the form of a false rate, which softened the constraints, improving the performances of the method.The technique is asymptotically exact.However, when the population size grows, the probability to sample a trajectory x which satisfies the observation constraints P[O|x] dramatically decreases.Therefore, we expect the method to be quite slow for large population sizes in order to well-approximate the posterior distribution.
Appendix H: Models of contact networks used in the numerical simulations Proximity model -The proximity contact graph is obtained by distributing N individuals uniformly in a square of side √ N and generating the contacts as follows.A contact can be established between two individuals i and j with a probability e −dij / , where d ij is the Euclidean distance between the points i and j are located and is a length scale that can be tuned to change the density of the contact graph.
OpenABM-Covid19 -A synthetic contact network has been generated through OpenABM, a platform used to set up realistic epidemic instances on dynamic contact networks where a set of intervention measures have been applied and studied in Ref. [24].The underline contact pattern is a superposition of two static graphs, one a complete graph representing interaction within households and a second small-word network mirroring occupation relationships.A random time-varying network is instead used and rebuilt on a daily basis to model contacts in public transportation, transient social gatherings, etc.The number of interactions through the (I4) where the symbol ∝ is intended with respect to x s so all terms that do not depend on x s can be removed and we defined: Where again, any term that is constant with respect to x s was removed.As P[x(s)|x(s − 1), O] ∝ P[x(s)|x(s − 1), . . ., x(0), O] and both are normalized with respect to x s , they must be identical.Therefore, using equation (4) we conclude the proof.We can conclude that, at fixed observations, the posterior distribution at a given time s can be represented as a Markov process, i.e. depending only on the previous time s − 1.However, the transition probabilities between these two times will depend on the observations at future times, i.e. t ≥ s and a closed formula would require to recast f (x(s), x(s − 1) | O) = P[x(s)|x(s − 1), O], which involves an exponential number of terms to be computed (remember that x is the full state of all the degrees of freedom).In this perspective, the CVA further assumes that even these transition probabilities factorize, i.e. there is a spatial conditional independence over the posterior, a necessary ingredient for sampling efficiency.
T t=0 P[O t |x(t)], then it can be shown (See Appendix I) that Markovianity extends to the posterior distribution, P[x(t)|x(t − 1), . . ., x(0), O] = P[x(t)|x(t − 1), O].There are simple but instructive examples where CVA leads to the exact posterior, see for example the SI epidemic model for N = 2 individuals in Appendix A.

FIG. 1 .
FIG. 1. Panel (a) Unconditioned homogeneous random walk on a one-dimensional lattice.Time is reported on the vertical axis (up to T = 40) and the spatial coordinate x is on the horizontal axis.Panel (b) Some trajectories are sampled from the unconditioned homogeneous distribution.The black (red) ones (do not) satisfy the constraints, i.e. they (do not) avoid the black horizontal barriers.The fraction of feasible trajectories among a given pool can be numerically estimated, and it approaches 10 −6 .In other words, only one of a million trajectories sampled from the unconditioned distribution satisfies the constraint.Panel (c) The distribution of the trajectories sampled from the CVA distribution.The color of each pixel indicates the probability for a trajectory to visit the corresponding state at a specific time.

FIG. 2 .
FIG.2.Area under the ROC (AUC) as a function of the number of observations for the risk assessment problem, i.e. t = T , in panel (a), and for the patient-zero problem, t = 0, panel (b).The simulated contact graph is a proximity network with average connectivity 2.2/N .For both simulations in panels (a) and (b), the total number of individuals is N = 50, the probability of being the zero patient is set to γ = 1/N , and the infection rate is λ = 0.1.For each epidemic realization, the inference is performed for an increasing number of noiseless observations (here pFNR = 0) at time t obs = T .Thick lines and shaded areas indicate the averages and the standard errors computed over 40 different instances.

FIG. 3 .
FIG. 3. AUC associated with the prediction of the infected individuals, for the Causal Variational Approach (CVA), Belief Propagation (sib) and SoftMargin (soft), and MCMC (MC) as a function of time during the epidemic propagation of a SI model on several instances of dynamic contact network generated using the OpenABM model [14] (in panel (a) N = 2000, in (c) N = 1000) and the StEM in Ref. [30] (panels (b) and (d)) for N = 904.The infection rate is set to λ = 0.15 for the latter and λ = 0.02 for the former; observations are noiseless in both cases.For panels (c) and (d), observations are performed at the last time of the dynamics, i.e. t obs = T .For the results in panels (a) and (b) observation times are extracted uniformly in the range [1, T ]; at each observation time t obs , infected nodes are observed with a biased probability equal to 1.1 × NI (t obs ) /N where NI (t obs ) is the number of infected individuals at time t obs and N is the total number of individuals.The total number of observations is n obs = N • 0.1 for OpenABM and n obs = 100 for the StEM.

Figure 4 (
b) displays a scatter plot of the inferred values against the true ones.Results suggest a good agreement between the result of the inference and the generative process.

FIG. 4 .
FIG. 4. Panel (a) Heat map of the free energy (F := − log P(O)) computed at the convergence of CVA as a function of the assumed hyperparameters of the generative SI model.The experiment is performed on a proximity graph with N = 50 individuals and density ρ = 2/N ; the epidemic model is characterized by the zero-patient probability γ * = 1/N and the infection rate λ * = 0.1, shown here as a green star.We perform a large number of observations (n obs = 2N ) at uniformly randomly distributed times.As expected, the lowest values of this free energy are concentrated around the exact value (γ * , λ * ).The oriented paths (white arrows) represent the convergence towards the minimum of − log P[O] obtained by performing a gradient descent algorithm over the hyperparameters starting from three different initial points in the plane (γ, λ).Panel (b) Scatter plot of inferred values for the infection probability against the ground truth.In these experiments, we fix and assume to know the zero patient probability γ = 1/N while the infection parameter λ is varied.For each λ an epidemic simulation is performed and n obs = 10N observations are taken at uniformly randomly distributed times.

FIG. 5 .
FIG. 5. Effects of model reduction on inferential performances and generative capabilities.The numerical experiments are performed on a proximity graph with N = 100 individuals and density 2.2/N .The observed epidemic realizations are generated using an SEIR model with γ = 1/N , λ = 0.3 (panels (a) and (b)) and 0.15 (panel (c)), latency delay ν = 0.5 and recovery delay µ = 0.1.Panel (a) Values of the AUC as a function of time obtained using the CVA in two observation regimes (when the number of observations is n obs = N/10 and n obs = N/2), with the three different inferred posterior distributions: an SEIR model with known hyperparameters (green diamonds), an SEIR model with unknown hyperparameters (blue circles), and a SI model with unknown hyperparameters (red squares).Shaded areas represent the error around the average value, computed using 22 instances.Panels (b) and (c) The average fraction of infected individuals as a function of time estimated using the correct SEIR prior model (green diamonds), an SEIR prior with the inferred hyperparameters (blue circles), and a SI prior model with the inferred hyperparameters (red squares).The regimes shown correspond to unbiased observations (center, for λ = 0.3), and to observations preferentially sampled from large outbreaks (right, for λ = 0.15).The black curves represent the same quantity computed from the observed epidemic realizations.Shaded areas represent the standard error computed from 40 realizations of the dynamics.