MeteoModels.jl

struct AlgebraicModel{T,A<:AbstractMatrix{T}} <: LinearModel
  matrix::A
end

Standard implementation of a LinearModel. The field matrix represents the constant Jacobian of the model itself.

const DEnKF{A<:Model,B<:Model,E<:Law,F<:Law} = EnsembleKalmanFilter{A,B,<:Ensemble{DEnKFStrategy},<:Ensemble,E,F}

Implements the DEnKF. Simply requires a specialization of the update! function.

struct DEnKFStrategy <: EnsembleStyle end

Trait for ensembles mimicking the DEnKF (deterministic EnKF) method:

• run the forecast step on each ensemble member (see forecast!);
• compute the Kalman gain $K$ as usual (see kalman_gain!);
• compute the ensemble innovations $ỹ$ (see innovation!);
• update the ensemble mean according to the formula:

\[μ ← μ + K ⋅ mean(ỹ),\]

where $μ$ is the ensemble mean;

• update the ensemble anomaly according to the formula:

\[A ← (I - rac{1}{2} K⋅H) A \]

where $H$ is the jacobian of the observation model, evaluated in the forecasted ensemble mean, and $A$ is the ensemble anomaly. This is the so-called deterministic approximation of DEnKF;

• update the ensemble $E$ according to the formula:

\[E[:,i] = A[:,i] + μ \]

for every $i = 1,...,nₑ$. The covariance $P$ is then estimated via its ensemble counterpart:

\[P = A ⋅ Aᵀ / (nₑ - 1)\]

const EnKF{A<:Model,B<:Model,E<:Law,F<:Law} = EnsembleKalmanFilter{A,B,<:Ensemble{EnKFStrategy},<:Ensemble,E,F}

Implements the standard EnKF. Simply requires a specialization of the update! function.

struct EnKFStrategy <: EnsembleStyle end

Trait for ensembles mimicking the EnKF method:

• run the forecast step on each ensemble member (see forecast!);
• compute the Kalman gain $K$ as usual (see kalman_gain!);
• compute the ensemble innovations $ỹ$ (see innovation!);
• update the ensemble according to the formula:

\[E ← E + K ⋅ ỹ \]

where $E$ is the ensemble matrix. The mean $μ$ and covariance $P$ are then estimated via their ensemble counterparts:

\[μ = ∑ᵢ E[:,i] / nₑ P = ∑ᵢ (E[:,i] - μ)⋅(E[:,i] - μ)ᵀ / (nₑ - 1)\]

struct Ensemble{C<:EnsembleStyle,A<:AbstractMatrix,B<:AbstractVector,D<:AbstractMatrix} <: SecondMoment{B}
  values::A
  mean::B 
  covariance::D
  anomaly::A
  strategy::C
end

This SecondMoment distribution represents the ensemble needed to run the EnsembleKalmanFilter. Fields:

• values: $n × n_e$-dimensional matrix storing the ensemble values;
• mean: $n$-dimensional vector representing the ensemble mean;
• covariance: $n × n$-dimensional matrix representing the ensemble covariance;
• anomaly: $n × n_e$-dimensional matrix storing the ensemble anomaly;
• strategy: trait of type EnsembleStyle which determines the update type of the ensemble.

const EnsembleKalmanFilter{A<:Model,B<:Model,C<:Ensemble,D<:Ensemble,E<:Law,F<:Law} = GenericKalmanFilter{A,B,C,D,E,F}

Implements an Ensemble Kalman Filter. In particular:

• instead of propagating a single probability distribution as in a Kalman filter, we do so for several different (ensemble) distributions.
• the explicit update of the state covariance matrix is not required. Indeed, the variability of the state is implicitly encoded in the ensemble's spread.
• the remaining steps are equivalent to a standard Kalman Filter.

Subtypes:

• EnKF
• DEnKF

abstract type EnsembleStyle end

Trait specifying how the ensemble covariance of an Ensemble distribution should be updated. The reason why this is kept as a parameter is that, in ensemble filtering strategies, computing the ensemble covariance with the usual formula

\[P = ∑ᵢ (E[:,i] - μ)*(E[:,i] - μ)ᵀ / (nₑ - 1)\]

is generally expensive, and thus alternative strategies are sought. Here, $E$ are the ensemble members. Subtypes:

• EnKFStrategy
• DEnKFStrategy

