{'fraction_non_alphanumeric': 0.05028984789179012, 'fraction_numerical': 0.026358816826041718, 'mean_word_length': 3.9121440739839377, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A High Accuracy Stochastic Estimation of a Nonlinear Deterministic Model

In this paper, an approach to estimating a nonlinear deterministic model is presented. We introduce a stochastic model with extremely small variances so that the deterministic and stochastic models are essentially indistinguishable from each other. This point is explained in the paper. The estimation is then carried out using stochastic optimisation based on Markov chain Monte Carlo (MCMC) methods.


Background
Given a time series coming from a physical process of interest, the researcher is concerned with inferring unknown quantities, thus modeling the unknown dynamical system in such a way that both systems, the physical and the proposed one, share the same qualitative features.
In recent years, it has been recognized that purely deterministic systems, even with a small number of degrees of freedom, can lead to complicated and irregular behavior. For example, types of nonlinear-chaotic systems have been used in modeling a variety of physical phenomena. Cases include: a simple model for turbulence (Ruelle and Takens, 1971), understanding the behavior of plasmas (Lichtenberg and Lieberman, 1983), the formation of prices in the stock market (Day, 1994 and1997). These models can even incorporate dynamic noise as model error resulting in what is known as random dynamical systems (Lasota and Mackey, 1994).
In dynamical systems theory, when a functional form for the model has been chosen, accurate estimation of quantities such as the control parameters and initial conditions is crucial. For example, when stochastic components are present, the functional form of the deterministic part is polynomial and the information for the system is given in the form of a time series, a number of methodologies using Bayesian, parametric or nonparametric, probability models have been proposed (Meyer and Christensen, 2000and 2001, Hatjispyros et al. 2007and 2009. In this paper, we develop the methodology for obtaining high accuracy estimates of the control parameter and the initial condition when the proposed system is a discrete time, onedimensional nonlinear deterministic dynamical system. The time series information is given in the form of a highly censored sequence of observations, or equivalently a kneading sequence (Metropolis et al. 1973, Wang andKazarinoff, 1987).
In Wu et al. (2004), an alternative algorithm is presented for the estimation of the control parameter and the initial condition when a two-valued symbolic sequence generated by iterations of the quadratic map is given. This algorithm utilizes the relationship between the Gray Ordering Numbers and one-dimensional quadratic maps introduced in Alvarez et al. (1998).
Due to the use of MCMC methods, which may be not familiar to the readers, we provide a brief explanation here.

MCMC and the Gibbs sampler
The approach to estimation of the dynamical system will rely on stochastic models and the sampling of these models. These samples are then used for estimation purposes. Markov chain Monte Carlo methods (Hastings 1970;Smith and Roberts, 1993;Tierney, 1994) are now routinely used in statistical inference problems, notably for Bayesian inference, when the models are complicated and direct estimation via numerical methods is unfeasible.
Suppose the density function f X,Y (x, y) is to be sampled and direct sampling is not possible. A Markov chain is set up with stationary density f X,Y (x, y). Hence, as the chain proceeds, the samples taken, say (x (i) , y (i) ), converge to samples from this density and given a large number of these samples, say N G , estimates of desired functionals can be computed. For example, the marginal mean of x can be computed as where x (i) is the output from the chain for the x values.
A simple chain with the correct stationary density can be constructed in the following way. Start with a fixed point x (0) and sequentially sample, for i = 1, 2, . . ., from the conditional densities f Y |X (y|x (i−1) ), to get y (i) , and then f X|Y (x|y (i) ), to get x (i) . See Smith and Roberts (1993) for the details of the theory for such a chain, known specifically as the Gibbs sampler.
To remove the bias of the starting point x (0) , the first N b samples can be ignored and so inference relies only on the output Of course this sampler relies on the availability of the two conditional densities to be sampled directly. Often this can be the case. When it is not then a recourse is to introduce a latent variable, z say, (Damien et al., 1999) which does not change the marginal density of f X,Y (x, y).
and an augmented Gibbs sampler operates now on (x, y, z). This works by sequentially sampling the conditional densities, f X|Y,Z (x|y, z), f Y |X,Z (y|x, z) and f Z|X,Y (z|x, y). The (x, y) output can still be used for inference purposes.

Layout of paper
The layout is as follows. In Section 2 we provide the probability model for the observations which contain the latent variables. The joint probability model therefore contain the parameters to be estimated and the unobserved variables. In Section 2 we also describe the full conditional densities required to implement the Gibbs sampler. Details of how to sample these and latent variables useful for the implementation are provided in the Appendix.
Section 3 then elaborates on the MCMC method which necessarily needs to be non-standard. This is because the joint density we are actually sampling is highly modal and with sharp and highly peaked modes. In the end we are searching the largest mode and for estimation purposes only, so we can "zoom" in on the correct mode. Therefore the actual algorithm involves a number of steps, which are described in Section 3. Further elaboration is provided in Section 4 where we undertake an illustrative example involving the quadratic map.

The probability model
Suppose we have a one dimensional nonlinear dynamical model g : X → X, whereby for some initial condition y 0 , the trajectory y * i is formed from y * i−1 via the deterministic recurrence relation y * i = g(ϑ, y * i−1 ), i = 1, . . . , n. The map g is continuous in y * , ϑ ∈ Θ is the parametrization of the map and X is invariant under g(ϑ, ·) for all values of ϑ ∈ Θ. We choose g to be quadratic in y i−1 .
The underlying data y * = (y * i ) 1≤i≤n , are censored. So if y * i ∈ A we observe b i = 0 whereas if y * i / ∈ A then we observe b i = 1. Given the binary symbolic data sequence of b = (b i ) 1≤i≤n our aim is to estimate (ϑ, y 0 ).
To permit a stochastic approach to estimation we adopt a stochastic model for the unobserved sequence y * . So let the stochastic counterpart of y * i be written as y i . Therefore, for some specified sequence of variances (σ 2 i ), we assume that where the (ǫ i ) are independent and identically distributed as standard normal random variables. The stochastic model can be made arbitrarily close to the deterministic model by allowing the variances σ 2 i to collapse to zero. In practice, due to the nature of the data (b i ), we only need to choose the variances small enough to ensure that the stochastic model gives rise to the same data with very high probability.
Here we start to describe some of the joint densities of interest that we would wish to sample in order to estimate the model. These are intractable from a numerical point of view but we will resolve this by using the Gibbs sampler. The random variables y i , given (ϑ, y i−1 , σ i ) are conditionally independent normals with means g(ϑ, y i−1 ) and variances σ 2 i ; our notation for this will be . . , n, and the joint density of (y 1 , . . . , y n ) conditional on the values y 0 , ϑ and (σ i ) is The conditional probability of D i given y i is 1 whenever y i ∈ D i and 0 elsewhere. We express that by using the characteristic function notation Then the joint probability model becomes π(y 1 , . . . , y n , D 1 , . . . , D n | y 0 , ϑ, One way to make this model accurate is to take the specific sequence of standard deviation constant (σ i = σ) and to make σ very small. This follows since the (y i ) are deterministically generated, albeit with rounding errors, which can be used to justify the inclusion of some small σ.
However, in this case, the likelihood function is not only multi-modal, but incredibly spiked at the modes. This renders direct estimation of (ϑ, y 0 ) very difficult, if not impossible. This is essentially caused by the fact that the combination of: 1. The noninformative prior specification of ϑ and y 0 for model map g, namely here formulated in terms of prior distributions, where U(·|A) denotes the uniform distribution over an appropriate subset A of R.
2. The available sample of length n of censored observations of the form makes inefficient a direct (ϑ, y 0 |D)-inference via Markov chain Monte Carlo (MCMC) methods. Therefore we need to think about smart MCMC methods.
Definition 2. The estimating strength (ES) of a (g, θ, D)-candidate pointȳ i for i = 0, . . . , n is the number The cumulative estimating strength (CES) of a n-vector candidateȳ is the sum Our strategy will be to maximize the CES function S(θ|D) over Θ, via MCMC methods. In turn, the approximation of a global maximum, if it exists, of the CES function using a sample of censored observations of length n, over a Θ-grid of length δ, and the use of the inverse of the map g over D, will provide us with approximations (θ n,δ ,ȳ n,δ 0 ) to the true values (ϑ, y 0 ) responsible for the sequence of censored observations D. In other words, we will be able to specify numerically a subset Θ n,δ of the parameter space, containing the true value ϑ, along with an approximation to the (g, θ n,δ , D)-candidate vectorȳ n,δ . The latter will enable us to substitute the prior information for θ with π(θ) = U(θ|Θ n,δ ) and the initial data sets (D i ) with the refined sets D n,δ,ǫ for appropriate small ǫ, thus obtaining an efficient θ-augmented MCMC scheme that improves the accuracy of the previous estimates (θ n,δ ,ȳ n,δ 0 ).

The strength maximizing MCMC
We now describe the MCMC algorithm for the estimation of the truncation set Θ n,δ and the candidate approximationỹ n,δ .
The full conditional distributions for y 0 , . . . , y n−1 are truncated versions of nonstandard densities. The full conditional for y n is a D n -truncated normal distribution with mean g(θ, y n−1 ) and variance σ 2 . In all cases to sample, we use a slice sampling technique. More details on the particular sampling scheme are given in the appendix.

Polishing the MCMC
Having computed the truncating set Θ n,δ , and, the more informative, data sets D n,δ,ǫ i for i = 0, 1, . . . , n we describe the MCMC scheme for polishing the (θ n,δ ,ȳ n,δ 0 ) approximations to higher accuracy.
2. The full conditionals for y i , 0 ≤ i ≤ n, are π(y i | · · ·) ∝ I y i ∈ D n,δ,ǫ It is not difficult to verify that when g(θ, y) is affine in θ, say g(θ, y) = α(y) + β(y) θ, the full conditional for θ is a normal truncated on the interval Θ n,δ , with mean In the next section we will provide an application of the augmented MCMC algorithm for the quadratic map.

Grid maximization results
We begin by maximizing the CES function S(θ|D). We set σ = 10 −3 , this value of σ is to balance the noninformative prior specifications in (10) and the small ranges of the posterior marginal distributions for (y i ) as smaller σ's will cause the corresponding MCMC to converge extremely slow. At this point we do not have to use the full length sample, a smaller sample D ′ of length m = 600 is enough to estimate both θ and y 0 to 4 decimal places and to obtain a Θ m,δ truncating set of length δ ≈ 3 × 10 −3 . In our simulated example we consider the Θ-grid sequence in such way that the true value ϑ is not a grid point at any level. Initially we consider the grid then for each value θ ∈ Θ 1 we compute iteratively via the strength maximizing MCMC scheme the vector of running averagesȳ (j) and the associated CES function S (j) (θ r | D ′ ). After a burn-in period of about N b = 4 × 10 4 iterations and N G = 2 × 10 5 in total Gibbs iterations we compute an approximation to the CES function as The result is shown in Fig. 1(a). The CES function has a maximum over Θ 1 at θ 6,1 = −1.7250 which is the center of the interval (−1.77, −1.68). The associated value for the estimate of the initial condition isȳ 0 = 0.80256. This gives us the next refinement grid Θ 2 = {θ 0,2 = −1.6800, θ s,2 = θ 0,2 − 0.0082 s, 1 ≤ s ≤ 11}.
Reiterating the whole process, we find a maximum of the S (N G ) (θ| D ′ ) function over Θ 2 at θ 5,2 = −1.71273 andȳ 0 = 0.80073, as shown in Fig. 1(b), which is the center of the interval (−1.72091, −1.70455). Our final refinement with the D ′ data set is where we find a maximum at θ 5,3 = −1.70995 andȳ 0 = 0.80026, see Fig. 1(c). Both approximations agree to the true values to 4 decimal places (see table 1). As a θ-prior information to the polishing MCMC we are using the interval (θ 6,3 , θ 4,3 ) = (−1.71132, −1.70859). A plot of the sequence of strengths at the maximum θ 5,3 for D ′ can be found in Fig. 2. The sequence is highly oscillatory and a number of maximums are descernible. The worst estimatedȳ i values are on the tail of the sequence.

Increasing the accuracy of theȳ i estimates
In the sequel we are making use of the full sample D of size n = 1000. We compute the sequence of strengths at the closest maximum θ * to the true value ϑ -in our example is θ 5,3 -and we pickȳ κ as the candidate point that is closer to the end of the sequence and has a strength value greater than 20. We discard estimates of order higher than κ, and we apply g −1 | D (θ * , ·), the inverse of g over D at θ * , to the remaining candidate points. For example, for the quadratic map in (8), we havẽ and in our example κ = 984. We plot the sequence (S i (θ * |D)) 0≤i≤80 of the first 80 strengths, including the strength of y 0 , in Fig. 3(a). The absolute error of the estimatesȳ i at θ * for 0 ≤ i ≤ 80, is in Fig. 3(b). Finally, in Fig. 3(c), we plot the absolute error of the estimatesỹ i , obtained by (14). We can see that the order of the errors after the application of g −1 | D to the candidate points (ȳ i ) are considerably reduced from order 10 −3 to order 10 −5 (see table 2

4.3
Obtaining the high accuracy estimates of (θ, y 0 ) We first use a 100 length sub-sample of D with ǫ = 5 × 10 −5 . We iterate the augmented MCMC using as an initial condition the vector (θ * ,ỹ 0 ,ỹ 1 , . . . ,ỹ 100 ) and σ = 10 −8 . The approximations to the mean of the θ and y 0 posteriors are shown in Fig. 4(a) and 4(b) respectively. We can see that after a large number of iterations are still oscillating in the intervals θ ∈ [−1.7100024, −1.7100013] and y 0 ∈ [0.7999952, 0.7999951] and are correct to 5 decimal places.
Finally we increase the sample size to κ = 984. The posterior ergodic averages for θ and y 0 are shown in Fig. 5(a) and 5(b) respectively. Here we approximate the true values to 7 decimal places. The values are: θ = −1.70999996 and y 0 = 0.79999995.

Discussion
We have described an algorithm for the estimation of the parameters of a one-dimensional quadratic nonlinear dynamical system. The stochastic model provides us with an objective function to be maximized. However, the function is multimodal with sharp spikes and hence standard application of stochastic search procedures based on MCMC methods will fail. We therefore introduce a zooming or polishing process which allows us to target the largest mode. Estimating maximum of functions using MCMC methods is by now a standard practice and dates back to the work of Geyer and Thompson (1992) and Douset et al. (2002) for a Bayesian setting. The maximum can be obtained to an arbitrarily level of accuracy.
Our proposed algorithm besides its practical use on estimating (θ, x 0 ) from a set of incomplete measurements can be applied to a chaotic cryptography setting Baptista 1998, Alvarez et al 2003 andreferences therein. In such systems (θ, x 0 ) is used as a shared secret key, from both the encoder of the plain text and the decoder of the cipher text. Reusing the same secret key to encrypt messages weakens the security of the cryptosystem. We can think of a chaotic cryptosystem transmitting together with the coded message a long m-valued symbolic sequence which is the projection of a (θ, x 0 )-orbit on an uneven m-partition of the invariant set of the map, being a refinement of some initial generating partition. The new partition then becomes a shared secret key for the recovery of the control parameters and initial condition from the decoder.
Our method in principle can be applied for estimating more complex nonlinear recurrences. Although the complexity of the algorithm will increase, the associated computational workload will not increase dramatically, for example: 1. When g(θ, y) has degree m > 2 in y. Then there are m − 1 critical points in the invariant set X of g and each data set D i will be one of the m sets of the partition of X induced by the critical points of the map. Such a partition is generating, so that there is a one to one correspondence between m-valued symbolic sequences and dynamical orbits. The only fact that adds to the complexity of our algorithm is that sampling from the full conditional distributions for y i in (7) will now involve solving numerically inequalities of degree m of the form where u i2 is an auxiliary random variable, see appendix equation (15).
2. Higher dimensional dynamical systems of the form y i = g(θ, y i−1 , . . . , y i−d ). In this case there is no systematic way for the construction of a generating partition Eisele (1999). Nevertheless it is possible to design a partition that maximizes the entropy content of the associated symbolic sequences Strelioff and Cruchfield (2007), Ebeling and Nicolis (1991). Then the conditional distributions for y i for i = d, . . . , n − d are given by Sampling from such full conditionals will involve embedded MCMC schemes with d + 1 auxiliary random variables.
It is clear that the y i marginal of the latter density is the density we want to sample from. We now describe the associated Gibbs sampler.
Hence, u i1 and u i2 are both exponential variables, with mean 1/λ, truncated in the intervals (h θ (y i+j−2 , y i+j−1 ), ∞) for j = 1, 2, and thus they can be sampled as where q ij ∼ U(0, 1) independently. The full conditional for y i is which is in general a mixture of uniforms. For example when g(θ, y) is a polynomial of degree m, in y, the full conditional for y i is a mixture of at most m uniforms. In the special case when g(θ, y) = 1 + θy 2 we have where we finally obtain the uniform mixture    Fig. 3(a) we plot the sequence S i (θ * |D) of individual strengths for 0 ≤ i ≤ 80 based on a data set of length 600. In Fig. 3(b) we plot the absolute error of the estimatesȳ i at θ * . The sequence of strengths and absolute errors are negatively correlated. Finally in Fig.  3(c) we plot the absolute error of the estimatesỹ i , obtained with the help of (14).  Figure 4: Here the sample size is 100 with ǫ = 5 × 10 −5 and σ = 10 −8 . The initial conditions for the MCMC scheme in (6)-(7) is the vector (θ * ,ỹ 0 ,ỹ 1 , . . . ,ỹ 100 ) withỹ i 's coming from (14). The approximations of the θ and y 0 posterior means -shown in Fig. 4(a) and 4(b) respectively -are converging to θ = −1.7100024 and y 0 = 0.7999952 and are correct to 5 decimal places.  Figure 5: Here the sample size is κ = 984 with ǫ = 5×10 −5 and σ = 10 −8 . The initial conditions for the MCMC scheme in (6)-(7) is the vector (θ * ,ỹ 0 ,ỹ 1 , . . . ,ỹ κ ) withỹ i 's coming from (14). The approximations of the θ and y 0 posterior means -shown in Fig. 5(a) and 5(b) respectively -are converging to θ = −1.70999996 and y 0 = 0.79999995 and are correct to 7 decimal places.  Table 2: The true values y * i = g (i) (ϑ, y * 0 ) for 0 ≤ i ≤ 4, ϑ = −1.71 and y 0 = 0.8 are on the 2nd column. The approximations to the candidate valuesȳ i are on the 3rd column. After the application of g −1 | D using (14) with κ = 984, we obtainỹ i on the 4th column. The absolute errors of the previous two approximations are in the last two columns.