Skip to content

particle⚓︎

Weight- & resampling-based DA methods.

OptPF ⚓︎

'Optimal proposal' particle filter, also known as 'Implicit particle filter'.

Ref: bocquet2010a.

!!! note Regularization (Qs) is here added BEFORE Bayes' rule. If Qs==0: OptPF should be equal to the bootstrap filter PartFilt.

PFa ⚓︎

PF with weight adjustment withOUT compensating for the bias it introduces.

'alpha' sets wroot before resampling such that N_effective becomes >alpha*N.

Using alpha≈NER usually works well.

Explanation: Recall that the bootstrap particle filter has "no" bias, but significant variance (which is reflected in the weights). The EnKF is quite the opposite. Similarly, by adjusting the weights we play on the bias-variance spectrum.

NB: This does not mean that we make a PF-EnKF hybrid -- we're only playing on the weights.

Hybridization with xN did not show much promise.

PFxN ⚓︎

Particle filter with buckshot duplication during analysis.

Idea: sample xN duplicates from each of the N kernels. Let resampling reduce it to N.

Additional idea: employ w-adjustment to obtain N unique particles, without jittering.

PFxN_EnKF ⚓︎

Particle filter with EnKF-based proposal, q.

Also employs xN duplication, as in PFxN.

Recall that the proposals: Opt.: q_n(x) = c_n·N(x|x_n,Q )·N(y|Hx,R) (1) EnKF: q_n(x) = c_n·N(x|x_n,bar{B})·N(y|Hx,R) (2) with c_n = p(y|x^{k-1}_n) being the composite proposal-analysis weight, and with Q possibly from regularization (rather than actual model noise).

Here, we will use the posterior mean of (2) and cov of (1). Or maybe we should use x_a^n distributed according to a sqrt update?

PartFilt ⚓︎

Particle filter ≡ Sequential importance (re)sampling SIS (SIR).

Refs: wikle2007, van2009, chen2003

This is the bootstrap version: the proposal density is just

\[ q(x_{0:t} \mid y_{1:t}) = p(x_{0:t}) = p(x_t \mid x_{t-1}) p(x_{0:t-1}) \]

Tuning settings:

  • NER: Trigger resampling whenever N_eff <= N*NER. If resampling with some variant of 'Multinomial', no systematic bias is introduced.
  • qroot: "Inflate" (anneal) the proposal noise kernels by this root to increase diversity. The weights are updated to maintain un-biased-ness. See chen2003, section VI-M.2

particle_method ⚓︎

Declare default particle arguments.

all_but_1_is_None(*args) ⚓︎

Check if only 1 of the items in list are Truthy.

auto_bandw(N, M) ⚓︎

Optimal bandwidth (not bandwidth^2), as per Scott's rule-of-thumb.

Refs: doucet2001sequential section 12.2.2, [Wik17]_ section "Rule_of_thumb"

mask_unique_of_sorted(idx) ⚓︎

Find unique values assuming idx is sorted.

NB: returns a mask which is True at [i] iff idx[i] is not unique.

raw_C12(E, w) ⚓︎

Compute the 'raw' matrix-square-root of the ensemble' covariance.

The weights are used both for the mean and anomalies (raw sqrt).

Note: anomalies (and thus cov) are weighted, and also computed based on a weighted mean.

regularize(C12, E, idx, no_uniq_jitter) ⚓︎

Jitter (add noise).

After resampling some of the particles will be identical. Therefore, if noise.is_deterministic: some noise must be added. This is adjusted by the regularization 'reg' factor (so-named because Dirac-deltas are approximated Gaussian kernels), which controls the strength of the jitter. This causes a bias. But, as N→∞, the reg. bandwidth→0, i.e. bias→0. Ref: doucet2001sequential, section 12.2.2.

resample(w, kind='Systematic', N=None, wroot=1.0) ⚓︎

Multinomial resampling.

Refs: doucet2009, van2009, liu2001theoretical.

  • kind: 'Systematic', 'Residual' or 'Stochastic'. 'Stochastic' corresponds to rng.choice or rng.multinomial. 'Systematic' and 'Residual' are more systematic (less stochastic) varaitions of 'Stochastic' sampling. Among the three, 'Systematic' is fastest, introduces the least noise, and brings continuity benefits for localized particle filters, and is therefore generally prefered. Example: see docs/images/snippets/ex_resample.py.

  • N can be different from len(w) (e.g. in case some particles have been elimintated).

  • wroot: Adjust weights before resampling by this root to promote particle diversity and mitigate thinning. The outcomes of the resampling are then weighted to maintain un-biased-ness. Ref: liu2001theoretical, section 3.1

Note: (a) resampling methods are beneficial because they discard low-weight ("doomed") particles and reduce the variance of the weights. However, (b) even unbiased/rigorous resampling methods introduce noise; (increases the var of any empirical estimator, see [1], section 3.4). How to unify the seemingly contrary statements of (a) and (b) ? By recognizing that we're in the sequential/dynamical setting, and that future variance may be expected to be lower by focusing on the high-weight particles which we anticipate will have more informative (and less variable) future likelihoods.

reweight(w, lklhd=None, logL=None, innovs=None) ⚓︎

Do Bayes' rule (for the empirical distribution of an importance sample).

Do computations in log-space, for at least 2 reasons:

  • Normalization: will fail if sum==0 (if all innov's are large).
  • Num. precision: lklhd*w should have better precision in log space.

Output is non-log, for the purpose of assessment and resampling.

If input is 'innovs', then \(\(\text{likelihood} = \mathcal{N}(\text{innovs}|0,I)\)\).

sample_quickly_with(C12, N=None) ⚓︎

Gaussian sampling in the quickest fashion.

Method depends on the size of the colouring matrix C12.

trigger_resampling(w, NER, stat_args) ⚓︎

Return boolean: N_effective <= threshold. Also write self.stats.