stats⚓︎

Statistics for the assessment of DA methods.

Stats records per-timestep statistics during assimilation and exposes them as DACycleSeries attributes. The default statistics come in two shapes and sizes:

Vector stats (err, spread, mu, …) — one value per state dim. Each has child scalar field-summary (.rms, .m, .ma, …)
Scalar stats (mad, skew, kurt, trHK, …) — inherently scalar.

Avrgs holds the same hierarchy after time-averaging.

Default statistics⚓︎

For a tabular overview of all default statistics, see the Statistics section of the README.

`DACycleSeries` subscripts⚓︎

The time series of statistics (the attributes of .stats) are DACycleSeries objects. They may have sub-arrays .f, .a, .s, .i, referring to whether the statistic is for a "forecast", "analysis", "smoothed" (smoothers only), or "integrational" (forecast including intermediate, non-obs-time step) estimate, as determined by the fais argument passed to Stats.assess.

The same subscripts apply for the time-averages stored in .avrgs.

Field summaries⚓︎

The statistics are also averaged in space. This is done according to the methods listed in rc.field_summaries of the dpr_config.

Note

Although sometimes pretty close, rmv (a.k.a. spread.rms) is not (supposed to be) an un-biased estimator of rmse (a.k.a. err.rms). This is because of the square roots involved in the field summary. Instead, spread.ms (i.e. the mean variance) is the unbiased estimator of err.ms.

Regional field summaries⚓︎

If the HiddenMarkovModel has the attribute .sectors with value, e.g.,

HMM.sectors = {
    "ocean": inds_of_state_of_the_ocean,
    "atmos": inds_of_state_of_the_atmosphere,
}

then .stats.rms and .avrgs.rms will also have attributes named after the keys of HMM.sectors, e.g. stats.err.rms.ocean. This also goes for any other (than rms) type of field summary method.

Declaring new, custom statistics⚓︎

Only the time series created with Stats.new_series will be in the format operated on by Stats.average_in_time. For example, to add a full DA-cycle series for estimated inflation values:

self.stats.new_series('infl', 1)

For a statistic that only exists at analysis times (no .f or .i sub-arrays):

self.stats.new_series('infl', 1, analysis_only=True)

Alternatively you can overwrite a default statistic; for example:

error_time_series_a = xx - ensemble_time_series_a.mean(axis=1)
self.stats.err.rms.a = np.sqrt(np.mean(error_time_series_a**2, axis=-1))

To register a scalar (non-series) custom statistic so it appears in xp.avrgs, use self.stats.register:

self.stats.register("your_custom_stat", value)

This calls register internally, which sets the attribute and records "your_custom_stat" so that xp_launch.run_experiment (which deletes .stats after assimilation to save memory) can still copy it into xp.avrgs.

`Avrgs` ⚓︎

Bases: StatPrint, DotDict

Time-averaged statistics produced by stats.Stats.average_in_time.

Attributes get populated by Stats.average_in_time().

Vector stats time series (e.g. err, spread) first get reduced to scalar time series (stats.FieldAvrgs objects). Scalar time series (including those inherently so, e.g. mad, trHK) get reduced to stats.DACycleAvrgs objects.

Supports:

Deep attribute access: avrgs.err.rms.a
Dotted-key getattr: getattr(avrgs, "err.rms.a")
Named abbreviations (as properties): avrgs.rmse → avrgs.err.rms
Tabulation: avrgs.tabulate(["rmse.a", "rmv.a"])
Dict-like iteration and in tests (inherited from DotDict)

`N_eff` ⚓︎

Time-averaged effective ensemble size.

`crps` ⚓︎

Time-averaged CRPS field (and its scalar summaries).

`duration` ⚓︎

Wall-clock time (seconds) for the full assimilate() call.

`err` ⚓︎

Time-averaged error field (and its scalar summaries).

`gscore` ⚓︎

Time-averaged Gaussian score field (and its scalar summaries).

`infl` ⚓︎

Time-averaged inflation factor.

`iters` ⚓︎

Time-averaged iteration count.

`kurt` ⚓︎

Time-averaged excess kurtosis.

`mad` ⚓︎

Time-averaged mean absolute deviation.

`mu` ⚓︎

Time-averaged mean-estimate field (and its scalar summaries).

`resmpl` ⚓︎

Time-averaged resampling count.

`skew` ⚓︎

Time-averaged skewness.

`spread` ⚓︎

Time-averaged spread field (and its scalar summaries).

`svals` ⚓︎

Time-averaged singular values (and their scalar summaries).

`trHK` ⚓︎

Time-averaged trace of HK.

