Identification of Stochastically Perturbed Autonomous Systems from Temporal Sequences of Probability Density Functions

The paper introduces a method for reconstructing one-dimensional iterated maps that are driven by an external control input and subjected to an additive stochastic perturbation, from sequences of probability density functions that are generated by the stochastic dynamical systems and observed experimentally.


Introduction
There is considerable interest in modeling and analyzing dynamical systems that generate densities of states. Examples of such systems include chaotic systems (Boyarsky and Góra 1997;Lasota and Mackey 1994) and stochastically perturbed dynamical sys-Communicated by George Haller.
In many practical situations, the system that generates the density of states is unknown and only the densities of states generated by the system or the invariant density associated with the system can be observed, while the individual point trajectories are not measurable. Conventional solutions (Maguire et al. 1998;Han et al. 2004;Príncipe and Kuo 1995;Lai et al. 1999;Lai and Tél 2011;Bollt et al. 2001) rely on time series observations, but for such situations they become unsuitable. The problem of inferring the unknown dynamical system from the observed densities is known as the inverse Frobenius-Perron problem (Boyarsky and Góra 1997;Ershov and Malinetskii 1988). The problem of reconstructing an unknown onedimensional autonomous chaotic map given only knowledge of the invariant density function of the system has been considered by a number of authors (Ershov and Malinetskii 1988;Góra and Boyarsky 1993;Diakonos and Schmelcher 1996;Pingel et al. 1999), while there are special cases in which this problem has a unique solution. Given the invariant symmetric beta density functions, methods were introduced to construct a class of symmetric maps (Diakonos and Schmelcher 1996) and a broader class of continuous unimodal maps whose each brand covers the complete interval . Given arbitrary invariant densities other similar approaches were proposed for identifying the maps with specified forms: two types of one-dimensional symmetric maps (Koga 1991), smooth chaotic map with closed form (Huang 2006, multi-branches complete chaotic map ). Problems of synthesizing one-dimensional maps with prescribed invariant density function or autocorrelation function were tackled in Baranovsky and Daems (1995) and Diakonos et al. (1999). Using positive matrix theory an approach to synthesizing chaotic maps with arbitrary piecewise constant invariant densities and arbitrary mixing properties was developed in Rogers et al. (2004). This method was further extended to synthesizing dynamical systems with desired statistical properties (Rogers et al. 2008a), developing communication networks (Berman et al. 2004) and designing randomly switched chaotic maps and two-dimensional chaotic maps used for image generation (Rogers et al. 2008b). In Bollt (2000) and Bollt and Santitissadeekorn (2013), a global and open-loop strategy of controlling chaos was presented to solve the inverse problem. The problem was reduced to that of finding a perturbation of the original Frobenius-Perron matrix to achieve the target invariant density function. In general, given only invariant density function, the solution to the inverse problem is not unique, as different maps exhibiting remarkably different dynamics may possess a same invariant density function. Therefore, additional assumptions or constraints are required to ensure the uniqueness of the identification results. A more recent approach (Nie and Coca 2015) addresses the uniqueness issue by considering sequences of density functions generated by the system rather than just the invariant density function of the system. This method allows inferring the map that exhibits the same transient and asymptotic dynamics as the underlying system that generated the data. Although it is shown that the method is robust to noise, the approach does not exploit any a priori knowledge of the noise distribution. In addition, to our knowledge, all existing methods consider only autonomous maps.
In this context, this paper introduces for the first time a method to infer a onedimensional map that is driven by an external control input while being subjected to an additive stochastic perturbation from sequences of observed density functions generated by the unknown system. We formulate the operator transferring the state density function of the stochastic dynamical system in terms of the Frobenius-Perron operator associated with the unperturbed underlying system that we aim to estimate, and derive the matrix representation of the transfer operator in terms of the Frobenius-Perron matrix. Based on this representation the new algorithm is developed to estimate the Frobenius-Perron matrix using temporal sequences of probability density functions generated by the stochastic dynamical system given the density functions of the control input and noise. The approach also determines the monotonicity of general nonlinear transformations over each interval of the partition, which is a crucial step to reconstruct the true dynamical system.
The paper is structured as follows. Section 2 introduces the inverse problem. The stochastic Frobenius-Perron operator associated with stochastically perturbed autonomous systems is derived in Sect. 3. A matrix approximation of the operator is given in Sect. 4. Section 5 introduces the methodology of reconstructing general nonlinear maps from sequences of density functions. Section 6 presents a numerical simulation example to demonstrate the effectiveness of the developed algorithm for the stochastically perturbed autonomous systems. Conclusions are given in Sect. 7. Let (R [0,b], B, μ) be a normalized measure space, where μ is a measure on (R, B) and B is a Borel σ -algebra of subsets in R. Consider the following discrete-time stochastic dynamical system

