Stochastic filters (ensemble Kalman filters)

Summary

File: src/estimators/ensembles.py. EnsembleEstimator is the concrete intermediate base; leaf classes only implement _analysis_kernel(Af, d, Cdd).

Key attributes: m (ensemble size), std_phi, std_alpha, regularization_factor

Class	Filter type	Reference
EnKF	Stochastic EnKF (perturbed observations)	Evensen (2003)
EnSRKF	Deterministic square-root EnKF	Tippett et al. (2003)
rBA_EnKF	Regularized bias-aware EnKF, weight γ	Nóvoa, Racca & Magri (2023)

Covariance inflation (src/estimators/inflation.py, Evensen 2009 Chap. 15): fixed multiplicative inflation via inflation_factor. Applied factors are logged in estimator.inflation_history.

---

romda.estimators.EnsembleEstimator(parent_model, parent_bias=None, **kwargs)

Bases: Estimator

Abstract base for ensemble-based estimators (EnKF, EnSRKF, rBA-EnKF).

Owns the model, ensemble, and bias. Subclasses only need to implement _analysis_kernel(Af, d, Cdd, **kwargs) -> Aa.

Attributes:

Name	Type	Description
m	int	Number of ensemble members.
std_phi	float	Fractional std for initial state perturbations.
std_alpha	float or dict	Std (or {name: (lo, hi)}) for initial parameter perturbations.
distribution_phi	str	Sampling distribution for state members ("normal" or "uniform").
distribution_alpha	str	Sampling distribution for parameter members.
ensure_mean_at_init	bool	Force one member to equal the ensemble mean at initialisation.
ensemble_psi0	ndarray or None	Pre-built initial ensemble; bypasses generation if provided.
activate_parameter_estimation	bool	Whether to include parameter rows in the analysis update.
regularization_factor	float	Bias-regularisation weight (used by rBA-EnKF; ignored otherwise).

Initialise ensemble estimator.

Parameters:

Name	Type	Description	Default
parent_model	Model instance or Model subclass	If a class is passed, remaining kwargs are forwarded to its constructor (e.g. dt, psi0, model-specific parameters).	required
parent_bias	Bias instance, Bias subclass, or None		None
**kwargs		Ensemble config keys (m, std_phi, std_alpha, …) and/or model constructor keys are accepted here.	{}

Source code in src/estimators/ensembles.py

@typechecked
def __init__(
    self,
    parent_model: Model | type[Model],
    parent_bias: Bias | type[Bias] | None = None,
    **kwargs,
):
    """Initialise ensemble estimator.

    Parameters
    ----------
    parent_model : Model instance or Model subclass
        If a class is passed, remaining *kwargs* are forwarded to its
        constructor (e.g. ``dt``, ``psi0``, model-specific parameters).
    parent_bias : Bias instance, Bias subclass, or None
    **kwargs
        Ensemble config keys (``m``, ``std_phi``, ``std_alpha``, …) and/or
        model constructor keys are accepted here.
    """
    #  Apply general DA attributes from Estimator base.
    super().__init__(**kwargs)

    # Split the remaining kwargs: ensemble-generation keys go to
    # init_ensemble; everything else (dt, psi0, model parameters, ...) to
    # the model constructor. Keys consumed by Estimator.__init__ are
    # forwarded to neither.
    ensemble_kwargs = allowed_kwargs_for_func(parent_model.init_ensemble, kwargs)
    model_kwargs = {k: v for k, v in kwargs.items()
                    if k not in self._consumed_kwargs and k not in ensemble_kwargs}

    # Instantiate or copy the model.
    if isinstance(parent_model, Model):
        parent_model = parent_model.copy()
    else:
        parent_model = parent_model(**model_kwargs)

    # Generate the initial ensemble inside the model, falling back to the
    # estimator's config for keys not passed explicitly.
    for attr in ('std_phi', 'std_alpha', 'distribution_phi', 'distribution_alpha'):
        if attr not in ensemble_kwargs:
            ensemble_kwargs[attr] = getattr(self, attr)
    parent_model.init_ensemble(**ensemble_kwargs)

    self._model = parent_model

    self._init_bias(parent_bias)

m property

Ensemble size.

inflation_history property

Namedtuple with fields times and factors of each applied inflation.

get_observables(Nt=1, **kwargs)

Return ensemble observables, bias-corrected if applicable.

Source code in src/estimators/ensembles.py

