"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 \n",
"\n",
"- the observed-minus-computed observations <code>Delta_y_i</code> $=\\Delta y_{[i]}$,\n",
"- the observed-minus-computed observations <code>Delta_y_i_i</code> $=\\Delta y_{[i]}$,\n",
"- the estimate <code>Delta_x_hat_i</code>$=\\Delta\\hat{\\mathrm{x}}_{[i]}$,\n",
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/2024/book/observation_theory/07_NLSQ.html). In this case our model is non-linear in the unknown parameters of interest:
In other words, we can no longer use a (linear) A matrix to evaluate the functional model. However, as you have seen in the reading, we still formulate a _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 functional 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:markdown id: tags:
## Part 0: Set Up and Understand Functions
As described at the top of the notebook, our functional model is no longer a linear system of equations. This means we can no longer use a simple A matrix to evaluate it. Instead, we will use a Python function. Recall that the general form of the functional model is $\mathrm{Y}=q(\mathrm{x}) + \epsilon$, where $\mathrm{x}$ represents the vector of model parameters. Note also that to implement this model in practice often additional (deterministic!) parameters are required. This can be written in a Python function with the following form of input arguments:
```
y_comp = compute_y(x, <auxiliary_arguments>)
```
Where the name `compute_y()` is used to indicate that this is not the functional model, as the Python function does not include the random errors. The <auxiliary_arguments>` will be different in type and/or number on a case-by-case basis. Your code will generally be more compact and adaptable to other cases if the parameters are specified in a list, array or tuple.
Read the code below to understand the functional model. Be sure to understand which objects are the parameters of interest, and what the auxiliary arguments are (type and number).
</p>
</div>
%% Cell type:code id: tags:
``` python
def compute_y(x, times, rain):
'''Functional model, q: response due to rain event.
Inputs:
x: tuple, list or array of parameters (d, a, r)
times: array of times
rain: tuple, list or array describing rain event
(rain (m), start day, stop day)
Returns: ndarray of groundwater levels due to rain event
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
test_n_days = 25
test_times = np.arange(1, test_n_days+1, 0.1)
test_h_t = compute_y((d, a, r), test_times, (0.05, 4, 7))
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_init, a_init, r_init</code>. Hint: for the choice of <code>d_init</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 functional 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>Delta_y_i</code> $=\Delta y_{[i]}$,
- the observed-minus-computed observations <code>Delta_y_i_i</code> $=\Delta y_{[i]}$,
- the estimate <code>Delta_x_hat_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_init = YOUR_CODE_HERE
a_init = YOUR_CODE_HERE
r_init = YOUR_CODE_HERE
rain_event = (0.05, 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)
# x_hat_i: array to store the estimated parameters at each iteration
param_init = np.array([d_init, a_init, r_init])
x_hat_i = np.zeros((3, max_iter))
x_hat_i[:] = np.nan
x_hat_i[:, 0] = param_init
iteration = 0
while x_norm >= 1e-12 and iteration < max_iter - 1:
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?