Inverse Problem Formulation
where S: R → R is a measurable and nonsingular transformation [i.e., μ(S −1 ( A)) B for any A ∈ B and if μ(S −1 ( A)) 0 for all A ∈ B, then μ( A) 0], x n ∈ R is the state variable having the probability density function f n ∈ D (R, B, μ), is the control input of the system with a probability density function f u ∈ L 1 (R) that can be assigned, and ξ n is an independent random variable with a known probability density function g that has compact support on [− ε, ε], that is, , K be random vectors of initial and final state observations, respectively, such that where i 1, . . . , K . Assuming that for practical reasons it is not possible to track individual point trajectories during the experiment, that is to associate an initial state x i 0, j with its image x i 1, j under the transformation, the inverse problem considered in this paper is to infer the point transformation S in (1) from the probability density functions f 0, j and f 1, j of the initial and final states X 0, i {x 0, j } θ j 1 and X 1, i {x 1, j } θ j 1 , i = 1, …, K.

The Stochastic Frobenius-Perron Operator Associated with the Stochastically Perturbed Transformation
In this section the transfer of density function at n to n + 1 is derived given the input and noise density functions f u and g. For a dynamical system with a constantly applied random perturbation written in the following general form (Lasota and Mackey 1994) x n+1 where S: R → R is a given transformation and ξ n is an independent random variable having a density function g. The operator transferring state density functions of the perturbed dynamical system are called the stochastic Frobenius-Perron operator, denoted byP,P where τ (x, y) g(x − S(y)) is a stochastic kernel, satisfying τ (x, y) > 0, and R τ (x, y) 1. For a nonsingular unperturbed transformation S, the Frobenius-Perron operator (Boyarsky and Góra 1997) corresponding to S exists, denoted by P S , and (4) is further written asP Let G: R × R → R be defined by such that (1) can be written as Letx n+1 G(x n , u n ) ∈ R. From (5) it follows that the probability density function ofx n+1 is given bȳ where χ Δ (x) is the indicator function defined by Equation (7) becomes where the probability density function of x n+1 is given by Substituting (8) into (11) leads to the following formulation of the stochastic Frobenius-Perron operator, denoted byP, associated with stochastic dynamical system (1) (12) Equation (12) relates the operatorP, corresponding to the stochastic system, to the Frobenius-Perron operator P S associated with the map S. This equation forms the basis for the new approach to reconstruct the map S based on sequences of density functions.

Remark 1
The additive noise ξ n is an i.i.d. random variable that normally satisfies that max(|ξ n |) < b in general practical measurements. For the unusual case ε > b, (10) can be rewritten as where k 1 . Since k 1 and k 2 have infinite results given only x andx, f n+1 cannot be uniquely generated fromf n+1 in (11). Hence, ξ n is treated as a variable bounded in [− ε, ε], ε ≤ b.
Remark 2 An alternative compact way of formulating the stochastic Frobenius-Perron operator is to apply the joint density function denoted by f α ∈ L 1 (R) for the control input and noise to (4). Let α n u n + ξ n (mod b). Thus, f α can be given in terms of f u and g by It follows that x n+1 S(x n ) + α n (mod b), and from (5) we have that Substituting (14) into (15) gives that . It follows that (12) is obtained from (16).
In the first instance, we assume that S belongs to a special class of nonlinear transformations called piecewise linear semi-Markov transformations and develop the algorithm to reconstruct it. We then show how the reconstruction approach can be applied to approximate more general one-dimensional maps.