def get_observables(self, Nt: int = 1, **kwargs) -> np.ndarray:
    """Return ensemble observables, bias-corrected if applicable."""
    y_model = self.model.get_observables(Nt=Nt, **kwargs)
    if self.bias is None or isinstance(self.bias, NoBias):
        return y_model
    if Nt != 1:
        raise NotImplementedError(
            "Bias correction for get_observables with Nt > 1 is not yet implemented."
        )
    # unbiased = y + b, cf. _recover_unbiased_solution
    return y_model + self.current_bias_estimate

get_observable_hist(Nt=0)

Return (y_unbiased, y_model) history, interpolating bias if needed.

Source code in src/estimators/ensembles.py

def get_observable_hist(
    self, Nt: int = 0
) -> tuple[np.ndarray | None, np.ndarray]:
    """Return ``(y_unbiased, y_model)`` history, interpolating bias if needed."""
    pb = self.bias
    y_model = self.model.get_observable_hist(Nt=Nt)

    if pb is None or isinstance(pb, NoBias):
        return None, y_model

    t_model = self.model.hist_t[-Nt:] if Nt else self.model.hist_t
    y_unbiased = self._recover_unbiased_solution(pb.hist_t, pb.hist, t_model, y_model)
    return y_unbiased, y_model

analysis_step(d, Cdd, return_analysis=False)

Bayesian analysis step.

Builds the augmented forecast ensemble, delegates the filter update to _analysis_kernel, applies inflation, validates parameters, and updates the model history in-place.

Parameters:

Name	Type	Description	Default
d	ndarray(Nq)	Observation vector.	required
Cdd	ndarray(Nq, Nq)	Observation noise covariance.	required
return_analysis	bool	If True, return the analysed state array.	False

Source code in src/estimators/ensembles.py

def analysis_step(
    self,
    d: np.ndarray,
    Cdd: np.ndarray,
    return_analysis: bool = False,
) -> np.ndarray | None:
    """Bayesian analysis step.

    Builds the augmented forecast ensemble, delegates the filter update to
    ``_analysis_kernel``, applies inflation, validates parameters, and
    updates the model history in-place.

    Parameters
    ----------
    d : ndarray (Nq,)
        Observation vector.
    Cdd : ndarray (Nq, Nq)
        Observation noise covariance.
    return_analysis : bool
        If True, return the analysed state array.
    """
    assert self.rng is not None, "RNG not set; ensure __init__ completed."

    Af_state = self.current_state          # (Nphi+Na, m)

    # State-estimation-only warm-up, as in Ensemble.analysis_step.
    if self.start_param > 0:
        self.activate_parameter_estimation = (
            len(self.assimilated_data.times) >= self.start_param
        )

    # Optionally exclude parameters from the analysis update. The forecast
    # parameters are kept aside and re-appended below, so that the analysis
    # written to history always has the model's full Nphi+Na rows.
    Af_params = None
    if self.Na > 0 and not self.activate_parameter_estimation:
        Af_params = Af_state[self.Nphi:self.Nphi + self.Na, :].copy()
        Af_state = Af_state[:self.Nphi, :]

    # Append observables to form augmented state.
    y = self.model.get_observables()       # (Nq, m)
    Af_aug = np.vstack((Af_state, y))      # (Nphi+Na+Nq, m)  or  (Nphi+Nq, m)

    # ── Call the subclass filter kernel ──────────────────────────────────


    Aa = self._analysis_kernel(Af_aug, d, Cdd)


    # ── Covariance inflation (Evensen 2009, Chap. 15) ────────────────────
    rho = self.inflation_factor
    if rho != 1.0:
        Aa = multiplicative_inflation(Aa, rho)
        self.inflation_history.times.append(self.current_time)
        self.inflation_history.factors.append(rho)

    # ── Spread and parameter validity checks ─────────────────────────────
    if not self.has_valid_spread(Aa[:self.model.Nphi, :]):
        self.rejected_analysis = (self.current_time, 'Invalid analysis spread')

    if self.Na > 0 and self.alpha_limits_matrix is not None:
        Aa_alpha = Aa[self.Nphi:self.Nphi + self.Na, :]
        is_physical, idx_alpha, _ = self.has_valid_params(
            Aa_alpha, self.alpha_limits_matrix, get_deltas=False
        )
        if not is_physical:
            self.rejected_analysis = (
                self.current_time,
                f'Non-physical parameters at indices {idx_alpha}',
            )
            Aa = multiplicative_inflation(Af_aug, self.inflation_factor_rejection)

    # An EnKF-type gain corrects each state row from its own covariance with y,
    # independent of what other rows are present -- so leaving the reservoir out
    # of the analysis is the same as running it and discarding its correction.
    if not self.update_reservoir:
        N_units = getattr(self.model, 'N_units', 0)
        if N_units > 0:
            N_dim = self.model.N_dim
            Aa[N_dim:N_dim + N_units, :] = Af_aug[N_dim:N_dim + N_units, :]

    # ── Update model history ─────────────────────────────────────────────
    # Only the state rows go to the model; the observables are derived.
    Aa_psi = Aa[:Af_state.shape[0], :]
    if Af_params is not None:
        Aa_psi = np.vstack((Aa_psi, Af_params))   # parameters left at their forecast

    self.update_history(
        Aa_psi,
        self.current_time,
        modify_saved_states=True,
    )
    self.assimilated_data = (d, self.current_time)

    # ── Update the bias estimator with the analysis innovation ───────────
    # i^a = d - y^a, evaluated on the *analysed* observables. Ported from
    # Ensemble.update_history_analysis; smooth's EnsembleEstimator dropped
    # this, which left the bias estimator frozen across analysis steps.
    if self.bias is not None:
        Ya = self.model.get_observables()               # (Nq, m)
        ia = d[:, np.newaxis] - Ya                      # (Nq, m)
        updated_state = self.bias.update_state_from_innovation(ia)
        self.bias.update_history(
            updated_state, t=self.bias.current_time, modify_saved_states=True
        )

    if return_analysis:
        # The analysed *state* rows, matching model.hist and the convention in
        # DeterministicEstimator.analysis_step. The trailing rows of Aa are the
        # augmented observables, which are derived rather than stored.
        return Aa_psi

