{'fraction_non_alphanumeric': 0.040238055877004464, 'fraction_numerical': 0.02159034551165482, 'mean_word_length': 4.126440677966102, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solving Inverse Stochastic Problems from Discrete Particle Observations Using the Fokker-Planck Equation and Physics-informed Neural Networks

The Fokker-Planck (FP) equation governing the evolution of the probability density function (PDF) is applicable to many disciplines but it requires specification of the coefficients for each case, which can be functions of space-time and not just constants, hence requiring the development of a data-driven modeling approach. When the data available is directly on the PDF, then there exist methods for inverse problems that can be employed to infer the coefficients and thus determine the FP equation and subsequently obtain its solution. Herein, we address a more realistic scenario, where only sparse data are given on the particles' positions at a few time instants, which are not sufficient to accurately construct directly the PDF even at those times from existing methods, e.g., kernel estimation algorithms. To this end, we develop a general framework based on physics-informed neural networks (PINNs) that introduces a new loss function using the Kullback-Leibler divergence to connect the stochastic samples with the FP equation, to simultaneously learn the equation and infer the multi-dimensional PDF at all times. In particular, we consider two types of inverse problems, type I where the FP equation is known but the initial PDF is unknown, and type II in which, in addition to unknown initial PDF, the drift and diffusion terms are also unknown. In both cases, we investigate problems with either Brownian or Levy noise or a combination of both. We demonstrate the new PINN framework in detail in the one-dimensional case (1D) but we also provide results for up to 5D demonstrating that we can infer both the FP equation and} dynamics simultaneously at all times with high accuracy using only very few discrete observations of the particles.


Introduction
Stochastic dynamical systems, defined as deterministic systems of differential equations perturbed by random disturbances that are not necessarily small, have been used effectively across disciplines, including gene regulatory networks, physics, chemistry, population dynamics and medicine. The general form of equations may be similar, however the specific noise type and the constant or variable coefficients in these equations are case-dependent. Hence, it is important to use data to infer these coefficients and infer the stochastic dynamical system and forecast its dynamics.
There exist a few methods to determine stochastic differential equations from observations. Ruttor et al [15] used a Gaussian process prior over the drift as a function of the state vector, and developed an approximate expectation maximization algorithm to deal with the unobserved latent dynamics between observations. Archambeau et al [2] presented a variational approach for the approximate inference of stochastic differential equations (SDEs) from a finite set of noisy observations. Subsequently, Mapper [11] used the variational inference approach to approximate the distribution over the unknown path of the SDE conditioned on the observations and provided approximations for the intractable likelihood of the drift. Also, in [3], the authors used the sparse learning method to learn the stochastic dynamical systems. In [17], the authors used small samples from just a few snapshots of unpaired data to infer the drift and diffusion terms in SDEs with Brownian and Lévy noise in high-dimensions. Here, we follow the recent work of [17] but we replace the SDE with the Fokker-Planck (FP) equation, which governs the dynamics of the probability density function (PDF).
For the forward problem, there are several methods to solve the FP directly [6,8,18,5] assuming that the initial condition is known. Recently, Xu et al [16] used PINNs to solve the one-dimensional stationary FP equation (for Brownian motion) with homogeneous Dirichlet boundary conditions, where the PINN loss was augmented with an additional loss for the constraint of numerical integration. More broadly, there are several papers by now using deep learning to solve general partial differential equations [14,13,9,12,4]. In these works, to solve the forward problem, the initial condition and the boundary condition is needed, while for the inverse problem, observation of the solution in some points in space-time is required.
In the current paper, we assume that we only have the observation data of the sample paths at just a few discrete times but we do not know the initial PDF. We consider problems with both Brownian and Lévy noise as well as a combination of the two. Our aim is to combine the observation data and the corresponding parametrized FP equation with unknown coefficients to infer its precise form and hence obtain the PDF at all times, including the initial condition. In particular, we employ the variational form of the Kullback-Leibler divergence [10] to connect and combine the observations of particle positions and the corresponding FP equation. To this end, we employ a physics-informed neural network (PINN) as shown in Figure 1 that enables a seamless integration of the data and the FP equation using automatic differentiation, hence avoiding the need for grid generation, which is computationally prohibitive in high-dimensions. Figure 1: Schematic of the PINNs for solving inverse stochastic problems from discrete particle observations. The total loss consists of a loss from residual of the Fokker-Planck equation and another loss from particle observations. The paper is organized as follows. In Section 2, we introduce the methodology to solve the FP equation and the connection between the trajectories and the PDF. In Section 3, we present the results of PINNs for solving the two types of inverse problems of the FP equation for up to five dimensions (5D). Finally, we conclude with a short summary. We have included more detailed results in the Appendices (A-F), studying both Brownian noise and Lévy noise cases.

Problem setup
Consider the following stochastic differential equation: where a(·) is the vector drift term, σ and ε are n × n matrices, B t is a Brownian motion in R n , and L α t is a symmetric α-stable Lévy process in R n . The generating triplet of the Lévy process is (0, 0, ν α ). The processes B t and L α t are taken to be independent.
The generator of the SDE 1 is [1,7]: where ν α (dy) is the α-stable Lévy measure and ν α (dy) = C n,α ||y|| −n−α dy, C n,α = αΓ((n+α)/2) 2 1−α π n/2 Γ(1−α/2) . The FP equation describing the time evolving PDF of the stochastic process X t , starting from p 0 at T 0 , is written as follows: where A * is the adjoint operator of the generator A and has the following form: For simplicity, we only consider the case where the α-stable Lévy noises are independent across dimensions, which means the Lévy measure has the formula [5]: where δ 0 is the Dirac delta measure concentrated at 0. Consequently, We consider the scenario where the data available are the discrete observations of the particles at different time instants ti are i.i.d. random samples of X ti drawn from p(·, t i ). We refer to D i as a snapshot for each i. Our objective is to infer the possibly unknown terms in the FP equation as well as the PDF at all times by utilizing the equation and the observation data in a unified learning system. More specifically, we consider the following two types of problems.
1. Problem I: We have full knowledge of the operator A * but we do not know the initial condition p 0 . Instead, we have a few snapshots consisting of limited samples, and wish to learn the PDF as a function of time, including the initial condition.
2. Problem II: While we know the general form of A * , we do not have full knowledge of some terms in A * , e.g., the drift term a, or the entries in the coefficient matrices σ and ε. We wish to use the snapshots to infer the unknown terms in the FP equation.

PDF estimation using neural networks
Given two distributions P and Q on R n with PDF p(x) and q(x), respectively, with p, q > 0 for simplicity, the Kullback-Leibler (KL) divergence from Q to P is formulated as: From [10] we have the following variational form of the KL divergence: where the supremum is attained at In other words, Note that Eqn. 8 is a key to our formulation and provides an effective means of estimating the PDF from particle observations. We denote the estimator of p(x) asp(x). Suppose that we are given a finite number of i.i.d. samples from the unknown distribution P on R n with PDF p(x), for which the empirical distribution is denoted asP . We define Q as a known distribution with PDF q(x) > 0. Then While other distributions, e.g., the Gaussian distribution also applies in principle, for simplicity, we set Q as an uniform distribution on D, i.e., q(x) = 1/µ(D) for x ∈ D, where µ denotes the Lebesgue measure and D is a hypercube that contains the support ofP . We also set the estimatorp(x) = 0 for x outside D. Then, from 9 we havê where the integral over D can be calculated via Monte Carlo sampling or quadrature. G is a specific class of positive real value functions. Intuitively, we wish G to be rich enough so that it can contain p(x)/q(x), while also exclude "morbid" functions to avoid overfitting. A natural idea is to use a family of neural networks as G due to the huge capacity of neural networks as well as our empirical knowledge that the post-training neural networks tend to be "regularized".

Solving the FP equation using discrete particle observations
We use a deep neural networkp(x, t) with the concatenation of x and t as input and a real number as output to approximate p(x, t), the unknown solution of the FP equation. Different from vanilla neural networks, we add an additional softplus activation function before the output layer to makep(x, t) strictly positive. Based on the FP equation and observations of discrete particle observations, we have two requirements for the post-trainingp(x, t), as follows: where {(x i , t i )} Nr i=1 are N r randomly sampled so-called "residual points" on the spatio-temporal domain D × [T 0 , T 1 ]. The differentiation applied to the neural networkp(x, t) could be performed via automatic differentiation, while the integral operator in A * can be estimated by a quadrature as in the fractional PINNs [13].
More specifically, in 1D problem, where we truncate the solution on the domain (−L, L), and use the following discretization scheme [8] to approximate the non-local part at point (x i , t j ): where C h =−C α ζ(α−1)h 2−α and means that the quantities corresponding to the two end summation indices are multiplied by 1 2 . In problems with dimension larger than one, the scheme of non-local part in 4 in each dimension is the same as in the 1D case.
For the second requirement, recall thatp(x, t i ) as an estimator of p(x, t i ) should satisfỹ for i = 1, 2...N , where represents equality up to a constant. While the first term in the right-hand side is averaged over all the samples in ti } ni j=1 , in practice we can take the average over a minibatch of size b smaller than n i . The integral over D in the second term can be estimated by quadrature or Monte Carlo sampling.
The total loss is then a combination of PDE loss and data loss: where τ is the weight to scale and balance the losses from the FP equation and data.
We note that when we have observations of X t , we can also use the kernel density estimation to obtain the PDF, denote as as p KD . Correspondingly, we can replace the data loss in 15 with another loss from the estimated density, as follows: where and