A Matrix Representation of the Transfer OperatorP
Let S be a piecewise linear and expanding semi-Markov transformation over the Ninterval partition, (Góra and Boyarsky 1993) The restriction S| R i is a homeomorphism from R i to a union of intervals of where and p(i) denotes the number of disjoint subintervals Q (i) k corresponding to R i . Let f n be a piecewise constant function over the partition such that f n (x) According to the property of semi-Markov map (Boyarsky and Góra 1997), its image under transformation P S f n is also a piecewise constant function over such that P S f n (x) . In this case, the Frobenius-Perron operator can be represented by a finite-dimensional matrix such that where M (m i, j ) 1≤i, j≤N is the Frobenius-Perron matrix induced by S with entries given by From (18) it follows that be the coefficient vectors of the piecewise constant density functions f n and P S f n over the partition , respectively. We have ϕ P S f N n w f N n M. By integrating both sides of (12) over R i ∈ , it follows that For f n ∈ L 1 we define where λ denotes the Lebesgue measure and f N n+1 denotes the piecewise constant approximation of f n+1 over the partition . We have the following result (Li 1976).
Substituting (17) in (21) gives Let H (h i, j ) 1≤i, j≤N be a matrix with entries given by It follows from (20) and (24) that Let Q M H . The evolution of density functions is formulated as w f N n+1 w f N n Q. Q is the matrix representation of the transfer operatorP. Formula (26) yields the final density function estimated over the N-interval partition , mapping from the initial piecewise constant density function over the N-interval partition . This establishes the basis of the new algorithm of reconstructing the unknown transformation S from sequences of probability density functions.
Remark 3 Given the nonsingular transformation S: R → R that induces the Frobenius-Perron matrix M with respect to the partition , input density function f u ∈ L 1 and noise density function g ∈ L 1 , from (26) the estimated state density function over of stochastic dynamical system (1) can be predicted from a piecewise constant initial density function Remark 4 Let Q (q i, j ) 1≤i, j≤N , where from (26) q i, j is given by It follows that This implies that matrix Q is a stochastic matrix that has 1 as the eigenvalue of maximum modulus, of which the algebraic and geometric multiplicities are 1. Since Q and Q have the same eigenvalues, we then have represents the equilibrium density vector of Q.

Remark 5 Remark 4 suggests that there exists a stationary density function
for the transfer operatorP N . It follows from Lemma 1 that f N * (x) converges to f * (x) of the stochastic dynamical system as N → +∞.

Solving the Stochastic Inverse Frobenius-Perron Problem for Continuous Nonlinear Transformations
This section introduces a method to reconstruct the underlying map S in Eq. (1) based on a sequence of probability density functions estimated from data, under the general assumption that S is a continuous nonlinear map. Specifically, the method infers a piecewise linear semi-Markov mapŜ with respect to a uniform partition a 2 ], . . . , (a N −1 , a N ]}, a N b, given K random vectors of initial states X 0, i {x 0, i j } θ j 1 , from K initial state densities f 0, i , i = 1, …, K, the corresponding final state vectors X 1, i {x 1, i j } θ j 1 , i = 1, …, K under transformation (1) and the density of the noise and of the control input, g and f u , respectively. The matrix M associated with P S can be approximated arbitrarily well, and thus,Ŝ approximates S to an arbitrary accuracy as N → +∞. While g is fixed, f u can be defined by the user when the experiment is conducted. It is assumed that the correspondence between an initial state measurement x 0, j and its image x 1, j under the transformation is not known and hence the point transformation S in (1) has to be inferred based on the probability density functions f 0, j K j 1 , f 1, j K j 1 , g and f u . The proposed reconstruction algorithm for general nonlinear and continuous maps is summarized below. However, it is worth emphasizing that this method can also be used in cases when S is piecewise semi-Markov.
Step 1: corresponding to the piecewise constant density functions f N t, i (x) that approximate the new state density functions f t, i (x) over the regular partition . Compute the matrix H; Step 3: Identify a trial Frobenius-Perron matrixM firstly to determine the indices of consecutive positive entries of the matrix M that represents the Frobenius-Perron operator P S associated with the optimal approximate mapŜ and subsequently a refined matrix M; Step 4: Construct the approximate piecewise linear semi-Markov transformation on , and smooth it to obtain the continuous nonlinear map.
These steps are described below in more detail.