`umisf` ⚓︎

Time-averaged projected error (and its scalar summaries).

`w` ⚓︎

Time-averaged importance weights (and their scalar summaries).

`wroot` ⚓︎

Time-averaged weight-tempering root.

`getattr(key)` ⚓︎

Forward dotted-key lookups such as getattr(avrgs, "err.rms.a").

Called only when normal lookup fails. If key contains a dot, resolve it as a chain of attribute accesses (which lets the abbreviation properties — rmse, etc. — participate naturally). Plain missing keys raise AttributeError immediately to avoid infinite recursion.

`tabulate(statkeys=(), decimals=None)` ⚓︎

Tabulate using tabulate_avrgs.

`DACycleAvrgs` ⚓︎

Bases: StatPrint, DotDict

Time-averaged scalar stat: one tools.rounding.UncertainQtty per subscript.

Only subscripts for which finite data exist are set as attributes. Inherits DotDict so the object is iterable and supports in tests.

`a` ⚓︎

Analysis (posterior) time average.

`f` ⚓︎

Forecast (prior) time average.

`i` ⚓︎

Integrational (between-obs) time average.

`s` ⚓︎

Smoothed time average (smoothers only).

`FieldAvrgs` ⚓︎

Bases: StatPrint, DotDict

Time-averaged vector stat: one stats.DACycleAvrgs per field summary.

Only the field summaries that were actually computed are set as attributes. Inherits DotDict so the object is iterable and supports in tests.

`gm` ⚓︎

Geometric-mean average.

`m` ⚓︎

Mean-field average.

`ma` ⚓︎

Mean-absolute average.

`ms` ⚓︎

Mean-square average.

`rms` ⚓︎

Root-mean-square average.

`Stats` ⚓︎

Bases: StatPrint

Records and computes per-timestep statistics for a DA method.

DA methods may register additional stats via self.stats.register(name, value) for the purpose of automatic plotting and averaging.

`N_eff` ⚓︎

Effective ensemble size 1 / Σwᵢ².

`crps` ⚓︎

Continuous ranked probability score.

`duration` ⚓︎

Wall-clock time (seconds) for the full assimilate() call.

`err` ⚓︎

Error of the mean estimate vs. truth (mu - x).

`gscore` ⚓︎

Gaussian (log) score: 2 log(spread) + (err/spread)².

`infl` ⚓︎

Inflation factor applied at this analysis step.

`iters` ⚓︎

Number of iterations (iterative methods only).

`kurt` ⚓︎

Excess kurtosis of the ensemble (0 for Gaussians).

`mad` ⚓︎

Mean absolute deviation of the ensemble from its mean.

`mu` ⚓︎

Mean estimate (ensemble mean, or Kalman/variational estimate).

`resmpl` ⚓︎

Resampling flag/count at this analysis step.

`rh` ⚓︎

Rank histogram: rank of the truth among sorted ensemble members.

`skew` ⚓︎

Skewness of the ensemble.

`spread` ⚓︎

Ensemble spread (std dev), or Kalman posterior std dev.

`svals` ⚓︎

Singular values (principal-component scores) of the ensemble anomalies.

`trHK` ⚓︎

Trace of the observation-space gain matrix HK.

`umisf` ⚓︎

Error projected onto the leading ensemble/covariance directions.

`w` ⚓︎

Importance weights (particle filter; absent when weights are uniform).

`wroot` ⚓︎

Root-finding output for optimal weight tempering (particle filter).

`init(xp, HMM, xx, yy, liveplots=False, store_i=rc.store_i, LP_kwargs=None)` ⚓︎

Init the default statistics.

LP_kwargs overrides LivePlot pause timings (seconds; 0 disables). Per-experiment xp.LP_kwargs takes priority over this batch value.

`assess(k, ko=None, fais=None, E=None, w=None, mu=None, Cov=None)` ⚓︎

Common interface for both Stats.assess_ens and Stats.assess_ext.

The _ens assessment function gets called if E is not None, and _ext if mu is not None.

fais: One or more of ['f', 'a', 'i', 's'], indicating that the result should be stored in (respectively) the forecast/analysis/integrational/smoothed attribute. Default: 'i' if ko is None else 'ai' ('a' and 'i').

`assess_ens(now, x, E, w)` ⚓︎

Ensemble and Particle filter (weighted/importance) assessment.

`assess_ext(now, x, mu, P)` ⚓︎

Kalman filter (Gaussian) assessment.

`average_in_time(kk=None, kko=None, free=False)` ⚓︎

Average all scalar time series, producing xp.avrgs.

Parameters:

Name	Type	Description	Default
`kk`	`ndarray \| None`	Model-step indices to include when averaging the `i` (integrational) subscript. Defaults to `tseq.mask` (post-burnin).	`None`
`kko`	`ndarray \| None`	Obs-time indices for the `f`/`a`/`s` subscripts. Defaults to `tseq.masko`.	`None`
`free`	`bool`	If `True`, delete `xp.stats` after averaging to free memory.	`False`

`derivative_stats(now)` ⚓︎

Stats that derive from others, and are not specific for _ens or _ext).

`new_series(name, shape, analysis_only=False, field_mean=False, **kws)` ⚓︎

Create (and register) a statistics time series, initialized with nans.

Creates a tools.series.DACycleSeries. If analysis_only=True, the series has no .f sub-array and .i is a ring buffer of size 1.

NB: The sliding_diagnostics liveplotting relies on detecting nan's to avoid plotting stats that are not being used. Thus, you cannot use dtype=bool or int for stats that get plotted.

`replay(figlist='default', speed=np.inf, t1=0, t2=None, **kwargs)` ⚓︎

Replay LivePlot with what's been stored in 'self'.

t1, t2: time window to plot.
'figlist' and 'speed': See LivePlot's doc.

Note

store_i (whether to store non-obs-time stats) must have been True to have smooth graphs as in the actual LivePlot.

Note

Ensembles are generally not stored in the stats and so cannot be replayed.

`summarize_marginals(now)` ⚓︎

Compute Mean-field and RMS values.

`write(stat_dict, k, ko, sub)` ⚓︎

Write stat_dict to series at (k, ko, sub).

`align_col(col, pad='␣', missingval='', just='>')` ⚓︎

Align column.

Treats ints and fixed-point float/str especially, aligning on the point.

Examples:

>>> xx = [1, 1., 1.234, 12.34, 123.4, "1.2e-3", None, np.nan, "inf", (1, 2)]
>>> print(*align_col(xx), sep="\n")
␣␣1␣␣␣␣
␣␣1.0␣␣
␣␣1.234
␣12.34␣
123.4␣␣
␣1.2e-3
␣␣␣␣␣␣␣
␣␣␣␣nan
␣␣␣␣inf
␣(1, 2)

`center(E, axis=0, rescale=False)` ⚓︎

Center ensemble.

Makes use of np features: keepdims and broadcasting.

Parameters:

Name	Type	Description	Default
`E`	`ndarray`	Ensemble which going to be inflated	required
`axis`	`int`	The axis to be centered. Default: 0	`0`
`rescale`	`bool`	If True, inflate to compensate for reduction in the expected variance. The inflation factor is $\sqrt{\frac{N}{N - 1}}$ where N is the ensemble size. Default: False	`False`

Returns:

Name	Type	Description
`X`	`ndarray`	Ensemble anomaly
`x`	`ndarray`	Mean of the ensemble

`crps_ens(x, ensemble, weights=None)` ⚓︎

Compute CRPS for ensemble given obs/truth x.

Tested to reproduce values from properscoring.crps_ensemble().

The 0^th axis of ensemble is taken to enumerate the members, and must have same length as the 1d weights (if any) If x.ndim == ensemble.ndim, the 0^th axis of x is taken to enumerate multiple x. The CRPS is computed independently (and efficiently) for any/all other dimensions. Thus ensemble.shape[1:] must be matched by x.shape[1:] or x.shape, and the output gets the shape of x.

Examples:

In 1D:

>>> ens = np.array([-1.5, -1.0, 1.0, 1.5])
>>> crps_ens(0, ens)
array(0.5625)

>>> crps_ens([0], ens)
array([0.5625])

>>> crps_ens([-2, -1.5, 0, 1, 1.5, 2], ens)
array([1.3125, 0.8125, 0.5625, 0.5625, 0.8125, 1.3125])

In 2D:

>>> ens2 = np.vstack([ens, ens]).T
>>> crps_ens([1, 2], ens2)
array([0.5625, 1.3125])

>>> crps_ens([[1, 2]], ens2)
array([[0.5625, 1.3125]])

>>> crps_ens([[1, 2], [-2, -1.5]], ens2)
array([[0.5625, 1.3125],
       [1.3125, 0.8125]])

Try weighting:

>>> from scipy.stats import norm
>>> rng = np.random.default_rng(3000)
>>> ens = rng.standard_normal(10**3)
>>> grd = np.linspace(-5, 5, num=len(ens))
>>> a1 = crps_ens(0, ens)
>>> a2 = crps_ens(0, ens, weights=norm.pdf(grd))
>>> np.allclose(a1, a2, atol=1e-2)
True

