Output error behavior for discretizations of ergodic, chaotic systems of ordinary differential equations

Progress in computational processing capabilities and numerical methods has led to increasing use of unsteady, scale resolving simulation methods for turbulent flows for applied research.^1–3 In spite of advances, the cost of these simulations, particularly large-eddy simulation (LES), remains an obstacle to industrial use. Cost reduction remains a key research area for these systems, including more aggressive modeling approaches.^1,4 On the other hand, advancements in numerical methods have demonstrated that a stable, high-order solution of turbulent systems is possible⁵ with ongoing developments enabling new efficiencies.⁶ As these advances approach maturity, a key outstanding question is how to use them to extract estimates as efficiently as possible when simulations are unsteady and chaotic and discretization methods are not significantly constrained by explicit stability limits.

The resolution of the turbulent—and therefore chaotic—nature of the underlying flows is a primary complication for estimation using these more complex simulations. For chaotic systems, estimation of long-time behavior is challenging, because the ordinary differential equations (ODEs) to which they can be reduced have limited predictability.⁷ Of the general class of chaotic systems, a subset is ergodic systems, whose long-term states are drawn from a stationary distribution, independent of initial conditions (ICs).⁸ For ergodic chaotic problems, we frequently want to quantify the unique infinite-time average of some instantaneous quantity of interest of the system

where g is the instantaneous output functional and the state $u(t)$ is governed by a dynamical system of the form

with a given initial condition (IC), $u(0)=uIC$⁠.

Often, the complexity of chaotic systems of interest is high, and accordingly, the cost of an accurate computational estimate of $J∞$ becomes formidable.^9–12 As the scale, cost, and complexity of a computational simulation increase, efficient discretization methods become critical for accurately estimating quantities of interest.

Understanding the error in approximations of $J∞$ is nontrivial because statistical errors (errors due to a finite-time approximation) and discretization error (error due to a numerical approximation of solutions) are always simultaneously present. In the largest direct numerical simulation (DNS) and large eddy simulation (LES) cases, for example, it is typical to fix the period of simulation time over which an average is performed at some large number of characteristic times and validate that discretization error converges as expected, assuming negligible sampling errors.^13–16 A recent work has sought to quantify the effect of statistical error more robustly, using a turbulent flow theory,¹⁷ advanced spatiotemporal statistical post-processing methods,¹⁸ statistical windowing techniques,¹⁹ or by extending the concept of Richardson extrapolation to chaotic flows using auto-regressive models and Bayesian methods.²⁰ The latter work is notable for its use to estimate the statistical errors in the DNS of a high-$ℜ$ turbulent channel flow.²¹ 

The objective of this paper is to investigate the behavior of statistical and discretization errors as a function of the computational cost for ergodic systems. Following a similar approach of Oliver et al.,²⁰ we propose a simple error model for finite-time, discrete approximations of infinite-time averages on attractors. Using the Lorenz system as an example, we demonstrate that the discretization error converges as the time step size decreases. However, if we denote by “sampling time” the amount of simulation time on which the integral is taken (after discarding spin-up), we show that error does not increase exponentially with sampling time as might be expected from classical numerical analysis but rather asymptotes to a constant, time step-dependent value with respect to the sampling time. Furthermore, for a given computational cost (e.g., the number of timsteps), an optimal choice of discretization (i.e., time step) exists that minimizes the expected error in simulation, when accounting for both the effects of discretization error and sampling error. We show that this optimal choice results in a convergence rate with respect to computational cost that is bounded by the sampling convergence rate with a minor impact from the discretization order of accuracy. Finally, we consider the implications of spin-up time (i.e., unsampled time needed to arrive at the stationary distribution) and parallelism on the optimal error. We develop a method for estimating transient-related errors and then evaluate optimal choices incorporating the results.

To approximate $J∞$⁠, we compute finite-time, discrete estimates of the outputs of interest of the true system as

where the notation $Iab(·)$ here represents the quadrature approximation of the integral $∫ab(·)dt$ of a quantity $(·)$ between a and b. Here, we have made a discrete approximation of the state using an order-p discretization with a temporal grid with characteristic size $h=Δt$⁠, where an order-p discretization is one for which the discretization error behaves as

when the discretization is applied to a well-posed (non-chaotic) system. Then we sample a discrete state over a finite sampling period, T_s, starting at some initial time t₀. We can define the error that is incurred as

By introducing a third value,

we can re-write the error using an identity

