{ "cells": [ { "cell_type": "markdown", "metadata": { "colab_type": "text", "execution": {}, "id": "view-in-github" }, "source": [ "\"Open   \"Open" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "# Tutorial 4: Autoregressive models\n", "\n", "**Week 2, Day 2: Linear Systems**\n", "\n", "**By Neuromatch Academy**\n", "\n", "**Content Creators**: Bing Wen Brunton, Biraj Pandey\n", "\n", "**Content Reviewers**: Norma Kuhn, John Butler, Matthew Krause, Ella Batty, Richard Gao, Michael Waskom\n", "\n", "**Post-Production Team:** Gagana B, Spiros Chavlis" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "# Tutorial Objectives\n", "\n", "*Estimated timing of tutorial: 30 minutes*\n", "\n", "The goal of this tutorial is to use the modeling tools and intuitions developed in the previous few tutorials and use them to _fit data_. The concept is to flip the previous tutorial -- instead of generating synthetic data points from a known underlying process, what if we are given data points measured in time and have to learn the underlying process?\n", "\n", "This tutorial is in two sections.\n", "\n", "**Section 1** walks through using regression of data to solve for the coefficient of an OU process from Tutorial 3. Next, **Section 2** generalizes this auto-regression framework to high-order autoregressive models, and we will try to fit data from monkeys at typewriters." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Tutorial slides\n", "# @markdown These are the slides for the videos in all tutorials today\n", "from IPython.display import IFrame\n", "link_id = \"snv4m\"\n", "print(f\"If you want to download the slides: https://osf.io/download/{link_id}/\")\n", "IFrame(src=f\"https://mfr.ca-1.osf.io/render?url=https://osf.io/{link_id}/?direct%26mode=render%26action=download%26mode=render\", width=854, height=480)" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "# Setup" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Install and import feedback gadget\n", "\n", "!pip3 install vibecheck datatops --quiet\n", "\n", "from vibecheck import DatatopsContentReviewContainer\n", "def content_review(notebook_section: str):\n", " return DatatopsContentReviewContainer(\n", " \"\", # No text prompt\n", " notebook_section,\n", " {\n", " \"url\": \"https://pmyvdlilci.execute-api.us-east-1.amazonaws.com/klab\",\n", " \"name\": \"neuromatch_cn\",\n", " \"user_key\": \"y1x3mpx5\",\n", " },\n", " ).render()\n", "\n", "\n", "feedback_prefix = \"W2D2_T4\"" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# Imports\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Figure settings\n", "import logging\n", "logging.getLogger('matplotlib.font_manager').disabled = True\n", "import ipywidgets as widgets # interactive display\n", "%config InlineBackend.figure_format = 'retina'\n", "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/main/nma.mplstyle\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Plotting Functions\n", "\n", "def plot_residual_histogram(res):\n", " \"\"\"Helper function for Exercise 4A\"\"\"\n", " plt.figure()\n", " plt.hist(res)\n", " plt.xlabel('error in linear model')\n", " plt.title(f'stdev of errors = {res.std():.4f}')\n", " plt.show()\n", "\n", "\n", "def plot_training_fit(x1, x2, p):\n", " \"\"\"Helper function for Exercise 4B\"\"\"\n", " plt.figure()\n", " plt.scatter(x2 + np.random.standard_normal(len(x2))*0.02,\n", " np.dot(x1.T, p), alpha=0.2)\n", " plt.title(f'Training fit, order {r} AR model')\n", " plt.xlabel('x')\n", " plt.ylabel('estimated x')\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Helper Functions\n", "\n", "def ddm(T, x0, xinfty, lam, sig):\n", " '''\n", " Samples a trajectory of a drift-diffusion model.\n", "\n", " args:\n", " T (integer): length of time of the trajectory\n", " x0 (float): position at time 0\n", " xinfty (float): equilibrium position\n", " lam (float): process param\n", " sig: standard deviation of the normal distribution\n", "\n", " returns:\n", " t (numpy array of floats): time steps from 0 to T sampled every 1 unit\n", " x (numpy array of floats): position at every time step\n", " '''\n", " t = np.arange(0, T, 1.)\n", " x = np.zeros_like(t)\n", " x[0] = x0\n", "\n", " for k in range(len(t)-1):\n", " x[k+1] = xinfty + lam * (x[k] - xinfty) + sig * np.random.standard_normal(size=1)\n", "\n", " return t, x\n", "\n", "def build_time_delay_matrices(x, r):\n", " \"\"\"\n", " Builds x1 and x2 for regression\n", "\n", " Args:\n", " x (numpy array of floats): data to be auto regressed\n", " r (scalar): order of Autoregression model\n", "\n", " Returns:\n", " (numpy array of floats) : to predict \"x2\"\n", " (numpy array of floats) : predictors of size [r,n-r], \"x1\"\n", "\n", " \"\"\"\n", " # construct the time-delayed data matrices for order-r AR model\n", " x1 = np.ones(len(x)-r)\n", " x1 = np.vstack((x1, x[0:-r]))\n", " xprime = x\n", " for i in range(r-1):\n", " xprime = np.roll(xprime, -1)\n", " x1 = np.vstack((x1, xprime[0:-r]))\n", "\n", " x2 = x[r:]\n", "\n", " return x1, x2\n", "\n", "\n", "def AR_prediction(x_test, p):\n", " \"\"\"\n", " Returns the prediction for test data \"x_test\" with the regression\n", " coefficients p\n", "\n", " Args:\n", " x_test (numpy array of floats): test data to be predicted\n", " p (numpy array of floats): regression coefficients of size [r] after\n", " solving the autoregression (order r) problem on train data\n", "\n", " Returns:\n", " (numpy array of floats): Predictions for test data. +1 if positive and -1\n", " if negative.\n", " \"\"\"\n", " x1, x2 = build_time_delay_matrices(x_test, len(p)-1)\n", "\n", " # Evaluating the AR_model function fit returns a number.\n", " # We take the sign (- or +) of this number as the model's guess.\n", " return np.sign(np.dot(x1.T, p))\n", "\n", "def error_rate(x_test, p):\n", " \"\"\"\n", " Returns the error of the Autoregression model. Error is the number of\n", " mismatched predictions divided by total number of test points.\n", "\n", " Args:\n", " x_test (numpy array of floats): data to be predicted\n", " p (numpy array of floats): regression coefficients of size [r] after\n", " solving the autoregression (order r) problem on train data\n", "\n", " Returns:\n", " (float): Error (percentage).\n", " \"\"\"\n", " x1, x2 = build_time_delay_matrices(x_test, len(p)-1)\n", "\n", " return np.count_nonzero(x2 - AR_prediction(x_test, p)) / len(x2)" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "# Section 1: Fitting data to the OU process" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Video 1: Autoregressive models\n", "from ipywidgets import widgets\n", "from IPython.display import YouTubeVideo\n", "from IPython.display import IFrame\n", "from IPython.display import display\n", "\n", "\n", "class PlayVideo(IFrame):\n", " def __init__(self, id, source, page=1, width=400, height=300, **kwargs):\n", " self.id = id\n", " if source == 'Bilibili':\n", " src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'\n", " elif source == 'Osf':\n", " src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'\n", " super(PlayVideo, self).__init__(src, width, height, **kwargs)\n", "\n", "\n", "def display_videos(video_ids, W=400, H=300, fs=1):\n", " tab_contents = []\n", " for i, video_id in enumerate(video_ids):\n", " out = widgets.Output()\n", " with out:\n", " if video_ids[i][0] == 'Youtube':\n", " video = YouTubeVideo(id=video_ids[i][1], width=W,\n", " height=H, fs=fs, rel=0)\n", " print(f'Video available at https://youtube.com/watch?v={video.id}')\n", " else:\n", " video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,\n", " height=H, fs=fs, autoplay=False)\n", " if video_ids[i][0] == 'Bilibili':\n", " print(f'Video available at https://www.bilibili.com/video/{video.id}')\n", " elif video_ids[i][0] == 'Osf':\n", " print(f'Video available at https://osf.io/{video.id}')\n", " display(video)\n", " tab_contents.append(out)\n", " return tab_contents\n", "\n", "\n", "video_ids = [('Youtube', 'VdiVSTPbJ7I'), ('Bilibili', 'BV1fK4y1s7AQ')]\n", "tab_contents = display_videos(video_ids, W=854, H=480)\n", "tabs = widgets.Tab()\n", "tabs.children = tab_contents\n", "for i in range(len(tab_contents)):\n", " tabs.set_title(i, video_ids[i][0])\n", "display(tabs)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Autoregressive_models_Video\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "To see how this works, let's continue the previous example with the drift-diffusion (OU) process. Our process had the following form:\n", "\n", "\\begin{equation}\n", "x_{k+1} = x_{\\infty} + \\lambda(x_k - x_{\\infty}) + \\sigma \\eta\n", "\\end{equation}\n", "\n", "where $\\eta$ is sampled from a standard normal distribution.\n", "\n", "For simplicity, we set $x_\\infty = 0$. Let's plot a trajectory for this process again below. Take note of the parameters of the process because they will be important later." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @markdown Execute to simulate the drift diffusion model\n", "np.random.seed(2020) # set random seed\n", "\n", "# parameters\n", "T = 200\n", "x0 = 10\n", "xinfty = 0\n", "lam = 0.9\n", "sig = 0.2\n", "\n", "# drift-diffusion model from tutorial 3\n", "t, x = ddm(T, x0, xinfty, lam, sig)\n", "\n", "plt.figure()\n", "plt.title('$x_0=%d, x_{\\infty}=%d, \\lambda=%0.1f, \\sigma=%0.1f$' % (x0, xinfty, lam, sig))\n", "plt.plot(t, x, 'k.')\n", "plt.xlabel('time')\n", "plt.ylabel('position x')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "What if we were given these positions $x$ as they evolve in time as data, how would we get back out the dynamics of the system $\\lambda$?\n", "\n", "Since a little bird told us that this system takes on the form\n", "\n", "\\begin{equation}\n", "x_{k+1} = \\lambda x_k + \\eta \\,,\n", "\\end{equation}\n", "\n", "where $\\eta$ is noise from a normal distribution, our approach is to solve for $\\lambda$ as a **regression problem**.\n", "\n", "As a check, let's plot every pair of points adjacent in time ($x_{k+1}$ vs. $x_k$) against each other to see if there is a linear relationship between them." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @markdown Execute to visualize X(k) vs. X(k+1)\n", "# make a scatter plot of every data point in x\n", "# at time k versus time k+1\n", "plt.figure()\n", "plt.scatter(x[0:-2], x[1:-1], color='k')\n", "plt.plot([0, 10], [0, 10], 'k--', label='$x_{k+1} = x_k$ line')\n", "plt.xlabel('$x_k$')\n", "plt.ylabel('$x_{k+1}$')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Hooray, it's a line! This is evidence that the _dynamics that generated the data_ is **linear**. We can now reformulate this task as a regression problem.\n", "\n", "Let $\\mathbf{x_1} = x_{0:T-1}$ and $\\mathbf{x_2} = x_{1:T}$ be vectors of the data indexed so that they are shifted in time by one. Then, our regression problem is\n", "\n", "\\begin{equation}\n", "\\mathbf{x}_2 = \\lambda \\mathbf{x}_1\n", "\\end{equation}\n", "\n", "This model is **autoregressive**, where _auto_ means self. In other words, it's a regression of the time series on itself from the past. The equation as written above is only a function of itself from _one step_ in the past, so we can call it a _first order_ autoregressive model.\n", "\n", "Now, let's set up the regression problem below and solve for $\\lambda.$ We will plot our data with the regression line to see if they agree." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @markdown Execute to solve for lambda through autoregression\n", "# build the two data vectors from x\n", "x1 = x[0:-2]\n", "x1 = x1[:, np.newaxis]**[0, 1]\n", "\n", "x2 = x[1:-1]\n", "\n", "# solve for an estimate of lambda as a linear regression problem\n", "p, res, rnk, s = np.linalg.lstsq(x1, x2, rcond=None)\n", "\n", "# here we've artificially added a vector of 1's to the x1 array,\n", "# so that our linear regression problem has an intercept term to fit.\n", "# we expect this coefficient to be close to 0.\n", "# the second coefficient in the regression is the linear term:\n", "# that's the one we're after!\n", "lam_hat = p[1]\n", "\n", "# plot the data points\n", "fig = plt.figure()\n", "plt.scatter(x[0:-2], x[1:-1], color='k')\n", "plt.xlabel('$x_k$')\n", "plt.ylabel('$x_{k+1}$')\n", "\n", "# plot the 45 degree line\n", "plt.plot([0, 10], [0, 10], 'k--', label='$x_{k+1} = x_k$ line')\n", "\n", "\n", "# plot the regression line on top\n", "xx = np.linspace(-sig*10, max(x), 100)\n", "yy = p[0] + lam_hat * xx\n", "plt.plot(xx, yy, 'r', linewidth=2, label='regression line')\n", "\n", "mytitle = 'True $\\lambda$ = {lam:.4f}, Estimate $\\lambda$ = {lam_hat:.4f}'\n", "plt.title(mytitle.format(lam=lam, lam_hat=lam_hat))\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Pretty cool! So now we have a way to predict $x_{k+1}$ if given any data point $x_k$. Let's take a look at how accurate this one-step prediction might be by plotting the residuals.\n" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Coding Exercise 1: Residuals of the autoregressive model\n", "\n", "*Referred to as exercise 4A in video*\n", "\n", "Plot a histogram of residuals of our autoregressive model, by taking the difference between the _data_ $\\mathbf{x_2}$ and the _model_ prediction. Do you notice anything about the standard deviation of these residuals and the equations that generated this synthetic dataset?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "code", "execution": {} }, "outputs": [], "source": [ "##############################################################################\n", "## Insert your code here take to compute the residual (error)\n", "raise NotImplementedError('student exercise: compute the residual error')\n", "##############################################################################\n", "# compute the predicted values using the autoregressive model (lam_hat), and\n", "# the residual is the difference between x2 and the prediction\n", "res = ...\n", "\n", "# Visualize\n", "plot_residual_histogram(res)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# to_remove solution\n", "\n", "# compute the predicted values using the autoregressive model (lam_hat), and\n", "# the residual is the difference between x2 and the prediction\n", "res = x2 - (lam_hat * x[0:-2])\n", "\n", "# Visualize\n", "with plt.xkcd():\n", " plot_residual_histogram(res)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Residuals_of_the_autoregressive_model_Exercise\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "# Section 2: Higher order autoregressive models\n", "\n", "*Estimated timing to here from start of tutorial: 15 min*" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Video 2: Monkey at a typewriter\n", "from ipywidgets import widgets\n", "from IPython.display import YouTubeVideo\n", "from IPython.display import IFrame\n", "from IPython.display import display\n", "\n", "\n", "class PlayVideo(IFrame):\n", " def __init__(self, id, source, page=1, width=400, height=300, **kwargs):\n", " self.id = id\n", " if source == 'Bilibili':\n", " src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'\n", " elif source == 'Osf':\n", " src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'\n", " super(PlayVideo, self).__init__(src, width, height, **kwargs)\n", "\n", "\n", "def display_videos(video_ids, W=400, H=300, fs=1):\n", " tab_contents = []\n", " for i, video_id in enumerate(video_ids):\n", " out = widgets.Output()\n", " with out:\n", " if video_ids[i][0] == 'Youtube':\n", " video = YouTubeVideo(id=video_ids[i][1], width=W,\n", " height=H, fs=fs, rel=0)\n", " print(f'Video available at https://youtube.com/watch?v={video.id}')\n", " else:\n", " video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,\n", " height=H, fs=fs, autoplay=False)\n", " if video_ids[i][0] == 'Bilibili':\n", " print(f'Video available at https://www.bilibili.com/video/{video.id}')\n", " elif video_ids[i][0] == 'Osf':\n", " print(f'Video available at https://osf.io/{video.id}')\n", " display(video)\n", " tab_contents.append(out)\n", " return tab_contents\n", "\n", "\n", "video_ids = [('Youtube', 'f2z0eopWB8Y'), ('Bilibili', 'BV1si4y1V7Ru')]\n", "tab_contents = display_videos(video_ids, W=854, H=480)\n", "tabs = widgets.Tab()\n", "tabs.children = tab_contents\n", "for i in range(len(tab_contents)):\n", " tabs.set_title(i, video_ids[i][0])\n", "display(tabs)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Monkey_at_a_typewriter_Video\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Now that we have established the autoregressive framework, generalizing for dependence on data points from the past is straightforward. **Higher order** autoregression models a future time point based on _more than one point in the past_.\n", "\n", "In one dimension, we can write such an order-$r$ model as\n", "\n", "\\begin{equation}\n", "x_{k+1} = \\alpha_0 + \\alpha_1 x_k + \\alpha_2 x_{k-1} + \\alpha_3 x_{k-2} + \\dots + \\alpha_{r+1} x_{k-r} \\, ,\n", "\\end{equation}\n", "\n", "where the $\\alpha$'s are the $r+1$ coefficients to be fit to the data available." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "These models are useful to account for some **history dependence** in the trajectory of time series. This next part of the tutorial will explore one such time series, and you can do an experiment on yourself!\n", "\n", "In particular, we will explore a binary random sequence of 0's and 1's that would occur if you flipped a coin and jotted down the flips.\n", "\n", "The difference is that, instead of actually flipping a coin (or using code to generate such a sequence), you -- yes you, human -- are going to generate such a random Bernoulli sequence as best as you can by typing in 0's and 1's. We will then build higher-order AR models to see if we can identify predictable patterns in the time-history of digits you generate." ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "**But first**, let's try this on a sequence with a simple pattern, just to make sure the framework is functional. Below, we generate an entirely predictable sequence and plot it." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "both", "execution": {} }, "outputs": [], "source": [ "# this sequence is entirely predictable, so an AR model should work\n", "monkey_at_typewriter = '1010101010101010101010101010101010101010101010101'\n", "\n", "# Bonus: this sequence is also predictable, but does an order-1 AR model work?\n", "#monkey_at_typewriter = '100100100100100100100100100100100100100'\n", "\n", "# function to turn chars to numpy array,\n", "# coding it this way makes the math easier\n", "# '0' -> -1\n", "# '1' -> +1\n", "def char2array(s):\n", " m = [int(c) for c in s]\n", " x = np.array(m)\n", " return x*2 - 1\n", "\n", "\n", "x = char2array(monkey_at_typewriter)\n", "\n", "plt.figure()\n", "plt.step(x, '.-')\n", "plt.xlabel('time')\n", "plt.ylabel('random variable')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Now, let's set up our regression problem (order 1 autoregression like above) by defining $\\mathbf{x_1}$ and $\\mathbf{x_2}$ and solve it." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# build the two data vectors from x\n", "x1 = x[0:-2]\n", "x1 = x1[:, np.newaxis]**[0, 1]\n", "\n", "x2 = x[1:-1]\n", "\n", "# solve for an estimate of lambda as a linear regression problem\n", "p, res, rnk, s = np.linalg.lstsq(x1, x2, rcond=None)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# take a look at the resulting regression coefficients\n", "print(f'alpha_0 = {p[0]:.2f}, alpha_1 = {p[1]:.2f}')" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Think! 2: Understanding autoregressive parameters\n", "\n", "Do the values we got for $\\alpha_0$ and $\\alpha_1$ make sense? Write down the corresponding autoregressive model and convince yourself that it gives the alternating 0's and 1's we asked it to fit as data." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# to_remove explanation\n", "\"\"\"\n", "The corresponding autoregressive model is:\n", "\n", "x_{k+1} = 0 - x_{k}\n", "\"\"\";" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Understanding_autoregressive_parameters_Discussion\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Truly random sequences of numbers have no structure and should not be predictable by an AR or any other models.\n", "\n", "However, humans are notoriously terrible at generating random sequences of numbers! (Other animals are no better...)\n", "\n", "To test out an application of higher-order AR models, let's use them to **model a sequence of 0's and 1's that a human tried to produce at random**. In particular, I convinced my 9-yr-old monkey to sit at a typewriter (my laptop) and enter some digits as randomly as he is able. The digits he typed in are in the code, and we can plot them as a time series of digits here.\n", "\n", "If the digits really have no structure, then we expect our model to do about as well as guessing, producing an error rate of 0.5. Let's see how well we can do!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "both", "execution": {} }, "outputs": [], "source": [ "# data generated by 9-yr-old JAB:\n", "# we will be using this sequence to train the data\n", "monkey_at_typewriter = '10010101001101000111001010110001100101000101101001010010101010001101101001101000011110100011011010010011001101000011101001110000011111011101000011110000111101001010101000111100000011111000001010100110101001011010010100101101000110010001100011100011100011100010110010111000101'\n", "\n", "# we will be using this sequence to test the data\n", "test_monkey = '00100101100001101001100111100101011100101011101001010101000010110101001010100011110'\n", "\n", "x = char2array(monkey_at_typewriter)\n", "test = char2array(test_monkey)\n", "\n", "## testing: machine generated randint should be entirely unpredictable\n", "## uncomment the lines below to try random numbers instead\n", "# np.random.seed(2020) # set random seed\n", "# x = char2array(np.random.randint(2, size=500))\n", "# test = char2array(np.random.randint(2, size=500))\n", "\n", "plt.figure()\n", "plt.step(x, '.-')\n", "plt.xlabel('time')\n", "plt.ylabel('random variable')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "## Coding Exercise 2: Fitting AR models\n", "\n", "*Referred to in video as exercise 4B*\n", "\n", "Fit a order-5 ($r=5$) AR model to the data vector $x$. To do this, we have included some helper functions, including ``AR_model``.\n", "\n", "We will then plot the observations against the trained model. Note that this means we are using a sequence of the previous 5 digits to predict the next one.\n", "\n", "Additionally, output from our regression model is continuous (real numbers) whereas our data are scalar (+1/-1). So, we will take the sign of our continuous outputs (+1 if positive and -1 if negative) as our predictions to make them comparable with data. Our error rate will simply be the number of mismatched predictions divided by the total number of predictions." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @markdown Execute this cell to get helper function `AR_model`\n", "\n", "def AR_model(x, r):\n", " \"\"\"\n", " Solves Autoregression problem of order (r) for x\n", "\n", " Args:\n", " x (numpy array of floats): data to be auto regressed\n", " r (scalar): order of Autoregression model\n", "\n", " Returns:\n", " (numpy array of floats) : to predict \"x2\"\n", " (numpy array of floats) : predictors of size [r,n-r], \"x1\"\n", " (numpy array of floats): coefficients of length [r] for prediction after\n", " solving the regression problem \"p\"\n", " \"\"\"\n", " x1, x2 = build_time_delay_matrices(x, r)\n", "\n", " # solve for an estimate of lambda as a linear regression problem\n", " p, res, rnk, s = np.linalg.lstsq(x1.T, x2, rcond=None)\n", "\n", " return x1, x2, p\n", "\n", "help(AR_model)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "code", "execution": {} }, "outputs": [], "source": [ "##############################################################################\n", "## TODO: Insert your code here for fitting the AR model\n", "raise NotImplementedError('student exercise: fit AR model')\n", "##############################################################################\n", "# define the model order, and use AR_model() to generate the model and prediction\n", "r = ...\n", "x1, x2, p = AR_model(...)\n", "\n", "# Plot the Training data fit\n", "# Note that this adds a small amount of jitter to horizontal axis for visualization purposes\n", "plot_training_fit(x1, x2, p)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": {} }, "outputs": [], "source": [ "# to_remove solution\n", "\n", "# define the model order, and use AR_model() to generate the model and prediction\n", "r = 5 # remove later\n", "x1, x2, p = AR_model(x, r)\n", "\n", "# Plot the Training data fit\n", "# Note that this adds a small amount of jitter to horizontal axis for visualization purposes\n", "with plt.xkcd():\n", " plot_training_fit(x1, x2, p)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @title Submit your feedback\n", "content_review(f\"{feedback_prefix}_Fitting_AR_models_Exercise\")" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Let's check out how the model does on the test data that it's never seen before!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @markdown Execute to see model performance on test data\n", "\n", "x1_test, x2_test = build_time_delay_matrices(test, r)\n", "plt.figure()\n", "plt.scatter(x2_test+np.random.standard_normal(len(x2_test))*0.02,\n", " np.dot(x1_test.T, p), alpha=0.5)\n", "mytitle = f'Testing fit, order {r} AR model, err = {error_rate(test, p):.3f}'\n", "plt.title(mytitle)\n", "plt.xlabel('test x')\n", "plt.ylabel('estimated x')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Not bad! We're getting errors that are smaller than 0.5 (what we would have gotten by chance).\n", "\n", "Let's now try **AR models of different orders** systematically, and plot the test error of each.\n", "\n", "_Remember_: The model has never seen the test data before, and random guessing would produce an error of $0.5$." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "cellView": "form", "execution": {} }, "outputs": [], "source": [ "# @markdown Execute to visualize errors for different orders\n", "\n", "# range of r's to try\n", "r = np.arange(1, 21)\n", "err = np.ones_like(r) * 1.0\n", "\n", "for i, rr in enumerate(r):\n", " # fitting the model on training data\n", " x1, x2, p = AR_model(x, rr)\n", " # computing and storing the test error\n", " test_error = error_rate(test, p)\n", " err[i] = test_error\n", "\n", "\n", "plt.figure()\n", "plt.plot(r, err, '.-')\n", "plt.plot([1, r[-1]], [0.5, 0.5], 'r--', label='random chance')\n", "plt.xlabel('Order r of AR model')\n", "plt.ylabel('Test error')\n", "plt.xticks(np.arange(0,25,5))\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "Notice that there's a sweet spot in the test error! The 6th order AR model does a really good job here, and for larger $r$'s, the model starts to overfit the training data and does not do well on the test data.\n", "\n", "In summary:\n", "\n", "\"**I can't believe I'm so predictable!**\" - JAB" ] }, { "cell_type": "markdown", "metadata": { "execution": {} }, "source": [ "---\n", "# Summary\n", "\n", "*Estimated timing of tutorial: 30 minutes*\n", "\n", "In this tutorial, we learned:\n", "\n", "* How learning the parameters of a linear dynamical system can be formulated as a regression problem from data.\n", "* Time-history dependence can be incorporated into the regression framework as a multiple regression problem.\n", "* That humans are no good at generating random (not predictable) sequences. Try it on yourself!" ] } ], "metadata": { "colab": { "collapsed_sections": [], "include_colab_link": true, "name": "W2D2_Tutorial4", "provenance": [], "toc_visible": true }, "kernel": { "display_name": "Python 3", "language": "python", "name": "python3" }, "kernelspec": { "display_name": "Python 3", "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.9.17" } }, "nbformat": 4, "nbformat_minor": 0 }