Reconstruction of one-dimensional chaotic maps from sequences of probability density functions

In many practical situations, it is impossible to measure the individual trajectories generated by an unknown chaotic system, but we can observe the evolution of probability density functions generated by such a system. The paper proposes for the first time a matrix-based approach to solve the generalized inverse Frobenius–Perron problem, that is, to reconstruct an unknown one-dimensional chaotic transformation, based on a temporal sequence of probability density functions generated by the transformation. Numerical examples are used to demonstrate the applicability of the proposed approach and evaluate its robustness with respect to constantly applied stochastic perturbations.

physics and economics [1], which generate density of states. Examples include particle formation in emulsion polymerization [2], papermaking systems [3], bursty packet traffic in computer networks [4,5], cellular uplink load in WCDMA systems [6]. A major challenge is that of inferring the chaotic map that describes the evolution of the unknown chaotic system, solely based on experimental observations. Starting with seminal papers of Farmer and Sidorovich [7], Casadgli [8] and Abarbanel et al. [9], the problem of inferring dynamical models of chaotic systems directly from time series data has been addressed by many authors using neural networks [10], polynomial [11] or wavelet models [12].
In many practical applications, it is more convenient to observe experimentally the evolution of the probability density functions, instead of individual point trajectories, generated by such systems. For example, the particle image velocimetry (PIV) method of flow visualization [13] allows identifying individual tracer particles in each image, but not to track these between images. In biology, flow cytometry is routinely used to measure the expression of membrane proteins of individual cells in a population [14]. However, it is impossible to track the cells between subsequent analyses.
The problem of inferring an unknown chaotic map given the invariant density generated by the map is known as the inverse Frobenius-Perron problem (IFPP). This inverse problem has been investigated by Friedman and Boyarsky [15] who treated the inverse problem for a very restrictive class of piecewise constant density functions, using graph-theoretical methods. Ershov and Malinetskii [16] proposed a numerical algorithm for constructing a one-dimensional unimodal transformation which has a given invariant density. The results were generalized in Góra and Boyarsky [17], who introduced a matrix method for constructing a 3-band transformation such that an arbitrary given piecewise constant density is invariant under the transformation. Diakonos and Schmelcher [18] considered the inverse problem for a class of symmetric maps that have invariant symmetric Beta density functions. For the given symmetry constraints, they show that this problem has a unique solution. A generalization of this approach, which deals with a broader class of continuous unimodal maps for which each branch of the map covers the complete interval and considers asymmetric beta density functions, is proposed in [19]. Huang presented approaches to constructing smooth chaotic transformation with closed form [20,21] and multibranches complete chaotic map [22] given invariant densities. Boyarsky and Góra [23] studied the problem of representing the dynamics of chaotic maps, which is irreversible by a reversible deterministic process. Baranovsky and Daems [24] considered the problem of synthesizing one-dimensional piecewise linear Markov maps with prescribed autocorrelation function. The desired invariant density is then obtained by performing a suitable coordinate transformation. An alternative stochastic optimization approach is proposed in [25] to synthesize smooth unimodal maps with given invariant density and autocorrelation function. An analytical approach to solving the IFPP for two specific types of one-dimensional symmetric maps, given an analytic form of the invariant density, was introduced in [26]. A method for constructing chaotic maps with arbitrary piecewise constant invariant densities and arbitrary mixing properties using positive matrix theory was proposed in [5]. The approach has been exploited to synthesize dynamical systems with desired characteristics, i.e. Lyapunov exponent and mixing properties that share the same invariant density [27] and to analyse and design the communication networks based on TCPlike congestion control mechanisms [28]. An extension of this work to randomly switched chaotic maps is studied in [29]. It is also shown how the method can be extended to higher dimensions and how the approach can be used to encode images. In [30], the inverse problem is formulated as the problem of stabilizing a target distribution using an open-loop perturbation approach. In order to characterize the patterns of activity in the olfactory bulb, an optimization approach was proposed to infer the elements of the Frobenius-Perron matrices that encode the invariant density functions of interspike intervals corresponding to different odours [31].
All existing methods can be used to construct a map with a given invariant density. A limitation of these approaches is that the solution to the inverse problem is not unique. Typically, there are many transformations, exhibiting a wide variety of dynamical behaviours, which share the same invariant density. Therefore, the reconstructed map does not necessarily exhibit the same dynamics as the underlying systems even though it preserves the required invariant density. Additional constraints and model validity tests have to be used to ensure that the reconstructed map captures the dynamical properties of the underlying system (Lyapunov exponents, fixed points, etc.) and predicts its evolution. This is of paramount importance in a many practical applications ranging from modelling and control of particulate processes [48], characterizing the formation and evolution of the persistent spatial structures in chaotic fluid mixing [32], characterizing the chaotic behaviour of electrical circuits [33], chaotic signal processing [34,35], analysing and interpreting cellular heterogeneity [36,37] and identification of molecular conformations [38]. Potthast and Roland [39] solved the Frobenius-Perron equation describing the evolution of uniform probability distributions under the action of nonlinear dynamical automata that implement Turing machines. In this context, the realization of nonlinear dynamical automata that describe cognitive processes based on experimental data can be formulated as an inverse Frobenius-Perron problem.
Another major limitation of the existing matrixbased reconstruction algorithms is the assumption that the Markov partition is known in advance. In general, no a priori information about the unknown map is available, so the partition identification problem has to be solved as part of the reconstruction method.
This paper proposes a systematic method for determining an unknown chaotic map given sequences of density functions estimated from data. In other words, the inverse problem studied in this paper is that of determining the map that exhibits the same transient as well as asymptotic dynamics as the underlying system that generated the data. The proposed methodology involves the identification of the Markov partition, esti-mation of the Frobenius-Perron matrix and the reconstruction of the underlying map.
To our knowledge, our approach provides for the first time a solution to the problem of inferring, from sequences of density functions, a broad class of onedimensional transformations that admit an invariant density when the Markov partition is not known in advance.
This paper is organized as follows: in Sect, 2, we give a brief introduction to the inverse Frobenius-Perron problem. The methodology for reconstructing piecewise-linear semi-Markov transformations from sequences of densities is presented and demonstrated using a numerical simulation example in Sect. 3. An extension of the method to continuous nonlinear transformations is introduced and demonstrated numerically in Sect. 4. Conclusions are presented in Sect. 5.