Here, we define the “discretization error” and “sampling error,” respectively, as

We can take an absolute value of both sides of (7), followed by a manipulation using the triangle inequality

Thus, the total error incurred by the approximation is bounded by the sum of the absolute discretization and sampling errors. Next, we define the attractor of the operator f, $A$⁠, as the set of long-term states toward which all trajectories converge independently of the initial condition.²² We can define the expectation $EA[ϕ(u0)]$ for a generic function $ϕ$ as the expectation taken over all the trajectories that can result from starting from points on the attractor, $A$⁠,

For the case in question, we will be considering either

or

with, for these examples, $u(t0)=u0∈A$⁠. Given these definitions, we can now take the expectation of (10), giving

by linearity.

From here, we propose asymptotic forms for the two right-hand side terms in (12). Consider the definition of e_T in (9)

Assuming that we choose t₀ such that each $u0$ is effectively an independent sample from the attractor's stationary distribution, then the quantity $g(u(t))$ is a random variable drawn from a stationary distribution. The states of ergodic systems, in general, are not independent in time, but as long as the system has satisfactorily strong mixing properties, the central limit theorem (CLT) can be applied to finite time averages of its outputs. This is the case whenever the condition of α-mixing is met,^23,24 which has been shown for the Lorenz system.²⁵ Thus, we can write e_T as

with a constant coefficient A₀, where $N(μ,σ2)$ gives the normal distribution with mean μ and variance $σ2$⁠. If we take the absolute value of this random variable, the result is a halfnormal distribution

where $H(σ2)$ gives a halfnormal distribution such that $|X|∼H(σ2)$ when $X∼N(0,σ2)$⁠. The expectation of the half-normal distribution is well defined, allowing

as T_s goes to infinity.

Now consider the use of a time-stepping method to give a discrete approximation $uhp(tn)$ of $u(tn)$ for each $tn=n(Δt)$⁠. Following the classical analysis,²⁶ we might expect that the discretization error should take a form

This analysis is based on bounding the growth of local truncation error at each time step by the Lipschitz constant, Λ, of the underlying system, with C_p being a constant parameter that depends on the choice of method. However, Viswanath showed²⁷ that, the global error could be modeled by a form

where $E(Ts;p)$ could be bounded by a constant for some nonlinear but non-chaotic systems that are exponentially stable. While this result has not been extended to an ergodic system, the expected convergence onto the attracting set suggests a bound of the form

As our results in Sec. III will show, (19) is a good description of the expected discretization error.

Thus, taking (12), (16), and (19), we assume a bound of the form

that bounds $EA[|eT,hp|]$ when $Δt$ is small enough and T_s is large enough to satisfy the asymptotic assumptions. Here, we have used $Ts=NsΔt$⁠, where N_s is the number of timesteps computed for sampling; additionally, q refers to the observed discretization convergence behavior, which, in practice, may differ from p due to numerical cancelations or if the solutions of the system are insufficiently regular. Similarly, r is an observed sampling convergence rate coefficient, which we expect to be 1/2 asymptotically under the CLT.

In this section (Sec. III), we will fit numerical results for the Lorenz system to determine q, r, C_q, and A₀ and show that this model is representative of the observed behavior. The Lorenz system is given by:²⁸ 

where $u=[u0,u1,u2]⊤$ and $α=[α0,α1,α2]⊤$⁠. The Lorenz system is known to be chaotic for the classic Lorenz parametrization:²⁹ $α=[10,28,8/3]$⁠, which is used everywhere in this text. For the output, we choose $g(u)=u2$⁠. We consider a set of explicit methods: forward Euler (FE, p = 1), third-order Runge–Kutta (RK3, p = 3), and fourth-order Runge–Kutta (RK4, p = 4). In all of these methods, we expect asymptotic convergence of $JT,hp$ to J_T to be at least $O(Δtp)$ for non-chaotic systems.³⁰ The stability limit for each of these methods is $Δt≳10−1$⁠, which is above the limits of asymptotic convergence for all of the methods here; tests with implicit methods indicated similar error trends at higher cost and are, therefore, omitted.

For any given discrete instance, we will start the simulation at an initial state at t = 0 that is sampled randomly from a normal distribution

To guarantee that the initial sampling state $u0$ at t₀ is on the attractor (as well as further guaranteeing the independence from the other Monte Carlo instances), we evolve the state of any given Lorenz system discretization from its starting state $uinit$ for $t0=100$ before proceeding to the sample; we refer to the process of evolving the solution until it is on the attractor as “spin-up.” Then, we evolve the state over the next T_s, during which we integrate and compute (3) using the same numerical integration scheme that was used for the state itself.