`crps_ext(x, mu, var)` ⚓︎

Adapted from properscoring.crps_gaussian().

The shapes of x, mu, and var must match, but can be anything, and yields output of same shape.

Ref http://cran.nexr.com/web/packages/scoringRules/vignettes/crpsformulas.html

`inflate_ens(E, factor)` ⚓︎

Inflate the ensemble (center, inflate, re-combine).

Parameters:

Name	Type	Description	Default
`E`	`ndarray`	Ensemble which going to be inflated	required
`factor`	`float`	Inflation factor	required

Returns:

Type	Description
`ndarray`	Inflated ensemble

`mean0(E, axis=0, rescale=True)` ⚓︎

Like center, but only return the anomalies (not the mean).

Uses rescale=True by default, which is beneficial when used to center observation perturbations.

`np2builtin(v)` ⚓︎

Sometimes necessary since NEP-50

`register(self, name, value)` ⚓︎

Do self.name = value and register name in self._stat_names.

Note: self is not always a Stats object, but could be a "child" of it (e.g. a DACycleSeries). Child objects get _stat_names lazily.

`tabulate_avrgs(avrgs_list, statkeys=(), decimals=None)` ⚓︎

Tabulate avrgs (val±prec).

`unbias_var(w=None, N_eff=None, avoid_pathological=False)` ⚓︎

Compute unbias-ing factor for variance estimation.

Parameters:

Name	Type	Description	Default
`w`	`ndarray`	Importance weights. Must sum to 1. Only one of `w` and `N_eff` can be `None`. Default: `None`	`None`
`N_eff`	`float`	The "effective" size of the weighted ensemble. If not provided, it is computed from the weights. The unbiasing factor is $$ N_{eff} / (N_{eff} - 1) $$.	`None`
`avoid_pathological`	`bool`	Avoid weight collapse. Default: `False`	`False`

Returns:

Name	Type	Description
`ub`	`float`	factor used to unbiasing variance

Reference

Wikipedia

`unpack_uqs(uq_list, decimals=None)` ⚓︎

Convert list of uqs into dict of lists (of equal-length) of attributes.

The attributes are obtained by vars(uq), and may get formatted somehow (e.g. cast to strings) in the output.

If uq is None, then None is inserted in each list. Else, uq must be an instance of tools.rounding.UncertainQtty.

Parameters:

Name	Type	Description	Default
`uq_list`	`list`	List of `uq`s.	required
`decimals`	`int`	Desired number of decimals. Used for (only) the columns "val" and "prec". Default: `None`. In this case, the formatting is left to the `uq`s.	`None`

`weight_degeneracy(w, prec=1e-10)` ⚓︎

Check if the weights are degenerate.

If it is degenerate, the maximum weight should be nearly one since sum(w) = 1

Parameters:

Name	Type	Description	Default
`w`	`ndarray`	Importance weights. Must sum to 1.	required
`prec`	`float`	Tolerance of the distance between w and one. Default:1e-10	`1e-10`

Returns:

Type	Description
`bool`	If weight is degenerate True, else False

stats⚓︎

Default statistics⚓︎

DACycleSeries subscripts⚓︎

Field summaries⚓︎

Regional field summaries⚓︎

Declaring new, custom statistics⚓︎

Avrgs ⚓︎

N_eff ⚓︎

crps ⚓︎

duration ⚓︎

err ⚓︎

gscore ⚓︎

infl ⚓︎

iters ⚓︎

kurt ⚓︎

mad ⚓︎

mu ⚓︎

resmpl ⚓︎

skew ⚓︎

spread ⚓︎

svals ⚓︎

trHK ⚓︎

umisf ⚓︎

w ⚓︎

wroot ⚓︎

__getattr__(key) ⚓︎

tabulate(statkeys=(), decimals=None) ⚓︎

DACycleAvrgs ⚓︎

a ⚓︎

f ⚓︎

i ⚓︎

s ⚓︎

FieldAvrgs ⚓︎

gm ⚓︎

m ⚓︎

ma ⚓︎

ms ⚓︎

rms ⚓︎

Stats ⚓︎

N_eff ⚓︎

crps ⚓︎

duration ⚓︎

err ⚓︎

gscore ⚓︎

infl ⚓︎

iters ⚓︎

kurt ⚓︎

mad ⚓︎

mu ⚓︎

resmpl ⚓︎

rh ⚓︎

skew ⚓︎

spread ⚓︎

svals ⚓︎

trHK ⚓︎

umisf ⚓︎

w ⚓︎

wroot ⚓︎