The inverse Frobenius-Perron problem
Let I = [a, b], B be a Borel σ -algebra of subsets in I , and μ denote the normalized Lebesgue measure on I . Let S : I → I be a measurable, non-singular transformation, that is, μ(S −1 (A)) ∈ B for any A ∈ B and μ(S −1 (A)) = 0 for all A ∈ B with μ(A) = 0. If x n is a random variable on I having the probability density function f n ∈ D(I, B, μ), then x n+1 given by is distributed according to the probability density function f n+1 = P S f n , where P S : is the Frobenius-Perron operator [40] associated with the transformation S. In this case, P S can be written explicitly as The Frobenius-Perron operator of a non-singular transformation S is a Markov operator [41].
Definition 1 A linear operator P S : L 1 → L 1 satisfying (a) P S f n ≥ 0 for f n ≥ 0, f n ∈ L 1 ; (b) P S L 1 < 1, and P S f n L 1 = f n L 1 , for f n ≥ 0, f n ∈ L 1 , is called a Markov operator.
Definition 2 A Markov operator P S : L 1 → L 1 is called strongly constrictive if there exists a compact set F ⊂ L 1 such that for any f n ∈ D, If P S is strongly constrictive, according to the spectral decomposition theorem [40], there exist a sequence of densities f 1 , . . . , f r and a sequence of bounded linear functionals g 1 , . . . , g r such that where P n S is the nth iteration of P, the densities f 1 , . . . , f r have mutually disjoint supports ( f i f j = 0 for i = j), and P S f i = f α(i) , i = 1, . . . , r and {α(1), . . . , α(r )} is a permutation of the integers {1, . . . , r }. Furthermore, P n S f converges to an invariant density f * , which satisfies f * = P S f * .
Let R = {R 1 , R 2 , . . . , R N } be a partition of I into intervals, and int( where S i is the monotonic segment of S on each interval R i . The inverse Frobenius-Perron problem is usually formulated [17] as the problem of determining the point transformation S such that the dynamical system x n+1 = S(x n ) has a given invariant probability density function f * . In general, the problem does not have a unique solution.
The generalized inverse problem addressed in this paper is that of inferring the point transformation which generated a sequence of density functions and has a given invariant density function. Specifically, let be two sets of initial and final states observed in K separate experiments, where x j 1,i = S(x j 0,i ), i = 1, . . . , θ, j = 1, . . . , K , and S : I → I is an unknown, non-singular point transformation. We assume that for practical reasons, we cannot associate with an initial state x j 0,i the corresponding image x j 1,i , but we can estimate the probability density functions f j 0 and f j 1 associated with the initial and final states, Moreover, let f * be the observed invariant density of the system. The inverse problem is to determine S :

A solution to the IFPP for piecewise-linear semi-Markov transformations
This section presents a method for solving the IFPP for a class of piecewise monotonic and expanding trans- is The restriction S| Q (i) k is a homeomorphism from R i to a union of intervals of R where The following theorem [17] establishes an important property of such transformation, namely that its invariant density is piecewise constant over R.
denotes the space of piecewise constant functions defined over the partition R, then the Frobenius-Perron operator P S associated with the piecewise-linear Rsemi-Markov transformation satisfies and be a sequence of probability density functions generated by the unknown map S, given a set of initial density func- Assuming that the invariant density function f * of the Frobenius-Perron operator associated with the unknown transformation S can be estimated from experimental data, the proposed identification approach can be summarized as follows: a. Given the samples, construct a uniform partition C and an initial piecewise constant density estimate f * C of the true invariant density f * which maximizes a penalized log-likelihood function. b. Select a sub-partition C d (l j ) of C. c. Estimate the matrix representation of the Frobenius-Perron operator over the partition C d (l j ) based on the observed sequences of densities generated by S. d. Construct the piecewise-linear mapŜ (l j ) corresponding to the matrix representation. e. Compute the piecewise constant invariant density f * C d (l j ) associated with the identified transforma-tionŜ (l j ) and evaluate performance criterion. f. Repeat steps (b) to (e) to identify the partition and map which minimize the performance criterion.

Identification of the Markov partition
Let f * ∈ H (R) be the invariant density associated with a R-semi-Markov transformation S. Let {x * i } θ i=1 be a finite number of independent observations of f * . The aim is to determine an orthogonal basis set where χ R i (x) is the indicator function and h i are the expansion coefficients given by λ(R i ) denotes the length of the interval R i . We start by constructing a uniform partition with intervals N that maximizes the following penalized log-likelihood function [42] The coefficients h i for the regular histogram are given by Let C = {c 1 , . . . , c N −1 } be the strictly increasing sequence of cut points corresponding to the resulting uniform partition The final Markov partition R is determined by solving where denotes the piecewise constant invariant density associated with the transfor-mationŜ (l j ) identified over the partition 3.2 Identification of the Frobenius-Perron matrix be the piecewise constant densities on R, which are estimated from the samples.
Let f 0 (x) be an initial density function that is piecewise constant on the partition R where the coefficients satisfy be the set of initial conditions obtained by sampling f 0 (x) and be the set of states obtained by applying t times the transformation S such that The density function associated with the states X t is given by where the coefficients w t, j = .
In practice, the observed f t (x), t = 0, . . . , T , are approximations of the true density functions, which are inferred from experimental observations. It follows that where and The matrix M is obtained as a solution to a constrained optimization problem where || · || F denotes the Frobenius norm. The matrix = W T 0 W 0 has to be non-singular for the solution to be unique.
Using Cauchy-Binet formula, the determinant of can be written as Since are also linearly independent. Moreover, given that S is a R-semi-Markov, piecewise linear and expanding, we have where . . , N . Alternatively, this can be written as where and  Since W 0 is non-singular, the Frobenius-Perron matrix M is given by The derivative of S| Q (i) k is 1/m i, j , and the length of which allows computing iteratively q (i) k for each interval R i starting with q (i) 0 = c i−1 . By assuming each branch S| R i is monotonically increasing, the piecewise-linear semi-Markov mapping is given by for k = 1, . . . , p(i), and j is the index of image R j of The map is constructed as depicted in Fig. 1. In practice, we can choose the piecewise constant probability density functions f j 0 (x) = 1 λ(R j ) χ R j (x). These are sampled in order to generate N sets of initial conditions that will be used in the experiments. For each set of initial conditions X i 1 , we measure a corresponding set of final states The density function f i 1 associated with the set X i 1 of final states is given by where Remark We only need to generate initial conditions for the densities that correspond to the finest uniform partition N = N . Coarser partitions are obtained by merging adjacent intervals, for example R j and R j+1 , leading to the new partition {R 1 , . . . , R N −1 }. It follows that the initial and final states corresponding to the merged interval R j = R j ∪ R j+1 are given by X In general, initial density functions are not piecewise constant over the partition R. Let f ∈ L 2 ⊃ H (R Q ), P N Q : L 2 → H (R Q ) be the orthogonal projector operator and p(i), and their images Proof The Frobenius-Perron operator associated with S is given by It follows that Then, Hence, Alternatively, (27) can be written as where is the Frobenius-Perron matrix that corresponds to a unique piecewise linear and expanding transformation S given by

Numerical example 1
The applicability of the proposed algorithm is demonstrated using numerical simulation. Consider the following piecewise-linear and expanding transformation S The graph of S is shown in Fig. 2.
A set of initial states X 0 = {x 0, j } θ j=1 , θ = 5 × 10 3 , generated by sampling from a uniform probability density function f 0 (x) = χ [0,1] (x), were iterated using S to generate a corresponding set of final states X T = {x T, j } θ j=1 where T = 20,000. The data set X T was used to determine the uniform partition with N intervals, 1 ≤ N ≤ θ/ log θ = 587, which maximizes the penalized log-likelihood function in equation (10). In this example, N = 10; i.e. C = {0.1, . . . , 0.9} and the estimated invariant density f * C (x) with respect to the 10-interval partition is shown in Fig. 3. The is obtained for l 7 = 4.72, as shown in Fig. 5. This corresponds to the final Markov partition R =   The corresponding identified mappingŜ is shown in Fig. 7.
The estimated coefficients of the identified piecewise-linear semi-Markov transformationŜ  To evaluate the performance of the reconstruction algorithm, we computed the absolute percentage error for x ∈ X = {0.01, 0.02, . . . , 0.99}. As shown in Fig. 8, the relative error between the identified and original map is less than 2 %. Furthermore, Fig. 9 shows the true invariant density f * associated with S superimposed on the invariant densityf * associated with the identified mapŜ (Fig. 7). In practical situations, measurements are corrupted by noise. Given the process where S : R → R is a measurable transformation and {ζ n } is a sequence of independent random variables with density g, it can be shown [41] that the evolution of densities for this transformation is described by the Markov operatorP :  Furthermore, ifP is constrictive thenP has a unique invariant density f * and the sequence {P n f } is asymptotically stable for every f ∈ D [41].
To study how noise affects the performance of our algorithm, we considered the process where S : [0, 1] → [0, 1] is a measurable transformation that has a unique invariant density f * , {ζ n } is i.i.d. with density g and ε is a known noise level. This leads to an integral operator P ε which has a unique invariant density f * ε [41]. It can be shown that lim ε→0 P ε f − P f = 0 for all f ∈ D and that, for 0 < ε < ε 0 , if lim ε→0 f * ε exists, then the limit is f * . To evaluate the performance of the proposed algorithm in the presence of noise, we assumed ζ ∼ N(0, 1) and reconstructed the map for different values of ε. We computed the mean absolute percentage error (MAPE) between S andŜ.
where {x i } θ δS i=1 = {0.01, . . . , 0.99}, θ δS = 99. The results, summarized in (Table 1), demonstrate that the algorithm is robust with respect to constantly applied stochastic perturbations. Remarkably, the approximation errors remain relatively small even for noise levels that would make it almost impossible to reconstruct the map based on time series data [11,43].

Extension to general nonlinear transformations
The approach to reconstructing piecewise-linear and expanding transformations from densities can be extended to more general nonlinear maps. Ulam [44] conjectured that for one-dimensional systems, the infinite-dimensional Frobenius-Perron operator can be approximated arbitrarily well by a finite-dimensional Markov transformation defined over a uniform partition of the interval of interest. The conjecture was proven by Li [45] who also provided a rigorous numerical algorithm for constructing the finite-dimensional operator when the one-dimensional transformation S is known. Here, the aim is to construct from data a piecewiselinear semi-Markov transformationŜ which approximates the original map S.
The main assumptions are that (a) S : I → I is continuous, I = [a, b]; (b) the Frobenius-Perron operator P S : L 1 → L 1 associated with the transformation has a unique stationary density f * and c) P n S f → f * for every f ∈ D; i.e. the sequence {P n S } is asymptotically stable.
Asymptotic stability of {P n S } has been established for certain classes of piecewise C 2 maps. For example, we have the following result [41].
There is a finite constant ψ such that then {P n S } is asymptotically stable.
By using a change of variables, it is sometimes possible to extend the applicability of the above theorem to more general transformations, such as the logistic map [41], which do not satisfy the restrictive conditions on the derivatives of S.

