https://colab.research.google.com/drive/1Gdhg31uDsLyZsbpy5FXtvx925ogxGATr

{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "Copy of hw4_yourName.ipynb",
"provenance": [],
"collapsed_sections": []
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "eKIz_ttKiG_1",
"colab_type": "text"
},
"source": [
"# IFT6269 - Homework 4 - Hidden Markov Models\n",
"**Due:** Tuesday, November 26, 2019"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "whO8SD53Y9Vl",
"colab_type": "text"
},
"source": [
"#### Name:\n",
"#### Student ID: \n",
"#### Collaborators: \n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "02fg3bxiZOMv",
"colab_type": "text"
},
"source": [
"## Introduction\n",
"\n",
"The file `EMGaussian.train` contains samples of data $\\{x_t\\}_{t=1}^T$ where $x_t \\in \\math
{R}^2$, with one datapoint per row. `EMGaussian.test` is structured similarly. This is the same data we used in Homework 3, but this time we use an HMM model to account for the possible temporal structure of the data. This means that we now consider each row of the dataset to be a point $x_t \\in \\math
{R}^2$ co
esponding to some temporal process, rather than thinking of them as *independent* samples as we did in the last homework. \n",
"\n",
"We consider the following HMM model: the chain $(z_t)_{t=1}^T$ has $K=4$ possible states, with an initial probability distribution $\\pi\\in\\Delta_4$ and a probability transition matrix $A \\in \\math
{R}^{4 \\times 4}$ where $A_{ij} = p(z_t=i | z_{t-1} = j),$ and conditionally on the cu
ent state $z_t$, we have observations obtained from Gaussian emission probabilities $x_t| (z_t=k) \\sim \\mathcal{N}(x_t | \\mu_k, \\Sigma_k)$. This is thus a generalization of a GMM since we now allow for time dependencie across the latent states $z_t$.\n",
"\n",
"This exercise has several implementation objectives:\n",
"* **Sum-product**: probabilistic inference on the HMM\n",
"* **Expectation-Maximization**: parameter estimation\n",
"* **Vite
i**: decoding.\n",
"\n",
"**Note:** You may use the (*possibly co
ected*) code you created for the previous assignment. Furthermore, notice there are some math questions in this notebook: do not forget to solve them! \n",
"\n",
"### Tasks\n",
"0. Get your own copy of this file via \"File > Save a copy in Drive...\",\n",
"1. Fill your personal information and collaborators at the top of this assignment, and rename the notebook accordingly, e.g., `hw4_thomasBayes.ipynb`\n",
"2. Read the instructions provided on each section and cell carefully,\n",
"4. Complete the exercises in the sections **Sum-product**, **Expectation-Maximization**, **Vite
i**, **Comparing methods** and **What about K?**.\n",
" \n",
"**Important**: You are allowed to collaborate with other students in both the math and coding parts of this assignment. However, the answers provided here must reflect your individual work. For that reason, you are not allowed to share this notebook, except for your submission to the TA for grading. **Don't forget to pin and save the version of the notebook you want to be graded on!**"
]
},
{
"cell_type": "code",
"metadata": {
"id": "S3gk4JHNY9yW",
"colab_type": "code",
"outputId": "85d33c40-d77d-46a2-c73b-d36f8d15e228",
"colab": {
"base_uri": "https:
localhost:8080/",
"height": 252
}
},
"source": [
"!wget http:
www.iro.umontreal.ca/~slacoste/teaching/ift6269/A19/notes/hwk3data.zip\n",
"!unzip hwk3data.zip\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('seaborn-white')\n",
"\n",
"X_train = np.loadtxt(\"/content/hwk3data/EMGaussian.train\")\n",
"X_test = np.loadtxt(\"/content/hwk3data/EMGaussian.test\")"
],
"execution_count": 0,
"outputs": [
{
"output_type": "stream",
"text": [
"--2019-11-26 06:01:46-- http:
www.iro.umontreal.ca/~slacoste/teaching/ift6269/A19/notes/hwk3data.zip\n",
"Resolving www.iro.umontreal.ca (www.iro.umontreal.ca)... 132.204.26.36\n",
"Connecting to www.iro.umontreal.ca (www.iro.umontreal.ca)|132.204.26.36|:80... connected.\n",
"HTTP request sent, awaiting response... 200 OK\n",
"Length: 7269 (7.1K) [application/zip]\n",
"Saving to: ‘hwk3data.zip’\n",
"\n",
"\rhwk3data.zip 0%[ ] 0 --.-KB/s \rhwk3data.zip 100%[===================>] 7.10K --.-KB/s in 0s \n",
"\n",
"2019-11-26 06:01:46 (585 MB/s) - ‘hwk3data.zip’ saved [7269/7269]\n",
"\n",
"Archive: hwk3data.zip\n",
" inflating: hwk3data/EMGaussian.test \n",
" inflating: hwk3data/EMGaussian.train \n"
],
"name": "stdout"
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "cBJGYsmug-ZD",
"colab_type": "text"
},
"source": [
"## Playground\n",
"\n",
"You are allowed to add as many cells and functions as you wish in this section, but not allowed to change the signature (name and inputs) of the functions we provided!"
]
},
{
"cell_type": "code",
"metadata": {
"id": "A2sgnsCmhG9p",
"colab_type": "code",
"colab": {}
},
"source": [
"# ---------------------------------------------------------------------------\n",
"# Code for plotting the results \n",
"# ! DO NOT MODIFY THESE FUNCTIONS !\n",
"# ---------------------------------------------------------------------------\n",
"\n",
"def plot_smoothing(gamma, K=4, time_limit=100):\n",
" plt.figure(figsize=(14, 2*K))\n",
" plt.suptitle('Smoothing probabilities $p(z_t|x_1, ..., x_T)$', fontsize=16)\n",
" for k in range(K):\n",
" plt.subplot(K, 1, 1+k)\n",
" plt.plot(range(1, time_limit+1), gamma[:time_limit, k] )\n",
" plt.ylabel(r'$p(z_t = ' + str(k+1) + ' | x_{1:T})$')\n",
" plt.ylim(0, 1)\n",
" plt.grid(True)\n",
" plt.xlabel('t')\n",
" plt.show()\n",
"\n",
"def plot_labelling(X, labels, mus, title=\"\"):\n",
" shapes = ['o', '*', 'v', '+'] \n",
" colors = [[31, 119, 180], [255, 127, 14], [44, 160, 44], [148, 103, 189],\n",
" [140, 86, 75], [227, 119, 194], [127, 127, 127], [188, 189, 34]]\n",
"\n",
" fig = plt.figure(figsize=(5, 5))\n",
" cs = [colors[int(_) % len(colors)] for _ in labels]\n",
" plt.scatter(X[:, 0], X[:, 1], c=np.a
ay(cs)/255.)\n",
" plt.scatter(mus[:, 0], mus[:, 1], marker='o', c='#d62728')\n",
" plt.xlim(-12, 12), plt.ylim(-12, 12)\n",
" plt.title(title, fontsize=16) \n",
" plt.show()\n",
"\n",
"def plot_dominoes(data):\n",
" # Pick max from data per timestep\n",
" data_maxhot = (data == data.max(axis=1, keepdims=True))\n",
" \n",
" fig, ax = plt.subplots()\n",
" fig.set_size_inches(12, 3)\n",
" ax.pcolor(1 - data_maxhot[:100,::-1].T, cmap=plt.cm.gray, alpha=0.6)\n",
" ax.set_yticks(np.arange(4) + 0.5, minor=False)\n",
" ax.set_yticklabels([4,3,2,1], minor=False)\n",
" plt.grid(True)\n",
" plt.tight_layout()\n",
" \n",
" plt.show()"
],
"execution_count": 0,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "oYffvc08munF",
"colab_type": "text"
},
"source": [
"## Sum-product [15 pts]"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "nFM6NpRKYvKr",
"colab_type": "text"
},
"source": [
"### $\\alpha$, $\\beta$ recursions\n",
"\n",
"Implement the $\\alpha$ and $\\beta$-recursions seen in class (and that can be found in chapter 12 of Mike's book with slightly different notation). Recall that $\\alpha(z_t) := p(z_t, x_{1:t})$, $\\beta(z_t) := p(x_{(t+1):T} | z_t)$. Implement also a function to compute the emission probabilties $\\epsilon_k(x_t) := p(x_t|z_t=k) = \\mathcal{N}(x_t|\\mu_k, \\Sigma_k)$.\n",
"\n",
"For numerical stability reasons, you are expected to implement your algorithms using **log probabilities** unless noted explicitly!\n",
"\n"
]
},
{
"cell_type": "code",
"metadata": {
"id": "4njzQRaJa5f3",
"colab_type": "code",
"colab": {}
},
"source": [
"\n",
"def log_emission(X, mus, sigmas):\n",
" sigma_inverse = None\n",
" p_x = np.nda
ay(K)\n",
" \n",
" diff = X - mus\n",
" \n",
" if len(sigmas.shape) == 3: # general covariance matrix\n",
" \n",
" if sigma_inverse is None:\n",
" sigma_inverse = np.linalg.pinv(sigmas)\n",
" \n",
" for k in range(K):\n",
" p_x[k] = np.exp(-1/2 * np.matmul(diff[k].T, sigma_inverse[k].dot(diff[k])))\n",
" p_x[k] /= (2 * np.pi)**(d/2) * np.linalg.det(sigmas[k])**(1/2)\n",
" \n",
" elif len(sigmas.shape) == 1: # proportional to identity covariance matrix\n",
" \n",
" for k in range(K):\n",
" p_x[k] = np.exp(-1/2 * np.dot(diff[k], diff[k]) / sigmas[k])\n",
" p_x[k] /= (2 * np.pi * sigmas[k])**(d/2)\n",
" \n",
" else:\n",
" raise Exception('Covariance must be of dim 3 (general covariance) or dim 1 (proportional to identity)!')\n",
" \n",
" return p_x\n",
"\n",
"\n",
"\n",
"def log_alpha_recursion(X, A, log_eps , pi):\n",
" pi = π\n",
" T = X.shape[0]\n",
" α = np.nda
ay((T,K))\n",
...