To approximate $eT,hp$⁠, we must first estimate $J∞$ by a reference value $Jref$⁠. $Jref$ is calculated using an ensemble mean of $JT,hp$ over $Mens=5122$ instances of the Lorenz system. Each instance is started from a different $uinit$ as given in (22) and simulated using RK4 with $Δt=17.7×10−6$ and $Ts=6646.9$⁠. The resulting $Jref$ is

with a 95% confidence estimate based on the ensemble mean estimator.

The computation of $Jref$ allows us to estimate errors $eT,hp≈JT,hp−Jref$⁠. For a given $Δt$⁠, T_s pair, we then approximate $EA[|eT,hp|]$ using a Monte Carlo method over $M=10 000$ independent instances of the discrete system, each started from initial states drawn from (22) and spun-up to independent sampling starting points on the attractor $u0(m)$⁠,

In Figs. 1–3, we compare the results of simulations with the FE, RK3, and RK4 discretizations with different values of N_s. In these figures, T_s scales with $Δt$ for a given N_s, so the T_s values on the x-axis will vary between lines on the plot. The fits shown are computed with truncated data in order to attempt to eliminate non-convergent data at small T_s or large $Δt$⁠; the limits used for truncation are found in Table I. The results of the nonlinear least squares fits for $Ns=104$⁠, 10⁵, and 10⁶ are given in Table II. In the table, we observe that $r→1/2$ as the discretization error is reduced, either by increasing N_s or by pushing p higher. These figures demonstrate that (19) has explanatory value, as the errors in the discretization-dominated region collapse independently of T_s. It is also worth noting that Table II demonstrates higher-than-expected discretization error convergence rates for FE and RK4. In Fig. 4, we can examine the sampling error behavior between discretization methods for a single shared choice of N_s. Here, we can see that the sampling error effects on the left-hand side of the plot collapse independently of the discretization method. This indicates that the statistical effects are properties of the dynamical system, not artifacts of the discretization, as we might expect in the limit as $Δt→0$⁠. We note that Table II shows some N_s-dependence in the fit values of r; we expect that this is an artifact of the fit process, in which discretization effects and statistical effects are intermingled. As N_s or q become small, the region in which the convergent CLT behavior is the dominant effect becomes smaller, and for these cases, the value of r detected is subject to error and deviation from the expected rate r = 1/2.

Finally, we attempt to compare the computational costs across the various discretizations. In this case, the number of timesteps N_s is not a good proxy for fixed cost, since the computation time for a time step will vary between methods. Instead, we now fix U_s, the total number of evaluations of the right-hand side f used in sampling timesteps. For the explicit schemes used in this work, we will have p right-hand side evaluations (e.g., forward Euler has p = 1 right-hand side evaluations) and, thus, $Us=pNs$⁠. In Fig. 5, we can see the effect of changing $Δt$ at fixed sampling cost U_s across discretizations. The error that can be achieved with the Runge–Kutta methods is lower than that of the forward Euler scheme, a factor of 4.8 improvement in the error from FE to RK4. However, the best-case improvement for going from third-order to fourth-order Runge–Kutta schemes is the only factor of about 1.4. Moreover, the results show that to achieve the lowest possible error, the optimal time step will be discretization dependent. We investigate this further in Sec. IV.

We now study the implications of the error model (20), specifically seeking to understand the convergence of the error with respect to computational effort. In this analysis, we will assume that r = 1/2.

Consider a non-dimensional form of the error model in which the error is normalized by the standard deviation of the instantaneous output σ_g, and the timescales $Δt$ and T_s are normalized by decorrelation time T_d. The decorrelation time relates the amount of variance from independent draws from the distribution on the attractor and the amount of variance in the finite-time mean estimators based on the correlated output signal, given by the relation:³¹ 

Furthermore, combining (15) and (25) allows us to write

In general, T_d is hard to estimate accurately; this is a crux of the work of Oliver et al.²⁰ In our formulation of the error model, we identify A₀, which avoids outright estimation of T_d. However, for the purpose of understanding the behavior of the error, T_d is an intrinsic timescale, which can be used to normalize $Δt$ and T_s.

The resulting non-dimensional form of the error model is

We can also write the optimizers and optimal value of (27) in terms of the non-dimensional variables. These are given by