Computational Results
All the neural networks in our tests have 4 hidden layers and 20 neurons per layer, with tanh activation function. We use the Adam optimizer with a learning rate of 10 −4 . We have two methods to compute the inverse problem II. In the first one, we use a polynomial to approximate the drift term, while in the second one, we use directly one or more neural networks to approximate the drift term, which could be any general function. The training takes 100, 000 iterations for problem I, while 200, 000 iterations for problem II, except 3D, 4D and 5D cases where the training iterations are shown in the corresponding figures.

One-dimensional Case
Consider the following 1D SDE: the corresponding FP equation is:

Problem I: Brownian noise
In this section, we use the deep neural network to compute problem I for the FP equation driven by Brownian noise. We consider the drift term a(x) = x − x 3 , σ = 1 and ε = 0. We choose the observation data of X t at several times, and use the the Monte-Carlo method with 5,000,000 samples to obtain the reference PDF for estimating the accuracy of our method.
For the inverse problem I, we assume that the observation data of samples of X t are available only at one snapshot. To assess the accuracy of our predictions, we consider 100, 1000 and 10000 samples at t = 0.3, i.e., N = 1. The minibatch size is set as the training data size. When we compute (14), we use the compound trapezoid formula with 301 points to approximate the integral part. The number of residual points for evaluating the FP residual is N r = 1000. The predicted results are shown in Figure 2 at times t = 0.1, 0.3, 0.7, 1. We can see that the predicted PDFs are better if we have more observation data, as expected; however, even with only 100 samples, the two peaks of the PDF can be captured accurately. The results here shows that physics-informed learning can capture the dynamics accurately at all times given data of one snapshot only. Next, we also compute the inverse problem I of Brownian noise with more observation data at t = 0.1, 0.3, 0.5. For each time, we have 100 independent samples. We use the kernel density estimation to approximate the density using the observation data (p KD ). When we obtain the p KD at t = 0.1, 0.3, 0.5, we use this approximate density as observation data for the PDF and use the neural network to train the model to obtain the density at other times; we call p N N KD the output of the PINN, whose loss function is defined in (16). We use N p = 301 equidistant points as the collocation points in (16). We use the compound trapezoid formula with 301 point to approximate the integral part of (14) and 1000 residual points for the FP residual in (12). We compare these results with the PINN results obtained directly from the observation data in (15). For each case we test different values of weights τ in (15), and for each τ we run the code for 3 times with different samples, initialization and random seeds. The errors are shown in the top line of the Figure 3. For the second and third row, we choose one of the three cases with τ = 100. The comparison is shown in Figure 3 in terms of the L 2 error of the two methods. We can see that the results obtained with the PINN using directly the particle observations are better than the PINN results that uses the estimated PDF as input data. Also, our results are superior to p KD from the kernel density estimation for individual time instants, showing that physics-informed learning is very important in capturing the dynamics and integrating the data at multiple time instants.