Step 1: Observe Sets of States to Assemble Sequences of Densities
Let f 0, i be a set of different initial density functions that is piecewise constant on the partition where the coefficients satisfy N j 1 w 0, i j N b , i = 1, …, K. Let X 0, i {x 0, i j } θ j 1 be the set of initial conditions obtained by sampling f 0, i (x), and X t, i {x t, i j } θ j 1 be the set of states obtained by applying t times Eq.
(1) such that are generated by sampling f u and g, respectively.
From Remark 4, given the input and noise density functions, the generated densities converge to a stationary density function regardless of the initial conditions. Therefore, finite densities characterizing the transient dynamics and evolving from an initial density function can be observed. For K sequences of densities, the most dynamical behavior exhibited by the perturbed underlying system can be observed from T m iterations given by T m min{t m }, representing the minimum steps taken to approach the stationary density, where t m is an integral set given by Thus, 1 ≤ T ≤ T m . Typically, the interval number of is set by 1 < N ≤ K T .

Step 2: Estimate the Coefficients w and Compute the Matrix H
The piecewise constant density function f N 1, i (x) on the partition is given by These following matrices are then derived.
and (35) Given the input and noise density functions f u and g, the matrix H is computed from (25).

Step 3: Identify the Frobenius-Perron Matrix M
For a continuous nonlinear map, the corresponding Frobenius-Perron matrix M must satisfy that the positive entries in each row are contiguous. Without enough constraints to optimize the matrix, it is generally difficult to identify a very fine matrix straightforwardly. Therefore, initially, a trail Frobenius-Perron matrix is derived to determine the indices of contiguous positive entries in each row, which are then used to refine the matrix. This is carried out in two stages. In the first instance, given (21) the coordinate vectors ϕ P S f n , n 0, . . . , T − 1, are obtained by solving the following constrained optimization problem min 0≤{ϕ n j } n 0,...,T −1 and ||·|| F denotes the Frobenius norm. Subsequently, the trial matrix denoted byM (m i, j ) 1≤i, j≤N is obtained as a solution to the following constrained optimization problem for i 1, . . . , N , k 1, . . . ,p(i) − 1. The approximation to the continuous map may have an infinite number of pieces of monotonicity, and each piece S| R i can be linearly approximated. Thus, for a piecewise linear semi-Markov approximationŜ, the maximum and minimum column indices of positive entries on two contiguous rows of M are further refined by r (i, p(i)) [r(i,p(i)) +r (i + 1, 1)] 2 , r (i + 1, 1) [r (i,p(i)) +r (i + 1, 1)] 2 , and S Q (i) [r (i, 1) +r (i + 1,p(i + 1))] 2 , r (i + 1, p(i + 1)) [r (i, 1) +r (i + 1,p(i + 1))] 2 and S Q (i) k 1 r (i + 1, k) and r (i, 1) −r (i + 1,p(i + 1)) > 1, and that S Q (i) k is the newly formed subinterval, and {r (i, 1), . . . , r (i, p(i))} are the identified column indices of positive entries in the ith row of the matrix M.
The refined Frobenius-Perron matrix M is then obtained by solving the following optimization problem k 1 m i, r (i, 1)+k−1 1 and m i, r (i, k) > 0, for i 1, . . . , N and m i, j 0, if j r (i, k), for k 1, . . . , p(i).

