Learning Stochastic Behaviour from Aggregate Data
Authors: Shaojun Ma
Shu Liu
Hongyuan Zha
Hao-Min Zhou
Venue:  (2021)
Keywords: Fokker Planck Equation
weak form
Wasserstein GAN

Abstract: Learning nonlinear dynamics from aggregate data is a challenging problem since the full trajectory of each individual is not available, namely, the individual observed at one time point may not be observed at next time point, or the identity of individual is unavailable. This is in sharp contrast to learning dynamics with trajectory data, on which the majority of existing methods are based. We propose a novel method using the weak form of Fokker Planck Equation (FPE) to describe density evolution of data in a sampling form, which is then combined with Wasserstein generative adversarial network (WGAN) in training process. In such a sample-based framework we are able to study nonlinear dynamics from aggregate data without solving the partial differential equation (PDE). The model can also handle high dimensional cases with the help of deep neural networks. We demonstrate our approach in the context of a series of synthetic and real-world data sets.

Under review as a conference paper at ICLR 2021
LEARNING STOCHASTIC BEHAVIOUR FROM AGGREGATE DATA
Anonymous authors Paper under double-blind review
ABSTRACT
Learning nonlinear dynamics from aggregate data is a challenging problem since the full trajectory of each individual is not available, namely, the individual ob- served one time point may not be observed next time point, the identity of individual is unavailable. This is in sharp contrast to learning dynamics with trajectory data, on which the majority of existing methods are based. We pro- pose a novel method using the weak form of Fokker Planck Equation (FPE) to describe density evolution of data in a sampling form, which is then combined with Wasserstein generative adversarial network (WGAN) in training process. In such a sample-based framework we are able to study nonlinear dynamics from aggregate data without solving the partial differential equation (PDE). The model can also handle high dimensional cases with the help of deep neural networks. We demonstrate our approach in the context of series of synthetic and real-world data sets.
INTRODUCTION
In the context of dynamic system, Aggregate data refers to the data sets that full trajectory of each individual is not available, meaning that there is no known individual level correspondence. Typical examples include data sets collected for DNA evolution, social gathering, density in control problems, bird migration during which it is impossible to identify individual bird, and many more. In those applications, some observed individuals one time point may be unobserved at the next time spot, or when the individual identities are blocked or unavailable due to various technical and ethical reasons. Rather than inferring the exact information for each individual, the main objective of learning dynamics in aggregate data is to recover and predict the evolution of distribution of all individuals together. Trajectory data, in contrast, is a kind of data that we are able to acquire the information of each individual all the time, although some studies considered the case that some individual trajectories are partially missing. However, the identities of those individuals, whenever they are observable, is always assumed available. For example, stock price, weather, customer behaviors and most training data sets for computer vision and natural language processing.
There are many popular models to learn dynamics of full- trajectejerr data. Typical ones include Hidden Markov Model (EMM))Alshamae et al. 2019 Eddy 1996), Kalman Filter (FFraaarahi & Yazdi
2020 Harvey 1990, Kalman 1 1960) and Partiele Filter (PF) (Santos et al., 2019: Djuric et al. [2003] as well as the models built upon HMM, KF and PF(Deriche et al. 2020 Fang et al.] [0199] Hefny al., 2015 Langford et al. 2009) they all require full trajectories of each individual, which may not be applicable in the aggregate data situations. On the other side, only a few methods are focused on aggregated data in the recent learning literature. In the work of Hashimoto al. (2016), authors assumed that the hidden dynamic of particles follows stochastic differential equation(SDE) in particular, they use a recurrent neural network to parameterize the drift term. Furthermore, Wang et al. (2018) improved traditional HMM model by using an SDE to describe the evolving process of hidden states. To the best of our knowledge, there is no method directly learning the evolution of the density of objects from aggregate data yet.
We propose to learn the dynamics of density through the weak form of Fokker Planck Equation (FPE), which is parabolic partial differential equation (PDE) governing many dynamical systems subject to random noise perturbations, including the typical SDE models existing studies. Our learning is accomplished by minimizing the Wasserstein distance between predicted distribution given by FPE and the empirical distribution from data samples. Meanwhile we utilize neural networks to handle higher dimensional cases. More importantly, by leveraging the framework of Wasserstein Generative

---

