In this notebook you will fit a model to a time series of height observations of a point on a glacier, to assess whether it is melting.
Learning objectives:
- apply least-squares (LS) and best linear unbiased (BLU) estimation
- evaluate the precision of the estimated parameters
- discuss the differences between least-squares and best linear unbiased estimation
You have 12 monthly measurements of the height of a point on a glacier. The measurements are obtained from a satellite laser altimeter. The observations are given in the code below.
We will fit a model assuming a linear trend (constant velocity) to account for the melting, plus an annual signal due to the seasonal impact on ice. Formulating a generalized form of this model is part of your task below. Note, however, that the periodic signal can be implemented in many ways. For this assignment you can assume that the periodic term is implemented as follows:
$$ \cos(2 \pi t / 12 + \phi) $$where $\phi$ is the phase angle, which we will assume is 0 for this workshop and will leave out of the model formulation in order to keep our model linear (i.e., leave the $\phi$ out of the A-matrix).
The precision (1 $\sigma$) of the first six observations is 0.7 meter, the last six observations have a precision of 2 meter due to a change in the settings of the instrument. All observations are mutually independent.
$\textbf{Solution (note)}$
One goal of the assignment was to formulate the model yourself (via the A-matrix). The complete form can also be described using the following analytic expression:
$$ y = x_1 + x_2 t + x_3cos(2 \pi t / 12 + \phi) $$Although it may initially seem that the periodic term leaves out the amplitude, you should recognize that the amplitude is in fact implicitly included in the model due to the way $x_3$ is included in the A matrix (i.e., as a linear parameter with the cosine term).
It is important to confirm that your own A-matrix (found during the session) includes all three terms as illustrated above.
Note also that this formulation was part of the reading in the textbook in the exercise page of the Least Squares chapter (3.2). Here is a link.
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
# this will print all float arrays with 3 decimal places
float_formatter = "{:.3f}".format
np.set_printoptions(formatter={'float_kind':float_formatter})
Part 1: Observation model¶
First we will construct the observation model, based on the following information that is provided: times of observations, observed heights and number of observations.
$\textbf{Task 1.1}:$
Complete the definition of the observation time array times and define the other variables such that the print statements correctly describe the design matrix A.
Use times=0 as the first point.
y = [59.82, 57.20, 59.09, 59.49, 59.68, 59.34, 60.95, 59.34, 55.38, 54.33, 48.71, 48.47]
# times = YOUR_CODE_HERE
# number_of_observations = YOUR_CODE_HERE
# number_of_parameters = YOUR_CODE_HERE
# print(f'Dimensions of the design matrix A:')
# print(f' {YOUR_CODE_HERE} rows')
# print(f' {YOUR_CODE_HERE} columns')
# SOLUTION:
times = np.arange(12)
number_of_observations = len(times)
number_of_parameters = 3
# SOLUTION:
print(f'Dimensions of the design matrix A:')
print(f' {number_of_observations:3d} rows')
print(f' {number_of_parameters:3d} columns')
Dimensions of the design matrix A:
12 rows
3 columns
Next, we need to assemble the design matrix in Python. Rather than enter each value manually, we have a few tricks to make this easier. Here is an example Python cell that illustrates how one-dimensional arrays (i.e., columns or rows) can be assembled into a matrix in two ways:
- As diagonal elements, and
- As parallel elements (columns, in this case)
Note also the use of np.ones() to quickly make a row or column of any size with the same values in each element.
column_1 = np.array([33333, 0, 0])
column_2 = 99999*np.ones(3)
example_1 = np.diagflat([column_1, column_2])
example_2 = np.column_stack((column_1, column_2))
print(example_1, '\n\n', example_2)
[[33333.000 0.000 0.000 0.000 0.000 0.000] [0.000 0.000 0.000 0.000 0.000 0.000] [0.000 0.000 0.000 0.000 0.000 0.000] [0.000 0.000 0.000 99999.000 0.000 0.000] [0.000 0.000 0.000 0.000 99999.000 0.000] [0.000 0.000 0.000 0.000 0.000 99999.000]] [[33333.000 99999.000] [0.000 99999.000] [0.000 99999.000]]
See documentation on np.diagflat and np.column_stack if you would like to understand more about these Numpy methods.
$\textbf{Task 1.2}:$
Complete the $\mathrm{A}$-matrix (design matrix) and covariance matrix $\Sigma_Y$ (stochastic model) as a linear trend with an annual signal. The assert statement is used to confirm that the dimensions of your design matrix are correct (an error will occur if not).
# A = YOUR_CODE_HERE
# Sigma_Y = YOUR_CODE_HERE
# SOLUTION:
A = np.column_stack((np.ones(number_of_observations),
times,
np.cos(2*np.pi*times/12)))
diag_1 = 0.7**2 * np.ones((1,6))
diag_2 = 2**2 * np.ones((1,6))
Sigma_Y = np.diagflat([diag_1,
diag_2])
assert A.shape == (number_of_observations, number_of_parameters)
2. Least-squares and Best linear unbiased (BLU) estimation¶
The BLU estimator is a linear and unbiased estimator, which provides the best precision, where:
- linear estimation: $\hat{X}$ is a linear function of the observables $Y$,
- unbiased estimation: $\mathbb{E}(\hat{X})=\mathrm{x}$,
- best precision: sum of variances, $\sum_i \sigma_{\hat{x}_i}^2$, is minimized.
The solution $\hat{X}$ is obtained as: $$ \hat{X} = \left(\mathrm{A}^T\Sigma_Y^{-1} \mathrm{A} \right)^{-1} \mathrm{A}^T\Sigma_Y^{-1}Y $$ Note that here we are looking at the estimator, which is random, since it is expressed as a function of the observable vector $Y$. Once we have a realization $y$ of $Y$, the estimate $\hat{\mathrm{x}}$ can be computed.
It can be shown that the covariance matrix of $\hat{X}$ is given by: $$ \Sigma_{\hat{X}} = \left(\mathrm{A}^T\Sigma_Y^{-1} \mathrm{A} \right)^{-1} $$
$\textbf{Task 2.1}$
Apply least-squares, and best linear unbiased estimation to estimate $\mathrm{x}$. The code cell below outlines how you should compute the inverse of $\Sigma_Y$. Compute the least squares estimates then the best linear unbiased estimate.
Remember that the @ operator (matmul operator) can be used to easily carry out matrix multiplication for Numpy arrays. The transpose method of an array, my_array.T, is also useful.
inv_Sigma_Y = np.linalg.inv(Sigma_Y)
# xhat_LS = YOUR_CODE_HERE
# xhat_BLU = YOUR_CODE_HERE
# SOLUTION:
xhat_LS = np.linalg.inv(A.T @ A) @ A.T @ y
xhat_BLU = np.linalg.inv(A.T @ inv_Sigma_Y @ A) @ A.T @ inv_Sigma_Y @ y
print('LS estimates in [m], [m/month], [m], resp.:\t', xhat_LS)
print('BLU estimates in [m], [m/month], [m], resp.:\t', xhat_BLU)
LS estimates in [m], [m/month], [m], resp.: [62.569 -1.046 -3.717] BLU estimates in [m], [m/month], [m], resp.: [62.318 -1.043 -3.235]
Task 2.2: The covariance matrix of least-squares can be obtained as well by applying the covariance propagation law. Calculate the covariance matrix and vector with standard deviations of both the least-squares and BLU estimates. What is the precision of the estimated parameters? The diagonal of a matrix can be extracted with np.diagonal.
Hint: define an intermediate variable first to collect some of the repetitive matrix terms.
# LT = YOUR_CODE_HERE
# Sigma_xhat_LS = LT YOUR_CODE_HERE
# std_xhat_LS = YOUR_CODE_HERE
# SOLUTION:
LT = np.linalg.inv(A.T @ A) @ A.T
Sigma_xhat_LS = LT @ Sigma_Y @ LT.T
std_xhat_LS = np.sqrt(np.diag(Sigma_xhat_LS))
# Sigma_xhat_BLU = YOUR_CODE_HERE
# std_xhat_BLU = YOUR_CODE_HERE
# SOLUTION:
Sigma_xhat_BLU = np.linalg.inv(A.T @ inv_Sigma_Y @ A)
std_xhat_BLU = np.sqrt(np.diagonal(Sigma_xhat_BLU))
print(f'Precision of LS estimates in [m], [m/month], [m], resp.:', std_xhat_LS)
print(f'Precision of BLU estimates in [m], [m/month], [m], resp.:', std_xhat_BLU)
Precision of LS estimates in [m], [m/month], [m], resp.: [0.602 0.138 0.673] Precision of BLU estimates in [m], [m/month], [m], resp.: [0.545 0.137 0.498]
Task 2.3: We are mostly interested in the melting rate. Discuss the estimated melting rate with respect to its precision.
Write your answer here.
The melting rate is approximately -1.03 m/month (estimated with BLU), and its precision is around 0.14 m/month, so the estimated signal is large compared to the uncertainty (signal-to-noise ratio is about 1.03/0.14 = 7.5).
3. Residuals and plot with observations and fitted models.¶
Run the code below (no changes needed) to calculate the weighted squared norm of residuals with both estimation methods, and create plots of the observations and fitted models.
eTe_LS = (y - A @ xhat_LS).T @ (y - A @ xhat_LS)
eTe_BLU = (y - A @ xhat_BLU).T @ inv_Sigma_Y @ (y - A @ xhat_BLU)
print(f'Weighted squared norm of residuals with LS estimation: {eTe_LS:.3f}')
print(f'Weighted squared norm of residuals with BLU estimation: {eTe_BLU:.3f}')
plt.figure()
plt.rc('font', size=14)
plt.plot(times, y, 'kx', label='observations')
plt.plot(times, A @ xhat_LS, color='r', label='LS')
plt.plot(times, A @ xhat_BLU, color='b', label='BLU')
plt.xlim(-0.2, (number_of_observations - 1) + 0.2)
plt.xlabel('time [months]')
plt.ylabel('height [meters]')
plt.legend(loc='best');
Weighted squared norm of residuals with LS estimation: 10.438 Weighted squared norm of residuals with BLU estimation: 8.147
Task 3.1:
- Explain the difference between the fitted models.
- Can we see from this figure which model fits better (without information about the stochastic model)?
Write your answer here.
- With LS we give equal weight to all observations, while with BLU the first six observations get larger weight and therefore the fitted line is closer to those observations such that the corresponding residuals are smaller.
- By only looking at figure, we cannot conclude which model fits better, since it does not show that some observations are more precise than others. Therefore we often add error bars.
4. Confidence bounds¶
In the code below you need to calculate the confidence bounds, e.g., CI_yhat_BLU is a vector with the values $k\cdot\sigma_{\hat{y}_i}$ for BLUE. This will then be used to plot the confidence intervals:
$$
\hat{y}_i \pm k\cdot\sigma_{\hat{y}_i}
$$
Recall that $k$ can be calculated from $P(Z < k) = 1-\frac{1}{2}\alpha$.
$\textbf{Task 4.1}$
Complete the code below to calculate the 98% confidence intervals of both the observations $y$ and the adjusted observations $\hat{y}$ (model values for both LS and BLU).
Use norm.ppf to compute $k_{98}$. Also try different percentages.
# yhat_LS = YOUR_CODE_HERE
# Sigma_Yhat_LS = YOUR_CODE_HERE
# yhat_BLU = YOUR_CODE_HERE
# Sigma_Yhat_BLU = YOUR_CODE_HERE
# SOLUTION:
yhat_LS = A @ xhat_LS
Sigma_Yhat_LS = A @ Sigma_xhat_LS @ A.T
yhat_BLU = A @ xhat_BLU
Sigma_Yhat_BLU = A @ Sigma_xhat_BLU @ A.T
# alpha = YOUR_CODE_HERE
# k98 = YOUR_CODE_HERE
# SOLUTION:
alpha = 0.02
k98 = norm.ppf(1 - 0.5*alpha)
# CI_y = YOUR_CODE_HERE
# CI_yhat_LS = YOUR_CODE_HERE
# CI_yhat_BLU = YOUR_CODE_HERE
# SOLUTION:
CI_y = k98 * np.sqrt(np.diagonal(Sigma_Y))
CI_yhat_LS = k98 * np.sqrt(np.diagonal(Sigma_Yhat_LS))
CI_yhat_BLU = k98 * np.sqrt(np.diagonal(Sigma_Yhat_BLU))
Tip: don't forget to include a factor 0.5 (or 2) in your calculation of k98 when using the inverse CDF norm.ppf.
You can directly run the code below to create the plots.
plt.figure(figsize = (10,4))
plt.rc('font', size=14)
plt.subplot(121)
plt.plot(times, y, 'kx', label='observations')
plt.plot(times, yhat_LS, color='r', label='LS')
plt.plot(times, yhat_LS + CI_yhat_LS, 'r:', label=f'{100*(1-alpha):.1f}% conf.')
plt.plot(times, yhat_LS - CI_yhat_LS, 'r:')
plt.xlabel('time [months]')
plt.ylabel('height [meters]')
plt.legend(loc='best')
plt.subplot(122)
plt.plot(times, y, 'kx', label='observations')
plt.errorbar(times, y, yerr = CI_y, fmt='', capsize=5, linestyle='', color='blue', alpha=0.6)
plt.plot(times, yhat_BLU, color='b', label='BLU')
plt.plot(times, yhat_BLU + CI_yhat_BLU, 'b:', label=f'{100*(1-alpha):.1f}% conf.')
plt.plot(times, yhat_BLU - CI_yhat_BLU, 'b:')
plt.xlim(-0.2, (number_of_observations - 1) + 0.2)
plt.xlabel('time [months]')
plt.legend(loc='best');
Task 4.2: Discuss the shape of the confidence bounds. Do you think the model (linear trend + annual signal) is a good choice?
Write your answer here.
The LS method does not take into account the precision of the observations; this is why there are no error bars plotted on the left figure.
By including the precision of the observations (the error bars in the right figure), as done in the BLU method, you can better assess whether the fit is good.
With BLU: higher precision of first 6 observations, results in tighter confidence intervals as expected.
Similarly as with a linear trend, you can observe a 'widening' towards to end, the uncertainty in the fitted model is increasing since there is less 'support' from observations on both sides.
Overall the fitted model seems a good choice: it follows the observations nicely and confidence region seems not too wide compared to the fitted model trend.
Note that for the LS case we use the covariance matrix of the observations; if we had this information to begin with in practice we would not want to use the LS method, as BLU can take it into account in a more complete way. We do it here to compare LS and BLU and to illustrate how the computed model and confidence intervals are improved (more narrow CI) when using BLU.
Task 4.3:
What is the BLU-estimated melting rate and the amplitude of the annual signal and their 98% confidence interval? Hint: extract the estimated values and standard deviations from xhat_BLU and std_xhat_BLU, respectively.
# rate =YOUR_CODE_HERE
# CI_rate = YOUR_CODE_HERE
# SOLUTION:
rate = xhat_BLU[1]
CI_rate = k98 * std_xhat_BLU[1]
# amplitude = YOUR_CODE_HERE
# CI_amplitude = YOUR_CODE_HERE
# SOLUTION:
amplitude = xhat_BLU[2]
CI_amplitude = k98 * std_xhat_BLU[2]
print(f'The melting rate is {rate:.3f} ± {CI_rate:.3f} m/month (98% confidence level)')
print(f'The amplitude of the annual signal is {amplitude:.3f} ± {CI_amplitude:.3f} m (98% confidence level)')
The melting rate is -1.043 ± 0.318 m/month (98% confidence level) The amplitude of the annual signal is -3.235 ± 1.158 m (98% confidence level)
Task 4.4: Can we conclude the glacier is melting due to climate change?
Write your answer here.
- the melting rate is signficant
- however, the seasonal signal is quite large compared to the rate and has large uncertainty
- need much longer time span to draw firm conclusions