Identification of the Markov partition
Although, for a nonlinear transformation, the invariant density f * ∈ D is not piecewise constant, the approach used to determine the Markov partition for piecewiselinear transformation in Sect. 3 is also used to determine the optimal partition for the piecewise-linear approximation of the unknown nonlinear map.

Identification of the Frobenius-Perron matrix
For the identified Markov partition R, a tentative Frobenius-Perron matrix can be identified using the approaches described in Sect. 3.
Let the obtained Frobenius-Perron matrix be denoted byM = (m i, j ) 1≤i, j≤N . The indices of the contiguous nonzero entries on the i-the row are denoted by for i = 1, . . . , N , j = 1, . . . , p(i) − 1. This means that m i,r (i,k) > 0 for k = 1, . . . , p(i) such that the solution to the optimization problem satisfies and m i, j = 0 if j = r (i, k), k = 1, . . . , p(i).

Reconstruction of the transformation from the Frobenius-Perron matrix
The method for constructing a piecewise-linear approx-imationŜ(x) over the partition R is augmented to take into account the fact that the underlying transformation is continuous and that on each interval of the partition, S R i is either monotonically increasing or decreasing. The entries of the positive Frobenius-Perron matrix are used to calculate the absolute value of the slope of Let for i = 2, . . . , N and σ (1) = σ (2).
Given that the derivative of S| Q (i) where k = 1, . . . , p(i) − 1 and q (i) p(i) = c i . The piecewise-linear semi-Markov transformation for each subinterval Q (i) j is given bŷ   The construction of the piecewise-linear semi-Markov transformation to approximate the original continuous nonlinear map is depicted in Fig. 10.
A smooth version of the estimated transformation can be obtained by fitting a polynomial smoothing spline.  The extended reconstruction algorithm is demonstrated using the quadratic (logistic) map (Fig. 11) It can be shown that {P n S } associated with this transformation is asymptotically stable [41].
A set of initial states X 0 = {x 0, j } θ j=1 , θ = 5 × 10 3 , generated by sampling from a uniform probability density function f 0 (x) = χ [0,1] (x), were iterated using S to generate a corresponding set of final states X T = where T = 30,000. The data set X T was used to search for an uniform partition with N intervals, 1 ≤ N ≤ θ/ log θ = 587, which maximizes the penalized log-likelihood function in Eq. (12), which in this case corresponds to N = 145. The estimated invariant density f * C (x) with respect to the 145interval partition is shown in Fig. 12. In this example, the longest strictly monotone subsequence L of L = {l j } 144 j=1 , l j = 145|h j+1 − h j | has 52 elements and the minimization of is achieved for l 20 = 0.1560, as shown in Fig. 13. This corresponds to a final Markov partition with 72 intervals. The invariant density on the irregular partition R with 72 intervals is shown in Fig. 14.
To identify the Frobenius-Perron matrix, 100 densities (see Appendix) were randomly sampled to gen- . . , 100, θ = 5 × 10 3 . The initial states X i 0 and their images X i 1 under the transformation S were used to estimate the initial and final density functions on R. Examples of initial and final densities are shown in Fig. 15.
The constructed piecewise-linear semi-Markov transformation with respect to the partition R is shown in Fig. 16. The smoothed map, obtained by fitting a cubic spline (smoothing parameter: 0.999), is shown in Fig. 17. The relative approximation error is shown in Fig. 18.  The estimated invariant density on R, obtained by iterating the smoothed map 20,000 times with the initial states X 0 , is shown in Fig. 19, compared with the true invariant density [41] f * (x) = 1/ π √ x(1 − x) . The performance of the algorithm for different noise levels was also evaluated, and the results are summarized in Table 2. As it can be seen, the approximation error remains relatively low (<5 %) for levels of noise (>50 %) that normally cause severe problems to reconstruction algorithms that use time series data.