Problem I: Lévy noise
Similar to the results above, in Appendix A we present a systematic study for problem I assuming that the SDE is driven by Lévy noise.

Problem II: Brownian and Lévy noise
We now consider problem II, where we want to infer in addition to the initial PDF, the drift and diffusion terms from snapshots. In Appendix B, we present a systematic study for the case of Brownian noise, and in Appendix C we present results for the case of the Lévy noise. Instead, here in this section, we consider the SDE driven by both Brownian noise and Lévy noise. Specifically, we consider the drift term a(x) = x − x 3 , α = 1.5, σ = 1 and ε = 1. We use the compound trapezoid formula with 201 points to approximate the integral part of (14) and 800 residual points to compute the FP residual in (12).
We compute the following cases assuming different snapshots of available data and different methods of representing the drift term, which is an unknown function. For case (a1), the observation data is available at t = 0.2, 0.5, 1 and we use a polynomial to approximate the drift term. For case (a2), the observation data is available at t = 0.2, 0.5, 1 and we use a neural network to approximate the drift term. For case (b1), the observation data is available at t = 0.1, 0.4, 0.7, 1 and we use a polynomial to approximate the drift term. For case (b2), the observation data is available at t = 0.1, 0.4, 0.7, 1 and we use a neural network to approximate the drift term. For each case, we assume that we have 10,000 samples at each time instant. We set the minibatch size b = 1000. The results for the approximation of the drift term and the evolution of the inferred PDF for all four cases are shown in Figure 4. The corresponding results for the parameter estimation for the polynomial approximation as well as the magnitudes of the diffusion coefficients of the Brownian term and the Lévy term are shown in Table 1. It seems that the drift term is better represented using the polynomial approximation but the diffusion coefficients are predicted more accurately using the extra neural network to represent the drift function. Overall, the accuracy is higher when we consider separately the Brownian noise and the Lévy noise as shown in Appendices B and C, respectively, compared to the combined case shown here.
The corresponding FP equation is: In the next sections we present results for the 2D Brownian and Lévy noise. In Appendices D and E we also show results of easier problems, where X t and Y t are independent of each other.
In sub-case (a), we construct the NN for the PDF and learn the parameters λ 0 , λ 1 , λ 2 , λ 3 , σ 1 , σ 2 and σ 3 . The results are shown in Figure 5(a). In sub-case (b), we construct the NN for the PDF and two different NNs to approximate the two drift terms, and learn the parameters σ 1 , σ 2 and σ 3 . The results are shown in Figure 5(b). In sub-case (c), we construct the NN for the PDF and use one NN to approximate the potential of the drift terms and take the gradient to obtain the drift term, and learn the parameters σ 1 , σ 2 and σ 3 . The results are shown in Figure 5(c). The corresponding results of parameter estimation are shown in Table 2. We can see that if we use the PINN to approximate the PDF and learn the parameter in the drift term and noise intensity, the results will be better even if we only have five snapshots of the observation data. However, if we also use a neural network to approximate the drift term, the errors in estimating the parameters are larger than 10%. We also give the phase diagram of the drift term in Figure 6. For the drift term results, we can see that the best results are obtained by learning the parameters and fitting the drift term explicitly.
In sub-case (a), we construct the NN for the PDF, then learn the parameter λ 0 , λ 1 , λ 2 , λ 3 , ε x and ε y . The results are shown in Figure 7(a). In sub-case (b), we construct the NN for the PDF and two different NNs to approximate the two drift terms, and learn the parameters ε x and ε y . The results are shown in Figure 7(b). In sub-case (c), we construct the NN for the PDF and one NN to approximate the potential Φ of the drift terms, and learn the parameters ε x and ε y . The results are shown in Figure 7(c). The phase diagram of the drift term is shown in Figure 8. The results of parameter estimation are shown in Table 3;

