session9/sparse_bayesian_LR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variational sparse Bayesian linear regression\n",
    "\n",
    "This notebook is a follow-up of the one on Bayesian regression with linear basis function models. It is assumed that you have already done this previous practical work (actually, you will need the EM algorithm implemented during the previous lab session).\n",
    "\n",
    "To help you answering the questions in this notebook, you can have a look to the following references: \n",
    "\n",
    "- D. G. Tzikas et al., [The variational approximation for Bayesian inference](http://www.ssp.ece.upatras.gr/downloads/articles/galatsanos2008_4.pdf), IEEE Signal Processing Magazine, 2008 (especially pages 9 and 10).\n",
    "\n",
    "- C. M. Bishop & M. E. Tipping, [Variational relevance vector machines](https://arxiv.org/pdf/1301.3838.pdf), in Proceedings of the sixteenth conference on Uncertainty in Artificial Intelligence (UAI), 2000.\n",
    "\n",
    "- Section 2.2 of Thomas Buchgraber [Ph.D. Thesis](https://www2.spsc.tugraz.at/www-archive/downloads/phdthesis-buchgraber_online.pdf) from Graz University of Technology, Austria.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Likelihood model\n",
    "\n",
    "We first recall the observation model for Bayesian linear regression (see the previous notebook for more details). \n",
    "\n",
    "We consider a dataset of $N$ training samples $\\{(x_i, t_i) \\in \\mathbb{R}^2\\}_{i=1}^N$ where\n",
    "\n",
    "$$t_i = \\sum_{j=0}^{M-1}{w_j \\phi_j(x_i)} + \\epsilon_i = \\mathbf{w}^\\top \\boldsymbol\\phi(x_i) + \\epsilon_i, \\qquad i \\in \\{1,...,N\\},$$\n",
    "\n",
    "with \n",
    "\n",
    "- $\\phi_0(x_i) = 1$,\n",
    "\n",
    "\n",
    "- $\\epsilon_i \\overset{i.i.d}{\\sim} \\mathcal{N}(0, \\beta^{-1})$, \n",
    "\n",
    "\n",
    "- $\\boldsymbol\\phi(x_i) = [\\phi_0(x_i), ..., \\phi_M(x_i)]^\\top \\in \\mathbb{R}^M$, \n",
    "\n",
    "\n",
    "- $\\mathbf{w} = [w_0, ..., w_M]^\\top \\in \\mathbb{R}^M$.\n",
    "\n",
    "This linear regression model with basis functions can be rewritten in matrix form as follows:\n",
    "\n",
    "$$ \\mathbf{t} = \\boldsymbol\\Phi \\mathbf{w} + \\boldsymbol\\epsilon, $$\n",
    "\n",
    "where\n",
    "\n",
    "- $\\mathbf{t} = [t_1, ..., t_N]^\\top \\in \\mathbb{R}^N$,\n",
    "\n",
    "\n",
    "- $\\boldsymbol\\epsilon = [\\epsilon_1, ..., \\epsilon_N]^\\top \\in \\mathbb{R}^N$ and $\\boldsymbol\\epsilon \\sim \\mathcal{N}(0, \\beta^{-1}\\mathbf{I})$,\n",
    "\n",
    "\n",
    "- and the *design matrix* $\\boldsymbol\\Phi \\in \\mathbb{R}^{N \\times M}$ is defined as\n",
    "\n",
    "$$\n",
    "\\boldsymbol\\Phi = \n",
    "\\begin{pmatrix}\n",
    "\\phi_0(x_1) &  \\phi_1(x_1) & \\cdots & \\phi_{M-1}(x_1) \\\\ \n",
    "\\phi_0(x_2) &  \\phi_1(x_2) & \\cdots & \\phi_{M-1}(x_2) \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "\\phi_0(x_N) &  \\phi_1(x_N) & \\cdots & \\phi_{M-1}(x_N)\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "Due to the presence of the additive Gaussian noise $\\boldsymbol\\epsilon$, the likelihood model is given by:\n",
    "\n",
    "$$ p(\\mathbf{t} \\lvert \\mathbf{w}, \\beta) = \\mathcal{N}(\\mathbf{t}; \\boldsymbol\\Phi \\mathbf{w}, \\beta^{-1} \\mathbf{I}), $$ \n",
    "\n",
    "where for better readibility we write $p(\\mathbf{t} \\lvert \\mathbf{w}, \\beta)$ instead of $p(\\mathbf{t} \\lvert \\mathbf{x}, \\mathbf{w}, \\beta)$, i.e. we omit to denote the input vector $\\mathbf{x} = [x_1, ..., x_N]^\\top \\in \\mathbb{R}^N$.\n",
    "\n",
    "As suggested by the variables after the conditioning bar in this likelihood model, we consider both $\\mathbf{w}$ and $\\beta$ as *latent random variables*, for which we need to define priors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Priors\n",
    "\n",
    "**Prior over $\\mathbf{w}$** \n",
    "\n",
    "For a Bayesian treatment of linear regression we need a prior probability distribution over model parameters $\\mathbf{w}$. In the lab session, we were using a stationary Gaussian prior distribution for the weights of the linear model. By stationary we mean that all the weights $w_j$, $j \\in \\{1,...,M\\}$ were independent and had the same Gaussian prior $\\mathcal{N}(0, \\alpha^{-1})$, we say that they are identically distributed.\n",
    "\n",
    "However, in many situations, it is important to allow the model to be more flexible. For this reason, a non-stationary Gaussian prior distribution with a distinct precision $\\alpha_j$ for each weight will be considered here:\n",
    "\n",
    "$$ p(\\mathbf{w} \\lvert \\boldsymbol\\alpha) = \\prod_{j=1}^J \\mathcal{N}(w_j; 0, \\alpha_j^{-1}) = \\mathcal{N}(\\mathbf{w}; 0, \\text{diag}(\\boldsymbol\\alpha)^{-1}).$$\n",
    "\n",
    "Notice that the weights $w_j$ are still independent, but they are not identically distributed anymore as each weight has its own precision $\\alpha_j$. In other words, the covariance matrix is diagonal but not proportional to the identity anymore.\n",
    "\n",
    "**Prior over $\\boldsymbol{\\alpha}$**\n",
    "\n",
    "We treat the above precision parameters $\\alpha_j$ as i.i.d latent random variables following a [Gamma distribution](https://en.wikipedia.org/wiki/Gamma_distribution) (the conjugate prior for the precision of the Gaussian):\n",
    "\n",
    "$$ p(\\boldsymbol\\alpha; a, b) = \\prod_{j=1}^J \\mathcal{G}(\\alpha_j; a, b).$$\n",
    "\n",
    "**Prior over $\\beta$**\n",
    "\n",
    "Finally, the noise precision $\\beta$ from the likelihood model is also considered as a latent random variable with a Gamma distribution:\n",
    "\n",
    "$$ p(\\beta; c, d) = \\mathcal{G}(\\beta; c, d).$$\n",
    "\n",
    "The prior hyperparameters $a,b,c,d$ are assumed to be deterministic and known, so that they will be omitted in the following."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Gamma distribution** \n",
    "\n",
    "We recall that the probability density function of a Gamma random variable $x > 0$ with shape and rate parameters $a >0$ and $b > 0$, respectively, is given by:\n",
    "\n",
    "$$ \\mathcal{G}(x; a, b) = \\frac{ b^a x^{a-1} e^{-b x}}{\\Gamma(a)}, $$\n",
    "\n",
    "where $\\Gamma(\\cdot)$ is the [Gamma function](https://en.wikipedia.org/wiki/Gamma_function). Moreover, we have the following useful properties:\n",
    "\n",
    "$$ \\mathbb{E}[x] = \\frac{a}{b}, \\qquad\\qquad \\mathbb{E}[\\ln(x)] = \\psi(a) - \\ln(b), $$\n",
    "\n",
    "where $\\psi(\\cdot)$ is the [digamma function](https://en.wikipedia.org/wiki/Digamma_function)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Complete-data likelihood**\n",
    "\n",
    "The previous model correpsonds to the following factorization of the complete-data likelihood:\n",
    "\n",
    "$$ \n",
    "\\begin{aligned}\n",
    "p(\\mathbf{t}, \\mathbf{w}, \\boldsymbol\\alpha, \\beta) &= p(\\mathbf{t} \\lvert \\mathbf{w}, \\beta) p(\\mathbf{w} \\lvert \\boldsymbol\\alpha) p(\\boldsymbol\\alpha) p(\\beta) \\\\\\\\\n",
    "&= \\mathcal{N}(\\mathbf{t}; \\boldsymbol\\Phi \\mathbf{w}, \\beta^{-1} \\mathbf{I}) \\mathcal{N}(\\mathbf{w}; 0, \\text{diag}(\\boldsymbol\\alpha)^{-1}) \\left( \\prod_{j=1}^J \\mathcal{G}(\\alpha_j; a, b) \\right) \\mathcal{G}(\\beta; c, d) \\\\\\\\\n",
    "&= \\left( \\prod_{i=1}^N\\mathcal{N}(t_i; \\mathbf{w}^\\top \\boldsymbol\\phi(x_i), \\beta^{-1} ) \\right) \\left(\\prod_{j=1}^J \\mathcal{N}(w_j; 0, \\alpha_j^{-1})\\right) \\left( \\prod_{j=1}^J \\mathcal{G}(\\alpha_j; a, b) \\right) \\mathcal{G}(\\beta; c, d).\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "--- \n",
    "\n",
    "**Exercise 1**\n",
    "\n",
    "**Q1.1**  \n",
    "\n",
    "Summarize in one or two sentences the key differences between this model and the one in the previous lab session (in terms of what we treat as latent variables or deterministic parameters, and in terms of the priors)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q1.2**  \n",
    "\n",
    "Draw the Bayesian network associated with this model (using the plate notation)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Posterior\n",
    "\n",
    "Bayesian inference requires the computation of the posterior distribution of all latent variables given the observations:\n",
    "\n",
    "$$ p(\\mathbf{w}, \\boldsymbol\\alpha, \\beta | \\mathbf{t}) = \\frac{p(\\mathbf{t} \\lvert \\mathbf{w}, \\beta) p(\\mathbf{w} \\lvert \\boldsymbol\\alpha) p(\\boldsymbol\\alpha) p(\\beta) }{p(\\mathbf{t})}. $$\n",
    "\n",
    "Unfortunately, the marginal likelihood $p(\\mathbf{t}) = \\int p(\\mathbf{t} \\lvert \\mathbf{w}, \\beta) p(\\mathbf{w} \\lvert \\boldsymbol\\alpha) p(\\boldsymbol\\alpha) p(\\beta) d\\mathbf{w} d\\boldsymbol\\alpha d\\beta$ cannot be computed analytically, and therefore the posterior either. \n",
    "\n",
    "Therefore, we will resort to approximate Bayesian inference, and more precisely variational inference.\n",
    "\n",
    "We define the following mean field approximation of the intractable posterior:\n",
    "\n",
    "$$ q(\\mathbf{w}, \\boldsymbol\\alpha, \\beta) = q(\\mathbf{w}) q(\\boldsymbol\\alpha) q(\\beta). $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "**Exercise 2 - Finding the optimal variational distributions**\n",
    "\n",
    "**Q2.1** - Applying the variational inference recipe seen during the lecture, give the expression of the parameters of $q^\\star(\\mathbf{w}) = \\mathcal{N}\\left(\\mathbf{w} ; \\boldsymbol\\mu_N, \\boldsymbol\\Sigma_N \\right)$. \n",
    "\n",
    "You \"simply\" have to develop and identify\n",
    "\n",
    "$$ \n",
    "\\begin{aligned}\n",
    "\\ln q^\\star(\\mathbf{w}) &= \\mathbb{E}_{q(\\boldsymbol\\alpha)q(\\beta)}[\\ln p(\\mathbf{t}, \\mathbf{w}, \\boldsymbol\\alpha, \\beta)] + cst \\\\\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Hint: $\\boldsymbol\\Sigma_N$ should depend on $\\mathbb{E}_{q(\\boldsymbol\\alpha)}[\\text{diag}\\{\\boldsymbol\\alpha\\}]$ and $\\mathbb{E}_{q(\\beta)}[\\beta]$, and $\\boldsymbol\\mu_N$ should depend on $\\mathbb{E}_{q(\\beta)}[\\beta]$. These expectations will be computed later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q2.2** - Similarly, give the expression for the parameters of $q^\\star(\\boldsymbol\\alpha) = \\prod_{j=1}^M \\mathcal{G}(\\alpha_j; \\tilde{a}, \\tilde{b}_j)$. Develop and identify \n",
    "\n",
    "$$ \n",
    "\\begin{aligned}\n",
    "\\ln q^\\star(\\boldsymbol\\alpha) &= \\mathbb{E}_{q(\\mathbf{w})q(\\beta)}[\\ln p(\\mathbf{t}, \\mathbf{w}, \\boldsymbol\\alpha, \\beta)] + cst \\\\\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Hint: $\\tilde{b}_j$ should depend on $\\mathbb{E}_{q(\\mathbf{w})}[w_j^2]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q2.3** - Similarly, give the expression for the parameters of $q^\\star(\\beta) = \\mathcal{G}(\\beta; \\tilde{c}, \\tilde{d})$. Develop and identify \n",
    "\n",
    "$$ \n",
    "\\begin{aligned}\n",
    "\\ln q^\\star(\\beta) &= \\mathbb{E}_{q(\\mathbf{w})q(\\boldsymbol\\alpha)}[\\ln p(\\mathbf{t}, \\mathbf{w}, \\boldsymbol\\alpha, \\beta)] + cst \\\\\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Hint: $\\tilde{d}$ should depend on $\\mathbb{E}_{q(\\mathbf{w})}[ \\parallel \\mathbf{t} - \\mathbf{\\Phi} \\mathbf{w} \\parallel_2^2]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Computing the required expectations**\n",
    "\n",
    "All the moments that you need to compute for the above variational parameters are obtained by using properties of the Gaussian and Gamma distributions:\n",
    "\n",
    "$\\mathbb{E}_{q(\\boldsymbol\\alpha)}[\\text{diag}\\{\\boldsymbol\\alpha\\}] = \\text{diag}\\{\\tilde{a}/\\tilde{b}_1, ..., \\tilde{a}/\\tilde{b}_M\\} \\tag{1}$\n",
    "$\\mathbb{E}_{q(\\beta)}[\\beta] = \\tilde{c}/\\tilde{d} \\tag{2}$\n",
    "$\\mathbb{E}_{q(\\mathbf{w})}[w_j^2] = [\\boldsymbol\\mu_N]_j^2 + [\\boldsymbol\\Sigma_N]_{j,j} \\tag{3}$\n",
    "$\\mathbb{E}_{q(\\mathbf{w})}[ \\parallel \\mathbf{t} - \\mathbf{\\Phi} \\mathbf{w} \\parallel_2^2] = \\parallel \\mathbf{t} - \\mathbf{\\Phi} \\boldsymbol\\mu_N \\parallel_2^2 + tr\\left(\\mathbf{\\Phi} \\boldsymbol\\Sigma_N \\mathbf{\\Phi}^\\top\\right) \\tag{4}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Evidence lower bound**\n",
    "\n",
    "The evidence lower bound (ELBO) is here given by:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathcal{L}(q(\\mathbf{w}), q(\\boldsymbol\\alpha), q(\\beta)) &= \\mathbb{E}_{q(\\mathbf{w}) q(\\boldsymbol\\alpha) \n",
    "q(\\beta)}\\left[\\ln p(\\mathbf{t} \\lvert \\mathbf{w}, \\beta) + \\ln p(\\mathbf{w} \\lvert \\boldsymbol\\alpha) + \\ln p(\\boldsymbol\\alpha) + \\ln p(\\beta)  \\right] \\\\\\\\\n",
    "&- \\mathbb{E}_{q(\\mathbf{w})} \\left[ \\ln q(\\mathbf{w}) \\right] - \\mathbb{E}_{q(\\boldsymbol\\alpha)}\\left[ \\ln q(\\boldsymbol\\alpha) \\right] - \\mathbb{E}_{q(\\beta)}\\left[ \\ln q(\\beta)) \\right] . \n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "This is the objective function which is cyclically maximized when computing $q^\\star(\\mathbf{w})$, $q^\\star(\\boldsymbol\\alpha)$ and $q^\\star(\\beta)$ (cf. the lecture on variational inference). Developing the above expression is not very difficult, but it is a bit painful, so we will give you a function to evaluate the ELBO.\n",
    "\n",
    "As seen during the lecture, the ELBO should be monotonically increasing after each update of the variational distribution factors, which is very useful to check that your implementation of variational inference does not contain bugs. \n",
    "\n",
    "Note that the hyperparameters $a, b, c, d$ are assumed to be known which is why they are not indicated as parameters of the ELBO. \n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3 - Implementation**\n",
    "\n",
    "Complete the function ```VI``` to cyclically update the variational parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy as sp\n",
    "import numpy as np\n",
    "np.random.seed(0)\n",
    "from utils import compute_ELBO\n",
    "\n",
    "def VI(Phi, t, a=1e-5, b=1e-5, c=1e-5, d=1e-5, max_iter=200, rtol=1e-5, verbose=False):\n",
    "    \"\"\"\n",
    "    \n",
    "    Args:\n",
    "    \n",
    "    Hyperparameters a, b, c and d\n",
    "    Maximum number of iterations\n",
    "    Tolerance to consider convergence\n",
    "    Verbose\n",
    "        \n",
    "    Returns:\n",
    "        a_tilde, b_tilde, c_tilde, d_tilde, mu_N, Sigma_N, ELBO.\n",
    "    \"\"\"\n",
    "    \n",
    "    N, M = Phi.shape\n",
    "\n",
    "    # Initialize variational parameters\n",
    "    a_tilde = ? # scalar\n",
    "    b_tilde = ? # vector of dimension M\n",
    "    c_tilde = ? # scalar\n",
    "    d_tilde = ? # scalar\n",
    "    mu_N = ? # vector of dimension M\n",
    "    Sigma_N = ? # matrix of dimension M x M\n",
    "    \n",
    "    # Intialize required expectations using the above formulas.\n",
    "    # Please avoid crazy loops:\n",
    "    # Use @ for matrix multiplication, np.trace to compute the trace, np.sum to compute the sum,\n",
    "    # np.diag to extract the diagonal of a matrix, .T to transpose, \n",
    "    # np.newaxis to add a dimension to a numpy array (see the 'numpy new axis trick.ipynb' in the GMM lab files).\n",
    "    exp_alpha = ? # vector of dimension M (diagonal of the matrix), equation (1)\n",
    "    exp_beta = ? # scalar, equation (2)\n",
    "    exp_w = ? # vector of dimension M, equation (3)\n",
    "    exp_w_epsilon = # scalar, equation (4)\n",
    "    \n",
    "    \n",
    "    ELBO = []\n",
    "        \n",
    "    for i in range(max_iter):\n",
    "        \n",
    "        # store current value of parameters\n",
    "        b_tilde_prev = b_tilde\n",
    "        d_tilde_prev = d_tilde\n",
    "        mu_N_prev = mu_N\n",
    "        Sigma_N_prev = Sigma_N\n",
    "        \n",
    "        # E-w-step\n",
    "        Sigma_N_inv = ?\n",
    "        Sigma_N = np.linalg.inv(Sigma_N_inv)\n",
    "        mu_N = ?\n",
    "        \n",
    "        # update expectations wrt q(w)\n",
    "        exp_w_epsilon = ?\n",
    "        exp_w = ?\n",
    "        \n",
    "        # E-alpha-step\n",
    "        b_tilde = ?\n",
    "        \n",
    "        # update expectation wrt q(alpha)\n",
    "        exp_alpha = ?\n",
    "        \n",
    "        # E-beta-step\n",
    "        d_tilde = ?\n",
    "        \n",
    "        # update expectation wrt q(beta)\n",
    "        exp_beta = ?\n",
    "        \n",
    "        # Compute the ELBO\n",
    "        ELBO.append(compute_ELBO(a, b, c, d, a_tilde, b_tilde, c_tilde, d_tilde, mu_N, Sigma_N, t, Phi))\n",
    "\n",
    "        if (np.isclose(b_tilde_prev, b_tilde, rtol=rtol).all() and \n",
    "            np.isclose(d_tilde_prev, d_tilde, rtol=rtol) and\n",
    "            np.isclose(mu_N_prev, mu_N, rtol=rtol).all() and\n",
    "            np.isclose(Sigma_N_prev, Sigma_N, rtol=rtol).all()):\n",
    "            if verbose:\n",
    "                print(f'Convergence after {i + 1} iterations.')\n",
    "            return a_tilde, b_tilde, c_tilde, d_tilde, mu_N, Sigma_N, ELBO\n",
    "\n",
    "    if verbose:\n",
    "        print(f'Stopped after {max_iter} iterations.')\n",
    "    \n",
    "    return a_tilde, b_tilde, c_tilde, d_tilde, mu_N, Sigma_N, ELBO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generate a sinusoidal training dataset of size 30 with variance $\\beta^{-1} = 0.3^2$ and then use `VI` to obtain the approximate posterior over $\\mathbf{w}$, $\\alpha$ and $\\beta$. The used regression model is a polynomial model of degree 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import g, expand, polynomial_basis_function\n",
    "\n",
    "N = 30\n",
    "\n",
    "degree = 4\n",
    "\n",
    "X = np.linspace(0, 1, N).reshape(-1, 1)\n",
    "t = g(X, noise_variance=0.3 ** 2)\n",
    "\n",
    "Phi = expand(X, bf=polynomial_basis_function, bf_args=range(1, degree + 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_tilde, b_tilde, c_tilde, d_tilde, mu_N, Sigma_N, ELBO = VI(Phi, t, a=1e-5, b=1e-5, c=1e-5, d=1e-5, \n",
    "                                                              max_iter=200, rtol=1e-2, verbose=True)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(ELBO)\n",
    "plt.title('ELBO')\n",
    "plt.xlabel('iterations')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ELBO should be monotonically increasing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Posterior predictive distribution**\n",
    "\n",
    "With this variational inference approach, the posterior predictive distribution is given by:\n",
    "\n",
    "$$\n",
    "p(t_{\\star} \\lvert \\mathbf{t}) = \\int{p(t_\\star \\lvert \\mathbf{w}, \\beta) q(\\mathbf{w}) q(\\beta) d\\mathbf{w} d\\beta}.\n",
    "$$\n",
    "\n",
    "Unfortunately, this integral is intractable, but it can be reasonably approximated by:\n",
    "\n",
    "$$\n",
    "p(t_{\\star} \\lvert \\mathbf{t}) \\approx \\int{p(t_\\star \\lvert \\mathbf{w}, \\hat{\\beta}) q(\\mathbf{w}) d\\mathbf{w}} = \\mathcal{N}\\left(t ; \\boldsymbol\\mu_N^\\top \\boldsymbol\\phi(\\mathbf{x}_{\\star}),\\hat{\\beta}^{-1} + \\boldsymbol\\phi(\\mathbf{x})^\\top \\mathbf{\\Sigma}_N \\boldsymbol\\phi(\\mathbf{x})\\right),\n",
    "$$\n",
    "\n",
    "where $\\hat{\\beta} = \\mathbb{E}_{q(\\beta)}[\\beta] = \\tilde{c}/\\tilde{d}.$\n",
    "\n",
    "This posterior predictive ends up to be very similar to equations (16)-(17) in the previous lab session, except that we replaced $\\mathbf{m}_N$ by $\\boldsymbol\\mu_N$, $\\mathbf{S}_N$ by $\\mathbf{\\Sigma}_N $ and $\\beta$ by $\\hat{\\beta}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test observations\n",
    "X_test = np.linspace(0, 1, 100).reshape(-1, 1)\n",
    "\n",
    "# Function values without noise \n",
    "y_true = g(X_test, noise_variance=0)\n",
    "    \n",
    "# Design matrix of test observations\n",
    "Phi_test = expand(X_test, bf=polynomial_basis_function, bf_args=range(1, degree + 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import posterior_predictive, plot_data, plot_truth, plot_predictive\n",
    "\n",
    "y_VI, y_var_VI = posterior_predictive(Phi_test, mu_N, Sigma_N, c_tilde/d_tilde)\n",
    "\n",
    "plt.figure()\n",
    "plot_data(X, t)\n",
    "plot_truth(X_test, y_true)\n",
    "plot_predictive(X_test, y_VI, np.sqrt(y_var_VI))\n",
    "plt.ylim(-1.0, 2.0)\n",
    "plt.legend()\n",
    "plt.title('Posterior predictive for the sparse Bayesian LR model')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare with the posterior predictive obtained using the Bayesian linear regression model of the previous lab session.\n",
    "\n",
    "You first have to complete the ```EM``` function in ```EM.py``` using the code you wrote in the previous lab session. Then, run the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from EM import EM\n",
    "\n",
    "alpha, beta, m_N, S_N, mlls = EM(Phi, t, rtol=1e-5, verbose=True)\n",
    "y_EM, y_var_EM = posterior_predictive(Phi_test, m_N, S_N, beta)\n",
    "\n",
    "plt.figure()\n",
    "plot_data(X, t)\n",
    "plot_truth(X_test, y_true)\n",
    "plot_predictive(X_test, y_EM, np.sqrt(y_var_EM))\n",
    "plt.ylim(-1.0, 2.0)\n",
    "plt.legend()\n",
    "plt.title('Posterior predictive for the Bayesian LR model')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now look at the estimated weights given by the posterior means ```m_N``` and ```mu_N```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(m_N, 'o')\n",
    "plt.plot(mu_N, 'x')\n",
    "plt.title(\"weights\")\n",
    "plt.xlabel(\"weight index\")\n",
    "plt.legend(['Bayesian LR', 'sparse Bayesian LR'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.hist(m_N)\n",
    "plt.title(\"weights histogram for Bayesian LR\")\n",
    "plt.figure()\n",
    "plt.hist(mu_N)\n",
    "plt.title(\"weights histogram for sparse Bayesian LR\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a polynomial model of degree 4, you should not observe much difference between the bayesian LR model (stationary prior for the weights resulting in tracatable posterior inference) and the sparse Bayesian LR model (non-stationary prior for the weights resulting in approximate posterior computation using variational inference). \n",
    "\n",
    "**Q3.1** Increase the polynomial degree to 40, re-run the algorithms, and compare the accuracy of the two models in terms of root mean square error. What do you obtain?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**A3.1**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q3.2** For the two models, compare the number of weights that have a posterior mean (or approximate posterior mean) greater that $10^{-3}$. Based on what you observe, explain why the model we consider in this notebook (with a non-stationary prior for the regression weights) is referred to as \"sparse Baysian linear regression\" or \"automatic relevance determination\"?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**A3.2**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q3.3** The prior over a weight $w_j$ is conditionally Gaussian, with a variance $\\alpha_j^{-1}$ following an Inverse-Gamma distribution (or equivalently a precision $\\alpha_j$ following a Gamma distribution). We already encountered this hierarchical prior during one of the lectures. What is the marginal distribution of $w_j$ and how does its density differ from the Gaussian distribution?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**A3.3**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q3.4** If the deterministic hyperparameters $a,b,c,d$ were unknown, how could we estimate them? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**A3.4**"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}