Step 4: Construct the Nonlinear Map
This step involves reconstructing the semi-Markov map that corresponds to the identified Frobenius-Perron matrix M. For a continuous map, it started with determining the monotonicity of each branch S| Q (i) ] be the image of the interval R i under the semi-Markov transformationŜ associated with the identified Frobenius-Perron matrix M. Denote a r (i, 1)−1 as the starting point of R r (i, 1) mapped from the subinterval Q (i) 1 , and a r (i, p(i)) as the end point of R r (i, p(i)) , the image of the subinterval Q for i 2, . . . , N and γ (1) γ (2). Given that the derivative of S| Q (i) k is 1 m i, j , the end point q within R i is given by where k 1, . . . , p(i)−1 and q for m i, j 0, i 1, . . . , N , j 1, . . . , N , k 1, . . . , p(i) − 1. A smooth nonlinear map is then obtained by fitting a polynomial smoothing spline.
In practice, there are no restrictions on the shape of this density function. The density functions of the input and noise, f u and g, are shown in Fig. 1. For the purpose of inferring the piecewise linear semi-Markov transformation that approximates the original logistic map S, we define a uniform partition of [0, 1] with N 40 intervals. To generate the data used in the reconstruction, K 40 piecewise constant initial density functions f 0, i χ R i (x), i 1, . . . , 40 were sampled to generate the initial states X 0, i {x 0, i j } θ 1 j 1 , θ 1 5 × 10 3 , i 1, . . . , 40. The input and noise densities were sampled to generate the input and noise data sets U {u i } θ 1 i 1 and Ξ {ξ i } θ 1 i 1 . In total, 40 sequences of new states X t, i {x t, i j } θ 1 j 1 , i 1, . . . , 40 were then observed by iterating t times system (44) and these were subsequently used to estimate the corresponding piecewise constant densities f N t, i , i 1, . . . , 40, t ≥ 1, over the uniform partition . Figure 2 shows the calculated results of the performance function J(t) in (32), 1 ≤ t ≤ 23, which represents the summing differences of each two successive densities of the K sequences of densities. As can be seen, the minimum can be found at t ≥ 4, which suggests that densities in all the sequences approach the equilibrium distribution. It follows that 1 ≤ T ≤ 4. Here we choose T = 1. Figure 3 shows the initial densities f 0, k and their image densities f N 1, k , which are used to reconstruct the approximate map, and also f N 2, k , f N 3, k and the equilibrium observed after 1 × 10 4 iterations. It can be seen that compared with f N 1, k , f N 2, k and f N 3, k are more close to the stationary density, and densities in each sequence rapidly converge to the same stationary density. This is also demonstrated in Fig. 2 that the derivative of f i, k at i = 3 is apparently lower than that of f 1, k and f 2, k . Using the lsqlin function in the optimization toolbox of MATLAB to solve optimization problems (36), (38) and (40) in the proposed algorithm, the Frobenius-Perron matrix M is obtained.
The piecewise linear semi-Markov mapŜ associated with the identified matrix M is shown in Fig. 4a. Finally, the continuous nonlinear mapS was estimated by fitting a cubic smoothing spline with the smoothing parameter 0.999, to a set of 10 3 data points obtained by uniformly sampling the piecewise linear mapŜ over [0,1]. The reconstructed continuous nonlinear map is shown in Fig. 4b. The performance of the reconstruction algorithm is evaluated by computing the relative percentage error (RPE) between the original and estimated mapsŜ and between the original and the smoothed mapS, respectively, which is illustrated in Fig. 5. As can be seen, the full error for S and 95% of the error forS is lower than 5%. With the increase of N, the estimated mapŜ is more close to S. To further evaluate the accuracy of the reconstruction, the constructed piecewise linear semi-Markov approximationŜ and the estimated continuous mapS were used to predict n-iteration ahead density functions, n 1, . . . , 60, respectively, using a Gaussian distribution function N (0.6, 0.4 2 ) truncated to [0, 1] as the initial state density function f 0 , the input density given in (45) and the noise density given in (46). With 100 sets of θ 2 1 × 10 4 input data U {u k, i } 100, θ 2 k 1, i 1 generated by sampling f u and 100 sets of the same number of noise data Ξ {ξ k, i } 100, θ 2 k 1, i 1 from g, 100 sets of θ 2 randomly distributed initial states X 0, k {x 0, k j } θ 2 j 1 , k = 1, …, 100, were separately iterated for 60 steps using stochastic model (1) with the original map S, the identified piecewise linear semi-Markov approximationŜ and the estimated continuous mapS, respectively. In each step piecewise constant density functions f 40 n, k (x) n, k j χ R j (x), k = 1, …, 100, n = 1, …, 60, are then estimated over from the generated states. The root-mean-square error (RMSE) between f 40 n, k andf 40 n, k and between f 40 n, k andf 40 n, k is calculated by RMSE(S,Ŝ) n,k 1 40 whereŵ n, k i andw n, k i are the coefficients of predicted density function usingŜ andS, respectively. The mean, 10 and 90% quantiles of the 100 RMSEs forŜ andS at each iteration are shown in Figs. 6 and 7. As can be seen from them, 90% quantiles of the error remain constantly less than 0.1 for bothŜ andS, and their mean values stabilize around 0.08 after 10 iterations.  Figure 8 shows the RMSE between S andŜ on 100 uniformly spaced points within for T = 1, …, 8. As can be seen, a small decrease of the error occurs from T = 1 till 4, and then the error maintains almost constant for T ≥ 4. This implies that all the sequences reach the equilibrium distribution after 4 iterations, which is in consistent with Fig. 2. From Fig. 2, distance between f 1, k and f * is remarkably larger than that between f n, k and f * for n = 2,3,4, and therefore the error is slightly diminished even though more new densities are added for the identification.

