# 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

Sep 6, 2022

### 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.

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.

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.

### 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:

$$sigma_{mu } = sqrt {frac{1}{T}mathop sum limits_{t = 1}^{T} left{ {mu left( t right) – x_{{{text{bias}}}} left( t right)} right}^{2} } ,$$

(3)

where μ(t) is the time series of the parameter ensemble mean for each time step, t; xbias(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.

### 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.

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 cR value was not symmetric to cR = 130. In addition, it is possible that the parameter likelihood itself was not symmetric around the true value.

### 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:

$$mu^{prime } left( t right) = phi mu^{prime } left( {t – 1} right) + varepsilon left( t right),$$

(4)

where t is each DA time step; μ′(t) ≡ μ(t) − xbias(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 σε2 = E[ε2], 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, σμ2 = E[μ2]. As the AR(1) model assumes that μ′ and ε have no correlation (i.e., E[με] = 0), the following relationship between σμ2 and σε2 is satisfied:

$$sigma_{mu }^{2} = phi^{2} sigma_{mu }^{2} + sigma_{varepsilon }^{2} .$$

(5)

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

$$rho left( 1 right) = frac{{Eleft[ {mu^{prime } left( t right)mu^{prime } left( {t – 1} right)} right]}}{{sigma_{mu }^{2} }} = phi .$$

(6)

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, σμ2 was estimated from Eq. (3); then, using Eq. (6), ϕ was estimated as:

$$phi = frac{1}{T – 1}mathop sum limits_{t = 2}^{T} mu^{prime } left( t right)mu^{prime } left( {t – 1} right)/sigma_{mu }^{2} .$$

(7)

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

$$sigma_{varepsilon } = sigma_{mu } sqrt {1 – phi^{2} } .$$

(8)

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 σ.

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

$$Sleft( f right) = frac{{4sigma_{varepsilon }^{2} /T}}{{1 + phi^{2} – 2phi cos left( {2pi f} right)}},$$

(9)

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 cR. 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).

### 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:

$$left[ {mu^{prime } left( t right)} right] = phi^{t} mu^{prime } left( 0 right),$$

(10)

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:

$$t = – frac{1}{ln left( phi right)}.$$

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:

$$t_{{{text{AR}}left( 1 right)}} = leftlceil- frac{1}{ln left( phi right)}rightrceil.$$

(11)

Here, (lceil a rceil) indicates the least integer greater than or equal to a. Figure 14 shows the tAR(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.

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

$$sigma_{{mu {text{AR}}left( 1 right)}} = frac{{sigma_{varepsilon } }}{{sqrt {1 – phi^{2} } }}.$$

(12)

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−2 for σ = 0.1, was only 4% larger than σμAR(1) = 9.27 × 10−2 for σ = 0.05. By contrast, σμAR(1) = 1.06 × 10−1 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.

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 cR 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 tAR(1) calculated using Eq. (11); the fourth row presents σμAR(1) obtained using Eq. (12); these values were calculated using the estimated ϕ and σε.