Predicting and estimating probability density functions of chaotic systems

In the present work, for the first time, we employ Ulam's method to estimate and to predict the existence of the probability density functions of single species populations with chaotic dynamics. In particular, given a chaotic map, we show that Ulam's method generates a sequence of density functions in L 1 -space that may converge weakly to a function in L 1 -space. In such a case we show that the limiting function generates an absolutely continuous (w.r.t. the Lebesgue measure) invariant measure (w.r.t. the given chaotic map) and therefore the limiting function is the probability density function of the chaotic map. This fact can be used to determine the existence and estimate the probability density functions of chaotic biological systems.


1.
Introduction. Several ecologists have studied the implications of chaotic dynamics both in field experiments [11,12,13,16,25,39,40] and experimental population research [7,14,20,22,36,45]. In wildlife, there are only a few proofs of chaos that seem convincing. This includes chaotic population dynamics of boreal rodents due to mustelid predators [15], dynamics of vole populations in northern Fennoscandia [17,32,41] and the epidemics of childhood diseases [34,33]. Although laboratory experiments are far different than the wildlife situations, they can be used for testing some ecological concepts. There are more convincing examples of chaotic behaviors in laboratory experiments. In particular, ecologists have observed experimental evidence of chaos in cultures of the flour beetle Tribolium [8,10], Nicholson's sheep blowflies [29,35,43] and the microbial food web with constant experimental conditions [3].
There are several mathematical and statistical models that can be used to explore chaotic behaviors. In the present work we consider discrete-time models of the form where f (x) is the density-dependent recruitment function.
In biology, the recruitment functions are often classified into overcompensating and compensating density dependent functions [2]. The former occurs when f has a local maximum and decreases for large population sizes, whereas the latter corresponds to a recruitment function that is bounded and increases monotonically. Compensating density dependence relates to contest competition, in which only a portion of the population obtains enough resources for survival and the rest will gradually die [5]. Whereas overcompensating density relates to scramble competition, in which resources are divided between all members and large population sizes will result in decreased survival rates [2,5].
It can be shown that there is no possibility of period doubling or chaotic behavior for populations with compensating density dependence. In particular, when the recruitment function f (x) is a bounded monotone increasing function, the slope of f (x) at the positive equilibrium is between zero and one. Therefore, the positive equilibrium remains stable at all times. On the contrary, populations that are ruled by overcompensating density dependence may exhibit period doubling and chaotic behaviors (see for example chapter 2 of [5]). When the model is chaotic, the solution x n behaves randomly and it is highly sensitive to the initial condition. Studies of discrete chaotic dynamical systems often lead to bifurcation diagrams, frequency plots and many other computational aspects of the chaotic systems [1,7,10,38]. In the present work, for the first time, we employ a mathematical method known as Ulam's method [42] to estimate the probability density function (PDF) of single species populations that exhibit chaotic behaviors. We will also establish the conditions, under which the PDF related to a chaotic system exists. Using the estimated PDF the ecologist will able to further investigate the behavior of the chaotic populations and compute the estimated probability of any interval of population sizes.
Mathematically, we study the existence of invariant measures of chaotic maps. The existence of a continuous invariant measure shows that the system in question is chaotic. The existence of an absolutely continuous invariant measure (acim) proves even more chaotic behavior and gives tools to statistically predict the future of the system.
The most popular method to practically approximate an acim is the so-called Ulam's method of approximating the Frobenius-Perron operator of the system by a sequence of finite-dimensional operators (matrices), [24]. We prove some facts about Ulam's approximations. In particular, we show that any *-weak limit of Ulam's approximations is an invariant measure, and give sufficient conditions for this limit measure to be absolutely continuous.
As an example, we show that the Smith-Slatkin model for some parameters has an acim. We use the Schwarzian derivative method and sketch the methodology which can give the existence of an acim for a large set of parameter values.

2.
Preliminaries. Throughout the paper we assume that [a, b] is a closed and bounded interval in the real number line, B is the Borel σ-algebra on [a, b], and µ is a measure on B. We recall some basic definitions. More details can be found in [4].
Definition 3. We say that the measurable map τ : Now we are ready to recall the Birkhoff Ergodic Theorem (see [4] for the proof).