reshape_ensemble(m=None, reset=True)

Re-perturb the current ensemble around its mean.

Source code in src/estimators/ensembles.py

def reshape_ensemble(self, m: int | None = None, reset: bool = True) -> None:
    """Re-perturb the current ensemble around its mean."""
    pm = self.model
    if m is None:
        m = self.m
    if m == 1:
        raise ValueError('Ensemble size m must be greater than 1.')
    current_psi = pm.current_state
    mean_psi = np.mean(current_psi, axis=-1)
    std_psi = np.std(current_psi, axis=-1)
    new_ensemble = mean_vector_to_ensemble(pm.rng, mean_psi, std_psi, m, method='normal')
    pm.update_history(psi=new_ensemble, t=pm.current_time, reset=reset)

inflate(A, rho, d=None, additive=True) staticmethod

Deprecated alias for romda.estimators.inflation.multiplicative_inflation.

Source code in src/estimators/ensembles.py

@staticmethod
def inflate(
    A: np.ndarray,
    rho: float,
    d: np.ndarray | None = None,
    additive: bool = True,
) -> np.ndarray:
    """Deprecated alias for `romda.estimators.inflation.multiplicative_inflation`."""
    return multiplicative_inflation(A, rho)

has_valid_spread(A, tol=1e-06) staticmethod

Return True if ensemble spread is non-degenerate.

Source code in src/estimators/ensembles.py

@staticmethod
def has_valid_spread(A: np.ndarray, tol: float = 1e-6) -> bool:
    """Return True if ensemble spread is non-degenerate."""
    return True  # placeholder — detailed check commented out in original

has_valid_params(A_alpha, alpha_limits_matrix, get_deltas=False) staticmethod

Check whether all parameter ensemble members lie within bounds.

Parameters:

Name	Type	Description	Default
A_alpha	(Na, m)		required
alpha_limits_matrix	(2, Na, 1) rows are [lower_bounds, upper_bounds]		required

Returns:

Name	Type	Description
is_physical	bool
idx_alpha	list of out-of-bounds parameter indices
d_alpha	allowed values array (only populated if get_deltas=True)

Source code in src/estimators/ensembles.py

