ACSE-7 Inversion & Optimisation – Part C (15 pts)
This question covers lecture 8 on Data Assimilation methods only.
Please read the general instructions in README.md first!
Probability
The exponential distribution is a one parameter distribution with a probability density
$$ p(x):=\begin{cases} \lambda\exp\left(-\lambda x\right), & x>0,\\ 0, & \mbox{otherwise} \end{cases}$$
for a known value of $\lambda$, which is the reciprocal of the mean of the distribution. Random number generators for this distribution are implementated in regular Python as random.expovariate and in numpy as numpy.random.exponential.
In [7]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
x = np.linspace(0, 20)
lmd = 1.0
plt.plot(x, lmd*np.exp(-lmd*x))
Out[7]:
[

Question C1
Calculate the formula for the cumulative distribution function $P(X\leq z) = \int_{-\infty}^z p(x)dx$ for the exponential PDF given above. Write python code to plot this CDF function against a histogram of at least 1000 samples from one of the implementations above using the scale factor $\lambda=1$ over the range $x\in[0,20]$.
In [ ]:
#write your code here.
# You can use the matplotlib hist function for the histogram, or any other method
#plot your function
plt.plot(x, …)
# take your samples
data = …
plt.hist(data, density=True)
Data Assimilation
A set of $N_1$ particles move in 2D space. The vector of their positions $\mathbf{x}$, where $$\mathbf{x}=\{\mbox{x-coordinate of particle 1}, \mbox{y-coordinate of particle 1}, \ldots, \mbox{y-coordinate of particle $N_1$}\} $$ can be modelled forward in time using an equation $$\mathbf{x}_{i+1}=\mathbf{M}\mathbf{x}_i+\mathbf{f}$$ for a known fixed linear update model $\mathbf{M}$ and forcing vector $\mathbf{f}$. You have been asked to forecast the future positions of the particles, and are provided with a photograph of their current state. It is believed that the vector, $\mathbf{y}$, of the intensities of the $N_2$ pixels of the photograph, satisfies a relation approximately proportional to the sum of the sum of the inverse squares of the distances of their pixel centres from the particle sources.
That is to say that for a pixel at position $\mathbf{q}$ its intensity is $$I= k\sum_{\mathbf{r}\in\mathbf{x}}\frac{1}{a+\|\mathbf{q}-\mathbf{r}\|^2}.$$ for known constants $k$ and $a$. This implicitly defines an observation operator, $\mathbf{h}(\mathbf{x})$.

You may assusme that $N_1\ll N_2$. You are also provided with estimates for the error covariance matrices for the background state, $\mathbf{B}$, the model, $\mathbf{Q}$, and the observations, $\mathbf{R}$, with all errors assumed Gaussian.
Question C2:
Write Python code to implement the observation operator described above in a format appropriate for data assimilation methods. You can assume that the input state vector $\mathbf{x}$ is a length $2N_1$ vector, giving the 2D spatial coordinates of the first point, then the 2nd and so on, and that that your function will also be provided with a $2\times N_2$ array containing the $\mathbf{q}$ positions of the pixels, as well a $k$ value.
In [1]:
def h(x, q, k=1.0, a):
“””Calculate expected observation for a given state vector.
Parameters
———-
x: ndarray
1d array of length 2N_1 containing point positions
q: ndarray
2d array of length [2, N_2] containing pixel centre positions
k: float
scaling coefficient
a: float
weighting coefficient
Result
——
1d array of length N_2 containing expected intensities
“””
y = …
return y
Question C3
Consider the following list of potential data assimilation strategies:
• A single OI assimilation event
• A single 3D-Var assimilation event
• A single strong 4D-Var assimilation event
• iterations of the ordinary Kalman filter
• iterations of the extended Kalman filter
Place the five possible strategies in order from most appropriate to least appropriate in order to efficiently generate an accurate forcast at the future positions of the particles, given the information you have. For each strategy include a short (i.e at least a sentence, but no more than a paragraph) explanation of why you have placed it where you have.
Write your answer to question C3 in this cell
Most appropriate
1. ..
2. ..
3. ..
4. ..
5. ..
Least appropriate
Question C4
How would your answer change if you were instead to be given a set of photographs at known times? Include a short explanation of any changes to your previous order that you make.
Write your answer to question C4 in this cell
Most appropriate
1. ..
2. ..
3. ..
4. ..
5. ..
Least appropriate
Question C5
Further analysis reveals that the particle intensities above only register on the photograph above a given cut off value,
$$ I_\mbox{registered} = \begin{cases} I-c,&I\geq c\\ 0,&I