The Birkhoff Ergodic Theorem
Furthermore, we havef In addition, if τ is ergodic w.r.t. the measure µ, then (3) implies thatf is constant µ-a.e., so from (4) we obtain So if µ([a, b]) = 1, that is, the measure µ is a probability measure, (2) becomes Thus, by letting f = χ B in (5) we obtain Here χ B denotes the characteristic function over the set B, that is, Remark 1. Let's assume the measure µ defined by where f * ∈ L 1 [a, b]. Then we say that the measure µ has density f * . Note that the measure µ defined by (7) is absolute continuous w.r.t. the Lebesgue measure m. So if τ is ergodic w.r.t. the Lebesgue measure m, then it is also ergodic w.r.t. the measure µ. Now we further assume that the measure µ defined by (7) is τ -invariant probability measure. If τ is also ergodic w.r.t. the Lebesgue measure µ, from (6) and (7) we obtain Note that if we regard (x, τ (x), τ 2 (x), . . .) as the occurrences of a random variable X at each discrete step, then (8) can be written So f * is the probability density function of the random variable X. In section 5, we will find conditions under which there exists a density function f * ∈ L 1 [a, b] such that the probability measure µ with density f * is τ -invariant using the concept of weak precompactness.
We also need to review the Frobenius-Perron operator [4] associated with a map τ .
is called the Frobenius-Perron operator associated with τ .

Remark 2. By choosing
So by differentiating both sides w.r.t. x we have which is the explicit formula of the Frobenius-Perron operator associated with τ . Now we have the following theorem that tells us under what condition on the density function f * ∈ L 1 [a, b], the probability measure µ with density f * is τinvariant.
] be a nonsingular map, and let P τ be the Frobenius-Perron operator associated with τ . Let f * ∈ L 1 [a, b] be a density function. Then the probability measure µ with density f * is τ -invariant if and only if f * is a fixed point of P τ , that is, P τ (f * ) = f * .
3. Partially invariant measures. In this section we show that Ulam's method gives an almost τ -invariant measure. We briefly explain Ulam's method. We let P n = {I i } n i=1 be a partition of [a, b] into n equal length subintervals. We define an n × n nonnegative matrix A = (a ij ), which is called Ulam's matrix associated with τ based on the partition P n of [a, b], as follows: Note that Ulam's matrix is a column stochastic matrix. So A can be viewed as the transition matrix of a finite state Markov chain so that the a ij can be interpreted as the probability that a point in I j is mapped into I i under τ . Since Ulam's matrix is a column stochastic matrix, 1 is its eignenvalue. Since A is nonnegative, there exists a nonnegative eigenvector y of A corresponding to the eigenvalue 1, that is, Ay = y and y ≥ 0. We may assume that n j=1 y j = n. We let f n = n j=1 y j χ Ij . Since y is a nonnegative vector such that n j=1 y j = n/(b − a), it follows that f n is a density function. We call f n Ulam's density function associated with τ based on the partition P n of [a, b]. Now we prove that if µ n is the measure with density f n , then µ n is τ -invariant on P n . Lemma 1. If we let µ n be the measure with density f n , where f n is Ulam's density function associated with τ based on the partition P n = {I i } n i=1 , then we have Proof. Let A = (a ij ) be Ulam's matrix associated with τ based on the partition P n , so that Recall that Ay = y and f n = n j=1 y j χ Ij . We first compute µ n (τ −1 (I i )). Observe that On the other hand, note that So we have proved that 4. Invariant measures. In this section we will prove that if a measure µ is a limit point of the measures {µ 2 n } n≥1 in * -weak topology, then µ is τ -invariant. We recall that the convergence of µ 2 n to µ in *-weak topology means that (See, for example, [4, Chapter 2, Definition 2.2.8]). We first prove the following lemma.
On the other hand, using Lemma 1 we also have So, we have shown that b a g dµ 2 n = b a g • τ dµ 2 n , provided that g is a step function on P 2 n = {I i } 2 n i=1 . We will need the following result in proving Proposition 1.
Theorem 2. Let {ν n } n≥1 be a sequence of probability measures. Then there exists a subsequence {ν n k } k≥1 such that ν n k → ν in *-weak topology as k → ∞ and ν is a probability measure, that is, the set of probability measures is compact in *-weak topology.
We also need the following result.
Let us consider a sequence of Ulam's densities {f 2 n } n≥1 associated with τ based on the sequence of partitions {P 2 n } n≥1 . Let µ 2 n be probability measure with density f 2 n . Since {µ 2 n } n≥1 is a family of probability measures, by Theorem 2 there exists a subsequence {µ 2 n k } k≥1 such that µ n k → µ as k → ∞ in *-weak topology, where µ is a probability measure. Then µ is τ -invariant.
Proof. For notational convenience, we will write {µ 2 n k } k≥1 as {µ 2 n } n≥1 . By The- where each I i ∈ P 2 n and c i is the infimum value of g on I i . Since g ∈ C[a, b], g 2 n → g uniformly on the interval [a, b]. In particular, we can find an N 1 ≥ 1 such that