Under review as a conference paper at ICLR 2021
Adversarial Network (WGAN) (Arjovsky et al., 2017), model is capable of approximating the distribution of samples at different time points without solving the SDE or FPE. More specifically, treat the drift coefficient, the goal of learning, in the FPE generator, and the test function in the weak form of FPE as a discriminator. In other words, our method can also be regarded as a data-driven method estimate transport coefficient in FPE, which corresponds the drift terms in SDEs. Additionally, though we treat diffusion term as constant in our model, it is straightforward generalize it to be neural network well, which can be extension of this work. We would like to mention that several methods of solving SDE and FPE (Weinan et al. 2017 Beck et al. 2018)] Li al. 2019) adopt opposite ways to our method, they utilize neural networks to estimate the distribution ((x,() with given drift and diffusion terms.
In conclusion, our contributions are:
We design algorithm that is able recover the density evolution of nonlinear dynamics via minimizing the Wasserstein discrepancy between real aggregate data and our generated data.
By leveraging the weak form of FPE, are able compute the Wasserstein distance directly without solving the FPE.
Finally, we demonstrate the accuracy and the effectiveness of our algorithm by several synthetic and real-world examples.
2 PROPOSED METHOD
2.1 FOKKER PLANCK EQUATION FOR THE DENSITY EVOLUTION
8(Xo. fo)
g(X,.))
(a-1 bn-1)
Xo
Xi
Xn
Figure 1: State model of the stochastic process X,
We assume the individuals evolve in pattern in the space RD shown in Figure One example satisfying such process is the stochastic differential eSDation(DDE), which is also known as the process (Øksendal] 2003): d.X, = g(X,, f)dt + odW, Here d',, represents an infinitesimal change of [X,] along with time increment dr, g(: f) = (g' (..1)....g (.. " is the drift term (drifting vector field) that drives the dynamics of the SDE, is the diffusion constant, (W,) is the standard Brownian Motion. Then the probability density of (X,) is governed by the Fokker Planck Equation(FPPF)
(Risken & Caugheyz 1991), stated below in Lemma
Lemma 1. Suppose (X, solves the SDE dX, g(X,, t)dt odW, denote pl,, 1) the probability density of the random variable X,, Then p(x, t) solves the following equation:
0p(x. 1)
Σ mm1
1 dx, (p(x,))(e.)) + + Σ lml
Ox p(..
(I)
Ō1
As linear evolution PDE, FPE describes the evolution of density function of the stochastic process driven by SDE, Due to this reason, FPE plays crucial role in stochastic calculus, statistical physics and modeling (Nelson 1985 Qi Majda 2016) Risken 1989). importance is also drawing more attention among statistic and machine learning communities (Liu & Wang 2016 Pavon et al. 2018 Rezende & Mohamed| 2015). In this paper, we utilize the weak form of FPE as a basis to study hidden dynamics of the time evolving aggregated data without solving FPE.
Our task can be described as: assume that the individuals evolve with the process indicated by Figure [ which can be simulated by Itô process. Then given observations æ, along time axis, aim to recover the drift coefficient g(x, t) in FPE, and thus we are able to recover and predict the density evolution of such dynamic. For simplicity we treat g(ær, f) as function uncorrelated to time namely,
2

---

Under review as a conference paper at ICLR 2021
g(z, g(z). Notice that though evolving process of individuals can be simulated by Itō process, in reality since we lose identity information of individuals, the observed data become aggregate data, thus need new way other than traditional methods to study the swarm's m' distribution.
2.2 WEAK FORM FOKKER PLANCK EQUATION
Given FPE stated in Lemma I if we multiply test function f (R.)) on both sides of the FPE, then the integration on both sides leads to:
ʃ [op f(x)dx
D --
a ōx,
((( f(e) + 0²
8 ūx, p(x. fa(dd.
Then integrating by parts on the right hand side leads to the weak form of FPE:
S dp 0r
S
a
1 _02 L 2
D Φ ax, f(z)p(x.1)dz.
f(x)ddx=
g g'(x) f(z)p(a, ndx 1=1 dx,
where H. () denote the Sobolev space. The first advantage of weak solution is that the solution of a PDE usually requires strong regularity and thus may not exist in the classical sense for a certain group of equations, however, the weak solution has fewer regularity requirements and thus their existence are guaranteed for a much larger classes of equations. The second advantage is that the weak formulation may provide new perspectives for numerically solving PDEs (Zienkiewicz & Cheung 1971), (Srrignano & Spiliopoulos) 2018) and (Zang et al. 2019) etc.
Suppose the observed samples time points 4m-1 and Im are following true densities p(: 1m-1) and p(,, fm). Consider the following SDE
dx, = (X,)dt + odW, 1m-1 S fm
X ~ P(=, fm-1)
(2)
Here 8u an approximation to the real drift term g. in our research, we treat &w as a neural network with parameters w. Let us now denote 1) as the density of X,.
Then it is natural to compute and minimize the discrepancy between our approximated density p(- Im) and true density (m), (m), within this process, we are optimizing ga, and thus will recover the true drift term g- In our research, we choose the Wasserstein-1 distance as our discrepancy function (Villani 2008) (Arjovsky et al. 2017). Now we apply Kantorovich- Rubinsttin duality (Villani 2008). this leads to:
W (p( 1m). ((..m)) > sup (Be_pptaa))((a..)) - //viis
The first term Se.-pte..ll in Equation B can be conveniently computed by Monte-Carde method since we are already provided with the real data points x, (m). To evaluate the second term, we first approximate p(:, (m) by trapezoidal rule (Atkinson 2008):
At 10p(x, 1m-1)
0p(x, fm) Ōt
p(z, fm) x p(x, +
here Ar = Im fm-1-
(4)
2 ōt
And thus we can compute:
S [(or)/(O, 4m- 1dd
Ar 2
At
(Sople
(ap(x,mm-)) 1m-1) f(x)d oppcx, fm)) M())d
Be, Be.
(5)
Ōt
In the above Equation (5), the second and the third term on the right-hand side can be reformulated via the weak form of FPE and we are able to derive a computable formulation for Wi (p(, ). pl, (m)). Furthermore, we can use Monte-Carly method to approximate the expectations in (5). the first and the second terms can be directly computed via data points from P(-, fm-).. For the third term, we need to generate samples from p(-, Im), achieve this, we apply Buler-Maruyan scheme (Kloeden Platen 2013) to SDE (2) in order to acquire our desired samples R..;
T t + (a ) + VAiz here N(0. J) Pinm ~ (c..m-1)..
Here N(0, /) is the standard Gaussian distribution on R Now we summarize these results in Proposition

---

Under review as a conference paper at ICLR 2021
Proposition 1. For set of points X (xx1) ..... x0 in R.. We denote:
J,()-
1
(z)li)) 0x,
a
f(olll))
+ 0²
1
Σ
0²
f(xl))
N
2
1=1
El
E=l
then at time point im the Wasserstein distance between (mm 1m) and pC, 1m) can be approximated by:
1
At
Wi(p(.. Im.. bm)) sup
fa
(F 1)
19s1s1
N
Here (a
(1) P(. 1m-1). (3 ~ fm) and we denote Xm-1-
= 13 (n-1
2
2 1. 8m =12 2 (00)
where each ( computed by Euler- Marer-Maayaama scheme 66.
2.3 WASSERSTEIN DISTANCE ON TIME SERIES
In real cases, it is not realistic to observe the data at arbitrary two consecutive time nodes, especially when is small. To make our model more flexible, we should extend our formulation so that we are able to plug in observed data at arbitrary time points. To be more precise, suppose observe data
set 8, (4: .... ( J different time points lo, I1 ........ And denote the generated data
set as X, = (T .... 1 (N)) here each a:!/ is derived from the n-step Fuler-Maruyama scheme:
R Fi + ga Ba ) o VAIz, with I) 0 Osján, n, So - Pc.10)
Let us denote p(..t) f) as the solution to FPE ① with replaced by & and with initial condition p(. to) p(., fo). then the approximation formula for evaluating the Wasserstein distance Wi(p(, (n). f.)) is provided in the following proposition:
Proposition 2. Suppose we keep all the notations defined as above, then we have the approximation:
e-1 (2, +F + Σ Fr(K,))
(((((..))P(.L)) = sup INSISI
Σ ((
N
N
R=
Minimizing the Objective Function: Based on Proposition 2 we obtain objective function by summing up the accumulated Wasserstein distances among observations along the time axis. Thus, our ultimate goal is to minimize the following objection function:
N
N
min
sup Will
Σ k=1
Σ
f
At 2 (F
Σ Fr (aa)))
N
N
Notice that since we have observations on J distinct time points(despite the initial point), for each time point we compute Wasserstein distance with the help of the dual function fu. thus we will involve test functions in total. In our actual implementation, we will choose these dual functions as neural networks. We call our algorithm Fokker Planck Process(FPP), the entire procedure is shown Algorithm |
We also provide an error analysis in Appendix.
3 EXPERIMENTS
In this section, we evaluate our model on various synthetic and realistic data sets by employing Algorithm [
We generate samples ã,, and make all predictions base on Equation (6) starting with dn-
We compare our model with two recently proposed methods. One model (NN) adopts recurrent neural netWork(RNN) to learn dynamics directly from observations of aggregate data (Hashimoto et al. 2016). The other one model (LEGEND) learns dynamics in a HMM framework (Wang et al. 2018)) The baselines in experiments are two typical representatives that have state- of-the-art performance on learning aggregate data. Furthermore, though we simulate the evolving process of the data as a SDE, which on the same track with NN, mentioned before, NN trains its RNN via optimizing Sinkhorn distance (Cuturi 2013). our model starts with a view of weak form of PDE, focuses more WGAN framework and easier computation.
Baselines:
4

---

Under review as a conference paper at ICLR 2021
Algorithm 1 Fokker Planck Process Algorithm
Require: Initialize fa, (1 S n S J).
Require: Set ef, as the inner loop learning rate for Sa, and ex as the outer loop learning rate for go
1: for # training iterations do
2: for steps do
3:
for observed time 1, in (f)..... ..... ty] do
4:
Compute the generated data set X, from Euler- Maller-a scheme 7 for S SSSJ
5:
Acquire data sets X, = 12 ..... end for For each dual function fa. compute:
( from real distribution p(-, r,) for 1 = S
6:
7:
8:
Fr Fr (8t) + Fr (X,,) 2 Σ F (X,)
Update each Sa, by 8, - 0g 5% V. EN So, (8 ) EN So. (2. F))
9:
10: end for
11: Update 8w by w- V EΣ fo (2.
* (2 (
12:
end for
3.1 SYNTHETIC DATA
We first evaluate our model on three synthetic data sets which are generated by three artificial dynamics: Synthetic-1, Synthetic-2 and Synthetic-3.
Experiment Setup: In all synthetic data experiments, we set the drift term g and the discriminator two simple fully-connected networks. The g network has one hidden layer and the f network has three hidden layers. Each layer has 32 nodes for both and f. The only one activation function n. the activation function of must be twice differentiable to avoid loss of weight gradient. In terms of training process, we use the Adam optimizer (Kingma & Ba. 2014) with learning rate 10-4, Furthermore, we use spectral normalization to realize |V fll s Miiiyto et al. 2018). We initialize the weights with Xavier initialization(Glome & Bengio 2010) and train our model by Algorithm
we choose is Tanh. Notice that since we need calculate #'f
Synthetic-1:
t0 N(O. Eo).
PuAsi f.,- -AA + 6) VAN(O.)
Synthetic-2:
④n N(O. En).
sun -, + VAIN(O.
Ni
(2)
m1)
(3) 1121)
G
0i Ni+No
0. Mi+Ni
N: N No (2? - m112)
1 No
(2? M22)
No+No
(R)
m11)
(o) - =12)²
(22 12211
(27
112 )2
Ni
exp
N₂ V2R02
exp
V2ROI
20²
20
20
20
Synthetic-3 (Van der Pol uscillator
Li 2018
to N(O, Ex)
1 $ ,A x/ 10 x -4-1())) (s.)) + + Ar VA,N(O,
2 = 2 3(1 2/)Ar + VAIN(O.1)
The data size at each time point N = 2000, the dimension of the data is = 2. We treat 1200 data points the training set and the other 800 data points as the test set. In Synthetic- the data is following a simple linear dynamic, where = [(4. 0), (0, 1)], b = [-12, - -12"". We let Ar = 0.01, Σn (1, (0, 1)]. We utilize true Da 2:20 and 2:200 in training process and predict the
5

---

Under review as a conference paper at ICLR 2021
(a) Syn-l: at 10A/
(b) Syn-1: at S0A/
(c) Syn-1: at 500A/
(d) Syn-2: at 2
(e) Syn-2: at 4A/
(f) Syn-2: at 10Ar
(g) Syn-3: at 10Ar
(h) Syn-3: at 30
(i) Syn-3: at 50Ar
Figure 2: Comparison of generated data(blue) and ground truth(red) of Synthetic- ((a) to (c)), Synthetic-2((d) to (f)) and Synthetic-3((g) to (i)). In each case, it finally converges to a stationary distribution.
distributions of 250 and Zsoo- In Synthetic-2, the data is following a complex nonlinear dynamic. We let A 0.01, o = 0 = 022 = μ4 = [15, 15J" and μ4 = [-15,- -15"" We utilize true æa, o3 and 26 in training process and predict 22, 2t and 20 In Synthetic-3, Ar = 0.01,
1, we utilize true Ø3. 2y and 222 in training process then predict the distributions of æ æ I and 2 We also consider cases in higher dimensions: D 6 and 10. In each high dimensional case, to be convenience, we set every two dimensions follow the dynamics of low dimensional case(D = 2) in each data set. Notice that in Syn-2 and Syn-3, 2/ represents the i-th dimension of at..
Results: We first show the capability of our model for learning hidden dynamics of low-dimensional (d data. As visualized in Figure [ the generated data(blue) covers all areas of ground truth(red), which demonstrates that our model is able to precisely learn the dynamics and correctly predict the future distribution of data. The samples we predict converge to stationary distributions finally (as the ground truths suggest). We then evaluate three models using Wasserstein distance error metric for both low- www-dimensional (d 2) and high- nigh-diienssonal (d = 6, 10) data. As reported in Appendix Table I our model achieves lower Wasserstein error than the two baseline models in all cases.
3.2 REALISTIC DATA RNA SEQUENCE OF SINGLE CELL
In this section, we evaluate our model on a realistic biology data set called Single-cell RNA-seq(Klein ett al. 2015), which is typically used for learning the evolvement of cell differentiation. The cell population begins to differentiate at day 0 (DO). Single-cell RNA-seq observations are then sampled day 0 (DO), day (DD2,, day 4 (D4) and day (D7). At each time point, the expression of 24,175 genes of several hundreds cells are measured (933, 303, 683 and 798 cells on DO, D2, D4 and D) respectively). Notice that there is only whole group's distribution but no trajectory information of each gene on different days. We pick 10 gene markers out of 24,175 to make a 10 dimensional data set. In first task we treat gene expression at DO D4 and D7 as training data learn the hidden dynamic and predict the distribution of gene expression at D2. In second task we train the model with gene expression at DO, D2 and D4, then predict the distribution of gene expression D7. We plot the prediction results of two out of ten markers, i.e. Mtl and Mt2 in Figure
Experiment Setup: We set both f and g as fully connected three-hidden tmee-hiddee-llyerr neural networks, each layer has 64 nodes. The only activation function we choose is Tanh. The other setups of neural networks and training process are the same with the ones we use in Synthetic data. Notice that in realistic cases, At and T/At become hyperparameters, here we choose A = 0.05, T/Ar= 35, which means the data evolves 10Ar from DO to D2, then 10Ar from D2 to D4 and finally 15Ar from D4 to D7. For preprocessing, we apply standard normalization procedures (Hicks et al. 2015) to correct batch effects and use non-neet negaaivee matrix factorization to impute missing expression Eevels(Hashimota ett al. 2016 Wang et al. 2018).

---

Under review as a conference paper at ICLR 2021
(a) D2 of Mt1
(b) D7 of Mt1
(c) D22 of Mr2
(d) D7 of Mr2
(e) Correlation on D2 (f) Correlation on D7 (g) W-loss of Mtl, D2(h) W-loss of Mtl, D7
Figure 3: (a) and (b): Wasserstein loss as iteration increases of Mtl. (c) to (f): The performance comparision among different models D)2 and D7 of Mtl and Mt2. (g) and (h): True (red) and predicted (blue) correlations between (x-(xsa and Mt2((-axis) on D2 (left) and D7 (right).
As shown Table in Appendix, when compared to other baselines, our model achieves lower Wasserstein error both Mtl and Mt2 data, which proves that our model is capable of learning the hidden dynamics of the two studied gene expressions. In Figure (a) to (d). we visualized the predicted distributions of the two genes. The distributions of Mtl and Mt2 predicted by model (curves in blue) are closer to the true distributions (curves in red) on both D2 and D7. Furthermore, model precisely indicates the correlations between Mtl and Mt2, shown in Figure (e) and (f), which also demonstrates the effectiveness of our model since closer to the true correlation represents better performance. In Figure (g) and (h), we see that with simpler structure, the training process of our model is easier with least computation time.
Results:
3.3 REALISTIC DATA DAILY TRADING VOLUME
In this section we would like to demonstrate the performance of our model in financial area. Trading volume the total quantity of shares or contracts traded for specified securities such as stocks, bonds, options contracts, future contracts and all types of commodities. It can be measured any type of security traded during trading day or specified time period. In our case, daily volume of trade measured on stocks. Predicting traded volume is an essential component in financial research since the traded volume, as a basic component or input of other financial algorithms, tells investors the market's activity and liquidity. The data set The data covers period from January 2018 to January 2020 and is obtained from Bloomberg. Each day from 14:30 to 20:55, we have 1 observation every 5 minutes, totally 78 observations everyday. Our task is described as follows: given first two years data, we use the traded volume at 14:30, 14:40, 15:05, 15:20 and 16:20 training data, namely, 25, 22, 307, 200. 2222 to train model, then for next 100 days we predict traded volume at 14:35, 15:15, 15:35 and 16:15, namely, I 29, 313, 221 One of baselines we choose is classical rolling means(RM) method, which predicts intraday volume of particular time interval by the average volume traded in the same interval over the past days. The other one baseline is a kalman filter based model (Chen et al. 2016) that outperforms all available models in predicting intrady trading volume.
is the historical traded volume of the stock "JPM".
Experiment Setup: Following similar setup as we did for RNA data set, we utilize the same structures for neural networks here. For hyperparameters we set Ar = 0.02, T/Ar = 22, which means takes one single from x; to 21+1- For preprocessing. we rescale data by taking natural logarithm of trading volume, which is common way in trading volume research. We conduct experiments on two groups to show advantages of our method, for first group we train our model on complete data set, in this case the data has full trajectory; for second group manually delete some trajectories of the data, for instance, we randomly kick out some samples of Zy, I2, Ly, Lo0, 22 then follow the same procedures of training and prediction.
As shown in first four figures, with full trajectory, prediction made by RM is almost straight line, the prediction value bouncing up and down within
Results: We present prediction results in Figure
7

---

Under review as a conference paper at ICLR 2021
(a) Comparison on 14:35
(b) Comparison on 15:15
(c) Comparison on 15:35
(d) Comparison on 16:15
(c) Comparison on 14:35 (f) Comparison on 15:15
(g) Comparison on 15:35 (h) Comparison on 16:15
Figure 4: (a) to (d): Group A: with full trajectory of training data, predictions of traded volume in next 100 days, RM(yellow) fails to capture the regularities of traded volume in time series, kalman filter based model(green) fails to capture noise information and make reasonable predictions, our model(blue) is able to seize the movements of traded volume and yield better predictions. (e) to (h): Group B: predictions of our model without full trajectory.
very small range, thus this model cannot capture the volume movements, namely, regularities existing in the time series; prediction made by the Kalman filter based model captures the regularities better than RM model, but fails to deal with noise component existing in the time series, thus some predictions are out of a reasonable range. Traded volume predicted by our model closer to the real case, moreover, our model captures regularities meanwhile gives stable predictions. Furthermore, without full trajectory, Kalman filter based model fails to be applied here and RM model still fails to capture the regularities, we randomly drop half of the training samples and display predictions made by our model in last four figures of Figure 4. we see our model still works well.
4 DISCUSSIONS
In this section we discuss the limitations and extension of our model.
An essential challenge for recovering the drift term: Mathematically it is impossible recover the exact drift term an SDE if we are only given the information of density evolution on certain time intervals, because there might be infinitely many drift functions induce the same density evolution. More precisely, suppose p(x, f) solves FPE (equation D. consider the following equation
0 =
ax; (x,xx + 2
Σ
D
0x ((x
82
(8)
One can prove, under mild assumptions. that there are infinitely many vector fields u(x, t) (u (x, ()..,,... ( f)) solving equation (8). One can check that the solution to FPE (1) with drift term 8(x,1) u(x, f) is still p(x, 1), i.e. the vector field u(x, 1) will never affect the density evolution of the dynamic. This illustrates that given the density evolution p(,,), the solution for drift term is not necessarily unique. This is clearly an essential difficulty of determining the exact drift term from the density. In this study, the main goal is to recover the entire density evolution (i.e. interpolate the density between observation time points) and predict how will the density evolve in the future. As result, although we cannot always acquire the exact drift term of the dynamic, can still accurately recover and predict the density evolution. This is still meaningful and may find its application in various scientific domains.
8

---

Under review as a conference paper at ICLR 2021
(a) Training data
(b) Prediction at
(c) Prediction at
(d) Vector field
Figure 5: Results of learning curl field:(a): the whole distribution is indicated by green, from which choose subsets for training purpose, training samples lo- 42 and indicated by pink, orange, red and brown respectively. (b) and (c):the prediction at f and f3 are indicated by blue points. (d): The vector field generated by our model.
(a) Prediction at
(b) Prediction at f30 (c) Prediction
Figure 6: Results of learning diffusion coefficient.
Curl Field: The drift function we showed in the numerical experiments will apparently cause the evolution of the distribution. If the drift function a curl, namely V x F, then the distribution does not change, it still possible to learn the density evolution? The answer is yes. To demonstrate this point of view, we simulate a curl field (y, xx) induced by A [(0, 1), (-1, 0)] on Gaussian distribution (indicated by green points in Figure 5 (a)). here we set noise part as Realizing that our learning is based samples but not on the density, suppose we observe samples at four time points which does not perfectly follow the distribution (indicated by four colors: pink(to), orange(11), red((2) and brown((3) in Figure 5 (a)), we learn and predict the distribution at fi and (indicated by blue points Figure 5 (b) and (c)). We also show the vector field learned by our model Figure (d). We that the predictions and vector field all satisfy the correct curl field pattern.
Our framework also works for learning unknown diffusion coeffi- cient in the Itô process. As extension of our work, if we approximate the diffusion coefficient with a neural network or (with parameters 17), we revise the operator of as:
Learning Diffusion Coefficient:
D
D
O² dx? (aae))
F,(0
&es (A)) Ox
((
+
which can derived by the same technique we used to derive Proposition We test this formulation on synthetic data set, where we only consider diffusion influence, namely, drift term in Equation (8) is ignored. We set the ground truth diffusion coefficient as o = [(1, 0). (0, 2)J. We design the neural network simple fully connected layer with 32 nodes, then show result in Figure we see that the predictions(blue) follow the same pattern as the ground truth(red) does.
CONCLUSION
In this paper, we formulate a novel method to recover the hidden dynamics from aggregate data. In particular, our work shows one can simulate the evolving process of aggregate data as an Itô process, in order to investigate aggregate data, we derive a new model that employs the weak form FPE well the framework of WGAN, Furthermore, in Appendix prove the theoretical guarantees of the error bound our model. Finally we demonstrate our model through experiments on three synthetic data sets and real-world data sets.

---

Under review as a conference paper at ICLR 2021
REFERENCES
D. Alshamaa, A. Chkeir, F. Mourad-Chehade, and P. Honeine. Hidden markov model for indoor trajectory tracking of elderly people. In IEEE Sensors Applications Symposium (SAS), 2019.
Martin Arjovsky, Soumith Chintala, and Léon Bottou. arXiv:1701.07875. 2017.
Wasserstein gan. In arXiv preprini
Kendall E Atkinson. An introduction to numerical analysis. John wiley & sons, 2008
C. Beck, S. Becker, P. Grohs, N. Jaafari, and A. Jentzen. Solving stochastic differential equations and kolmogorov equations by means of deep learning. 2018.
R. Chen, Y. Feng, and D. Palomar. Forecasting intraday trading volume: kalman filter approach. In Available at SSRN 3101695, 2016.
M. Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In Advances in Neural Information Processing Systems, 2013.
M. Deriche, A. Absa, Amin, and Liu. A novel approach for salt dome detection and tracking using hybrid hidden markov model with active contour model. Journal of Electrical Systems, 16(3):276-29,, 2020.
M. Djurie, J. H. Kotecha, J. Zhang, Y. Huang, Ghirmai, M. F. Bugallo, and J. Miguez. Particle filtering. IEEE Signal Processing Magazine, 20(5):19-38, 2003.
S. R. Eddy. Hidden markov models. Current Opinion in Structural Biology, 6(3):361-363 1996.
Y. Fang, C. Wang, W. Yao, X. Zhao, H. Zhao, and H. Zha. On-road vehicle tracking using part-based particle filter. IEEE Transactions on Intelligent Transportation Systems, 20 20(12):4538 4552, 2019.
Farahi and H. S. Yazdi. Probabilistic kalman filter for moving object tracking. Signal Processing: Image Communication, 82, 2020.
X. Glorot and Y. Bengio. Understanding the difficulty of training deep feedforward neural networks. In International Conference on Artifccill Intelligence and Statistics, 2010.
A. C. Harvey. Forecasting, structural time series models and the Kalman Cambridge University Press, 1990.
Hashimoto, D. Gifford, and Jaakkola. Learning population-level diffusions with generative mns. In International Conference on Machine Learning, pp. 2417-2426, 2016.
A. Hefny, C. Downey. and G. J. Gordon. Supervised learning for dynamical system learning. In Neural Information Processing Systems, 2015.
S. C. Hicks, M. Teng, and Irizarry. the widespread and critical impact of systematic bias and batch effects in single-cell ma-seq data. bioRxiv, 2015.
R. E. Kalman. A new approach to linear filtering and prediction problems. arXiv preprint arXiv:1805.04009 04099, 1960.
Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization, 2014.
A.M. Klein, L. Mazutis, Akartuna, Tallapragada, A. Veres, Li, L. Peshkin, D.A. Weitz, and M.W. Kirschner. Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells. Cell, 187-5)11177-1201, 2015.
Peter E Kloeden and Eckhard Platen. Numerical solution of stochastic differential equations, vol- ume 23. Springer Science & Business Media, 2013.
J. Langford, R. Salakhutdinov, and T. Zhang. Learning nonlinear dynamic models. In International Conference on Machine Learning, pp. 593-600, 2009.
W.C. Li, Liu, H. Zha, and H. Zhou. Parametric fokker tokker-planck equation. Geometry science of information, 2019.
10

---

Under review as a conference paper at ICLR 2021
Y. Li. A data-driven method for the steady state of randomly perturbed dynamics. arXiv preprint 0rXiv:180S.04099, 2018.
Q. Liu and D. Wang. Stein variational gradient descent: A general purpose bayesian inference algorithm. In Neural Information Processing Systems, pp. 2378-2386, 2016.
G. N. Milstein and M. V. Tretyakov (eds.). Stochastic numerics for mathematical physics. Springer Science & Business Media, 2013.
T. Miyato, T. Kataoka, M. Koyama, and Y. Yoshida. Spectral normalization for generative adversarial networks. arXiv preprint arXiv:1802.05957, 2018.
E. Nelson. Quantum fluctuations. Princeton University Press, 1985.
Bernt Øksendal. Stochastic differential equations. In Stochastic differential equations, pp. 65-84. Springer, 2003.
M. Pavon, E. G. Tabak, and G. Trigila. The data-driven schroedinger bridge. arXiv preprint arXiv:1806.0133 01364, 2018.
D. Qi and A.J. Majda. Low-dimensional reduced-order models for statistical response and uncertainty quantification: Two-layer baroclinic turbulence. Journal of the Atmospheric Sciences, 73(12): 4609 4639, 2016.
D.J. Rezende and S. Mohamed. Variational inference with normalizing flows. In arXiv preprint arXiv:1505.05770, 2015.
H. Risken. The fokker-planck equation. Springer Series in Synergetics, 18:4609-4639, 1989.
H. Risken and T.K. Caugheyz (eds.). The fokker-planck equation: Methods of solution and application. Springer, 1991.
N. P. Santos, V. Lobo, and A. Bernardino. Unmanned aerial vehicle tracking using a particle filter based approach. In JEEE Underwater Technology (UT), 2019.
Justin Sirignano and Konstantinos Spiliopoulos. Dgm: deep learning algorithm for solving partial differential equations. Journal of computational physics, 375: 1111-33993444 2018.
C. Villani. Optimal transport: old and new, volume 338. Springer Science & Business Media, 2008.
Y. Wang, B. Dai, L. Kong, S. M. Erfani, J. Bailey, and H. Zha. Learning deep hidden nonlinear dynamics from aggregate data. In Uncertainty in Artificial Intelligence, 2018.
E. Weinan, J. Han, and A. Jentzen. Deep learning based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations. In Commu- nications in Mathematics and Statistics, pp. 349-380, 2017.
Y. Zang, G. Bao, X. Ye, and H. Zhou. Weak adversarial networks for high-dimensional partial differential equations. In arXiv preprint arXiv:1907.08272 2019.
O.C. Zienkiewicz and I.K. Cheung. The Finite Element Method in Engineering Science. McGraw-
Hill European Publishing Programme. McGraw-Hill, 1971. ISBN 9780070941380. URL https
11

---

Under review as a conference paper at ICLR 2021
A
APPENDIX
Table 1: The Wasserstein error of different models on Synthetic-1/2/3 and -equeneee data sets.
Data
Task
Dimension
NN
LEGEND
Ours
2
1.37
0.44
0.05
6
4.79
2.32
0.06
9.13
2.89
0.10
Syn-1
2
0.84
0.18
0.03
a7500
3.28
0.30
0.03
10
8.05
1.79
0.09
2
5.53
1.27
0.05
44
6
9.97
4.20
0.08
10
16.51
6.08
0.12
Syn-2
2
1.68
0.53
0.05
6
4.84
1.08
0.07
9.23
2.85
0.11
2
4.13
1.29
0.08
0030
6,40
3.16
0.17
10
11.76
8.53
0.25
Syn-3
3.05
0.87
0.12
2550
6
6.72
1.52
0.16
10
9.81
3.55
0.23
D)2
10
33.86
10.28
4.23
RNA-Mt1
D7
10
12.69
7.21
2.92
D2
31.45
13.32
4.04
RNA-Mt2
D7
10
11.58
7.89
1.50
B ERROR ANALYSIS
In this section, we provide error analysis of our model. Suppose the hidden dynamics is driven by (x)),, the dynamics that we learn from data (x),,, then original Itô process, Euler processes computed by true g, and estimated g are:
dX 8 d( odW
/ x. S, + VAIN(O. 1)
// x/ + g)
) (x/))r++ V
where X is the ground truth, x' is computed by true 8, and æf is computed by estimated gy- Estimating the error between original Itô process and its Euler form can be very complex, hence we cite the conclusion from (Milstein & Tretyakov 2013) and focus more on the enror between original form and our model.
Lemma 2. With the same initial X 32 = Ca. if there is a global Lipschitz constani K which satisfies:
Ig(a, 1) g(y.t)s K| yl
then after steps, the expectation error between Itô process x,, and Euler forward process xf, is:
B|,, - 2 K(1 + 2 Ar
Lemma illustrates that the expectation error between original Itō process and its Euler form is not related to total steps but time step As.
12

---

Under review as a conference paper at ICLR 2021
Proposition 3. With the same initial zo, suppose the generalization error of neural network g is and existence of global Lipschitz constant
Ig(x)) g(y)ss KlØ yl
then after steps with step size A = T/n, the expectation error between Itô process a, and approximated forward process sd, is bounded by:
Ba - -2/15
8 (c - 1) K(1 Elzoo
(10)
K
Proposition 3 implies that besides time step size At, our expectation emror interacts with three factors, generalization error, Lipschitz constant of g and total time length. In our experiments, we find the best way to decrease the expectation error is reducing the value of K and n.
C
PROOFS
C.I PROOF OF PROPOSITION
Proof. Suppose 2/4) and 2 are our observed samples at In and fm-1 respectively, then expectations could be approximated by:
B Bo-Maas)(/))1 = f ((x))(a, mn)a ss
N ((
(11)
ʃ f(x)p(e, mm)x S f(x) pla, Im-1++
( Nn @p(ar, T)
Ef-ptusalffx))
dr dx
Öt
f f())p( bm 1)x S f(x))
Inn 0p(x, drdæ
Ōt
Nn f( ) ʃ [(x)
202
Σ
D Φ²
p(x, r)
didæ
Ox,
lam1
(12)
Then for the second term / above, it difficult to calculate directly, but we can use integration by parts rewrite
f(x) Ø
gl (x)p(t, +) o0 Σ
D
f(x)
8²
p(x, r)
dxdr
dx
dxp
D
Φ
& (x)
7)
f(x)
0² Σ
p(a,
T)
f(x)
dadr
dxi
dx²
Ea-pian
Σ
D
& (x)
0
f(x)
+ Eapp(())
D
0x f(x)
82
0x
N
D
(x
δ f(1
00²
2
Σ suati
D
Ø2
dr
(13)
Ox
To approximate the integral from fn-1 to Im- we adopt trapezoid rule, then we could rewrite the expectation in Equation 112
13

---

Under review as a conference paper at ICLR 2021
Bx-Mar [f()]
[f()]]
1
f(
Ar
(
+ 702
1
0
(k)
N
2
0x,
axp
k=
K=1
i=l
8²
8
(2
f(2%
0²
2
f(x
dx
0x,
1=
1=1
1
f(
At
[T
1))) (X
(14)
N
2
We subtract by 14 to finish the proof.
C.2 PROOF OF PROPOSITION 2
Proof. Given initial &. generate Tn T2y By C; sequentially by Equation 6 expectations can be rewritten as:
Then the
f (()) ta)d
B [[aa]]
1& (k))
(15)
N
N
D
D
f(
Σ
8
(1
f(1l))
50²
faalit dr
0x
dx,āx,
N
D
02
0
((
dr
N
Ox:
0x
-=1
N
1
+
(x
1(x)0)
+
=
0
((x10)
dr
(16)
Oxi
2
0xp
1=1
which is:
Eg Bl-
f(2
(
Al
Ar
2
2-1
A/ + 2
+ F811
(17)
Finally it comes
N
B Be- x
Σ r
At 2 F1180)+ + E
(18)
sml
We subtract
15
by
18 to finish the proof.
C.3 PROOF OF ERROR ANALYSIS
Proof. The proof process of Lemma 2 is quite long and out of the scope of this paper, for more details please see first two chapters in reference book (Milstein & Tretyakov 2013). While for the proof of
Proposition 3 with initial X and first one-step iteration:
x E
(19)
[ey w/ + g. ) VAIN(O.1)
(20)
x x + g )r VAIN(O,)
14

---

Under review as a conference paper at ICLR 2021
Then we have:
Baaa 2 = B/22, I/o1= 0
(21)
Baa 2,, = Bla", o/ - + VALNCO. 1) or VAINOO DI
I Blar,, (B (c4
= Blg.(ar))) ( + gr(og ) - 8, s,(2r IA
(ag 8y )(l + / - 8, (1
s &AN + Blg 8f 8 )Ai
&Ar + Blg (a Mayy - x |A
(o [ 2/111
&AI + KB 2 AAi
&Ai
(22)
Follow the pattern we have:
xy/y = Fy +g ) + VAN(OO 1)
(23)
1 I t A
1y1 = ب + tg ) + VA)
(24)
xy = Jhem + t A + VANN(O.)
Which leads to:
Bla% 2/ = Bla,, 2 gx Ar ( V I)- or VANN(O. 1)
I Bla, ,11+ ) (1) I/M
< Ble - x/I+ + KBJE,, - IAr
$ (1 + KAT)aA + EA1
(25)
Baa x S cAr
(1 KAD)
(26)
1=0
Now let En1 (1 + KAr), then consider followings:
S(KM) S(1 KAí) S
m-1
Σ (1+KA) KAr)
(1 + KAn)
7=0
= =( KAr)" -1
= =1+K->) kーx- -1
KT_ -1
(27)
Finally we have:
Bkc 2015 8 < 1)
(28)
K
kT_ B|C 21115 x [
1) K(1 + E|zo! 1/2
(29)
15


Reviews:
[{"review_id": "kPOz7UXZLR", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 8, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposes a new approach to learn the dynamics of density evolution of objects from aggregated data. The basic idea is to derive a closed-form Wasserstein distance between the empirical data distribution and the predicted distribution generated from the weak form of Fokker Planck Equation (FPE). Based on this measure, an objective function is developed to learn the underlying drift coefficients and the discriminator simultaneously using neural networks. Some theoretical results are provided and numerical experiments are carried out to illustrate the effectiveness of the proposed method.\n\nThe paper is overall well written and solves an important problem. I have only a few minor questions.\n\n1. How prediction is made based on the estimated model? One step ahead or multiple-step ahead prediction? Please provide more details on this issue.\n\n2. Some notations are unclear in Theorems 1-1. In Theorem 1, what are the definitions of $N$ and $x^{(k)}$? Do you generate data for N times? In Theorem 2, I  assume that $n$ stands for $n$ steps ahead from t_{m_0}?", "extra_scores": {}, "extra_text": {}}, {"review_id": "bVeh7yW-K0", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 4, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper is about an interesting setting, where observations correspond to lots of individuals evolving over time. The twist is that individual identifiers are not available.  Thus from one time point to the next, the observer does not know which individual is which. The problem is to learn a time series model from this, and models considered here are in the form of an SDE.\n\nIn my first reading of the paper, I thought the topic was very interesting but I struggled to follow some parts of the manuscript, which I thought was overly unclear. Some of the notation is undefined or ill-defined, the general reasoning is hard to follow, the notation does not help identifying the objects that are observed from the objects that are simulated, there is a lack of acknowledgement of limitations of the proposed method, the synthetic examples are not well-motivated, and various other fairly standard concerns.\n\nHowever, during my second reading I looked at the reference Wang et al, \"Learning Deep Hidden Nonlinear Dynamics from Aggregate Data\", published in UAI in 2018.  From this I got more serious concerns.\n\nThe authors' abstract starts with \"Learning nonlinear dynamics from aggregate data is a challenging problem since the full trajectory of each individual is not available, namely, the individual observed at one time point may not be observed at next time point, or the identity of individual is unavailable due to technical limitations, experimental cost and/or privacy issues.\" \n\nWang et al's abstract mentions: \"However, in most of the practical applications, these requirements are unrealistic: the evolving dynamics may be too complex to be modeled directly on observations, and individual-level trajectories may not be available due to technical limitations, experimental costs and/or privacy issues.\"\n\nThe authors write \"In the work of Hashimoto et al. (2016), a stochastic differential equation(SDE) was adopted to capture the dynamics of particles directly on observations, in particular the drift coefficient is parameterized by a recurrent neural network.\"\n\nWang et al. write: \"Modeling the dynamics on aggregate observations have been investigated recently in (Hashimoto et al., 2016), where a stochastic differential equation (SDE) has been used to capture the transition directly on observations Yt.\"\n\nThe authors write: \"There are many existing models to learn dynamics of full-trajectory data. Popular ones include Hidden Markov Model (HMM)(Alshamaa et al., 2019; Eddy, 1996), Kalman Filter (KF)(Farahi & Yazdi, 2020; Harvey, 1990; Kalman, 1960) and Particle Filter (PF) (Santos et al., 2019; Djuric et al., 2003). These models and their variants (Deriche et al., 2020; Fang et al., 2019; Hefny et al., 2015; Langford et al., 2009) require full trajectories of each individual, which may not be directly applicable to the aggregate data as we mentioned earlier. \"\n\nWang et al. write: \"Existing models such as Hidden Markov Model (HMM) (Eddy, 1996), Kalman Filter (KF) (Harvey, 1990) and Particle Filter (PF) (Djuric et al., 2003) are popular methods with hidden variables. However, these models and their variants (Langford et al., 2009; Hefny et al., 2015) require individual-level trajectories, which may not be available, as was mentioned earlier.\"\n\n--\n\nBased on these strong similarities I do not think that the proposed manuscript is suitable for publication.", "extra_scores": {}, "extra_text": {}}, {"review_id": "5XEz9k0gOr0", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 5, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper addresses an interesting problem: learning stochastic dynamics from aggregate data. The aggregate data refers to the setting when the data is anonymized. For example, one has data about the location of the birds at different instants in time, but the birds are not labeled and can not be distinguished from each other. For the purpose of learning, the paper proposes to use a generator, in the setting of WGAN, to learn the drift term in the dynamics. The novelty of the proposed approach is to consider the weak form of the Fokker-Planck equation to express the loss function. \n\nThe problem is interesting, but I found several issues about the proposed approach that need to be addressed.  \n\n1- The paper consider a special setting where the dynamics is given by a diffusion process. It is assumed that the diffusion coefficient is known and the objective is to learn the drift function. I think this is a strong assumption and learning diffusion coefficient is also very important and need to be included in this work.  \n\n2- Learning the drift function from aggregate data may have fundamental limitation that the paper needs to address. It seems impossible to learn mixing terms (the terms in the dynamics that do not effect the distribution) from aggregate data. In other words, even in and ideal setting with infinite capacity for neural network and infinite amount of data, the learned drift function differs from the actual one. This can be verified by running a simple experiment with dynamics [[0,1][-1,0]] so that the distribution does not change and the learned drift function should be zero. \n\n3- In the proposed formulation, eq (8), it is not clear which one is real data and which one is generated data. \n\n4- Theorem 1 is not stated correctly and it is missleading. It states that eq (7) is the Wasserstein distance. However, eq (7) is just an approximation. And I really don't think this should be stated as a result in a Theorem. This is basic definition of Wasserstein distance and eq (1) and it is more appropriate to be presented as part of text.  Also, Lemma 1 is well-known and a reference should be given. \n\n5- The numerical experiments are nice, but they should also consider the long-term behavior: If the learned dynamics results to the correct stationary distribution or not. Because the stationary distribution is one of the important characteristics of a stochastic dynamics.  ", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: While the reviewer's noted a number of strengths of your paper, the approach that you took, and agreed that you had tackled an important problem, concerns remained about presentation and clarity. I agree. (Here are just a few miscellaneous comments: the very first paragraph of the Introduction needs to be rewritten for clarity, in my opinion. Later on page 1, you use the term "the dynamics of density" but you should not assume that the reader knows what that means. There are typos as well, e.g. "make all predictions base[d] on Equation (6)" on page 4. It would be helpful to know something about why you chose the experimental setups in Synthetic-1, 2, and 3. )

Regarding the similarities between this paper and a previously published article I believe that the authors have addressed these concerns; I hope they are careful to avoid this situation in the future.
Decision: Reject