In terms of convergence with respect to sampling costs, the error model will scale at best as

In the limit as $q→∞$⁠, the rate $q/(2q+1)→1/2$⁠: the CLT limits the convergence rate. Table III gives the rates of convergence (28) for various values of q.

Using the reference simulation, we can also find

where $σ̂g$ is an estimate of the standard deviation of g. Together, these allow us to estimate

With these values, we can plot the non-dimensional error model with fixed r = 1/2, which is given for $Ns=105$ in Fig. 6.

We now consider the implications of these results for increasing N_s. To focus solely on control of the discretization error, increases in N_s can be used to refine $Δt=Ts/Ns$⁠, with T_s fixed. On the other hand, to focus solely on controlling the sampling error, $Ts=NsΔt$ can be increased, holding $Δt$ fixed. In Fig. 7, the two approaches are compared with the optimal use of resources. In orange is the discretization error control strategy. In this approach, the simulations converge at a high-order rate in N_s toward the optimal error behavior; once the error reaches this optimum, however, it asymptotes to a constant: statistical errors limit the estimation of $JT,hp$⁠. On the other hand, the sampling error control approach is shown in blue. In this approach, the central limit convergence rate of 1/2 is initially achieved until the error asymptotes to a constant: discretization errors limit the estimation of $JT,hp$⁠. In other words, this latter case represents having increasingly many significant figures of certainty about an answer that is significantly offset from the true quantity of interest, $J∞$⁠. In the literature for large simulations, discussed in the introduction, simulations tend to be planned using either the discretization or statistical error control approach. What (20) implies and Fig. 7 demonstrates is that, in fact, there is a particular optimal scheme in which $Δt$ and T_s are simultaneously varied that will extract the most accurate estimate of $J∞$ as N_s increases.

In this section, we show that our simulations of chaotic, ergodic ODEs are consistent with a bounded relationship between the local and global discretization errors. Consider an estimate of the global error based on N_s timesteps

where

and $u⋆(tn+1)$ is the exact solution integrated from $uhp(tn)$ through $Δt$

In (31), we have assumed that the error from any given local state perturbation is propagated forward in time by the dynamics, before being transformed into an error in the output; this process is captured by an operator $G$⁠. Because the effect of local error propagates forward and not backward in time, $G(t,tn)=0$ for $t<tn$⁠, and moreover, we assume that due to ergodicity $G(t,tn)=0$ when $t−tn≳Td$⁠, where T_d is the decorrelation time associated with the attractor. This allows us to write

Now, we assume that a constant $Gmax$ exists such that

for all $tn,tη∈ℝ$⁠, and $v∈B(u(tn))⊂ℝd$⁠, where $B(u)$ is the set of states possible by perturbation of u that remain in the basin of attraction of the attractor $A$ of f. When this is the case, we can create a bound on the magnitude of e_hp

We now attempt to bound the value of $Gmax$ for the Lorenz system by approximating the local truncation error. To make an estimate, we compute both the solution at the next time step as well as a surrogate for the true solution at each time step: $uhp(tn+1)$ and $ũ⋆(tn+1)$⁠, where the former is computed with one time step of the method of interest and the latter is always computed with the highest available accuracy method, RK4, and subdividing $t∈[tn,tn+1]$ into ten consecutive timesteps rather than one. Both $uhp(tn+1)$ and $ũ⋆(tn+1)$ are always advanced from $uhp(tn)$⁠. This allows us to estimate $eLT,p(n+1)$ locally

In Fig. 8, we characterize the convergence of local error estimates. Computations are run with T_s = 100 and $t0=100$ fixed, varying $Δt$⁠. At each time step, the local truncation error is estimated by computing (37). The figure shows the computed $maxn||ẽLT,p(n)||∞$ and demonstrates that the expected rate of $(p+1)$ is nearly exactly achieved.

Using (36), we can estimate a bounding value for $Gmax$ by

where c_p is the leading truncation error coefficient fit in Fig. 8 and C_q and q are taken from Table II. Of course when $q>qtheory=p$⁠, there will be $Δt$ dependence. However, as (38) requires that the discretization error has an asymptotic behavior, we will only consider $Δt$ in the asymptotic convergence regions given in Table I to compute $Gmax$⁠. In Fig. 9, we show the values of the right-hand side quantity in (38), which allow us to make an estimate

Next, we use classical truncation error estimates²⁶ to relate the discretization error to properties of the solution. We will assume that the local truncation error is bounded by a form

