Precision and convergence speed of the ensemble Kalman filter-based parameter estimation: setting parameter uncertainty for reliable and efficient estimation – Progress in Earth and Planetary Science

Overview

Figure 6 shows the time series of the parameter ensemble mean for all experiments listed in Table 1; Fig. 6a–d presents the experiments for μ_init = 0 (Exps. 1–4), where 50 trials with different initial parameter ensemble perturbations were carried out for each setting. This figure shows that the parameter ensemble mean successfully remains around the true value of the nature run (x = 0) with some fluctuation; these time series appeared to be statistically stationary. Figures 6a–d clearly illustrate that the fluctuation frequency varies depending on σ: the time series with a small ES fluctuated slowly, while a large ES fluctuated rapidly. The amplitude of the fluctuation was indicative of the precision; the time series with the largest ES (σ = 0.4) fluctuated with the largest amplitude, indicating low precision. The difference in the fluctuation amplitude between σ = 0.05 and 0.2 was unclear, while some extreme values appeared for the experiment of σ = 0.2.

Time series of parameter ensemble mean for all experiments: a–d experiments with an initial ensemble mean μ_init of 0 (Exps. 1–4). The parameter ensemble spread (σ) for each experiment was: a 0.05; b 0.1; c 0.2; and d 0.4. Black lines represent the time series for a total of 50 trials for each setting. Red lines represent averaged time series of all trials. The origin of the time axis corresponds to T = 4 h for the nature run; e–h experiments where μ_init = 0.55 (Exps. 5–8: red lines), and − 0.55 (Exps. 9–12: blue lines); e experiments where σ = 0.05 (Exps. 5 and 9); f σ = 0.1 (Exps. 6 and 10); g σ = 0.2 (Exps. 7 and 11); and h σ = 0.4 (Exps. 8 and 12). Initial times of the experiments were every 1 h from 3 h following the commencement of Exps.1–4. The estimation duration was 6 h for each experiment. The beginning and end points of estimation are denoted by circles and squares, respectively. Black lines represent the averaged time series of Exps.1–4, which are the same as the red lines in (a–d)

Figure 7 shows the power spectral density of the time series for Exps. 1–4. For the high-frequency range, the spectral amplitude increased with σ. This corresponded to the frequently fluctuating time series for the large σ, as shown in Fig. 6. In all cases, the spectral amplitude decreased with increasing frequency, showing the “red noise” spectrum characteristics. The red noise spectrum is a feature that indicates the first-order autoregression process. In the later section, we approximate the experimental time series as the first-order autoregression model (AR(1) model) in order to understand the behavior of the time series clearly.

Power spectral density of the time series of parameter ensemble mean for Exps. 1–4. The red line is for the parameter ensemble spread (σ) of 0.05, the purple line is for σ = 0.1, the cyan line is for σ = 0.2, and the yellow line is for σ = 0.4. Each spectrum was calculated by applying the discrete Fourier transformation to the time series of each trial and averaging these values

Figure 6e–h shows the results from experiments with μ_init = 0.55 (Exps. 5–8: red lines) and − 0.55 (Exps. 9–12: blue lines). For this set of experiments, there were 13 trials conducted for each setting. The initial times of the experiments were every 1 h from 3 h following the commencement of Exps. 1–4; the initial time of these experiments was shifted from Exps. 1–4 because the time series of the parameter ensemble mean for each trial of Exps. 1–4 (black lines in Fig. 6a–d) during the early simulation time did not sufficiently spread out from the initial value provided. As such, it was not appropriate to carry out a statistical comparison of the time series during this period. We also modified the initial times for each trial of Exps. 5–12, as the convergence speed may vary depending on the simulation time. The estimation duration was 6 h for each experiment. The parameters successfully converged toward the true value in the μ_init = 0.55 and − 0.55 series; the convergence speed accelerated as σ became larger, which supports the findings from Kotsuki et al. (2018).

Accuracy

Although all experiments appeared to attain the true value of the nature run successfully, it was not clear whether the true value was accurately estimated for the entire simulation time. Figure 8 shows the averaged time series for all 200 trials of Exps. 1–4. Red open circles on the black line indicate that the mean of the estimates differs from the true value (x = 0) with a significance level of 1% based on the Student’s t test. From 14 to 20 h following the commencement of parameter estimation, the estimated parameters tended to be higher than x = 0. During this period, time steps where the mean estimates significantly differed from the true value continued intermittently. This indicates that the estimation accuracy varies depending on the observations for each simulation time, despite carrying out the estimation under the twin experiment framework. In this study, we considered the 3 h moving average of the mean time series (blue line in Fig. 8) as the bias component common to all experiments. Sections 3.3 and 3.4 discuss the precision and convergence speed of the estimation, accounting for this estimation bias.