struct ExtendedKalmanFilter{A<:Law} <: KalmanFilter{A}
  filter::KalmanFilter{A}
  cache::ExtendedKalmanCache
end

Implements the Extended Kalman Filter (EKF). In particular:

• for the forecasting step: we linearise the transition operator around the analysed state of the

previous iteration;

• for the analysis step: we linearise the observation operator around the forecasted state of the

current iteration. The remaining scheme is equivalent to that of a standard Kalman Filter.

abstract type Filter end

Type reserved for filter operators for data processes, such as the Kalman filter (KF). A filter should be given as input an transition operator, a observation operator, and a prior distribution, which should be representative of the distribution of the initial sample. The prior should be updated iteratively, assimilating observations at multiple different time instants, as well as the last bayesian estimate of the (unknown) state variable. Check the function loop for suitable iterative updates of the prior. Check here for more details on Kalman filtering. Subtypes:

• KalmanFilter

const FirstMoment{V} = Law{1,V}

Type reserved for distributions characterised only by their first moment, i.e. the mean, accessed via the function mean.

struct FourDVarCache
  Bfact::Factorization
  Rfact::Factorization
end

Cache for FourDVarFilter holding precomputed inverse covariance matrices.

Fields:

• Bfact: inverse of the background error covariance matrix $B^{-1}$
• Rfact: inverse of the observation error covariance matrix $R^{-1}$

struct FourDVarFilter{A<:Model,B<:Model,C<:SecondMoment,D<:SecondMoment,E<:Law} <: Filter

Filter implementing the four-dimensional variational data assimilation (4DVar) method. Over each assimilation window of $T$ time steps, finds the initial state $x_0$ minimising the cost function

\[J(x_0) = \frac{1}{2}(x_0 - x^b)^\top B^{-1}(x_0 - x^b) + \frac{1}{2}\sum_{k=1}^{T}\bigl(y_k - H(M^k(x_0))\bigr)^\top R^{-1} \bigl(y_k - H(M^k(x_0))\bigr)\]

where $x^b$ is the background state with error covariance $B$,$M$ is the one-step transition model,$H$ is the observation operator,and $R$ is the observation error covariance. Gradients are computed via forward-mode automatic differentiation; $M$ and $H$ must therefore be written as generic Julia functions (no hard-coded Float64 types).

The background for each successive window is the analysed initial condition of the previous window (restart-from-analysis cycling).

Fields:

• transition: Model representing the one-step transition operator $M$;
• observation: Model representing the observation operator $H$;
• prior: SecondMoment distribution for the background state $(x^b,B)$;
• obs_prior: SecondMoment distribution in observation space (used for dimension queries);
• obs_noise: Law with zero mean and covariance $R$;
• cache: FourDVarCache holding precomputed $B^{-1}$ and $R^{-1}$.

abstract type FunctionFilter <: Filter end

Subtype reserved for filters whose transition and observation models are time-dependent functions. In practice, the only difference with a standard Filter is that the models of a FunctionFilter must be explicitly evaluated at each iteration before running the Kalman iterations (see loop for more details).

struct GenericFirstMoment{A<:AbstractVector} <: FirstMoment{A}
  mean::A 
end

Most basic implementation of a FirstMoment distribution.

struct GenericFunctionKalmanFilter{A<:Function,B<:Function,C<:Law,D<:Law,E<:Law,F<:Law} <: FunctionKalmanFilter{C}
  transition::A 
  observation::B
  prior::C
  obs_prior::D
  noise::E 
  obs_noise::F
  cache::KalmanCache
end

Filter subtype implementing a Kalman filter procedure. Fields:

• transition: Real -> Model function representing the transition operator. The real input it receives

could be, for example, the time instant of the current Kalman iteration. This field should be evaluated at each iteration to successfully run the Kalman iterations, e.g. via loop;

• observation: Real -> Model function representing the observation operator. The real input it receives

could be, for example, the time instant of the current Kalman iteration. This field should be evaluated at each iteration to successfully run the Kalman iterations, e.g. via loop;

• prior: Law representing the probability distribution for the state;
• obs_prior: Law representing the probability distribution for the observation;
• noise: Law representing the probability distribution for the process (state) noise;
• obs_noise: Law representing the probability distribution for the observation noise;
• cache: cached object allowing for efficient in-place operations.

