Notebook exercises: is my null hypothesis good enough?

Notebook exercises: is my null hypothesis good enough?#

Click –> Live Code on the top right corner of this screen and then wait until all cells are executed.

If you want to download the notebook, select the .zip after clicking on .

In the following exercise, we will look at the curve fitting problem to a synthetic sea level rise dataset. The dataset contains yearly sea level measurements over 20 years (in total 20 observations). The time of observations are given in years as [1, 2, 3, …, 20]. The observations are assumed to be normally distributed and have a precision of \(\sigma=2\) cm. The sea level at time zero is unknown and should also be estimated together with the other parameters.

In this exercise, we will:

  1. Import the data, create the design matrix \(\mathrm{A}\), and the covariance matrix \(\Sigma_Y\) (assuming the linear trend model).

  2. Apply BLUE and overall model test to check whether the null hypothesis (linear trend model) is sufficiently supported by the data.

We will look at 3 cases where the null-hypothesis is rejected:

  1. Dataset is corrupted by one large outlier.

  2. We apply an incorrect functional model.

  3. We use a too optimistic value for the precision.

If we detect that the overall model test (i.e., our null hypothesis) is rejected, we would in practice need to apply the Generalized Likelihood Ratio test to test an alternative hypothesis against the null hypothesis, to decide whether it is sufficiently more likely.

import numpy as np
import matplotlib.pyplot as plt

from scipy.stats import norm
from scipy.stats.distributions import chi2

%matplotlib inline
plt.rcParams.update({'font.size': 14})

def blue(A, y, Sigma_Y):
    """Calculate the estimates of the unknown parameters and the corresponding
    covariance matrix using the BLUE.
    
    Parameters
    ----------
    A: design matrix.
    y: observations.
    Sigma_Y: covariance matrix of the observables.
    
    Returns
    -------
    xhat: estimates of the unknown parameters.
    Sigma_Xhat: covariance matrix of the estimates.
    yhat: adjusted observations.
    Sigma_Yhat: covariance matrix of adjusted observations.
    e_hat: residuals.
    Sigma_Ehat: covariance matrix of residuals.
    """

    Sigma_Xhat = np.linalg.inv(A.T @ np.linalg.inv(Sigma_Y) @ A)
    xhat = Sigma_Xhat @ A.T @ np.linalg.inv(Sigma_Y) @ y
    
    yhat = A @ xhat  
    Sigma_Yhat = A @ Sigma_Xhat @ A.T
    
    e_hat = y - yhat
    Sigma_Ehat = Sigma_Y - Sigma_Yhat
    
    return xhat, Sigma_Xhat, yhat, Sigma_Yhat, e_hat, Sigma_Ehat

def confidence_interval(conf_level, Sigma):
    """Calculate lower and upper bounds of the confidence interval.
    
    Parameters
    ----------
    conf_level: confidence level (=1-alpha), e.g., 0.95.
    Sigma: covariance matrix of the variable for which to calculate the
             confidence interval.

    Returns
    -------
    cint: half of the width of the confidence interval.

    Examples
    --------
    The lower and upper bounds of the confidence interval can be calculated as:
    y_lower = y - cint
    y_upper = y + cint
    """
    
    alpha = 1 - conf_level
    k = norm.ppf(1 - alpha/2)
    cint =  k * (np.sqrt(np.diagonal(Sigma)))

    return cint

def overall_model_test(alpha, ehat, Sigma_Y, q):
    """Performs the overall model test and prints the outcome.
    
    Parameters
    ----------
    alpha: false alarm probability
    ehat: residuals.
    Sigma_Y: covariance matrix of the observables.
    q: degrees of freedom of Chi2-distribution.
    """

    k = chi2.ppf(1 - alpha, q)
    Tq = ehat.T @ np.linalg.inv(Sigma_Y) @ ehat
    
    if Tq < k:
        print(f"(T = {Tq:.1f}) < (K = {k:.1f}), OMT is accepted.")
    else:
        print(f"(T = {Tq:.1f}) > (K = {k:.1f}), OMT is rejected.")

