{ "cells": [ { "cell_type": "markdown", "id": "3a2df1ea-c2aa-49ce-8928-0c65397cc81b", "metadata": {}, "source": [ "DATA-DRIVEN MODEL LEARNING, LAB 2: Trend and Seasonality Learning with Least Squares\n", "\n", "\n", "1. Prepare your functions and data\n", "\n", "a. Store the functions that you created during Lab 1 in a file named \"DataDrivenML_fcts.py\". Import those functions in the preamble.\n", "\n", "b. Load the file containing the life expectancy of French women between 1968 and 2016, \"WomanLifeExpectancy.mat\" (e.g. using \"scipy.io\"), and store the time and measurements in 2 arrays (np): \"time\" and \"meas\"" ] }, { "cell_type": "raw", "id": "485dada9-b14b-4073-866a-01f5b405284e", "metadata": {}, "source": [ "import csv\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "from DataDrivenML_fcts import ???\n", "\n", "????\n", "\n", "mat = scipy.io.loadmat(file)\n", "time = ???\n", "meas = ???" ] }, { "cell_type": "markdown", "id": "96ff6371-41c3-476f-906b-40c3ac5c20f7", "metadata": {}, "source": [ "2. Linear least-squares from QR\n", "\n", "a. Build a Vandermonde matrix corresponding to a 2nd order polynomial function of time and perform QR decomposition of this matrix.\n", "\n", "b. Define a function \"back_substitution(.)\" that solves A x = b for x using the back substitution method. Use this function to estimate theta and compute the resulting $\\hat{y}$.\n", "\n", "c. Print the value of theta and plot the data, estimation and residuals.\n", "\n", "d. Increase the polynomial order and compare the estimated parameters and the residuals: what can you observe?" ] }, { "cell_type": "raw", "id": "46245ed8-ccb9-4523-90ad-70c605e8b9d5", "metadata": {}, "source": [ "???\n", "\n", "def back_substitution(A: np.ndarray, b: np.ndarray) -> np.ndarray:\n", " # Solves Ax = b for x using the back substitution method\n", "???\n", "\n", "???" ] }, { "cell_type": "markdown", "id": "71f7ac32-83ab-47e4-b137-31f9441e4524", "metadata": {}, "source": [ "3. Seasonal pattern estimation of the monthly mean carbon dioxyde measured at Mauna Loa Observatory, Hawaii, USA\n", "\n", "a. load the data contained in \"HawaiiCarbonDioxide.csv\" in vectors for time and meas, and visualize the measurements.\n", "\n", "b. Note that this time-series is more turbulent than the previous one: this suggests to associate a Fourier series for seasonality to the polynomial function of low order, i.e.\n", "$$\n", "f(t,\\theta) = a_0 + \\sum_{j=1}^3 a_j \\cos (j \\omega_p t) + \\sum_{j=1}^3 b_j \\sin (j \\omega_p t) + c_1 t + c_2 t^2, \\quad \\omega_p = \\frac{\\pi}{6 T_s}\n", "$$\n", "\n", "Build the $\\phi$ matrix corresponding to this $f(t,\\theta)$ and find $\\theta$ from QR factorization using back substitution.\n", "\n", "c. Plot the data, the estimated $\\hat{y}$ and the residuals.\n", "\n", "d. Interpret your results." ] }, { "cell_type": "markdown", "id": "ff7a8b14-b3a1-4f26-b74a-447f43b2ae86", "metadata": {}, "source": [ "???" ] }, { "cell_type": "markdown", "id": "06212512-9238-4e27-96f7-5023439996e8", "metadata": {}, "source": [ "4. SSA and SVD for AR models\n", "\n", "In the first lesson, we have seen that a nonparametric model $(y_t^{ssa})_{t\\in\\mathbb{T}}$ containing deterministic components of a time series can be constructed using the basic SSA algorithm. By construction, this signal should satisfy a recurrent equation, i.e. $\\exists \\,\\textbf{r} \\in \\mathbb{R}^{l+1}$, $l\\in \\mathbb{N}^*$ such that:\n", "$$\n", "\\textbf{r}^T \\left[ \\begin{align} y_t^{ssa} \\\\ y_{t-1}^{ssa} \\\\ \\ldots \\\\ y_{t-l}^{ssa} \\end{align} \\right] = 0.\n", "$$\n", "or, equivalently, a parameter vector $\\theta \\in \\mathbb{R}^{l\\times 1}$ exists such that\n", "$$\n", "y_t^{ssa} = \\theta_1 y_{t-1}^{ssa} + \\ldots + \\theta_l y_{t-l}^{ssa}\n", "$$\n", "which corresponds to the classical autoregressive (AR) model\n", "\n", "a. Define $\\textbf{y}^{ssa}$, ${\\theta}$ and ${\\phi}$ such that the problem of finding the AR parameters is equivalent to finding $\\theta$ that minimizes\n", "$$\n", "|| \\textbf{y}^{ssa} - \\phi \\theta ||^2_2\n", "$$\n", "\n", "b. Decompose the measurements signal to capture annual seasonality: build the Hankel matrix, perform the SVD, extract the principal components and reconstruct the signal $\\textbf{y}^{ssa}$.\n", "\n", "c. Build the matrix $\\phi$, perform SVD on this matrix and reconstruct the optimal parameters ${\\theta}$ and the estimate $\\hat{y}$.\n", "\n", "d. Plot the data, the estimated $\\hat{y}$ and the residuals. Print the parameters values.\n", "\n", "e. Compare the residuals with those resulting from QR decomposition.\n", "\n", "d. Interpret your results.\n" ] }, { "cell_type": "raw", "id": "ce3c23bc-2dba-4d14-a5bc-3629672f3fa8", "metadata": {}, "source": [ "???" ] }, { "cell_type": "markdown", "id": "8ebd9c03-5e3e-479a-a4f0-4fb84d388b41", "metadata": {}, "source": [ "5. Condition number\n", "\n", "Consider the hourly readings from a ceramic furnace contained in the file \"TemperaturesCeramicFurnace.mat\". The high variability of this signal motivates the use of high order polynomials. \n", "\n", "a. Load the data in the proper vectors, and build $\\phi$ using a Vandermonde matrix and using Legendre polynomials of 10th order. Note that the time should be normalized on the interval $[-1, 1]$. Visualize these polynomials on a figure.\n", "\n", "b. Compute the condition numbers (e.g. using \"np.linalg.cond\") for each case. What can you conclude?\n", "\n", "c. Compute the least squares estimate using QR factorization for both regressors and conclude." ] }, { "cell_type": "raw", "id": "4b2ee2d0-2313-4998-8466-c311c3e6834b", "metadata": {}, "source": [ "from scipy.special import legendre\n", "\n", "???\n", "\n", "####################################################\n", "def LegendreBasis(time,l) :\n", " ???\n", "\n", "print(\"Condition numbers for Vandermonde: \", np.linalg.cond(phi_v),\", Legendre: \",np.linalg.cond(phi_l))\n", "\n", "???\n", "\n", "####################################################\n", "def LSE_QR(meas,phi) :\n", " ???\n", "\n", "???\n", "\n", "plt.figure(figsize=(15, 7))\n", "plt.plot(time, meas, 'k-^', y_hat_l,'r', y_hat_v, 'g')\n", "plt.legend(['Measurements','Legendre','Vandermonde'])\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "72d2cab2-4c04-4581-b10d-8dc9c4f6c334", "metadata": {}, "source": [ "6. Nonlinear LS\n", "\n", "Consider the Box & Jenkins airline data contained in \"Airpassengers.mat\" (monthly total of international airline passengers), a data set known for exhibiting numerous nonlinear features.\n", "\n", "a. Choose the appropriate periodicity and use SSA for an efficient extraction of the principal components. Plot the singular values and lead singular vectors and conclude on the number of principal components that should be used.\n", "\n", "b. Show the initial and reconstructed time series, along with the residuals.\n", "\n", "c. Note that the first component's behavior is almost linear: approximate it with a 3rd order polynomial. Show the resulting approximate and residuals.\n", "\n", "d. The other significant components have a more nonlinear behavior and a nonlinear model should be used. The impulse response of a second order system (Oppenheim et al., 2014, ch 6):\n", "$$\n", "s_{k,t} = a_k \\frac{\\omega_k}{\\sqrt{1-\\xi_k^2}}e^{-\\xi_k \\omega_k t}\\sin \\left( \\omega_k \\sqrt{1-\\xi_k^2} + \\psi_k \\right),\\, k\\in\\{1,2\\},\n", "$$\n", "is proposed as an interesting candidate to capture the oscillatory behavior. The amplitudes $a_k$, the natural frequencies $\\omega_k$, the damping ratios $\\xi_k$, and the phase shifts $\\psi_k$ have to be determined from the data of the principal components 2 and 3, and 4 and 5. Find the optimal values (in the least squares sense) for these parameters by formulating the problem as done in classe. Use a Levenberg-Marquart algorithm to perform the optimization.\n", "Note that this reference may be usefull for that aim:\n", "https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html\n", "\n", "e. Illustrate the use of different performance measures and discuss your results.\n" ] }, { "cell_type": "raw", "id": "7a0f3646-ff52-4ce1-8bf4-d219d1a27d42", "metadata": {}, "source": [ "???" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.7" } }, "nbformat": 4, "nbformat_minor": 5 }