On the 17th of April from 18:00 to 20:00, We will be upgrading GitLab. During this time the service will be unavailable, apologies for any inconvenience
In this notebook we will apply the least-squares method to a non-linear model: [non-linear least-squares estimation](https://mude.citg.tudelft.nl/book/observation_theory/07_NLSQ.html). In this case our model is non-linear in the unknown parameters of interest, and as you have seen in the reading, we formulate the linearized model in terms of the difference between the initial guess and computed value of each parameter:
We will apply best linear unbiased estimation to estimate the parameters $\Delta \mathrm{x}_{[0]}$. We will then obtain a new 'estimate' for the unknown parameters $\mathrm{x}_{[1]}=\mathrm{x}_{[0]}+\Delta \mathrm{x}_{[0]}$ and repeat the procedure. This is the Gauss-Newton iteration.
## Background
In this exercise we will consider an aquifer: a subsurface reservoir where water is stored in the pores of a soil or rock formation. If this is unfamiliar to you, think of it as an underground sponge with water in it that responds slowly to change. A well has been installed with monitoring equipment such that a time series of water level in the aquifer has been obtained, and our objective is to create a mathematical model that captures the fluctuations of the water level with time.
Water enters or exits the aquifer due to certain events, such as precipitation, evaporation or groundwater pumping. We assume that the water level in the aquifer is a linear combination of the events, allowing us to write the observation equation as:
$$ L(t) = \sum_{i=1}^{n} h_i(t) + d + \epsilon(t) $$
where $L(t)$ is the observed water level, $h_{i}(t)$ is an event that adds or removes water from the aquifer, $d$ is the base water level of the aquifer and $\epsilon(t)$ is the noise (random error) term. In this exercise we will consider the response due to a single event, $h_k(t)$. Groundwater flows through pores in the aquifer much slower than the application of $h_k(t)$ at time $t_0$. Therefore, the aquifer responds to event $k$ according to a step-response function, $s(t-t_0)$, which has the general form illustrated in the figure below. If the event has a constant value and continues indefinitely after $t_0$ the aquifer will eventually reach a steady-state value $r$ (a scaling factor of $a=1$ is used).
<div>
<center>
<imgsrc="./figures/response.png",width=320/>.
</center>
</div>
For our exercise we will consider a single event, $k$, with a constant input of water, $p$ [m], from time $t_0$ until $t_{\text{end}}$, and zero input before and after. In this case $p$ is assumed to be a deterministic (known) parameter and constant in time (imagine that it represents water inflow to the aquifer due to a rain event). The steady-state aquifer response scales with $p$, resulting in the following step-response function for the aquifer:
$$
h_k(t) = p \cdot \underset{s(t-t_0)}{\underbrace{r\left(1-\exp\left(-\frac{t-t_0}{a}\right)\right)}}
$$
where $a$ [days] is a scaling parameter representing the memory of the system (determines how fast the level is rising) and $r$ [m/m] is the response of the aquifier to an input $p$. Once the event stops, the response is described by:
$$
h_k(t) = p\cdot r\left(\exp\left(-\frac{t-t_e}{a}\right)-\exp\left(-\frac{t-t_0}{a}\right)\right), \;\;\text{for}\;\; t > t_{\text{end}}
$$
%% Cell type:markdown id: tags:
## Functional model
For this example, we consider a single rainfall event. At $t_0=4$ days, it starts raining, and at $t_{\text{end}}=7$ days it stops raining (and we assume the amount of rainfall to be constant during these days, for the sake of the example...). The observation equations become:
$$
\begin{align}
\mathbb{E}[L(t_i)] &= d && \text{for} \; t_1\leq t_i < 4\\
From these equations, it follows that there are three unknown parameters that need to be estimated: $d$, $r$ and $a$. These equations can be combined into a single equation by making use of the <ahref="https://en.wikipedia.org/wiki/Heaviside_step_function"> Heaviside function</a> $H(\Delta t)$:
$$\mathbb{E}\left[L(t_i)\right] = d + H(t_i-t_{0}) \cdot p\cdot r \left(1-\exp\left(-\frac{t_i-t_0}{a}\right)\right) - H(t_i-t_{\text{end}}) \cdot p\cdot r \left(1-\exp\left(-\frac{t_i-t_{\text{end}}}{a}\right)\right)$$
with $H(\Delta t) = 1$ if $\Delta t \geq 0$ and $H(\Delta t) = 0$ otherwise. <b>Check yourself that this gives indeed the same observation equations as above.</b>
The forward model can be defined as follows, using <ahref="https://numpy.org/doc/stable/reference/generated/numpy.heaviside.html"> NumPy's heaviside</a> function:
%% Cell type:code id: tags:
``` python
importnumpyasnp
importmatplotlib.pyplotasplt
fromnumpy.linalgimportinv
fromscipy.stats.distributionsimportchi2
plt.rcParams.update({'font.size':14})
```
%% Cell type:code id: tags:
``` python
defrain_event(times,rainy_days,rain,d,a,r):
'''Forward model: response due to rain event.'''
h=(d
+(np.heaviside(times-rainy_days[0],1)
*rain*r*(1-np.exp(-(times-rainy_days[0])/a))
-np.heaviside(times-rainy_days[1],1)
*rain*r*(1-np.exp(-(times-rainy_days[1])/a))
)
)
returnh
```
%% Cell type:markdown id: tags:
### Forward Model
If we generate a vector of time steps, we can plot the (exact) response of a system with parameters using the ```rain_event``` function defined above.
Familiarize yourself with the model by assigning values to the unknown parameters below. Continue changing values and re-running the cell until you are comfortable with the effect that each parameter has on the model. Note that the scale of the water level response to 0.05 m water input (e.g., precipiation) should be meters (i.e., >0 m and <10 m).
<em>Hint: if you are having trouble visualizing the model, adjust the y-axis limits using <code>plt.ylim([ , ]</code>)
</p>
</div>
%% Cell type:code id: tags:
``` python
# d = YOUR_CODE_HERE
# a = YOUR_CODE_HERE
# r = YOUR_CODE_HERE
# SOLUTION:
d=0
a=1
r=100
test_n_days=25
test_times=np.arange(1,test_n_days+1,0.1)
test_h_t=rain_event(test_times,(4,7),0.05,d,a,r)
plt.figure(figsize=(10,6))
plt.plot(test_times,test_h_t,'r-',linewidth=4.0)
plt.xlabel('Time [days]')
plt.ylabel('Water level [m]')
plt.xlim([0,test_n_days])
plt.ylim([0,5]);
```
%% Output
%% Cell type:markdown id: tags:
## Reading in the observations
We collected observations of the water level in the aquifer for 25 consecutive days and are stored ```data``` folder. Observations start at $t=1$.
Derive the partial derivatives that you will need for the Jacobian (on paper), then write them in the code below.
<em>Hints: when implementing the partial derivatives consider the dimensions of the Jacobian matrix carefully (the first term should <b>not</b> be 1; think "vector"). Also, the heaviside function takes care of the piece-wise function, so you do not need to implement this "manually."</em>
In the code you will first need to choose the initial values <code>d_i, a_i, r_i</code>. Hint: for the choice of <code>d_i</code> = $d_{[0]}$ have a look at the plot with the observations (remember it is the base level of the aquifer before the rain starts). Using the intuition you gained from the plot from the forward model (Task 0), find appropriate initial values for $a$ and $r$ by comparing the graph you get with the observations.
Next we will set the stochastic model, and initialize some other variables needed for the <code>while</code>-loop, where the actual Gauss-Newton iteration takes place. Try to see if you can follow the steps above, and fill in the missing code, where you have to compute
- the observed-minus-computed observations <code>dy</code> $=\Delta y_{[i]}$,
- the estimate <code>dxhat_i</code>$=\Delta\hat{\mathrm{x}}_{[i]}$,
Note in particular the 2 stop criteria used for the while loop. You should reach a solution within 50 iterations, otherwise you should reconsider the initial values.
</p>
</div>
%% Cell type:code id: tags:
``` python
# d_i = YOUR_CODE_HERE
# a_i = YOUR_CODE_HERE
# r_i = YOUR_CODE_HERE
# SOLUTION:
d_i=y[0]
a_i=1
r_i=3
rain=0.05
rainy_days=(4,7)
sigma=0.01
var_y=sigma**2
inv_Sigma_Y=1/var_y*np.eye(len(y))
max_iter=50
x_norm=10000# initialize stop criteria (norm of x)
# suffix "_i" indicates initial values; xhat_i stores every iteration
Only a few iterations are needed! It is clear that the parameter values become stable quickly. Also important: the converged values are close to those we estimated in the previous Task using the figures and the equation.
</p>
</div>
%% Cell type:markdown id: tags:
## Overall model test
Looking at the final results, it seems the model provides a good fit to the data. Let's test it with a probability of false alarm of 0.05.
Play with the value of <code>alpha</code> and see how it changes the critical value <code>k</code>. Why does it become smaller/larger? What is the impact on the probability of rejecting the null hypothesis?
Larger $\alpha$ results in smaller critical value - the critical region is right-sided, so will be larger. Hence there will be a larger (smaller) probability of rejection with a larger (smaller) $\alpha$.
Note in this case we would have to choose $\alpha = 0.88$ to get a rejection, which is incredibly high. This means the model is quite good!