High-dimensional problems
In this section, we present results for 5D problem with Brownian motion. In Appendix F we present results for similar 3D and 4D cases. We consider the following 5D SDE: where the a i (x) = x − x 3 and n = 5. We assume that the observation data are available at t = 0.1, 0.2, 0.3, 0.5, 0.7, 0.9, 1 i.e., we have N = 7 snapshots. For each time, we have 100, 000 samples, and we set the minibatch size b = 50, 000 to train the first part in (14). We use the Monte-Carlo method with 200, 000 samples to approximate the integral part of (14) and 10, 000 residual points for the FP equation.

Summary
We have developed a new framework based on the physics-informed neural networks (PINNs) and the Fokker-Planck (FP) equation that enables us to learn the probability density function (PDF) at all times given a few snapshots of the particles. The straightforward approach to accomplishing this would have been to estimate the PDF from the given particle observations and then used this estimated PDF as data in PINNs to infer the solution at all times. We have indeed implemented this approach but we demonstrated that in most cases and certainly when the data available is small, it is more accurate to include directly the particle observations in the loss function. The key idea connecting the particle observations and the PDF is based on the variational form of the Kullback-Leibler (KL) divergence as shown in Eqn 8. Here, we have considered inverse problems from 1D to 5D, of first type I, where the only unknown is the initial distribution function and we aimed to infer the PDF dynamics at all times, and of type II, where, in addition, the drift terms and noise intensity (Brownian or Levy) are unknown. The drift terms are functions of the position, and hence we have attempted to either model them parametrically assuming a polynomial form or directly using a neural network; for example, for a 2D case we employed two neural networks to approximate the drift terms. Fig. 1 shows a schematic of the PINN that we employ in the current study.
While we did not attempt to investigate the minimum number of snapshots required to be able to solve the above two types of inverse problems, we found that for even a single snapshot and 100 samples we could obtain the initial distribution and infer the PDF dynamics with reasonable accuracy at all times (e.g., for problem I, see Fig. 2). For problem II, both for the Brownian and Lévy noise, the noise intensities are more difficult to learn, and hence more snapshots with higher number of samples are required to obtain accurate results. However, even in problem II, we employed a relatively few number of snapshots, e.g., 5 or 7. Overall, parametrizing the drift terms as polynomial functions leads to more accurate results for the drift terms but this requires prior knowledge of the form of the drift terms. With regards to the dimensionality, here we demonstrated that we can easily obtain accurate results in 5D but more work is needed on active learning to investigate smart sampling ways in order to achieve similar results in higher dimensions. In a companion work in [17], we employed generative neural networks (GANs) and the stochastic differential equation governing the particles' trajectories to obtain accurate results in 20D with a relatively small number of samples.