struct GenericKalmanFilter{A<:Model,B<:Model,C<:Law,D<:Law,E<:Law,F<:Law} <: KalmanFilter{C}
  transition::A 
  observation::B
  prior::C
  obs_prior::D
  noise::E 
  obs_noise::F
  cache::KalmanCache
end

Filter subtype implementing a Kalman filter procedure. Fields:

• transition: Model representing the transition operator;
• observation: Model representing the observation operator;
• prior: Law representing the probability distribution for the state;
• obs_prior: Law representing the probability distribution for the observation;
• noise: Law representing the probability distribution for the process (state) noise;
• obs_noise: Law representing the probability distribution for the observation noise;
• cache: cached object allowing for efficient in-place operations.

struct GenericModel{F<:FType} <: NonlinearModel
  form::F
end

Standard implementation of a NonlinearModel. The field form represents the function or Gridap Map characterising the model itself.

const InflationDEnKF = InflationKalmanFilter{<:Ensemble{DEnKFStrategy},<:MultInflationParam}

struct InflationKalmanFilter{A<:Ensemble,B<:InflationParameter} <: KalmanFilter{A}
  filter::KalmanFilter{A}
  inflation_param::B
  cache
end

const JointLaw{N,V<:BlockVector} = Law{N,V}

Type representing a joint probability distribution characterised by N moments.

abstract type Law{N,V} end

Type representing a probability distribution characterised by N moments, and values of type V. Subtypes:

• FirstMoment
• SecondMoment

struct Linear <: Linearity end

Trait used by LinearModel.

const LinearModel = Model{Linear}

Models that are fully characterised by an $m × n$-dimensional Jacobian matrix $J$, i.e.

\[x ↦ J⋅x\]

represents the action of a LinearModel on an $n$-dimensional vector $x$. If the input is a distribution $d$, then:

• if $d$ is a FirstMoment distribution with mean $μ$, then the output is a FirstMoment

with mean $J⋅μ$;

• $d$ is a SecondMoment distribution with mean $μ$ and covariance $P$, then the output

is a SecondMoment with mean $J⋅μ$, and covariance $J⋅P⋅Jᵀ$.

abstract type Linearity end

A Model trait that allows certain optimizations during the evaluation of the models themselves.

abstract type Model{A<:Linearity} <: Map end

Type used for operator-like quantities, such as functions or Gridap Maps. For performance reasons, we distinguish models depending on their Linearity trait. To evaluate a Model a in a point $x$ — usually an $n$-dimensional vector or a scalar — simply call

a(x)

which automatically triggers the function evaluate(a,x), originally defined in Gridap, and overwritten here a few times. For in-place evaluations, use instead the syntax

evaluate!(cache,a,x)

where cache = return_cache(a,x) is a suitable cached object.

The main characteristic of a Model is that it may also be evaluated in a probability distribution. Given an input Law prior, the output

posteriori = a(priori)

returns another distribution posteriori, which should be thought of the propagation of priori through the model a. The type of Model and input Law determine the expression of posterior.

const NLLInflationEnKF = NLLInflationKalmanFilter{<:Ensemble{EnKFStrategy}}

const NLLInflationKalmanFilter{A<:Ensemble} = InflationKalmanFilter{A,<:NLLInflationParam}

struct Nonlinear <: Linearity end

Trait used by NonlinearModel.

const NonlinearModel = Model{Nonlinear}

Models that are in general characterised by either a function, or a Gridap Map. Denoting such function/Map by by f, the action of a NonlinearModel on an $n$-dimensional vector $x$ is simply defined by

\[x ↦ f(x)\]

where the output is an $m$-dimensional vector, or a scalar. If the input is a distribution $d$, then:

• if $d$ is a FirstMoment distribution with mean $μ$, then the output is a FirstMoment

with mean f(μ);

• $d$ is a SecondMoment distribution with mean $μ$ and covariance $P$, then the output

is a SecondMoment with mean f(μ), and covariance whose definition depends on the types of boht a and $d$.

struct NormalLaw{A<:AbstractVector,B<:AbstractMatrix} <: SecondMoment{A}
  mean::A 
  covariance::B
end

Normal distribution of mean mean and covariance covariance; a SecondMoment distribution defaults to a NormalLaw.

const SecondMoment{V} = Law{2,V}

Type reserved for distributions characterised by their first two moments, i.e. mean and covariance, accessed via the functions mean and cov.