where $CLT$ is a local truncation constant term dependent on the numerical method and the $||·||∞$ in this context refers to the maximum value in time of the inf-norm of a vector-valued, time-dependent quantity $(·)$⁠. The derivatives of $u(t)$ can be computed by evaluating $f(u)$ and its derivatives using solutions from a reference RK4 solution of the Lorenz system with T_s = 1000, $t0=100$⁠, and $Δt=10−4$⁠. Norms of the derivatives are shown in Fig. 10. The resulting values of $CLT$ that can now be derived by fitting the asymptotic behavior in Fig. 8 can be found in Table IV. The result of these estimates is that we can reliably bound the global error of a dynamical system as an accumulation of the local errors over a region of correlation.

We now want to consider how the global error behavior demonstrated here might extrapolate to more complicated systems by evaluating the spectral behavior of the Lorenz system. Using a discrete Fourier transform (DFT) with a Hann window function,³² we perform a spectral analysis on the states of the Lorenz system with sampling time T_s = 1000, $t0=100$⁠, and $Δt=10−3$⁠. The resulting spectrum can be found in Fig. 11. The Lorenz system tends to have the most content in the frequencies with $f ≲ 101$ with a region of exponential decay in the range $1 ≲ f ≲ 300$⁠. On scales with $f≳300$⁠, machine precision plateaus are observed and omitted here.

The fact that the Lorenz spectrum is an exponentially decreasing function of frequency f makes the use of high-order methods theoretically appealing for the spectral convergence of hp-refinement strategies.³³ Unfortunately, the effect of statistical error in (28) limits the impact of this exponential decay, such that the benefits of higher-order discretization methods are limited compared to their steady-state and non-chaotic application. The convergence to the central limit rates can be seen in Fig. 12, which shows the convergence of (28) with the total sampling cost. The effect of increasing order improves the convergence rate in (28) toward the CLT-implied asymptotic rate of $−1/2$⁠, as well as decreasing the value of the leading constant and the error never achieves the spectral rates possible with hp-refinement in the steady case. Nevertheless, the cost to achieve a given amount of error in expectation—in terms of function evaluations—is significantly less with higher-order methods. Managing to achieve 1% non-dimensional error in expectation is possible with RK4 at a cost ten times less than would be possible using FE; this factor grows larger than 100 when the tolerance is tightened to $10−4$⁠.

In this section, we will consider how the error behaves when ensemble averaging (over multiple parallel instances) and when spin-up effects are present.

Sampling error can be reduced at a fixed wall clock time by ensemble averaging across multiple parallel processes.³⁴ Consider a Monte Carlo approach to approximate $J∞$ with a set of $Mens$ independent realizations

We can write a modified version of (20) to approximate the error that we expect in the Monte Carlo estimator in (41)

with an equivalent non-dimensional version, assuming $r→1/2$

and an optimum given by

at

and

Equation (44) shows that, for finite values of q, the Monte Carlo method will have a mitigated return compared to its purely stochastic application as in Ref. 34; the optimal error scales as $Mens−q/(2q+1)$ as opposed to $Mens−1/2$⁠. However, parallelization can achieve perfect scaling in the expected error, in the sense that the effect of running $Mens$ ensembles with N_s sampling timesteps each will have an equivalent error in expectation to simulate $MensNs$ timesteps in series. As $Mens$ is varied on the set of optimal solutions, (45) and (46) indicate that the time step and sampling time should be adjusted with the same factor $Mens−1/(2q+1)$ to achieve perfect scaling.

So far, we have considered the error and cost on the attractor, neglecting the impact of spin-up from t = 0 to t = t₀. This spin-up is necessary because simulations of ergodic systems invariably need some time for the state to proceed onto the attractor from the initial condition.

Consider $u(t)$⁠, a solution of the ergodic chaotic system f from an arbitrary initial condition $u(0)=uIC$ in the basin of attraction of an attractor, $A$⁠. The existence of the attractor implies the nonlinear stability of the system, such that all $uIC$ will converge to trajectories on the attractor $A$⁠. Denote by $uA(t)$ a trajectory that is on the attractor for all t and to which $u(t)$ collapses as $t→∞$⁠. The perturbation $δuA(t)≡u(t)−uA(t)$ that describes the IC, therefore, exists in a stable subspace of perturbations to $uA$ and can be associated with the negative Lyapunov exponents of the system. Thus, we can assume that such perturbations are governed asymptotically by