Averaged time series for all the trials of Exps. 1–4. The black line is the averaged time series. Color shade represents the standard deviation of all trials. Red open circles indicate that the mean of the estimates differs from the true value (x = 0) with a significance level of 1% based on the Student’s t test. The blue line shows the 3 h moving average of the black line

Precision

The estimation precision is defined as how close the estimated values for different trials are to one another. In this study, the metric of estimation precision was the sample standard deviation of estimation time series with the bias component removed. The sample standard deviation of the estimated parameter was calculated as:

where μ(t) is the time series of the parameter ensemble mean for each time step, t; x_bias(t) is the bias component (blue line in Fig. 8); and T is the number of time steps of the estimation. Here, a smaller σ_μ indicates higher precision; Fig. 9 presents the mean σ_μ, averaged for all 50 trials for each parameter ES (σ). The σ_μ became smallest when σ was 0.1 or less; this implied that σ = 0.1 was the threshold for the parameter ES. There were no further improvements to the estimation precision even if σ was further reduced from a value of 0.1. As the convergence speed accelerates as σ decreases (see Sect. 3.4), a σ of 0.1 was optimal to provide the best precision and better convergence speed in this experimental framework. Section 3.5 presents the quantification of the dependency between precision and the parameter ES using the AR(1) model.

Sample standard deviation of estimated parameters for Exps. 1–4. The standard deviation for each trial was calculated after removing the bias component (blue line in Fig. 8). Vertical bars represent a confidence interval of 99%

Convergence speed

Figure 10 illustrates the composite time series of the parameter ensemble mean for Exps.5–12. Note that the bias component shown in Fig. 8 was removed from each time series prior to compositing. The convergence speed accelerated as σ became larger. The open circles in Fig. 10 indicate that the set of time series did not differ from that for Exps. 1–4, with a significance level of 1% based on the Student’s t test. With the exception of Exp. 5 (red line in Fig. 10a), all experiments converged to statistically stationary time series within 6 h.

Composite time series of parameter ensemble mean for Exps. 5–12. Red (blue) lines are for the initial ensemble mean of μ_init = 0.55 (− 0.55); a experiments where parameter ensemble spread (σ) was 0.05 (Exps. 5 and 9); b σ = 0.1 (Exps. 6 and 10); c σ = 0.2 (Exps. 7 and 11); and d σ = 0.4 (Exps. 8 and 12). Note that the bias component (x_bias) shown as the blue line in Fig. 8 was removed from each time series prior to compositing. Color shade represents the standard deviation of all trials. The origin of the time axis corresponds to the initial times of each trial. Open circles on each line indicate that the set of time series does not differ from Exps.1–4, with a significance level of 1% based on the Student’s t test

Figure 11 shows the convergence times for each experiment; here, the convergence time was defined as the period from the initial time to the first open circle in Fig. 10. For experiments where σ = 0.4, the parameter converged within only a few DA cycles; this occurred at ~ 10 min. This means that if σ was not too small, successful parameter estimation occurred within a shorter time than the time scale of mesoscale convective systems (a few hours). The convergence speed is related to the persistence of the time series. This was discussed in quantitative terms by considering the AR(1) model in Sect. 3.5. Notably, when σ ≤ 0.2, the convergence speed for μ_init = 0.55 (Exps. 5–7) was slower than for μ_init = − 0.55 (Exps. 9–11). This indicates that the convergence speed differs depending on whether parameter estimation begins from a larger or smaller value than the true value. This is likely to imply that model sensitivity to the c_R value was not symmetric to c_R = 130. In addition, it is possible that the parameter likelihood itself was not symmetric around the true value.

Convergence time plotted against the parameter ensemble spread (σ). Red line shows Exps. 5–8 where the initial ensemble mean of μ_init = 0.55; blue line shows Exps. 9–12 where μ_init = − 0.55. The time series of Exp. 5 does not converge within 6 h, and as such, was not plotted

Mathematical model of parameter estimation time series

Thus far, this study has demonstrated that although the precision of parameter estimation may be maximized for σ ≤ 0.1 (Fig. 9), the convergence speed for estimation accelerates as σ becomes larger, and as such, the convergence time is shortened (Fig. 11). To quantify the precision and convergence speed of parameter estimation, we attempted to model the temporal variation of the estimated parameter. The spectral amplitude that decreases with increasing frequency in Fig. 7 suggests the estimation time series was close to the AR(1) time series. Here, we approximated the time series of the parameter ensemble mean using the AR(1) model:

where t is each DA time step; μ′(t) ≡ μ(t) − x_bias(t) is the time series of the parameter ensemble mean where the bias component is removed; ϕ is the autoregressive parameter; and ε(t) is a random perturbation. The expected value of ε(t) is assumed as zero (i.e., E[ε] = 0), and its variance is referred to as σ_ε² = E[ε²], where E[a] indicates the expected value of a. When the time series by Eq. (4) is in a statistically stationary state, E[μ′] = 0. Under this condition, the variance of μ′ is defined as, σ_μ² = E[μ′²]. As the AR(1) model assumes that μ′ and ε have no correlation (i.e., E[μ′ε] = 0), the following relationship between σ_μ² and σ_ε² is satisfied:

Then, the lag-1 autocorrelation coefficient of μ′, ρ(1), satisfies the following:

Note that the range of ϕ is (− 1, 1). In the AR(1) model, ϕ is a metric of the time series persistence, where if ϕ is close to 1, μ′(t) is highly correlated to the previous value μ′(t − 1). In terms of parameter estimation, this means that if a parameter value at one step deviates considerably from the true value, the parameter at the next step also tends to be highly deviating, resulting in a slow convergence speed.

For each trial, σ_μ² was estimated from Eq. (3); then, using Eq. (6), ϕ was estimated as:

Finally, using Eq. (5), we estimated the standard deviation of the random perturbation, σ_ε, as:

Figure 12 shows the ϕ and σ_ε values, estimated by calculating Eqs. (7) and (8) for all trials and averaged for each experimental setting; here, ϕ approaches 1 as σ decreases. This trend corresponds to a slow convergence speed of the estimated parameter for the experiment with a smaller σ. At the same time, σ_ε increases with σ; this trend corresponds to a larger amplitude of the time series fluctuation for experiments with a larger σ.

Two AR(1) model parameters estimated from the time series of Exps. 1–4: a estimated autoregressive parameter ϕ of Eq. (7); and b estimated standard deviation of the random perturbation σ_ε of Eq. (8). Each plot represents the mean value for all the trials of Exp. 1 (σ = 0.05), Exp. 2 (σ = 0.1), Exp. 3 (σ = 0.2), and Exp. 4 (σ = 0.4). The vertical bars represent a confidence interval of 99% for each estimated value

The power spectral density function of AR(1) model is expressed as:

where f is frequency (Wilks 2011). Here, f was discretized as f = n/T for n = 0, 1, …, T/2; f = 1/2 was the Nyquist frequency corresponding to a 10 min period. Figure 13 presents the power spectral density function of Eq. (9) using the estimates of ϕ and σ_ε in Fig. 12 (dotted lines), alongside the power spectral density directly derived from the experimental time series (solid lines). Note that the experimental power spectral density slightly differed from Fig. 7 because the bias component was removed from each time series. For σ = 0.05, the spectrum function was in good agreement with the experimental spectrum, except for the right end of the spectrum (red lines). We inferred that the AR(1) model could suitably approximate the parameter estimation time series for the target parameter c_R. As σ increased, the spectrum function deviated slightly from the experimental spectrum, particularly for low-frequency components. This disagreement may be related to the relatively large uncertainty of estimated ϕ, when σ is large (Fig. 12a). As for the disagreement in the low-frequency range, it is also possible that the frequency components having a period of several hours and more were not detected accurately since the simulation time was just 24 h. From the perspective of the AR(1) model, the rapidly fluctuating time series for a large σ is explained in terms of the dependence of power spectral density function on ϕ. As the parameter ES (σ) becomes larger, the autoregressive parameter ϕ value declines (Fig. 12a), resulting in the increase of power spectral density for the high-frequency (larger f) range as per Eq. (9).

Power spectral density function assuming the AR(1) model. The lines of spectral functions were based on Eq. (9) using estimated values of autoregressive parameter ϕ and standard deviation of the random perturbation σ_ε for experiments where σ = 0.05 (red), 0.1 (purple), 0.2 (cyan), and 0.4 (yellow) (dotted lines). Solid lines show power spectral density directly derived from the experimental time series. The experimental power spectral density differed slightly from Fig. 7 because the bias component (blue line in Fig. 8) was removed from each time series prior to calculating power spectral density

Quantitative understanding of precision and convergence speed