struct SigmaPoints{A<:AbstractVector,B<:AbstractMatrix,C<:AbstractMatrix,D<:AbstractVector,E<:Real} <: SecondMoment{A}
  mean::A 
  covariance::B
  points::C 
  weights_mean::D 
  weights_cov::D
  λ::E
end

This SecondMoment distribution represents the sigma points needed to run the UnscentedTransform. Fields:

• points: $n × (2*L + 1)$-dimensional matrix storing the values of the sigma points;
• mean: $n$-dimensional vector representing the weighted mean of points;
• covariance: $n × n$-dimensional matrix representing the weighted covariance of points;
• weights_mean: $(2*L + 1)$-dimensional vector storing weights for mean;
• weights_cov: $(2*L + 1)$-dimensional vector storing weights for covariance;
• λ: real value used in the update of the sigma points points.

In an Unscented Transformation, two things occur in an iterative fashion:

• the field points is updated in-place via a call to sigma_points!;
• the fields mean and covariance are updated in-place via a call to update!.

struct UniformLaw{A<:AbstractVector,B<:AbstractMatrix} <: SecondMoment{A}
  mean::A 
  covariance::B
  lower_bound::A
  upper_bound::A
end

Uniform distribution of mean mean and covariance covariance.

const UnscentedTransform{A<:Model,B<:Model,E<:Law,F<:Law} = GenericKalmanFilter{A,B,<:SigmaPoints,<:SigmaPoints,E,F}

Implements the unscented transform. This model is analogous to a Kalman filter procedure for nonlinear transition/observation models, and dealing with the nonlinearities by means of the so-called sigma points (see SigmaPoints). These are a set of interpolation points that are used to approximate the first and second moments of a transitioned/observed distribution, according to a weighted linear combination. In particular:

• for the forecasting step: we update the expression of the sigma points based on the analysed sigma points

from the previous iteration; we evaluate the state prior distribution in the sigma points, and update its mean and covariance according to weighted linear combinations (see SigmaPoints);

• for the analysis step: we evaluate the obvervation prior distribution in the sigma points, and

update its mean and covariance according to the same weighted linear combinations; we compute innovation, Kalman gain, and state posterior distribution as in a standard Kalman Filter scheme. From an implementation standpoint, an UnscentedTransform simply requires the transition and observation priors to both be SigmaPoints.

NLL(true_values::AbstractVector,d::Law) -> Real

Computes the Negative Log Likelihood between the true data true_values and the first moment of the distribution d.

NLL(true_values::AbstractMatrix,history::AbstractVector{<:Law}) -> AbstractVector

Computes the Negative Log Likelihood between the true data true_values and history, the historical distributions obtained by running the Kalman iterations.

NRMSE(true_values::AbstractVector,d::Law) -> Real

Computes the Normalised Root Mean Square Error between the true data true_values and the first moment of the distribution d. This is equal to the Root Mean Square Error divided by the standard deviation of d.

NRMSE(true_values::AbstractMatrix,history::AbstractVector{<:Law}) -> Real

Computes the Normalised Root Mean Square Error between the true data true_values and history, the historical distributions obtained by running the Kalman iterations.

RMSE(true_values::AbstractVector,d::Law) -> Real

Computes the Root Mean Square Error between the true data true_values and the first moment of the distribution d.

RMSE(true_values::AbstractMatrix,history::AbstractVector{<:Law}) -> Real

Computes the Root Mean Square Error between the true data true_values and history, the historical distributions obtained by running the Kalman iterations.

analyse!(posterior::Law,f::Filter,args...) -> Law

In-place execution of the analysis step of a Kalman filter algorithm. This step consists of the following operations:

• Running the observation model on the forecasted distribution (see observation!)
• Computing the innovation (see innovation!)
• Computing the Kalman gain (see kalman_gain!)
• Updating the forecasted distribution by accounting for the Kalman gain (see update!)

To run an iteration of the Kalman filter, one must run the forecasting step in forecast! prior to the analysis one. If no optional argument is provided, the analysis is not performed.

blockdiag(A::AbstractVector{<:AbstractMatrix{T}}) where T -> BlockMatrix{T}
blockdiag(A::AbstractMatrix{T}...) where T -> BlockMatrix{T}

Construct a BlockMatrix with A... on the diagonal and zeros elsewhere. All off-diagonal blocks exist and are filled with zeros (not empty).

dimension(d::Law) -> Int 
dimension(d::JointLaw) -> Vector{Int}