with $Tλ$ being a characteristic time associated with the stable Lyapunov modes. In practice, we are interested in averages of quantities on the attractor $g(uA(t))$⁠, but we can only calculate quantities $g(u(t))$ that will include some effect—if small—of the spin-up transient.

Next, we seek to quantify the effect of this gap on estimates $JT≈J∞$⁠. Consider the computation of J_T. In (7), we have effectively found an estimate of

by choosing t₀ sufficiently large. We now want to consider an error model of the form

where a new error $eλ$ is introduced, associated with the spin-up transient. The model for e_T in (9) will apply without modification, while the model for e_hp will be subject to slightly different assumptions, where in (8), C_q was bounded by the value on the attractor, $A$⁠, here we must assume that C_q is bounded from t = 0 to $t=t0+Ts$⁠, including both the attractor and the transient part of the trajectory. We only require that the transient part be in the basin of attraction of $A, B(A)$⁠. We assume that a model of the form used in (8) applies in expectation when the transient component is included.

Next, we concentrate on $eλ$⁠,

We now assume that, like u, g will decay exponentially in t as (47), such that

will apply for $t∈[0,∞)$⁠, with $Aλ$ being a constant that can be related to the deviation between $g(u(0))$ and $g(uA(0))$⁠.

From this assumption,

Taking the absolute value, we can find a bounding model

As before, manipulation of (49) allows

Now, we take an expectation of the absolute value of $eT,hp$⁠,

where $EB(A)$ gives the expectation on the basin of attraction of $A$⁠. Here, the expectation of $|eT,hp|$ does not reduce to an expectation on the attractor. The statistical term is handled on the attractor as before, and we have assumed that the discretization error is bounded by the same form in expectation on $B(A)$ as on $A$⁠. Finally, the expectation of $|eλ|$ is taken on the set of initial conditions used. This allows us to take the expectation of (53) to complete (56). Because we anticipate $Tλ$ will be bounded by a constant for a given system, this is given by

If a $Aλ$ and $Tλ$ can be identified by observation of $g(u(t))$ given an initial condition $uIC, |eλ|$ is no longer stochastic and the $E[|eT,hp|]→|eT,hp|$ as in (53).

Putting all the pieces together, we can now give an error model that incorporates the effects of spin-up and ensemble estimation

where $Ãλ$ can be either estimated on an instance-by-instance basis or by estimating the expectation on the family of initial conditions. Under this model, $eλ$ will scale with the exponent of a large negative value when $t0≫Tλ$⁠. Even when $t0≫̸Tλ$⁠, (53) suggests that the decay-induced error term will still scale with $Ts−1$⁠, faster than the expected CLT rate of $Ts−1/2$⁠, and thus, it will be dominated it as $Ts≫1$⁠. This implies two “paths” to control spin-up errors: either choosing t₀ long enough to shrink the mean offset error from t₀ or choosing T_s long enough so that the mean offset contribution to the simulation error is small in spite of the error at t = t₀.

We will now develop a method to fit the error model. In order to do so, consider observations $gn≡g(tn)$ and $gnA≡gA(tn)$ for t_n in ${t0,t0+NskipΔt,…,t0+Ts}$⁠. We will assume that $Nskip$ is large enough that the solution at each t_n is effectively independent. If this is the case, then we can assume that each $gnA$ will be an independent and identically distributed (i.i.d.) draw from a bounded, stationary distribution with mean $J∞$⁠. The distributions of $gA(u(t))$ and $g(u(t))$⁠, in general, are not known. In order to facilitate an estimate of the mean behavior, we will assume $gnA$ are i.i.d. drawn from a normal distribution with mean value $J∞$⁠. Then, we have

where the relationship between g_n and $gnA$ is taken from (51).

In order to understand the implications of this model, we can use a set of reference RK4 simulations of the Lorenz system with $Nt=105$ timesteps sampled without spin-up over a period T = 100 from initial conditions similar to those given in (22) with a scaled-up standard deviation of 100 in all three variables to highlight the initial transient. In order to treat each of $J∞$⁠, σ_g, $Aλ$⁠, and $Tλ$ in (59) as unknowns, we use Hamiltonian Monte Carlo with the likelihood function implied in (59). We discard from t = 0 to t = 5 and then take 10 000 equispaced samples from t = 5 to t = 100. For prior models, we start by computing naïve estimators of the mean and standard deviation of the trace, $J̃$ and $σ̃$ using the downsampled trace signal ${gn}$ and then use