dapper.mods.Lorenz96.bocquet2010

 1"""From Fig. 1 of `bib.bocquet2010beyond`."""
 2
 3import numpy as np
 4
 5import dapper.mods as modelling
 6from dapper.mods.Lorenz96 import step
 7
 8tseq = modelling.Chronology(0.05, dko=1, T=4**3, BurnIn=20)
 9
10Nx = 10
11Dyn = {
12    "M": Nx,
13    "model": step,
14    "noise": 0,
15}
16
17X0 = modelling.GaussRV(M=Nx, C=0.001)
18
19jj = np.arange(0, Nx, 2)
20Obs = modelling.partial_Id_Obs(Nx, jj)
21Obs["noise"] = 1.5
22
23HMM = modelling.HiddenMarkovModel(Dyn, Obs, tseq, X0)
24
25####################
26# Suggested tuning
27####################
28# Why are these benchmarks superior to those in the article?
29# We use, in the EnKF,
30# - inflation instead of additive noise ?
31# - Sqrt      instead of perturbed obs
32# - random orthogonal rotations.
33# The particle filters are also probably better tuned:
34# - jitter covariance proportional to ensemble (weighted) cov
35# - no jitter on unique particles after resampling
36#
37# For a better "picture" of the relative performances,
38# see benchmarks in presentation from SIAM_SEAS.
39# Note: They are slightly unrealiable (short runs).
40
41#                                                   rmse.a:
42# xps += EnKF_N(N=8,rot=True,xN=1.3)                # 0.31
43
44# xps += PartFilt(N=50 ,NER=0.3 ,reg=1.7)           # 1.0
45# xps += PartFilt(N=100,NER=0.2 ,reg=1.3)           # 0.36
46# xps += PartFilt(N=800,NER=0.2 ,reg=0.8)           # 0.25
47
48# xps += OptPF(   N=50 ,NER=0.25,reg=1.4,Qs=0.4)    # 0.61
49# xps += OptPF(   N=100,NER=0.2 ,reg=1.0,Qs=0.3)    # 0.37
50# xps += OptPF(   N=800,NER=0.2 ,reg=0.6,Qs=0.1)    # 0.25
51
52# xps += PFa(     N=50 ,alpha=0.4,NER=0.5,reg=1.0)  # 0.45
53# xps += PFa(     N=100,alpha=0.3,NER=0.4,reg=1.0)  # 0.38
54
55# xps += PFxN     (N=30, NER=0.4, Qs=1.0,xN=1000)   # 0.48
56# xps += PFxN     (N=50, NER=0.3, Qs=1.1,xN=100 )   # 0.43
57# xps += PFxN     (N=100,NER=0.2, Qs=1.0,xN=100 )   # 0.32
58# xps += PFxN     (N=400,NER=0.2, Qs=0.8,xN=100 )   # 0.27
59# xps += PFxN     (N=800,NER=0.2, Qs=0.6,xN=100 )   # 0.25
60
61# xps += PFxN_EnKF(N=25 ,NER=0.4 ,Qs=1.5,xN=100)    # 0.49
62# xps += PFxN_EnKF(N=50 ,NER=0.25,Qs=1.5,xN=100)    # 0.36
63# xps += PFxN_EnKF(N=100,NER=0.20,Qs=1.0,xN=100)    # 0.32
64# xps += PFxN_EnKF(N=300,NER=0.10,Qs=1.0,xN=100)    # 0.28