Dimension of the (vector) space on which the distribution is defined. For a JointLaw, the function returns a vector of integers corresponding to the dimensions of the each marginal.

draw(d::SecondMoment) -> AbstractVector 
draw(d::SecondMoment,nsamples::Int) -> AbstractMatrix

Draws a $n$-dimensional random vector from the distribution d, where $n$ represents the dimension of d (see distribution). If an integer nsamples is also provided, the output will be an $n × nsamples$ - dimensional matrix.

forecast!(posterior::Law,f::Filter) -> Law

In-place execution of the forecast step of a Kalman filter algorithm. This step consists of the following operations:

• Running the transition model on the prior distribution (see transition!)

To complete a single iteration of the Kalman filter, one must run the analysis step in analyse! following the forecast one.

get_noise(f::Filter) -> Law

Fetches the transition noise from the filter f. This noise, denoted by $θ$, is such that

\[xₙ₊₁ := F(xₙ) + θ\]

where ${xₖ}ₖ$ is the state process, and $F$ is the transition model.

get_observation_model(f::Filter) -> Model

Fetches the observation model from the filter f. This model, denoted by $H$, is such that

\[yₙ := H(xₙ) + η\]

where ${xₖ}ₖ$ is the state process, ${yₖ}ₖ$ is the observed process, and $η$ is a (usually Gaussian) random variable.

get_observation_noise(f::Filter) -> Law

Fetches the observation noise from the filter f. This noise, denoted by $η$, is such that

\[yₙ := H(xₙ) + η\]

where ${xₖ}ₖ$ is the state process, ${yₖ}ₖ$ is the observed process, and $H$ is the observation model.

get_observation_prior(f::Filter) -> Law

Fetches the distribution of the observation variable from the filter f.

get_prior(f::Filter) -> Law

Fetches the distribution of the state variable from the filter f.

get_transition_model(f::Filter) -> Model

Fetches the transition model from the filter f. This model, denoted by $F$, is such that

\[xₙ₊₁ := F(xₙ) + θ\]

where ${xₖ}ₖ$ is the state process, and $θ$ is a (usually Gaussian) random variable.

innovation!(f::Filter,posterior::Law,z::InType) -> InType

Given a distribution posterior and an observation z, returns the innovation $ỹ$ such that

\[ỹ = z - yₙ = z - [H(xᶠₙ) + η]\]

where $yₙ$ represents the observation forecasted by the filter f, and $xᶠₙ$ is distributed according to posterior.

jac(a::Model,x::InType) -> AbstractMatrix

Returns a's Jacobian matrix evaluated in $x$.

joint_law(d::AbstractVector{<:Law}) -> Law

Given a list of marginal distributions d = (d1,...,dn), returns their joint distribution.

kalman_gain!(f::Filter,posterior::Law) -> AbstractMatrix

In-place computation of the Kalman gain $K$ according to the formula

\[K = Pxy ⋅ Pyy⁻¹\]

where $Pxy$ and $Pyy$ are the state-observation and observation covariance matrices, respectively. $K$ is storede as a cached object in f, while the covariances $Pxy$ and $Pyy$ can be computed by suitably accessing the transition and observation distributions via get_prior and get_observation_prior, respectively.

linearise(a::Model,x::InType) -> LinearModel

Linearizes a model a around $x$. If a is a LinearModel, it returns a itself.

loop(f::Filter,obs::AbstractArray) -> AbstractVector{<:Law}
loop(f::Filter,obs::AbstractArray,grid::AbstractVector,fine_grid::AbstractVector) -> AbstractVector{<:Law}

Given a filter f and a list of observations obs, iteratively runs the forecast-analyse paradigm typical of a Kalman filter, producing a list of posterior distributions for the state variable. In practice, one iteration of the loop consists of one call to forecast!, followed by one to analyse!. The posterior resulting from each analysis is then fed as the prior distribution to the next forecast step.

Two additional vectors, grid and fine_grid, can also be provided and are interpreted as follows:

• grid: a grid of time steps at which the observations obs are recorded;
• fine_grid: a finer grid that fully contains grid.

In this case, the filter is run on the fine_grid, while analyse! is only executed at the time steps corresponding to grid. This setup is intended to simulate scenarios in which observations are available only at selected time steps.

mixed_cov!(P::AbstractMatrix,f::Filter,posterior::Law) -> AbstractMatrix