@staticmethod
def has_valid_params(
    A_alpha: np.ndarray,
    alpha_limits_matrix: np.ndarray,
    get_deltas: bool = False,
) -> tuple[bool, list[int], np.ndarray]:
    """Check whether all parameter ensemble members lie within bounds.

    Parameters
    ----------
    A_alpha : (Na, m)
    alpha_limits_matrix : (2, Na, 1)  rows are [lower_bounds, upper_bounds]

    Returns
    -------
    is_physical : bool
    idx_alpha   : list of out-of-bounds parameter indices
    d_alpha     : allowed values array (only populated if get_deltas=True)
    """
    if alpha_limits_matrix is None:
        return True, [], np.array([])

    low_limits, high_limits = alpha_limits_matrix
    below = A_alpha < low_limits
    above = A_alpha > high_limits
    oob = np.any(below | above, axis=1)
    is_physical = not np.any(oob)
    idx_alpha = np.where(oob)[0].tolist()
    d_alpha: list[np.ndarray] = []

    if get_deltas and not is_physical:
        if np.any(above) and np.any(below):
            raise ValueError(
                'Both above and below limits detected simultaneously; '
                'check alpha_limits and A_alpha.'
            )
        allowed = A_alpha[~above & ~below]
        if np.any(below):
            d_alpha.append(
                np.max(allowed, axis=1) if allowed.size > 0 else low_limits[:, 0]
            )
        elif np.any(above):
            d_alpha.append(
                np.min(allowed, axis=1) if allowed.size > 0 else high_limits[:, 0]
            )

    return is_physical, idx_alpha, np.array(d_alpha)

visualize_history(**kwargs)

Plot observable and parameter histories.

Source code in src/estimators/ensembles.py

def visualize_history(self, **kwargs) -> None:
    """Plot observable and parameter histories."""
    kwargs_obs = allowed_kwargs_for_func(plot_observable_history, kwargs)
    plot_observable_history(ensemble=self, **kwargs_obs)
    if self.Na > 0:
        kwargs_alpha = allowed_kwargs_for_func(plot_alpha_history, kwargs)
        plot_alpha_history(ensemble=self, **kwargs_alpha)

visualize_state(**kwargs)

Plot ensemble state distributions.

Source code in src/estimators/ensembles.py

def visualize_state(self, **kwargs) -> None:
    """Plot ensemble state distributions."""
    func = self.model.visualize_state
    kwargs_state = allowed_kwargs_for_func(func, kwargs)
    func(**kwargs_state)

romda.estimators.EnKF(parent_model, parent_bias=None, **kwargs)

Bases: EnsembleEstimator

Stochastic Ensemble Kalman Filter (perturbed-observation variant).

Each ensemble member assimilates a randomly perturbed copy of the observation, \(\mathbf{d}_j \sim \mathcal{N}(\mathbf{d}, \mathbf{C}_{dd})\). Writing \(\boldsymbol{\Psi}^\mathrm{f} = \mathbf{A}^\mathrm{f} - \overline{\mathbf{A}^\mathrm{f}}\) for the mean-subtracted forecast ensemble and \(\mathbf{S} = \mathbf{M}\boldsymbol{\Psi}^\mathrm{f}\):

\[ \mathbf{C} = (m-1)\,\mathbf{C}_{dd} + \mathbf{S}\mathbf{S}^\mathrm{T}, \qquad \mathbf{D} = \mathbf{d}\mathbf{1}^\mathrm{T} + \mathrm{chol}(\mathbf{C}_{dd})\,\boldsymbol{\varepsilon}, \quad \boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0}, \mathbb{I}), \]

\[ \mathbf{A}^\mathrm{a} = \mathbf{A}^\mathrm{f} + \boldsymbol{\Psi}^\mathrm{f}\mathbf{S}^\mathrm{T}\mathbf{C}^{-1} \left(\mathbf{D} - \mathbf{M}\mathbf{A}^\mathrm{f}\right). \]

Notes

Algebraically equivalent to the textbook Kalman-gain form \(\mathbf{K} = \mathbf{C}_{\psi\psi}\mathbf{C}_{yy}^{-1}\) with \(\mathbf{C}_{\psi\psi} = \boldsymbol{\Psi}^\mathrm{f}\mathbf{S}^\mathrm{T}/(m{-}1)\), \(\mathbf{C}_{yy} = \mathbf{S}\mathbf{S}^\mathrm{T}/(m{-}1) + \mathbf{C}_{dd}\), but rearranged to avoid forming the covariance matrices explicitly. The implementation equivalently expresses the update as a member-space transform \(\mathbf{A}^\mathrm{a} = \mathbf{A}^\mathrm{f}(\mathbb{I}_m + \mathbf{X})\), \(\mathbf{X} = \mathbf{S}^\mathrm{T}\mathbf{C}^{-1}(\mathbf{D}-\mathbf{M}\mathbf{A}^\mathrm{f})\), which coincides with the equation above because \(\mathbf{S}\) has zero column-sum.

References

Evensen (2003). The ensemble Kalman filter: theoretical formulation and practical implementation. Ocean Dynamics, 53, 343-367, Eq. (9.27).

Source code in src/estimators/ensembles.py

