Linear Generative Models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Generative Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook explains the background for generative models of data,\n",
    "explains the assumptions behind linear generative models of data,\n",
    "and then quickly describes several such models:\n",
    "\n",
    "- Probabilistic Principal Components Analysis\n",
    "- Factor Analysis\n",
    "- $\\ell_1$ Sparse Coding\n",
    "- Spike-and-Slab Sparse Coding\n",
    "\n",
    "The emphasis is on generating data according to these models,\n",
    "rather then on describing algorithms for fitting such models to data,\n",
    "since there is a surfeit of resources for the former\n",
    "and a surplus of same for the latter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generative Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A *generative model* of data is a probabilistic model that allows for the generation of samples.\n",
    "\n",
    "For example, one of the simplest generative models of English is that letters and spaces occur randomly with a frequency calculated by counting letter occurrences over a corpus of text. Our model can be represented by a 27-entry vector composed of the frequencies of each letter.\n",
    "\n",
    "Here's a sequence generated from such a model,\n",
    "as first done by Claude Shannon,\n",
    "father of Information Theory, shortly after World War II:\n",
    "\n",
    "OCRO HLO RGWR NMIELWIS EU LL NBNESEBYATH EEI ALHENHTTPA OOBTTVA NAH BRL\n",
    "\n",
    "This looks nothing like English. Anyone who has played Scrabble could've guessed as such: the seven letters one draws from the Scrabble bag look nothing like an English word, usually. As an even more concrete example, we are keenly aware that when a letter \"Q\" appears, the next letter is highly likely to be a \"U\".\n",
    "\n",
    "This might motivate a more complex model, of English as a first-order Markov chain, where the probability of observing a letter is dependent on only the previous letter. Now our model is a $27$ by $27$ matrix, with each row corresponding to the conditional frequency of each letter given the previous letter. This is a more complicated model, and the resulting samples are correspondingly better:\n",
    "\n",
    "ON IE ANTSOUTINYS ARE T INCTORE ST BE SDEAMY ACHIN D ILONASIVE TUCOOWE ATTEASONARE\n",
    "FUSO TIZIN ANDY TOBE SEACE CTISBE IN NO IST LAT WHEY CRATICT FROURE BIRS GROCID\n",
    "PONDENOME OF DEMONSTURES OF THE REPTAGIN IS REGOACTIONA OF CRE\n",
    "\n",
    "If we extend our model to track the last three letters (aka we model English as a third-order Markov Chain, resulting in a $27$ by $27^3\\approx 12000$ matrix), the samples become even better:\n",
    "\n",
    "THE GENERATED JOB PROVIDUAL BETTER TRAND THE DISPLAYED CODE, ABOVERY UPONDULTS WELL\n",
    "THE CODERST N THESTICAL IT DO HOCK BOTHE BERG. INSTATES CONS ERATIONS. NEVER ANY OF\n",
    "PUBLIE AND TO THEORY. ENVTIAL CALLENGAND TO ELAST BENERATED IN WITH PIES AS IS WITH\n",
    "THE\n",
    "\n",
    "Kind of looks like someone mixed in some Old English with a typo-ridden technical report.\n",
    "\n",
    "As we increase the complexity of our model, we generate increasingly plausible samples of the data at the cost of increased computational effort.\n",
    "\n",
    "In this notebook, we'll look at *factorized linear generative models*, some of the simplest models for continuous random variables -- in some sense, the continuous analogues of the marginal and 1st-order Markov Chains above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(All English-modeling examples come from [these lecture notes](http://www.cim.mcgill.ca/~langer/423/lecture14.pdf). See them for more details.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code in this cell is used to implement the linear generative models below.\n",
    "It is included for readers interested in implementation details,\n",
    "but needn't be read for understanding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class GenerativeModel(object):\n",
    "    \"\"\"\n",
    "    A generative model has a method for sampling the hidden variables\n",
    "    and a method for sampling the data, conditioned on the hidden variables.\n",
    "    \n",
    "    Also included are methods for plotting the hidden and data distributions\n",
    "    for the 1 and 2-dimensional cases.\n",
    "    \"\"\"\n",
    "    def __init__(self, conditional_data_sampler, hidden_sampler):\n",
    "        \n",
    "        self.hidden_sampler = hidden_sampler\n",
    "        self.conditional_data_sampler = conditional_data_sampler\n",
    "        self._has_data = False\n",
    "        self._has_plot = False\n",
    "        \n",
    "    def generate_hidden(self,N):\n",
    "        hidden_sample = self.hidden_sampler(N)\n",
    "        return hidden_sample\n",
    "    \n",
    "    def generate_data(self,N):\n",
    "        hidden_sample = self.generate_hidden(N)\n",
    "        data_sample = self.conditional_data_sampler(hidden_sample)\n",
    "        self.hidden_sample = hidden_sample\n",
    "        self.data_sample = data_sample\n",
    "        self._has_data = True\n",
    "        self.N = N\n",
    "    \n",
    "    def plot(self, which=\"both\"):\n",
    "\n",
    "        if not self._has_data:\n",
    "            self.generate_data(50)\n",
    "        \n",
    "        if which == \"both\":\n",
    "            self.plot(which=\"data\")\n",
    "            self.plot(which=\"hidden\")\n",
    "        elif which in [\"data\", \"hidden\"]:\n",
    "            self._plot(which)\n",
    "        else:\n",
    "            raise ValueError(\"which must be data or hidden, but got {0}\".format(which))\n",
    "            \n",
    "    def _plot(self, which):\n",
    "        \n",
    "        if which == \"data\":\n",
    "            values = self.data_sample\n",
    "        elif which == \"hidden\":\n",
    "            values = self.hidden_sample\n",
    "        else:\n",
    "            raise ValueError(\"which must be data or hidden, but got {0}\".format(which))\n",
    "        \n",
    "        if values.shape[0] == 2:\n",
    "            self._plotgrid = sns.jointplot(*values,\n",
    "                           kind='scatter', stat_func=None,\n",
    "                            joint_kws={\"alpha\":0.05, \"s\":48}\n",
    "                          )\n",
    "            self._plotgrid.ax_joint.set_aspect(\"equal\")\n",
    "                \n",
    "        elif values.shape[0] == 1:\n",
    "            plt.figure()\n",
    "            sns.distplot(np.squeeze(values), kde=False)\n",
    "        else:\n",
    "            raise ValueError(\"values have shape[0] of {0}, expected 1 or 2\".format(values.shape[0]))\n",
    "            \n",
    "        self._has_plot = True\n",
    "        title = which.capitalize()+\" Variable Distribution\"\n",
    "        plt.suptitle(title,\n",
    "                    y=1.1,\n",
    "                    weight=\"bold\", size=\"xx-large\")\n",
    "        \n",
    "    def add_solutions(self):\n",
    "        assert self._has_plot, \"plot not created yet\"\n",
    "        if (hasattr(self,\"W\")) & (self.W.shape == (2,1)):\n",
    "            values = self.data_sample\n",
    "            x_min, x_max = (np.min(values[0]), np.max(values[0]))\n",
    "\n",
    "            x_range = np.expand_dims(np.linspace(x_min, x_max, num=3),0)\n",
    "            self._plotgrid.ax_joint.plot(*np.dot(W,x_range),\n",
    "                                         linewidth=2, color='k',\n",
    "                                        label = r'$W\\cdot h$')\n",
    "\n",
    "            self._plotgrid.ax_joint.legend()\n",
    "        else:\n",
    "            raise NotImplementedError(\"couldn't determine how to plot solutions\")\n",
    "            \n",
    "class LinearGenerativeModel(GenerativeModel):\n",
    "    \"\"\"\n",
    "    A linear generative model has a hidden-conditional data sampler\n",
    "    that is parameterized by a matrix W and a noise-covariance matrix sigma.\n",
    "    \n",
    "    For convenience, sigma can also be a scalar, which produces isotropic noise.\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, W, sigma, hidden_sampler):\n",
    "        \n",
    "        conditional_data_sampler = self._make_linear_sampler(W,sigma)\n",
    "        GenerativeModel.__init__(self, conditional_data_sampler, hidden_sampler)\n",
    "        self.W = W\n",
    "        \n",
    "    def _make_linear_sampler(self, W, sigma):\n",
    "        return lambda h: np.dot(W,h) + \\\n",
    "                            np.dot(sigma,\n",
    "                                   np.random.standard_normal(size=(W.shape[0],h.shape[1])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Factorized Linear Generative Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modeling Random Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above, we simplified the challenge of modeling English, a sequence of discrete symbols with unkown structure, by considering models that capture the structure of much simpler sequences: sequences that had no dependence on their past (i.i.d sequences) and then sequences with increasing dependence on their past (Markov chains of increasing order).\n",
    "\n",
    "How do we simplify the challenge of modeling a continuous random variable?\n",
    "\n",
    "First, let's define what it means to model a random variable, as opposed to a sequence.\n",
    "\n",
    "A model of a random variable is a representation of its probability distribution. For example, I might model the random variable \"My Location\" (which takes values in latitude and longitude) by saying that it is approximately a mixture of two Gaussian distributions, one centered at my workplace and the other at my house.\n",
    "\n",
    "Written explicitly:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "p(\\text{lat},\\text{long}) \\approx 0.5 g_{home}(\\text{lat},\\text{long}) + 0.5 g_{work}(\\text{lat},\\text{long})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $g_{location}$ is a Gaussian distribution centered at $location$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If I think of \"I am at home\" and \"I am at work\" as values of a random variable, then we can rewrite the model above. Using $H$ for the random variable, with $h=1$ corresponding to \"I am at home\" and $h=0$ corresponding to \"I am at work\", we get: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "p(\\text{lat},\\text{long}\\lvert H=h) \\approx h \\cdot g_{home}(\\text{lat},\\text{long}) + (1-h)\\cdot g_{work}(\\text{lat},\\text{long})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, this model says that, if we know that I am at home, then"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "p(\\text{lat},\\text{long}\\lvert H=1) \\approx h \\cdot g_{home}(\\text{lat},\\text{long})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and therefore the chance that my latitude and longitude take on any given value is a decreasing function of the distance from my house."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The utility of this model is that I have simplified what could be a very complicated distribution into something extremely simple -- so simple that we can work with it mathematically, with only a few symbols.\n",
    "So long as the assumptions of the model aren't too incorrect, we can work directly with the model to guide our decision-making."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we really believed that this model did a good job of capturing the behavior of human beings,\n",
    "we could use it to uncover the location of their home and their workplace,\n",
    "given a dataset of their location data over time.\n",
    "\n",
    "A model that proposes that the data can be explained by means of unobserved variables is called a *hidden variable model* or a *latent variable model*, because the variables of interest are \"hidden\" from us, or equivalently are \"latent\" in the data. It is sometimes called a *factor model*, since it proposes that unknown factors explain the data. Many generative models are hidden variable models, including all the ones we will consider here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple Hidden Variable Models: Factorial Linear Generative Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we are allowed total freedom in creating a hidden variable model,\n",
    "then the problem becomes uninteresting:\n",
    "every dataset is perfectly explained by a hidden variable model\n",
    "where the hidden variable is the identity or index of the datapoint, $i$,\n",
    "and the conditional distribution is a delta distribution\n",
    "at the value of datapoint $i$.\n",
    "\n",
    "Hidden variable models are useful only when they are constrained.\n",
    "Furthermore, just like the Markov models of English above,\n",
    "simpler models are usually easier to learn, store, and sample from.\n",
    "They might even be so simple that we can make analytic, mathematical statements\n",
    "about their properties!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simplifying Assumption #1 - Factorial Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first major restriction we place on our hidden variables\n",
    "is that they are *independent* of one another.\n",
    "That is,\n",
    "\n",
    "$$\n",
    "p\\left(\\vec{h}\\right) = \\prod_i p\\left(h_i\\right)\n",
    "$$\n",
    "\n",
    "or, in terms of entropies:\n",
    "\n",
    "$$\n",
    "H\\left(\\vec{h}\\right) = \\sum_i H\\left(h_i\\right)\n",
    "$$\n",
    "\n",
    "In this case, sampling from the model becomes easier because\n",
    "the hidden variables can be sampled in parallel from their marginal distributions,\n",
    "rather than sequentially from their joint distribution.\n",
    "Such a model is called a *factorial* model,\n",
    "since the joint probability *factorizes* over,\n",
    "or turns into a product of,\n",
    "the marginal distributions.\n",
    "(This issue didn't come up in the \"location\" model above because\n",
    "the hidden variable was one dimensional)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Relatedly,\n",
    "we also must assume that,\n",
    "once the hidden variables are known,\n",
    "there is no more unknown structure in the observed variables.\n",
    "\n",
    "If there were,\n",
    "then it'd be impossible to know which relationships between observed variables\n",
    "were due to this unknown structure in the observed variables\n",
    "and which were due to the unknown structure in the dependence\n",
    "of the observed variables on the hidden variables.\n",
    "\n",
    "One way to formalize this is to say that the\n",
    "observed variables are *conditionally independent of each other*\n",
    "when conditioned on the hidden variables.\n",
    "In terms of the distribution,\n",
    "this condition becomes:\n",
    "\n",
    "$$\n",
    "p\\left(\\vec{x}\\lvert\\vec{h}\\right) = \\prod_i p\\left(x_i\\lvert\\vec{h}\\right)\n",
    "$$\n",
    "\n",
    "or, in terms of conditional entropies\n",
    "\n",
    "$$\n",
    "H\\left(\\vec{x} \\lvert \\vec{h}\\right) = \\sum_i H\\left(x_i\\lvert\\vec{h}\\right)\n",
    "$$\n",
    "\n",
    "Combining this conditional distribution of the data\n",
    "with the assumed distribution of the hidden variables,\n",
    "we arrive at an expression for the joint distribution:\n",
    "\n",
    "$$\n",
    "p\\left(\\vec{x},\\vec{h}\\right) = \\prod_i p\\left(x_i\\lvert\\vec{h}\\right) \\prod_i p(h_i)\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "H\\left(\\vec{x}, \\vec{h}\\right) = \\sum_i H\\left(x_i\\lvert \\vec{h}\\right) + \\sum_i H\\left(h_i\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Graphical models provide an elegant and intuitive visualization\n",
    "of joint distributions in terms of their conditional independences,\n",
    "or equivalently in terms of which information measures are 0.\n",
    "They are exceedingly useful for reasoning about algorithms\n",
    "for statistical inference and learning,\n",
    "just as\n",
    "[computational graphs](http://colah.github.io/posts/2015-08-Backprop/)\n",
    "are invaluable for reasoning about neural network architectures.\n",
    "See\n",
    "[this monograph](https://people.eecs.berkeley.edu/~wainwrig/Papers/WaiJor08_FTML.pdf)\n",
    "for a thorough introduction\n",
    "or [this blog post](http://charlesfrye.github.io/stats/2017/09/26/discrete-channel-graph-model.html)\n",
    "for a much shorter introduction in the context of information theory.\n",
    "\n",
    "Roughly and in short,\n",
    "we represent the joint probability distribution by a graph\n",
    "whose nodes represent individual random variables\n",
    "and whose edges represent direct dependences between random variables.\n",
    "\n",
    "For a factorial hidden variable model\n",
    "with conditionally independent observed variables,\n",
    "the graphical model is as below,\n",
    "reproduced from\n",
    "[Deep Learning](http://www.deeplearningbook.org/),\n",
    "by Goodfellow, Courville, and Bengio."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![factorized_hidden_variable_model](./factorized_hidden_variable_model.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our assumption that the $h$s are independent of each other\n",
    "is manifest in the lack of edges between $h$ nodes,\n",
    "while the assumption that the $x$s are conditionally independent\n",
    "of each other, given the $h$s,\n",
    "is manifest in the presence of edges from each $h$ to each $x$\n",
    "and the absence of edges between $x$s."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The nature of the dependence between the $x$s and $h$s,\n",
    "represented by the edges connecting those nodes,\n",
    "is left unspecified.\n",
    "\n",
    "Our next two simplifying assumptions will resolve this. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simplifying Assumption #2 - Conditionally Gaussian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the model of my location from above,\n",
    "we assumed that knowing the hidden variable\n",
    "would only reveal my location up to a Gaussian-distributed error.\n",
    "\n",
    "For a variety of reasons,\n",
    "the choice of conditional Gaussianity is a decent one.\n",
    "Briefly:\n",
    "- **From a Bayesian perspective**: \n",
    "we usually describe our errors as having an average value (ideally 0)\n",
    "and some finite spread.\n",
    "If we think of our errors that way,\n",
    "the principle of maximum entropy guides us to\n",
    "choose a Gaussian distribution to model them.\n",
    "- **From a classical perspective**:\n",
    "the effects that give rise to our errors are,\n",
    "we assume, legion but individually minute and unrelated to each other.\n",
    "If that's the case, then the Central Limit Theorem\n",
    "states that they must be distributed approximately as a Gaussian.\n",
    "\n",
    "Gaussians also possess of numerous fabulous properties\n",
    "that make calculating interesting quantities related to them\n",
    "a breeze, which is important for creating a simple model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mathematically, this assumption is:\n",
    "\n",
    "$$\n",
    "p(\\vec{x}\\lvert \\vec{h}) \\propto \\mathrm{e}^{-\\frac{1}{2} \\left(\\vec{x}-f\\left(\\vec{h}\\right)\\right)^T \\Sigma^{-1} \\left(\\vec{x}-f\\left(\\vec{h}\\right)\\right)}\n",
    "$$\n",
    "\n",
    "or, more \"simply\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "-\\log p\\left(\\vec{x}\\lvert \\vec{h}\\right) = -\\frac{1}{2}\n",
    "\\left(\\vec{x}-f\\left(\\vec{h}\\right)\\right)^T \\Sigma^{-1} \\left(\\vec{x}-f\\left(\\vec{h}\\right)\\right) + A(\\Sigma)\n",
    "$$\n",
    "\n",
    "where $f\\left(\\vec{h}\\right)$ is a function that takes in the hidden variables and predicts values of the observed variables, $\\Sigma$ is the covariance matrix of the errors induced by using $f\\left(\\vec{h}\\right)$ to predict $\\vec{x}$, and $A(\\Sigma)$ is a normalizing constant also known as the log-partition function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The expression $\\vec{z}^T P \\vec{z}$ is a common one. If $P$ is a positive semi-definite matrix, as all covariance matrices are, then it represents a kind of \"distance measure\" on the vector $\\vec{z}$. For example, if $P$ is $I$, then we have\n",
    "$$\\vec{z}^T I \\vec{z} = \\vec{z}^T \\vec{z} = \\sum_i z_i^2 $$\n",
    "\n",
    "or the squared length of $\\vec{z}$, aka the squared distance of the tip of the vector from the origin. Different choices of $P$ result in different \"tilted\" versions of the squared Euclidean distance, called *Mahalanobis distances*.\n",
    "See [Wikipedia](https://en.wikipedia.org/wiki/Mahalanobis_distance) for more."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interpreting $-\\log p(x)$\n",
    "as the *surprise*\n",
    "associated with observing the value $x$\n",
    "for a random variable with distribution $p$,\n",
    "as one should\n",
    "(see [this blog post](http://charlesfrye.github.io/stats/2016/03/29/info-theory-surprise-entropy.html) for details),\n",
    "we can read this equation:\n",
    "\n",
    "$$\n",
    "-\\log p\\left(\\vec{x}\\lvert \\vec{h}\\right) = -\\frac{1}{2}\n",
    "\\left(\\vec{x}-f\\left(\\vec{h}\\right)\\right)^T \\Sigma^{-1} \\left(\\vec{x}-f\\left(\\vec{h}\\right)\\right) + A(\\Sigma)\n",
    "$$\n",
    "\n",
    "in English as \"it is less suprising the find the variable closer to the value predicted from the hidden variable\",\n",
    "and the utility and simplicity of the Gaussian assumption is hopefully more clear."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Nota bene*:\n",
    "We have assumed already that the data values,\n",
    "$\\vec{x}$, are conditionally independent given the hidden variables.\n",
    "For Gaussians, independence is synonymous with decorrelation,\n",
    "meaning that the off-diagonal elements of $\\Sigma$\n",
    "(and hence, $\\Sigma^{-1}$)\n",
    "must be $0$.\n",
    "That is, $\\Sigma$ is a *diagonal matrix*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far,\n",
    "we have passed over the function $f$ with no comment.\n",
    "Our final simplifying assumption tackles this piece\n",
    "and completes our description of the relationship between $\\vec{h}$ and $\\vec{x}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simplifying Assumption #3 - Linearity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you're looking for a simple class of functions,\n",
    "one of the first places to look is at the class of *linear* functions.\n",
    "\n",
    "Every linear function of a vector that produces a scalar\n",
    "can be represented as a vector,\n",
    "and application of a linear function to a vector\n",
    "is just vector-vector multiplication, as in\n",
    "\n",
    "$$\n",
    "f\\left(\\vec{x}\\right) = \\vec{v}^T\\vec{x}\n",
    "$$\n",
    "\n",
    "for some $\\vec{v}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Linear functions of vectors that produce vectors can be represented with matrices,\n",
    "which are just stacks of linear functions that each produce scalars\n",
    "(aka row vectors),\n",
    "one for each element of the output vector,\n",
    "and the application of such a linear function\n",
    "is just matrix-vector multiplication.\n",
    "\n",
    "A matrix is, traditionally,\n",
    "represented with a capital letter,\n",
    "so a linear function $f$\n",
    "of a vector $\\vec{x}$\n",
    "that is represented by a matrix $W$\n",
    "can be written\n",
    "\n",
    "$$\n",
    "f(\\vec{x}) = W\\vec{x}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This simple relationship between functions and input to functions is one of the sources of the power and elegance of linear algebra.\n",
    "Entire treatises could be written on the meaning and purpose of linearity,\n",
    "but we will pass over them and just note that linear functions are simple and play nicely with Gaussians."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our linear assumption in hand,\n",
    "we can write the conditional distribution as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "-\\log p\\left(\\vec{x}\\lvert \\vec{h}\\right) = -\\frac{1}{2}\n",
    "\\left(\\vec{x}-W\\vec{h}\\right)^T \\Sigma^{-1} \\left(\\vec{x}-W\\vec{h}\\right) + A(\\Sigma)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $W$ is a matrix relating the hidden and observed variables.\n",
    "Becaues $W$ (linearly) combines together the hidden variables to produce the observed variables,\n",
    "it is also called the *mixing matrix*.\n",
    "In sparse coding, to be described below,\n",
    "it is called the *dictionary*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This expression,\n",
    "combined with our graphical model,\n",
    "is sufficient to describe a family of useful linear generative models.\n",
    "Despite the strictness of our assumptions,\n",
    "the resulting models are already sufficiently complex\n",
    "to necessitate careful analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quick Note on Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cells below will generate data from linear generative models.\n",
    "\n",
    "The code to generate data is organized into a class,\n",
    "`LinearGenerativeModel`,\n",
    "that requires the user to provide"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- a mixing matrix\n",
    "- a conditional variance, either a scalar or a (diagonal) matrix\n",
    "- a method for sampling the hidden variables\n",
    "that takes as an argument the number of samples to be generated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It will generates `N` samples from that model\n",
    "when the method `.generate_data` is called with argument `N`\n",
    "and plot the distribution of the sampled data variables\n",
    "(jointly and marginally)\n",
    "and hidden variable or variables (jointly and marginally)\n",
    "when the method `.plot(\"both\")` is called.\n",
    "\n",
    "Calling `.add_solutions`\n",
    "will plot $Wh$ for a one-dimensional hidden variable.\n",
    "This is the solution that a proper implementation of the algorithm should find,\n",
    "provided it has enough data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fully-Gaussian Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To develop our first set of models,\n",
    "we add an additional simplifying assumption:\n",
    "the hidden variables are also Gaussian distributed.\n",
    "\n",
    "Since linear combinations of Gaussians\n",
    "(aka the result of linear functions applied to Gaussian random variables)\n",
    "are themselves Gaussian,\n",
    "the marginal distributions of the $x$s are Gaussian,\n",
    "and so the linear generative models\n",
    "generate Gaussian-distributed data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Probabilistic PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we assume that all the elements along the diagonal are equal to some value $\\sigma^2$:\n",
    "\n",
    "$$\n",
    "\\Sigma = \\text{diag}(\\sigma^2)\n",
    "$$\n",
    "\n",
    "then the resulting generative model is called *Probabilistic Principal Components Analysis*.\n",
    "$\\text{diag}$, here, like `np.diag`, takes in a value or vector of values and returns the matrix\n",
    "with that value or values along the diagonal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "W = np.asarray([[0.5],\n",
    "                [1.2]])\n",
    "\n",
    "hidden_variable_dimension = W.shape[1]\n",
    "\n",
    "hidden_variance = np.eye(W.shape[1])\n",
    "\n",
    "hidden_sampler = lambda N: np.dot(hidden_variance,\n",
    "                                  np.random.standard_normal(size=(hidden_variable_dimension,N)))\n",
    "\n",
    "PPCA = LinearGenerativeModel(W, 0.5, hidden_sampler)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "PPCA.generate_data(5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "PPCA.plot(\"both\")\n",
    "PPCA.add_solutions()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you're familiar with Principal Components Analysis,\n",
    "the line indicated by $Wh$ should look familiar:\n",
    "it is exactly the principal eigenvector of the covariance matrix of the data.\n",
    "\n",
    "However, probabilistic PCA is a generative model,\n",
    "while PCA is simply a dimensionality reduction technique.\n",
    "The advantage of PCA over pPCA is that its agnosticism\n",
    "towards the distribution of the data makes it robust\n",
    "to variations in that distribution.\n",
    "pPCA, on the other hand,\n",
    "will provide incorrect answers for data distributed differently than expected.\n",
    "The disadvantage of PCA is that it will give incorrect answers\n",
    "for data actually generated according to the pPCA model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Factor Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A closely-related technique to pPCA is\n",
    "*Factor Analysis*.\n",
    "Factor analysis is also a linear generative model\n",
    "with Gaussian-distributed hidden variables,\n",
    "but in this model the conditional, or nuisance, variance\n",
    "$\\Sigma$\n",
    "has possibly different values along the diagonal.\n",
    "\n",
    "The result is that, for a fixed $h$,\n",
    "the values of $x$ are distributed\n",
    "in an ellipse that is axis-aligned.\n",
    "To see this,\n",
    "simply run the code below with the\n",
    "`hidden_variance` set to some small value, like `0.1`.\n",
    "It might be instructive to compare to the same\n",
    "for pPCA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "W = np.asarray([[0.5],\n",
    "                [1.2]])\n",
    "\n",
    "hidden_variable_dimension = W.shape[1]\n",
    "\n",
    "hidden_variance = np.eye(W.shape[1])\n",
    "\n",
    "hidden_sampler = lambda N: np.dot(hidden_variance,\n",
    "                                  np.random.standard_normal(size=(hidden_variable_dimension,N)))\n",
    "\n",
    "sigma = np.asarray([[0.5,0],[0,1]])\n",
    "\n",
    "FA = LinearGenerativeModel(W, sigma, hidden_sampler)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "FA.generate_data(10000)\n",
    "\n",
    "FA.plot(\"both\")\n",
    "FA.add_solutions()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because pPCA and Factor Analysis have different underlying models of the data,\n",
    "they will find different solutions on the same data.\n",
    "See [this set of lecture slides by Geoff Hinton](https://www.cs.toronto.edu/~hinton/csc2515/notes/lec7middle.pdf)\n",
    "for more."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In short,\n",
    "because Factor Analysis presumes that noise acts along the individual axes of $x$,\n",
    "it discounts variability that is only in the direction of one of the $x$ variables,\n",
    "i.e. high variability in one measurement that doesn't predict variability in other measurements.\n",
    "This can sometimes be desirable and sometimes indesirable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Non-Gaussian Models: Sparse Coding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In sparse coding models,\n",
    "we allow our distribution of\n",
    "hidden variables to be non-Gaussian --\n",
    "specifically,\n",
    "we aim for the hidden variables to be *sparse*,\n",
    "in that most of the time,\n",
    "a given hidden variable is exactly or approximately 0.\n",
    "\n",
    "Sparse coding was originally proposed\n",
    "as a model of the behavior of visual neurons.\n",
    "See [this blog post](http://charlesfrye.github.io/FoundationalNeuroscience/48/)\n",
    "for more."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $\\ell_1$ Sparse Coding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In $\\ell_1$ sparse coding, the hidden variables are Laplace-distributed.\n",
    "\n",
    "The Laplace distribution somewhat resembles the Gaussian distribution:\n",
    "\n",
    "$$\n",
    "p(x) \\propto \\mathrm{e}^{-\\lambda\\lvert x - \\mu\\rvert}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just as the log-probability view of the Gaussian revealed a connection to the Euclidean distance,\n",
    "the log-probability view of the Laplace distribution is illuminating:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "-\\log p(x) = \\lvert x-\\mu\\rvert+ A(\\lambda)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, it is more suprising to find $x$ further a way from some position $\\mu$,\n",
    "where now distance is measured using the\n",
    "[Manhattan distance](https://en.wikipedia.org/wiki/Taxicab_geometry),\n",
    "also known as the $\\ell_1$ norm,\n",
    "hence the name for $\\ell_1$ sparse coding."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Nota bene*: the similarity of form is due to the fact that both the Gaussian family of distributions and the Laplace family of distributions are examples of\n",
    "[exponential families](https://en.wikipedia.org/wiki/Exponential_family).\n",
    "These families,\n",
    "which are among the only ones which admit\n",
    "[sufficient statistics](https://en.wikipedia.org/wiki/Sufficient_statistic),\n",
    "are also known as *log-linear* families.\n",
    "They play very nicely with graphical models,\n",
    "as expounded upon in\n",
    "[Wainwright and Jordan's monograph](https://people.eecs.berkeley.edu/~wainwrig/Papers/WaiJor08_FTML.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "W = np.asarray([[1],[2.5]])\n",
    "hidden_sampler = lambda N: np.random.laplace(size=(1,N))\n",
    "\n",
    "SC = LinearGenerativeModel(W,1,hidden_sampler)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "SC.generate_data(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "SC.plot(\"both\")\n",
    "SC.add_solutions()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Spike-And-Slab Sparse Coding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Though convenient for numerous reasons related to its log-linearity,\n",
    "the Laplace distribution does not generate hidden variables that are sparse\n",
    "in the sense that their values are often $0$.\n",
    "In fact,\n",
    "since the Laplace distribution is non-atomic,\n",
    "the probability that a variable with that distribution\n",
    "takes on any particular value is $0$.\n",
    "\n",
    "To rectify this,\n",
    "we can adjust our distribution of each hidden variable\n",
    "to be a mixture distribution:\n",
    "with some probabiliy,\n",
    "the variable is exactly $0$,\n",
    "and with some probability,\n",
    "the variable is drawn from some other distribution,\n",
    "typically Gaussian.\n",
    "\n",
    "The resulting model is called\n",
    "*spike-and-slab sparse coding*,\n",
    "due to the characteristic appearance of its\n",
    "hidden variable distribution,\n",
    "as clearly visible below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "W = np.asarray([[1],[2.5]])\n",
    "\n",
    "def spike_and_slab_sampler(N, sparsity=0.3):\n",
    "    switches = np.random.uniform(low=0.5,high=1.5,size=N).astype(int)\n",
    "    samples = np.zeros((1,N))\n",
    "    \n",
    "    for idx,switch in enumerate(switches):\n",
    "        if switch:\n",
    "            samples[:,idx] = np.random.standard_normal(size=(1,))\n",
    "    \n",
    "    return samples\n",
    "\n",
    "SSSC = LinearGenerativeModel(W,1,spike_and_slab_sampler)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "SSSC.generate_data(1000)\n",
    "\n",
    "SSSC.plot(\"both\")\n",
    "SSSC.add_solutions()"
   ]
  }
 ],
 "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}