This implies also
Since µ 2 n → µ in *-weak topology and g, τ ∈ C[a, b], we can find N 2 and N 3 such that Let n ≥ N . Then using Lemma 2 and the fact that µ 2 n is a probability measure, we have Since > 0 was arbitrary, we have The following example shows that even though the measure µ in Proposition 1 is τ -invariant, it may not be absolutely continuous measure w.r.t. the Lebesgue measure m.
It is easy to show that f 2 n = 2 n χ [0,1/2 n ] . So for any g ∈ C[0, 1], we have If we define µ as then clearly 1 0 g dµ = g(0). So we see that µ n → µ as n → ∞ in *-weak topology. So by Proposition 1 µ is τ -invariant, but clearly µ is not absolutely continuous with respect to the Lebesgue measure m.

Absolutely continuous invariant measures.
In this section we consider constructing τ -invariant measures which are also absolutely continous measures w.r.t. the Lebesgue measure m.
We now prove the following: ] be nonsingular. Let us consider a sequence of Ulam's densities {f 2 n } n≥1 associated with τ base on the sequence of partitions {P 2 n } n≥1 . Let µ 2 n be the measure with density f 2 n for each n. If the family {f 2 n } n≥1 is weakly precompact in L 1 [a, b] and hence f 2 n k → f * weakly in L 1 [a, b] as k → ∞, then the measure µ with density f * is an absolutely continuous τ -invariant probability measure.
Proof. For notational convenience we may assume that f 2 n → f * weakly in L 1 [a, b] as n → ∞. We first show that µ is τ -invariant on ∪ n≥1 P 2 n . Suppose I ∈ ∪ n≥1 P 2 n . Then I ∈ P 2 N 1 for some N 1 . By Remark 4 we have µ 2 n (I) = µ 2 n (τ −1 (I)) for all n ≥ N 1 , that is, Let > 0 be given. Since f 2 n converges to f * weakly and χ I ∈ L ∞ [a, b], there exists Since > 0 was arbitrary, we have µ(I) = µ(τ −1 (I)). Since I was an arbitrary element in ∪ n≥1 P 2 n , µ is τ -invariant on ∪ n≥1 P 2 n . Now we show that µ is τ -invariant on B. Since τ is nonsingular, both µ and µ • τ −1 are absolutely continuous measures w.r.t. the Lebesgue measure m. So for any > 0 there exists δ > 0 such that Now suppose B ∈ B. Then there is a finite disjoint collection of finite intervals where I = m k=1 I k . Since the collection of the end points of the intervals in {P 2 n } n≥1 is dense in the interval [a, b], for each k we can choose an interval J k such that (a) J k ⊆ I k ; (b) J k is a finite disjoint union of the intervals in P 2 n k for some n k ; (c) m(I k − J k ) < δ/m.
Since µ = µ • τ −1 on {P 2 n } n≥1 , it follows that µ(J k ) = (µ • τ −1 )(J k ) for 1 ≤ k ≤ m, and hence it follows that So using (9), (11), and (12), we have  So, using (9), (13), and (14), we have Finally, using (9), (14), and (15), we have Since was arbitrary, we have µ(B) = (µ • τ −1 )(B). Since B was an arbitrary element of B, µ is τ -invariant on B. Since µ has its density function f * ∈ L 1 [a, b], it is also an absolutely continuous probability measure w.r.t. the Lebesgue measure.  1 3 ≤ x ≤ 1. Using Ulam's method, we found {f 2 n } 10 n=1 and obtained the data displayed in Table  1. It shows strong evidence that {f 2 n } n≥1 is bounded by 2. So, by Theorem 4, {f 2 n } n≥1 is weakly pre-compact, and hence there exists a sub-sequence {f 2 n k } n≥1 of {f 2 n } n≥1 which converges weakly to f * ∈ L 1 [a, b]. Thus by Proposition 2 the measure µ with density f * is an absolutely continuous τ -invariant probability measure µ. Then by Theorem 1, we have P τ (f * ) = f * . Since τ is ergodic, f * is the unique fixed density of P τ . Hence again by Theorem 1 we see that µ is the unique absolutely continuous τ -invariant probability measure on B. One can verify that Since f * is the unique fixed density of P τ , we may expect that {f 2 n } n≥1 converges strongly to f * (see [24]). See Figure 1 for the plot of f 2 5 and f * .  Figure 1. Graphs of the probability density function f * and its estimation f 2 5 for Example 2. Since {f 2 n } n≥1 is bounded by 2, by Theorem 4 {f 2 n } n≥1 is weakly pre-compact, and hence there exists a subsequence {f 2 n k } n≥1 of {f 2 n } n≥1 which converges weakly to f * ∈ L 1 [a, b]. In fact, we may expect that {f 2 n } n≥1 converges strongly to f * . n=5 . Each f n has two peak values. Those peak values happened near 0 and near 1 for each n. We compiled those peak values and displayed in Table 2. We also presented their ratios of two consecutive f n 's, and displayed them in Table 3, which show that the limit behavior of {f 2 n } n≥1 is like O(1/ √ x) near 0 and O(1/ √ 1 − x) near 1. Since both 1/ √ x and 1/ √ 1 − x are in L 1 [0, 1], the collection {f 2 n } n≥1 is bounded by some g ∈ L 1 [0, 1]. So, by Theorem 4, {f 2 n } n≥1 is weakly precompact, and hence by Proposition 2 there is a density f * ∈ L 1 [a, b] such that if the measure µ has f * as its density, then the measure µ is the absolutely continuous τ -invariant probability measure. Since τ is ergodic, f * the only density function which makes µ τ -invariant. In fact, it is well known that Again we may expect that {f 2 n } n≥1 converges strongly to f * . See Figure 2 for the plot of f 2 9 and f * .    Figure 2. Graphs of the probability density function f * and its estimation f 2 9 for Example 3. Again we may expect that {f 2 n } n≥1 converges strongly to f * .
hence Theorem 4 and Proposition 2 do not apply. In fact, in Example 1 we showed that there is no absolutely continuous invariant measure in this case.
6. Applications to biological systems. The recursive equation can be used for nonlinear discrete-time epidemic models. Here γ (0 ≤ γ < 1) denotes the constant probability of surviving per generation; f models the typically nonlinear birth or recruitment process; x(n) denotes the population size at time n; and time is measured in discrete units.
If f (x n ) = αx n 1 + (βx n ) p , the above recursive equation is called the Smith-Slatkin model. Hence the Smith-Slatkin model gives Here α is the maximal per-capita intrinsic growth rate of the population; p reflects the type and strength of intra-specific competition; and β scales the carrying capacity of the population. Details of this model and its applications in ecology and epidemiology have been mentioned in [6].
We write the above recursive equation as If we let x * denote the fixed population of the Smith-Slatkin model under τ , from the equation we obtain The following lemma establishes the condition under which |τ (x * )| > 1.
The critical point of τ to the left of x * is We will study the dynamics of τ on the invariant interval a = τ 2 (x c ), b = τ (x c ) , see Figure 4. We will use the Schwarzian derivative of τ defined by Its negativity can be used to prove the existence of absolutely continuous invariant measure [27,28,31,21,30]. First, we will prove that Sτ < 0 for a wide range of parameters.
A strong result giving sufficient condition for the existence of acim for such maps was given in [31]. We give a simplified version applicable to our map τ : Theorem 5. Suppose that τ is unimodal, C 3 and has a negative Schwarzian derivative. Moreover assume that ∞ n=0 |(τ n+1 ) (x c ))| −1/2 < +∞.
Then, τ has a unique absolutely continuous invariant probability measure µ which is ergodic and of positive entropy. Furthermore there exists a positive constant K such that µ(A) < K m(A) for any A ∈ B.
As it is shown in Figures 6 and 7, when α = 1.5, β = 1, and γ = 0.1, for p = 8 the third image of the critical point τ 3 (x c ) is above the fixed point x * and for p = 8.5 we have τ 3 (x c ) < x * . Thus, for some value p =p in between we have τ 3 (x c ) = x * . For α = 1.5 and γ = 0.1, the in Proposition 3 is fulfilled for all p > 2.
The famous theorem of Jakobson [27,18,19,26,37,44] states that for a one parameter family of unimodal maps τ α with negative Schwarzian derivative, similar to the family of parabolas x → αx(1 − x), the set of parameters for which τ α has an acim is "large" (e.g., it is of positive Lebesgue measure) and in particular it is of Lebesgue density 1 around each "Misiurewicz parameter" (a parameter for which some iterate of critical point hits a periodic point).
Unfortunately, Smith-Slatkin family (say, with α, β, γ fixed and p as a parameter) does not satisfy the assumptions of any version of Jakobson's theorem we know. On the other hand the bifurcation diagram of Figure 3 suggests that map τ has an acim for a large set of p's. The following Example reinforces this evidence. Now we want to study whether there is a density f * ∈ L 1 [a, b] such that if the measure µ has f * as its density function, then µ is a τ -invariant absolutely continuous probability measure on B. If we can find such a density f * ∈ L 1 [a, b], then by Remark 1 f * is the p.d.f. of the population random variable X whose occurrences are according to (x, τ (x), τ 2 (x), . . .). So even though the population at discrete times may exhibit a chaotic behavior, we can do statistical analysis of the population random variable X.
Using Ulam's method we computed {f 2 n } 11 n=8 , and we found out that there are six peak values. We compiled those peak values and displayed them in Tables 4-6. We also present their ratios of two consecutive f 2 n and displayed them in Tables 7-9. Note that at each singular pointx the limit behavior is like O(1/ √x − x) when x <x and x ≈x and O(1/ √ x −x) when x >x and x ≈x. Since 1/ √x − x * and 1/ √x − x * are in L 1 [a, b], the collection {f 2 n } n≥1 is bounded by some g ∈ L 1 [a, b]. So by Theorem 5, {f 2 n } n≥1 is weakly precompact and hence by proposition 2 there is a density f * ∈ L 1 [a, b] such that if the measure µ has f * as its density function, then the measure µ is the absolutely continuous τ -invariant probability measure. Since τ is ergodic, µ is the only τ -invariant measure. So we may expect that {f 2 n } n≥1  converges to f * . Figure 8 is the graph of f 2 8 , in which we can recognize six peak values. It is clear that the first peak value occurs at x = a and the last peak value at x = b.  Figure 8. Graph of f 2 8 for Example 5. There are six peak values, where the first peak value occurs at x = a and the last peak value at x = b. Although the closed form of the probability density function f * is unknown, the numerical simulations strongly show that {f 2 n } n≥1 is bounded by an integrable function. So, by Theorem 4, {f 2 n } n≥1 is weakly pre-compact, and hence there exists a sub-sequence {f 2 n k } n≥1 of {f 2 n } n≥1 which converges weakly to f * ∈ L 1 [a, b]. By Proposition 2 the measure µ with density function f * is an absolutely continuous τ -invariant measure. So f * is the PDF of the chaotic map τ . In fact, we may expect that {f 2 n } n≥1 converges strongly to f * . Note that f * is the p.d.f. of the population random variable X whose occurrences are according to (x, τ (x), τ 2 (x), . . .). Using (5) with f (x) = x, we have  Table 7. The consecutive ratio of the first two peak values of f 2 n for τ in Example 5 n 2 n f 2 n /f 2 n−1 at 1st peak value f 2 n /f 2 n−1 at 2nd peak value 9 512 1.4128 1.3938 10 1024 1.4735 1.4055 11 2048 1.3779 1.3203 Table 8. The consecutive ratio of the middle two peak values of f 2 n for τ in Example 5 n 2 n f 2 n /f 2 n−1 at 3rd peak value f 2 n /f 2 n−1 at 4th peak value 9 512 1.3753 1.3630 10 1024 1.3929 1.3668 11 2048 1.3146 1.3067 Table 9. The consecutive ratio of the last two peak values of f 2 n for τ in Example 5 n 2 n f 2 n /f 2 n−1 at 5th peak value f 2 n /f 2 n−1 at 6th peak value 9 512 1.3667 1.3367 10 1024 1.2577 1.4294 11 2048 1.2737 1.3439 Since the measure µ has f * as its density, the above equation becomes Using these two formulas with N = 10 5 , 10 6 , and 10 7 , we compiled the estimates of X and σ in Table 10.