@typechecked
def __init__(
    self,
    parent_model: Model | type[Model],
    parent_bias: Bias | type[Bias] | None = None,
    **kwargs,
):
    """Initialise ensemble estimator.

    Parameters
    ----------
    parent_model : Model instance or Model subclass
        If a class is passed, remaining *kwargs* are forwarded to its
        constructor (e.g. ``dt``, ``psi0``, model-specific parameters).
    parent_bias : Bias instance, Bias subclass, or None
    **kwargs
        Ensemble config keys (``m``, ``std_phi``, ``std_alpha``, …) and/or
        model constructor keys are accepted here.
    """
    #  Apply general DA attributes from Estimator base.
    super().__init__(**kwargs)

    # Split the remaining kwargs: ensemble-generation keys go to
    # init_ensemble; everything else (dt, psi0, model parameters, ...) to
    # the model constructor. Keys consumed by Estimator.__init__ are
    # forwarded to neither.
    ensemble_kwargs = allowed_kwargs_for_func(parent_model.init_ensemble, kwargs)
    model_kwargs = {k: v for k, v in kwargs.items()
                    if k not in self._consumed_kwargs and k not in ensemble_kwargs}

    # Instantiate or copy the model.
    if isinstance(parent_model, Model):
        parent_model = parent_model.copy()
    else:
        parent_model = parent_model(**model_kwargs)

    # Generate the initial ensemble inside the model, falling back to the
    # estimator's config for keys not passed explicitly.
    for attr in ('std_phi', 'std_alpha', 'distribution_phi', 'distribution_alpha'):
        if attr not in ensemble_kwargs:
            ensemble_kwargs[attr] = getattr(self, attr)
    parent_model.init_ensemble(**ensemble_kwargs)

    self._model = parent_model

    self._init_bias(parent_bias)

romda.estimators.EnSRKF(parent_model, parent_bias=None, **kwargs)

Bases: EnsembleEstimator

Ensemble Square-Root Kalman Filter (no observation perturbations).

Updates the ensemble mean with the Kalman gain and transforms the ensemble deviations with a symmetric square-root transform, so no stochastic observation perturbations are needed. Writing \(\boldsymbol{\Psi}^\mathrm{f}\) for the mean-subtracted forecast ensemble and \(\mathbf{S} = \mathbf{M}\boldsymbol{\Psi}^\mathrm{f}\):

\[ \mathbf{C} = (m-1)\,\mathbf{C}_{dd} + \mathbf{S}\mathbf{S}^\mathrm{T}, \qquad \mathbf{K} = \boldsymbol{\Psi}^\mathrm{f}\mathbf{S}^\mathrm{T}\mathbf{C}^{-1}, \]

\[ \overline{\boldsymbol{\psi}}^\mathrm{a} = \overline{\boldsymbol{\psi}}^\mathrm{f} + \mathbf{K}\left(\mathbf{d} - \mathbf{M}\overline{\boldsymbol{\psi}}^\mathrm{f}\right), \qquad \boldsymbol{\Psi}^\mathrm{a} = \boldsymbol{\Psi}^\mathrm{f}\,\mathbf{T}^{1/2}, \qquad \mathbf{T} = \mathbb{I}_m - \mathbf{S}^\mathrm{T}\mathbf{C}^{-1}\mathbf{S}, \]

with \(\mathbf{T}^{1/2}\) the symmetric square root of \(\mathbf{T}\) (computed via its eigendecomposition, since \(\mathbf{T}\) is symmetric).

References

Tippett, Anderson, Bishop, Hamill & Whitaker (2003). Ensemble square root filters. Mon. Wea. Rev., 131, 1485-1490.

Source code in src/estimators/ensembles.py