This section quantifies the precision and convergence speed of parameter estimation based on the AR(1) model. Figures 9 and 11 present the experimental results, illustrating the dependence of precision and convergence time on the parameter ES (σ). This section reconsiders this in terms of the dependence of ϕ and σ_ε in Eq. (4) for σ. Based on Eq. (4), the convergence from the initial wrong value to the true value is expressed as:

where μ′(0) is the initial parameter ensemble mean; we used E[ε] = 0. The e-folding time was estimated by substituting μ′(0)/e into the left-hand side of Eq. (10) and solving it for t:

Here, we assumed ϕ > 0, which is reasonable for this experiment (Fig. 12a). As the DA time step, t, should be an integer, the actual time step where the expected deviation from the true value is < 1/e of its initial deviation assuming the AR(1) model, is:

Here, (lceil a rceil) indicates the least integer greater than or equal to a. Figure 14 shows the t_AR(1) calculated using Eq. (11), where the estimated ϕ was adopted from Fig. 12a. The right vertical axis of the figure presents the actual time in the correspondence, where one time step is 5 min. The AR(1)-based convergence time in Fig. 14 was in good agreement with the convergence time shown in Fig. 11, indicating that the AR(1) model approximation is reasonable. The convergence speed accelerates as σ increases. In terms of convergence time, a larger σ was more advantageous for parameter estimation; however, the acceptable convergence time was dependent on the time scale of the phenomena (e.g., the time scale of deep convective clouds). If convergence time is within an acceptable range, the optimal parameter ES should be determined to maximize precision.

Time step until the expected deviation from the true value becomes < 1/e assuming the AR(1) model. Each value was calculated based on Eq. (11) using the estimated autoregressive parameter, ϕ, for experiments where σ = 0.05, 0.1, 0.2, and 0.4 (Fig. 12a). The right vertical axis shows the actual time in correspondence where one time step is 5 min

The standard deviation of the estimated σ_μ is a metric of estimation precision, which is expressed by rewriting Eq. (5) for σ_μ:

Figure 15 shows the σ_μAR(1) obtained using Eq. (12), where the estimated ϕ and σ_ε values are shown in Fig. 12; here, σ_μAR(1) was minimized at σ = 0.05. However, the σ_μAR(1) values for σ = 0.05 and 0.1 did not differ greatly; σ_μAR(1) = 9.63 × 10⁻² for σ = 0.1, was only 4% larger than σ_μAR(1) = 9.27 × 10⁻² for σ = 0.05. By contrast, σ_μAR(1) = 1.06 × 10⁻¹ for σ = 0.2, was 11% larger than for σ = 0.1. Thus, we suppressed the estimation uncertainty by setting σ to ≤ 0.1. With a decrease in σ, ϕ and σ_ε were asymptotic to 1 and 0, respectively (Fig. 12). This means that the numerator and denominator on the right-hand side of Eq. (12) becomes zero in the σ = 0 limit. By decreasing σ, the behavior of σ_μAR(1) is determined by how the numerator and denominator of Eq. (12) approach zero. For the current target parameter, the ratio of the speed at which the denominator and numerator were asymptotic to zero was almost constant for σ ≤ 0.1; σ_μAR(1) was minimized for this condition. This demonstrates that the parameter ES of σ ≤ 0.1 achieves the best estimation precision.

Standard deviation of estimated values assuming the AR(1) model. This standard deviation was based on Eq. (12) using the estimated ϕ and σ_ε for experiments where σ = 0.05, 0.1, 0.2, and 0.4 (Fig. 12)

Following this, we determined the minimum ES for an acceptable convergence time. The convergence time for σ = 0.1 was approximately 1 h, as shown in Fig. 14. This time scale was comparable or shorter than long-lasting convective systems. On the other hand, the convergence time for σ = 0.05 was approximately 3 h. This is likely to be excessive duration for practical parameter estimation to target deep convection. In terms of convergence speed, we accepted a value of σ ≥ 0.1; and therefore, concluded that σ = 0.1 was the optimal parameter ES to estimate true c_R in this twin experiment. Table 2 summarizes the estimated ϕ and σ_ε values for different σ, and the convergence time step and the standard deviation of the estimated σ_μ, on the basis of the AR(1) fitting.

The first and second rows show the estimated autoregressive parameter, ϕ, and standard deviation of the random perturbation, σ_ε, respectively. The third row presents the t_AR(1) calculated using Eq. (11); the fourth row presents σ_μAR(1) obtained using Eq. (12); these values were calculated using the estimated ϕ and σ_ε.