def plot_results(t, A, y, Sigma_Y, conf_level, alpha):#(t, y, yhat, yhat_int, ehat, ehat_int, conf_level, alpha):
    """Produce plots of the (adjusted) observations and the residuals,
    including the confidence bounds.
    
    Parameters
    ----------
    t: times.
    y: observations.
    A: design matrix.
    Sigma_Y: covariance matrix of the observables.
    conf_level: confidence level.
    """
    
    xhat, Sigma_Xhat, yhat, Sigma_Yhat, ehat, Sigma_Ehat = blue(A, y, Sigma_Y)
    
    yhat_int = confidence_interval(conf_level, Sigma_Yhat)
    ehat_int = confidence_interval(conf_level, Sigma_Ehat)
    
    overall_model_test(alpha, ehat, Sigma_Y, len(y)-len(xhat))

    plt.figure(figsize=(12,6))
    
    plt.subplot(121)
    plt.plot(t, y,'r.',markersize=15, label='observations')
    plt.plot(t, yhat, linewidth=3, label='fitted model')
    plt.plot(t, yhat-yhat_int, 'r--', linewidth=2,
             label='{} % conf.'.format(conf_level*100))
    plt.plot(t, yhat+yhat_int, 'r--', linewidth=2)
    plt.grid()
    plt.xlabel('time [years]')
    plt.ylabel('Sea level [cm]')
    plt.legend(loc='best')
    
    plt.subplot(122)
    plt.plot(t, ehat, 'r.',markersize=15, label='residuals')
    plt.plot(t, ehat_int, 'r--', linewidth=2,
             label='{} % conf.'.format(conf_level*100))
    plt.plot(t, - ehat_int, 'r--', linewidth=2)
    plt.grid()
    plt.xlabel('time [years]')
    plt.ylabel('Residuals [cm]')
    plt.legend(loc='best')

Case 0: linear trend model is correct#

First, import the data and compute the number of observations m, vector with observation times t, standard deviation of the measurements std_y, covariance matrix of observations Sigma_Y, design matrix A.

y_0 = np.genfromtxt('./data/Dataset_0.txt')
m = len(y_0)
t = np.arange(1,m+1)

std_y = 2 
Sigma_Y = std_y**2 * np.eye(m)
A = np.column_stack((np.ones(m), t))

The function plot_results will perform BLUE and the overall model test with a false alarm probability alpha, and then plot results with confidence intervals with confidence level conf_level:

plot_results(t, A, y_0, Sigma_Y, conf_level = 0.95, alpha = 0.05)

Case 1: linear trend model, outlier in data#

In Dataset_1 one of the observations is an outlier… We will apply BLUE and perform the overall model test with this new dataset.

y_1 = np.genfromtxt('./data/Dataset_1.txt')
plot_results(t, A, y_1, Sigma_Y, conf_level = 0.95, alpha = 0.05)
Question

Why are more residuals negative in this case?

Answer

The fitted line is pulled towards the large outlier, and therefore more observations will lie below.

CASE 2: linear trend model is too simplistic#

For Dataset_2 the linear trend model assumptions is not correct… We will apply BLUE and perform the overall model test with this new dataset.

y_2 = np.genfromtxt('./data/Dataset_2.txt')
plot_results(t, A, y_2, Sigma_Y, conf_level=0.95, alpha = 0.05)
Questions

How can you see from the residuals plot that the linear model is too simplistic?

What would be a better model? Try yourself by completing the code below.

Answers

There is a clear ‘non-random’ pattern over time.

A second-order polynomial. You can try by inserting the corresponding \(\mathrm{A}\)-matrix A_2 in the next code cell.

A_2 = ?
plot_results(t, A_2, y_2, Sigma_Y, conf_level=0.95, alpha = 0.05)

Solution

A_2 = np.column_stack((np.ones(m), t, t**2))

CASE 3: linear trend model, too optimistic about precision#

The observations in Dataset_3 are less precise than before… We will apply BLUE and perform the overall model test with this new dataset.

y_3 = np.genfromtxt('./data/Dataset_3.txt')
plot_results(t, A, y_3, Sigma_Y, conf_level=0.95, alpha = 0.05)
Question

How can you see that we were too optimistic about the precision of the observations? Play with the factor below to see if indeed the model is accepted if we choose a larger value for the precision.

Solution

Much more than 5% of the residuals is outside the 95% confidence interval, which is determined by the assumed precision of 2 cm. If we would use 2 * Sigma_Y we would see that the model is accepted (hence the precision should have been \(2\sqrt{2}\) cm).

factor = 1
plot_results(t, A, y_3, factor * Sigma_Y, conf_level=0.95, alpha = 0.05)

Attribution

This chapter was written by Sandra Verhagen. Find out more here.

Contents