Conclusions
This paper has addressed in a systematic manner the problem of inferring one-dimensional chaotic maps based on sequences of probability density functions. Compared with previous solutions to solving the inverse Frobenius-Perron, we have derived sufficient conditions and rigorously demonstrated that the proposed approach can uniquely identify the transfor-mation that describes the underlying chaotic dynamics. Specifically, the reconstructed maps exhibit the same dynamics as the original systems and therefore can be used to carry out stability analysis, determine invariant sets and manipulate the dynamical behaviour of the underlying system of interest. The applicability to the proposed methodology and its performance for different levels of noise was demonstrated using numerical simulations involving a piecewise-linear and expanding transformations as well as a continuous onedimensional nonlinear transformation.
One of the reasons for developing the method presented in this paper was to characterize the heterogeneity of human embryonic stem cell (hESC) cultures and to develop efficient protocols for controlling their differentiation. In essence, sorted sub-populations of stem cells expressing different levels of particular cellsurface markers can over time reconstitute the equilibrium distribution of the parent population [46]. Using flow cytometry, it is possible to follow the evolution of the initial density function of the sorted cell fraction, by sampling and re-plating cells, over a number of days. The sequence of density functions generated in this process can then be used to infer the underlying transformation that governs the process, which can help elucidate the existence of cellular substates [47] that, potentially, correspond to the unstable fixed points of the reconstructed map. We have designed the experiments and started generating data by using flow cytometry to sort out cells with different initial densities. These are re-plated and re-analysed in subsequent days to generate sequences of density functions that are required to solve the generalized inverse problem.
The method could be extended to higher-dimensional maps, but this is not necessarily straightforward. The main limitation is the lack of rigorous theoretical results for two-and higher-dimensional maps. While for one-dimensional maps, we have a complete and elegant theoretical framework, for higher-dimensional maps key results are only available for some special cases. A possible solution is to convert the Ndimensional problem to a 1-D problem [29] and estimate the corresponding F-P matrix using the approach introduced in this paper. The main challenge is solving the inverse Ulam problem, i.e. construct the transformation based on the estimated F-P matrix. As noted in [30], for higher-dimensional systems Ulam's conjecture has been proven only for some special cases [48][49][50][51]. The method we are interested to explore to construct the transformation which approximates the original high-dimensional map is that introduced by Bollt [30].