@typechecked
def __init__(
    self,
    parent_model: Model | type[Model],
    parent_bias: Bias | type[Bias] | None = None,
    **kwargs,
):
    """Initialise ensemble estimator.

    Parameters
    ----------
    parent_model : Model instance or Model subclass
        If a class is passed, remaining *kwargs* are forwarded to its
        constructor (e.g. ``dt``, ``psi0``, model-specific parameters).
    parent_bias : Bias instance, Bias subclass, or None
    **kwargs
        Ensemble config keys (``m``, ``std_phi``, ``std_alpha``, …) and/or
        model constructor keys are accepted here.
    """
    #  Apply general DA attributes from Estimator base.
    super().__init__(**kwargs)

    # Split the remaining kwargs: ensemble-generation keys go to
    # init_ensemble; everything else (dt, psi0, model parameters, ...) to
    # the model constructor. Keys consumed by Estimator.__init__ are
    # forwarded to neither.
    ensemble_kwargs = allowed_kwargs_for_func(parent_model.init_ensemble, kwargs)
    model_kwargs = {k: v for k, v in kwargs.items()
                    if k not in self._consumed_kwargs and k not in ensemble_kwargs}

    # Instantiate or copy the model.
    if isinstance(parent_model, Model):
        parent_model = parent_model.copy()
    else:
        parent_model = parent_model(**model_kwargs)

    # Generate the initial ensemble inside the model, falling back to the
    # estimator's config for keys not passed explicitly.
    for attr in ('std_phi', 'std_alpha', 'distribution_phi', 'distribution_alpha'):
        if attr not in ensemble_kwargs:
            ensemble_kwargs[attr] = getattr(self, attr)
    parent_model.init_ensemble(**ensemble_kwargs)

    self._model = parent_model

    self._init_bias(parent_bias)

romda.estimators.rBA_EnKF(parent_model, parent_bias=None, gamma=1.0, **kwargs)

Bases: EnsembleEstimator

Regularized bias-aware ensemble Kalman filter (r-EnKF).

Extends the stochastic EnKF with an explicit correction for the (estimated) observation bias \(\mathbf{b}\) and its Jacobian \(\mathbf{J} = \mathrm{d}\mathbf{b}/\mathrm{d}(\mathbf{M}\boldsymbol{\psi})\), weighted by the regularization factor \(\gamma \ge 0\) (\(\gamma=0\) recovers the standard EnKF). Writing \(\mathbf{Y} = \mathbf{M}\mathbf{A}^\mathrm{f} + \mathbf{B}\) for the bias-corrected forecast observables, and \(\mathbf{C}_{\psi q}\), \(\mathbf{C}_{qq}\) for the (sample) forecast cross- and auto-covariances of the state and the mapped observables:

\[ \mathbf{C}_{yy} = (m-1)\,\mathbf{C}_{dd} + (\mathbb{I}+\mathbf{J})^\mathrm{T}(\mathbb{I}+\mathbf{J})\,\mathbf{C}_{qq} + \gamma\, \mathbf{J}^\mathrm{T}\mathbf{J}\,\mathbf{C}_{qq}, \qquad \mathbf{K} = \mathbf{C}_{\psi q}\,\mathbf{C}_{yy}^{-1}, \]

\[ \mathbf{A}^\mathrm{a} = \mathbf{A}^\mathrm{f} + \mathbf{K} \left[(\mathbb{I}+\mathbf{J})^\mathrm{T}(\mathbf{D} - \mathbf{Y}) - \gamma\, \mathbf{J}^\mathrm{T}\mathbf{b}\right], \]

where \(\mathbf{D}\) is the perturbed-observation ensemble and \(\mathbf{b}\) the current (mean) bias estimate, broadcast over the ensemble. The bias covariance is fixed to \(\mathbf{C}_{bb} = \mathbf{C}_{dd}\), so the weight \(\mathbf{W} = \mathbf{C}_{dd}\mathbf{C}_{bb}^{-1}\) of the reference paper reduces to the identity and is omitted.

Notes

\(\mathbf{C}_{\psi q}\) and \(\mathbf{C}_{qq}\) here are the sample cross- and auto-covariances (divided by \(m-1\)), unlike the unnormalized cross-moments used by the legacy romda.legacy.data_assimilation.rBA_EnKF implementation of the same filter. If observations are biased (Bias.biased_observations), \(\mathbf{d}\) is shifted by the mean bias-innovation gap before the update. During the first start_bias analysis steps the plain EnKF update is applied instead (bias-blind warm-up window).

Parameters:

Name	Type	Description	Default
gamma	float	Bias-regularization factor (default 1.0), aliased as regularization_factor. Larger values give more weight to the bias norm in the cost function.	1.0

References

Nóvoa, Racca & Magri (2023). Inferring unknown unknowns: Regularized bias-aware ensemble Kalman filter. Comput. Methods Appl. Mech. Eng., 418, 116502. DOI: 10.1016/j.cma.2023.116502.

Source code in src/estimators/ensembles.py

def __init__(
    self,
    parent_model,
    parent_bias=None,
    gamma: float = 1.0,
    **kwargs,
):
    # regularization_factor is the canonical attribute; gamma is the alias.
    kwargs.setdefault('regularization_factor', gamma)
    super().__init__(parent_model, parent_bias, **kwargs)