Conclusions
This paper introduced a new algorithm for reconstructing the underlying onedimensional map for an autonomous dynamical system that is driven by an additive control input and also subjected to an additive stochastic perturbation, given the observed sequences of probability density functions generated by the unknown sys-tem, and the input and noise density functions. Evolution of densities was formulated and described by a stochastic Frobenius-Perron operator that has a matrix representation. This forms the basis for the algorithm to identify the Frobenius-Perron matrix associated with a piecewise linear semi-Markov approximation to the underlying nonlinear map. Based on the matrix representation of the stochastic Frobenius-Perron operator the densities generated by the dynamical system and evolving from a given initial condition can be predicted. Convergence of the evolving densities analyzed from the matrix representation reveals a fact that only a limited number of densities characterizing the transient dynamics is observable for arbitrary initial condition, and thus, this requires different initial conditions so as to generate as many as possible temporal sequences of densities to reconstruct the underlying map.
For the situations where only a limited number of initial conditions are available for generating the temporal sequences of densities that converges quickly to the equilibrium distribution, a potential effective solution is to apply multiple linearly independent input density functions to the stochastic dynamical systems so that the densities would diverge to different equilibrium distributions, which will be further explored. From a practical perspective, it is also worthwhile to extend the approach to higher-dimensional systems based on sequences of mixture densities generated by the more complex systems.
Furthermore, this paper provides a new insight into identification of stochastic dynamical systems given the density functions of control inputs. It triggers a new scheme to solve the control problem for such systems. Specifically, given the noise density function, the problem aims to determine the optimal input density function so that the dynamical system can have a desired equilibrium distribution that represents the targeted asymptotic dynamics.