Appendices
In the following, we present further tests for both the inverse problems I and II driven by either Brownian noise or Lévy noise in 1D and 2D. We also investigate problem II for Brownian noise in 3D and 4D.

A Problem I: Lévy noise
Here, we consider the drift term a(x) = x − x 3 , α = 1.5, σ = 0 and ε = 1. We use the compound trapezoid formula with 201 points to approximate the integral part of (14) and 800 points to compute the residual part (12). For the inverse problem I, we assume that the observation data is available only for one snapshot at t = 0.3 with 100 (case a), 1000 (case b) and 10000 (case c) samples. The results are shown in Figure 10. In the second column, we use the density estimation to obtain the PDF (blue line in this figure). Same as in the case with Brownian motion, we also observe convergence with respect to the number of samples, as expected. We also compute the inverse problem I of Lévy noise with observations at multiple snapshots, namely, at t = 0.1, 0.3, 0.5. For each time, we have 200 independent samples. We use the kernel density estimation to approximate the density p KD using the observation data. When we have the p KD , we use this approximate density as observation data of PDF and use the neural network to train the model to obtain the PDF at all times; we set p N N KD as the output. We use N p = 201 equidistant points as the collocation points in (16). For each case we test different values of weights τ in (15), and for each τ we run the code for 3 times with different samples, initialization and random seeds. In the first row of Figure 11, we present the errors for different τ . We just show one of the three cases with τ = 100 in the second and third row of Figure 11. We compare these results with the results by the PINN method using the observation data on the particles directly. The comparison is shown in Figure 11; we can see that the results from the latter method are more accurate.

