TUTORIAL Class Observations¶
This tutorial introduces the Observations class and shows how to create synthetic truth and observations
Key Observations arguments
| Argument | Meaning |
|---|---|
t_start, t_stop |
Assimilation window |
t_max |
How far the truth is integrated beyond t_stop |
Nt_obs |
Observation frequency (every Nt_obs model time-steps) |
noise_level |
Std of additive Gaussian noise (fraction of signal std) |
manual_bias |
Added bias to the state (model error) |
observe_dims |
Which state dimensions are observed |
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
t = np.linspace(0, 20, 200)
true_data = np.sin(t)
raw_data = true_data + 0.1 * np.random.randn(len(t)) # Add a noise to the true data
dt_obs = 10
observed_data = raw_data[::dt_obs]
plt.figure(figsize=(6, 3))
plt.plot(t, true_data, c='k', alpha=.2, lw=4, label='True Data')
plt.plot(t, raw_data, c='C0', lw=1, label='Raw Data')
plt.scatter(t[::dt_obs], observed_data, color='red', label='Observed Data', edgecolors='k')
plt.xlabel('Time'), plt.ylabel('Value'), plt.legend(loc='upper left', bbox_to_anchor=(1, 1));
1. Create timeseries from the low-order model¶
Note that this data is clean and biased. Both the noise and the bias will be added in the next steps.
from romda.observations import Observations
reference_data = Observations(y_true=true_data,
y_raw=raw_data,
t_true=t,
Nt_obs=10
)
Observations.plot_truth(reference_data, window=20)
Instead of providing the dataset, we can create the observations from a model. For instance, using the Lorenz 63 model.
from romda.models.physical import Lorenz63
from romda.utils import set_working_directories
# Compat shim: truth pickles cached before the dynamodels split store old flat
# module paths (model, history, integrator, models_physical.*); alias them so
# the cached Truth_Lorenz63_beta3.0 keeps loading.
import sys
import dynamodels.physical
for _old, _new in {'model': dynamodels.model,
'history': dynamodels.history,
'integrator': dynamodels.integrator,
'models_physical': dynamodels.physical,
'models_physical.lorenz63': dynamodels.physical.lorenz63}.items():
sys.modules.setdefault(_old, _new)
results_dir = set_working_directories('Lorenz63')[1]
truth = Observations(model=Lorenz63,
beta=3.,
t_start=20.0,
t_stop=40.0,
Nt_obs=10,
results_folder=results_dir)
Observations.plot_truth(truth, Nq=3, f_max=20)
Loaded true data model: Truth_Lorenz63_beta3.0 <dynamodels.physical.lorenz63.Lorenz63 object at 0x7f1810f63390>
3. Add noise to the truth to create the observations¶
The noise type can be selected between Gaussian or coloured noise (from colours white, pink, brown, blue, and violet). Further, one can define the noise to be either [additive] or multiplicative.
- The bias is added in the initialization of the Observations instance via the _apply_noise function. Note: 'gaussian, additive' is the pre-defined setting.
- The function utils.colour_noise allows us to add different types of noises. We can visualize the different options.
truth_noisy = Observations(model=Lorenz63,
t_start=20.0,
t_stop=40.0,
Nt_obs=10,
add_noise=True,
noise_type='gauss, add',
noise_level=0.05)
Observations.plot_truth(truth_noisy,Nq=1)
...Adding noise: gauss, add with level 0.05.
import matplotlib.pyplot as plt
from romda.utils import colour_noise
import numpy as np
rng = np.random.default_rng(6)
noise_level = 0.02
noise_type = 'pink' + 'additive'
N = 100000
noise_white = np.fft.rfft(rng.standard_normal(N))
freq = np.fft.rfftfreq(N)
NOISES, PSDS = [], []
COLOURS = ['white', 'pink', 'brown', 'blue', 'darkviolet']
for noise_c in COLOURS:
S = colour_noise(N, noise_colour=noise_c)
S = noise_white * S
noise = np.fft.irfft(S) # transform back into time domain
# Store
PSDS.append(abs(S))
NOISES.append(noise)
fig = plt.figure(figsize=(10, 5))
figs = fig.subfigures(1, 2)
ax = figs[1].subplots(1, 1)
alpha = 1
for psd, c in zip(PSDS, COLOURS):
c = c if c != 'white' else 'gray'
ax.loglog(freq, psd, color=c, alpha=alpha)
alpha -= 0.15
ax.legend(COLOURS, ncol=2)
ax.set(xlabel='Frequency', ylim=[1e-3, None], title='PSD')
axs = figs[0].subplots(len(COLOURS), 1, sharex='col', sharey='col')
alpha = 1
N_plot = 1000
for ax, noise, c in zip(axs, NOISES, COLOURS):
c = c if c != 'white' else 'gray'
ax.plot(np.arange(N), noise, color=c, alpha=alpha)
alpha -= 0.15
axs[0].set(title='Time domain', xlim=[0, 2*N_plot])
axs[-1].set(xlabel='$t$');
We can observe that * White noise has an almost-flat PSD in the frequency domain, i.e., all the frequencies are equally present * Brownian and pink noises strength is lowest at higher frequencies, which is visualized in the time domain with seamingly non-zero mean. * Blue and violet noise strenght increases with the frequency.
2. (optional) Add bias¶
If we do not add the bias, then the Model is unbiased. Alternatively,
\[
\mathbf{y}^\mathrm{true} = \mathbf{y}^\mathrm{Model} + \mathbf{b}^\mathrm{true}
\]
then the prediction from the Model becomes biased. The bias is added in the initialization of the Observations instance via the _set_bias function.
t = np.linspace(0, 20, 200)
true_data = np.sin(t)
bias = 0.05 * t # for intance, a time drift
biased_data = true_data + bias # this is noise-free but biased data
raw_data = biased_data + 0.1 * np.random.randn(len(t)) # Add a noise to the true data
dt_obs = 10
observed_data = raw_data[::dt_obs]
plt.figure(figsize=(6, 3))
plt.plot(t, true_data, c='k', alpha=.2, lw=4, label='True Data')
plt.plot(t, biased_data, c='b', alpha=.3, lw=4, label='Biased Data')
plt.plot(t, bias, c='C4', ls='--', lw=1, label='Bias')
plt.plot(t, raw_data, c='C0', lw=1, label='Raw Data')
plt.scatter(t[::dt_obs], observed_data, color='red', label='Observed Data', edgecolors='k')
plt.xlabel('Time'), plt.ylabel('Value'), plt.legend(loc='upper left', bbox_to_anchor=(1, 1));
help(Observations._set_bias)
Help on function _set_bias in module romda.observations:
_set_bias(self)
Compute and add the manual bias $\mathbf{b}^t$ on top of `y_true`.
Controlled by `manual_bias`:
- `None` (default): no bias, $\mathbf{b}^t = \mathbf{0}$.
- ``'time'``: $\mathbf{b}^t = 0.4\,\mathbf{y}^t \sin\!\big((2\pi t)^2\big)$.
- ``'periodic'``: $\mathbf{b}^t = 0.2\,\mathbf{y}^t_{\max} \cos\!\big(2\mathbf{y}^t / \mathbf{y}^t_{\max}\big)$,
where $\mathbf{y}^t_{\max}$ is the time-wise maximum of `y_true`.
- ``'linear'``: $\mathbf{b}^t = 0.1\,\mathbf{y}^t_{\max} + 0.3\,\mathbf{y}^t$.
- ``'cosine'``: $\mathbf{b}^t = \cos(\mathbf{y}^t)$.
- a callable ``f(y_true, t_true) -> (b_true, name)`` for a user-defined bias.
Updates `y_true` in place (adds the bias) and sets `b_true` and `name_bias`.
If `y_raw` was not provided at construction, it is also initialised here to
the (now biased) `y_true` — noise is added on top of it in `_apply_noise`.
truth_noisy_biased = Observations(model=Lorenz63,
t_start=20.0,
t_stop=40.0,
Nt_obs=10,
add_noise=True,
noise_type='gauss, add',
noise_level=0.05,
manual_bias='time')
Observations.plot_truth(truth_noisy_biased,Nq=1)
...Applying manual bias: time
...Adding noise: gauss, add with level 0.05.
# Alternative: user-defined bias function
def my_bias_function(y, t):
"""Example of user-defined bias function."""
b = 0.1 * t[:, np.newaxis, np.newaxis] + .2*abs(y)
return b, 'mybias'
truth_noisy_biased = Observations(model=Lorenz63,
t_start=20.0,
t_stop=40.0,
Nt_obs=10,
add_noise=True,
noise_type='gauss, add',
noise_level=0.05,
manual_bias=my_bias_function)
Observations.plot_truth(truth_noisy_biased, Nq=3)
...Applying user-defined manual bias
...Adding noise: gauss, add with level 0.05.
