The observation model

To define the relationship between the Lorenz-63 simulation model and observations of this system, we need to define an observation model. When the observation model can be described using a standard probability distribution, we only need to create a Univariate subclass that returns a SciPy distribution for the given state vectors \(\mathbf{x_t}\).

For simplicity, we assume that \(x(t)\), \(y(t)\), and \(z(t)\) can be directly observed (define the observed values as \(X_t\), \(Y_t\), and \(Z_t\), respectively) and that the observation error is distributed normally with zero mean and standard deviation \(\sigma = 1.5\):

\[\begin{split}X_t \sim \mathcal{N}(\mu = x(t), \sigma = 1.5) \\ Y_t \sim \mathcal{N}(\mu = y(t), \sigma = 1.5) \\ Z_t \sim \mathcal{N}(\mu = z(t), \sigma = 1.5)\end{split}\]

The implementation of these observation models is straightforward, and comprises two steps. First, extract the expected value for each particle from the state vectors, then construct a corresponding normal distribution for each particle:

class ObsLorenz63(Univariate):

    def distribution(self, ctx, snapshot):
        expect = snapshot.state_vec[self.unit]
        return scipy.stats.norm(loc=expect, scale=1.5)

Note

pypfilt support multiple observation models. Each observation model is associated with a unique identifier (an “observation unit”) that is used to identify the observations related to this model. Here, we use the observation unit (self.unit) to identify the field in the state vector that is being observed (see the highlighted line, above). This allows us to use three instances of the ObsLorenz63 class to observe \(x(t)\), \(y(t)\), and \(z(t)\).

We could define an observation model specifically for \(x(t)\) (which is named 'x' in the state vector) by replacing the highlighted line above with:

expect = snapshot.state_vec['x']