B Problem II: Brownian noise
For the inverse problem II, we compute two cases corresponding to different observation data. Considering the probability shown in Figure 2, we can see that for early times, there is only one peak, but later there are two peaks emerging. Hence, we want to verify if small observation data is sufficient in order to capture the two peaks. To this end, we consider two cases of observation data. In the first case, observation data is available at t = 0.1, 0.3, 0.5, while in the other case data is available at t = 0.2, 0.5, 1. We will use the polynomial (a(x) = λ 0 +λ 1 x+λ 2 x 2 +λ 3 x 3 ) or alternatively a neural network to approximate the drift term. For case (a1), the observation data is available at t = 0.1, 0.3, 0.5 and we use the polynomial fit to approximate the drift term. For case (a2), the observation data is available at t = 0.1, 0.3, 0.5 and we use the neural network to represent the drift term. For case (b1), the observation data is available at t = 0.2, 0.5, 1 and we use the polynomial to approximate the drift term. For case (b2), the observation data is available at t = 0.2, 0.5, 1 and we use the neural network to approximate the drift term.
The results are shown in Figure 12. we can see that when the observation data is available at t = 0.1, 0.3, 0.5, we can still learn the drift term well, but the error of the drift term is larger than the case when the observation data is available at t = 0.2, 0.5, 1. For the PDF, the small observation data can also capture the two peaks of the PDF. The corresponding parameter estimation is shown in Table 4. We can see that the error for case (a) is larger than the error of case (b).

C Problem II: Lévy noise
For the inverse problem II driven by Lévy noise, we compute four cases like in the Brownian noise. We assume we have 10, 000 available samples at each observation time, and set the minibatch size b = 1, 000 to train the first part in (14). We use the compound trapezoid formula with 201 points to approximate the integral part and randomly choose 800 residual points in (12). For case (a1), the observation data is available at t = 0.1, 0.3, 0.5 and we use a polynomial to approximate the drift term. For case (a2), the observation data is available at t = 0.1, 0.3, 0.5 and we use a neural network to approximate the drift term. For case (b1), the observation data is available at t = 0.2, 0.5, 1 and we use a polynomial to approximate the drift term. For case (b2), the observation data is available at t = 0.2, 0.5, 1 and we use a neural network to approximate the drift term. We infer the drift term and the PDF at different times in Figure 13. The parameter estimation results are shown in Table 5.
where a 1 (x, y) = x − x 3 , a 2 (x, y) = y − y 3 , and σ x = σ y = 1. we use the compound trapezoid formula with 301 × 301 points to approximate the integral part of (14) and 5000 points for the residual of the FP equation. We present four cases to compute the inverse problem. For case (a), the observation data is available at t = 0.1, 0.4, 0.7, 1 and the drift term is approximated by polynomials. For case (b), the observation data is available at t = 0.1, 0.4, 0.7, 1 and the drift term is approximated by neural networks. For case (c), the observation data is available at t = 0.1, 0.3, 0.5, 0.7, 1 and the drift term is approximated by polynomials. For case (d), the observation data is available at t = 0.1, 0.3, 0.5, 0.7, 1 and the drift term is approximated by neural networks. For each time snapshot, we have 100, 000 samples, and set the minibatch size as b = 10, 000 for the first part in (14).
The predictions for the drift terms are shown in Figure 14. We use the MC method to obtain the reference PDF and the results are shown in top of Figure  15. Specifically, we compare the error results of PDF at t = 1 for four cases in Figure 15. The parameters we learn are shown in Table 6. We can see that the results of case (a) are better than case (b), while the results of case (c) are better than case (d), especially for the drift term.

E Problem II: 2D Lévy noise
We consider another example of 2D Lévy noise: where a 1 (x, y) = x − x 3 , a 2 (x, y) = y − y 3 , and ε x = ε y = 1. We use the compound trapezoid formula with 201 × 201 points to approximate the integral part of (14) and 160, 000 residual points for the FP equation. For each time snapshot, we have 100, 000 samples, and set the minibatch size as b = 5, 000 for the first part in (14). For case (a), the observation data is available at t = 0.1, 0.4, 0.7, 1 and we use polynomials to approximate the drift term. For case (b), the observation data is available at t = 0.1, 0.4, 0.7, 1 and we use two neural networks to approximate the drift term. For case (c), the observation data is available at t = 0.1, 0.3, 0.5, 0.7, 1 and we use polynomials to approximate the drift term. For case (d), the observation data is available at t = 0.1, 0.3, 0.5, 0.7, 1 and we use two neural networks to approximate the drift term. The results are shown in Figure 16, 17 and Table 7.

F Problem II: 3D and 4D Brownian noise
We consider the 3D and 4D case of (22), i.e., n = 3 and n = 4. Similar as in case (b) in Section 3.3, we use a cubic polynomial to parameterize the drift.