__init__(xp, HMM, xx, yy, liveplots=False, store_i=rc.store_i, LP_kwargs=None) ⚓︎

assess(k, ko=None, fais=None, E=None, w=None, mu=None, Cov=None) ⚓︎

assess_ens(now, x, E, w) ⚓︎

assess_ext(now, x, mu, P) ⚓︎

average_in_time(kk=None, kko=None, free=False) ⚓︎

derivative_stats(now) ⚓︎

new_series(name, shape, analysis_only=False, field_mean=False, **kws) ⚓︎

replay(figlist='default', speed=np.inf, t1=0, t2=None, **kwargs) ⚓︎

summarize_marginals(now) ⚓︎

write(stat_dict, k, ko, sub) ⚓︎

align_col(col, pad='␣', missingval='', just='>') ⚓︎

center(E, axis=0, rescale=False) ⚓︎

crps_ens(x, ensemble, weights=None) ⚓︎

crps_ext(x, mu, var) ⚓︎

inflate_ens(E, factor) ⚓︎

mean0(E, axis=0, rescale=True) ⚓︎

np2builtin(v) ⚓︎

register(self, name, value) ⚓︎

tabulate_avrgs(avrgs_list, statkeys=(), decimals=None) ⚓︎

unbias_var(w=None, N_eff=None, avoid_pathological=False) ⚓︎

unpack_uqs(uq_list, decimals=None) ⚓︎

weight_degeneracy(w, prec=1e-10) ⚓︎

`DACycleSeries` subscripts⚓︎

`Avrgs` ⚓︎

`N_eff` ⚓︎

`crps` ⚓︎

`duration` ⚓︎

`err` ⚓︎

`gscore` ⚓︎

`infl` ⚓︎

`iters` ⚓︎

`kurt` ⚓︎

`mad` ⚓︎

`mu` ⚓︎

`resmpl` ⚓︎

`skew` ⚓︎

`spread` ⚓︎

`svals` ⚓︎

`trHK` ⚓︎

`umisf` ⚓︎

`w` ⚓︎

`wroot` ⚓︎

`getattr(key)` ⚓︎

`tabulate(statkeys=(), decimals=None)` ⚓︎

`DACycleAvrgs` ⚓︎

`a` ⚓︎

`f` ⚓︎

`i` ⚓︎

`s` ⚓︎

`FieldAvrgs` ⚓︎

`gm` ⚓︎

`m` ⚓︎

`ma` ⚓︎

`ms` ⚓︎

`rms` ⚓︎

`Stats` ⚓︎

`N_eff` ⚓︎

`crps` ⚓︎

`duration` ⚓︎

`err` ⚓︎

`gscore` ⚓︎

`infl` ⚓︎

`iters` ⚓︎

`kurt` ⚓︎

`mad` ⚓︎

`mu` ⚓︎

`resmpl` ⚓︎

`rh` ⚓︎

`skew` ⚓︎

`spread` ⚓︎

`svals` ⚓︎

`trHK` ⚓︎

`umisf` ⚓︎

`w` ⚓︎

`wroot` ⚓︎

`init(xp, HMM, xx, yy, liveplots=False, store_i=rc.store_i, LP_kwargs=None)` ⚓︎

`assess(k, ko=None, fais=None, E=None, w=None, mu=None, Cov=None)` ⚓︎

`assess_ens(now, x, E, w)` ⚓︎

`assess_ext(now, x, mu, P)` ⚓︎

`average_in_time(kk=None, kko=None, free=False)` ⚓︎

`derivative_stats(now)` ⚓︎

`new_series(name, shape, analysis_only=False, field_mean=False, **kws)` ⚓︎

`replay(figlist='default', speed=np.inf, t1=0, t2=None, **kwargs)` ⚓︎

`summarize_marginals(now)` ⚓︎

`write(stat_dict, k, ko, sub)` ⚓︎

`align_col(col, pad='␣', missingval='', just='>')` ⚓︎

`center(E, axis=0, rescale=False)` ⚓︎

`crps_ens(x, ensemble, weights=None)` ⚓︎

`crps_ext(x, mu, var)` ⚓︎

`inflate_ens(E, factor)` ⚓︎

`mean0(E, axis=0, rescale=True)` ⚓︎

`np2builtin(v)` ⚓︎

`register(self, name, value)` ⚓︎

`tabulate_avrgs(avrgs_list, statkeys=(), decimals=None)` ⚓︎

`unbias_var(w=None, N_eff=None, avoid_pathological=False)` ⚓︎

`unpack_uqs(uq_list, decimals=None)` ⚓︎

`weight_degeneracy(w, prec=1e-10)` ⚓︎