In-place computation of the state-observation "mixed" covariance P.

mixed_cov!(cache,a::SigmaPoints,b::SigmaPoints) -> AbstractMatrix

In-place computation of the covariance between the SigmaPoints distributions a and b. These two distributions should have the same $L$ (i.e. the same number of sigma points) and $λ$ parameters, which also implies that they share the same mean/covariance weights. The formula used here is:

\[P = ∑ᵢ₌₁²ᴸ⁺¹ w[i] ⋅ (χᵃ[:,i] - μᵃ) ⋅ (χᵇ[:,i] - μᵇ)\]

where $χᵃ$ and $μᵃ$ are the sigma points and their mean for a, $χᵇ$ and $μᵇ$ are the sigma points and their mean for b, and $w$ are the covariance weights of either a or b.

observation!(f::Filter,posterior::Law) -> Law

In-place application of the observation model stored in f on the distribution posterior, which represents the posterior distribution of the state variable. In essence, denoting by $F$ the transition model, if posterior represents the forecasted distribution of the state variable xᶠₙ at the $n$th iteration, this step runs the observation model

\[yₙ := H(xᶠₙ) + η\]

overwriting the result yₙ in the distribution of the observed variable stored in f, and accessible through get_observation_prior. This function should be run during the analysis step in a Kalman filter algorithm.

sigma_points!(dcache::SigmaPoints,d::SigmaPoints;kwargs...) -> AbstractMatrix

In-place update of the sigma points according to the formula:

\[χₖ₊₁[:,1] = μₖ χₖ₊₁[:,2:L+1] = μₖ + √((L + λ)Pₖ) χₖ₊₁[:,L+2:2L+1] = μₖ - √((L + λ)Pₖ)\]

where $μₖ$ and $Pₖ$ are the previous mean and covariance fields. The output $χₖ₊₁$ overwrites the field points.

sigma_points(d::SecondMoment;kwargs...) -> AbstractMatrix

Given an input distribution d, computes the sigma points $χ$ according to the formula:

\[χ[:,1] = μ χ[:,2:L+1] = μ + √((L + λ)P) χ[:,L+2:2L+1] = μ - √((L + λ)P)\]

where μ and P are the mean and covariance of d, respectively. The variables $L$ and $λ$ may be passed as keyword arguments, and assume the following default values:

\[L = dimension(d) λ = 3 - L\]

similar_law(d::Law,dim=dimension(d)) -> Law

Returns a distribution of same type as d, with a possibly different dimension specified by the optional argument dim.

transition!(posterior::Law,f::Filter) -> Law

In-place application of the transition model stored in f on the distribution posterior, which represents the posterior distribution of the state variable. In essence, denoting by $F$ the transition model, if posterior represents the distribution of the state variable $xₙ$ at the $n$th iteration, this step runs the transition model

\[xₙ₊₁ := F(xₙ) + θ\]

overwriting the result $xₙ₊₁$ in posterior. This function should be run during the forecast step in a Kalman filter algorithm.

update!(cache,d::Ensemble) -> Ensemble

Update the mean, anomaly and covariance of the ensemble d.

update!(cache,d::SigmaPoints) -> SigmaPoints

Update the mean and covariance of the sigma points d.

update!(posterior::Law,f::Filter) -> Law

In-place update of the distribution posterior through the action of Kalman gain matrix cached in f. Denoting by $K$ the Kalman gain computed by running kalman_gain!, and by $ỹ$ the innovation computed by running innovation!, if posterior represents the forecasted distribution of the state variable $xᶠₙ$ at the $n$th iteration, this step runs the formula

\[xᵃₙ := xᶠₙ + K ⋅ ỹ\]

and overwrites the analysed distribution of the state variable $xᵃₙ$ in posterior.

visualise(
  history::AbstractVector{<:Law},
  grid=eachindex(history);
  variable::Int=1,
  kwargs...
)

visualise(
  true_values::AbstractMatrix,
  history::AbstractVector{<:Law},
  grid=eachindex(history);
  variable::Int=1,
  kwargs...
)

Plot the historical distributions obtained by running the Kalman iterations, for e.g. via the function loop. The primary estimator is the mean of the distributions. If the distributions feature a second moment (i.e. they are equipped by a variance) then a confidence interval is also drawn. The true data true_values, if known, may be provided, and will be plotted on the same figure. The keyword variable – an integer – indicates the state member to be plotted.