<commit_before>{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os, time, warnings\n",
    "warnings.filterwarnings( \"ignore\" )\n",
    "# from IPython.core.display import HTML\n",
    "from IPython.display import display, HTML\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "## This is the correct way to use the random number generator,\n",
    "##  since it allows finer control.\n",
    "rand = np.random.RandomState( 321 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# _EM_ and NIST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before proceeding, it seems pedagogically necessary (at least for myself) to revise the **EM**-slgorithm and show its \"correctness\", so to say.\n",
    "\n",
    "The $\\TeX$ markup used here uses the \"align*\" environment and thus should not be viewed though nbViewer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A brief description of the **EM** algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The EM algorithm seeks to maximize the likelihood by means of successive application of two steps: the E-step and the M-step.\n",
    "\n",
    "For any probability measure $Q$ on the space of latent variables $Z$ with density $q$ the following holds:  \n",
    "\\begin{align*}\n",
    "\\log p(X|\\Theta)\n",
    "    &= \\int q(Z) \\log p(X|\\Theta) dZ\n",
    "     = \\mathbb{E}_q \\log p(X|\\Theta) \\\\\n",
    "   %% &= \\Bigl[p(X,Z|\\Theta) = p(Z|X,\\Theta) p(X|\\Theta) \\Bigr] \\\\\n",
    "    &= \\mathbb{E}_{Z\\sim q} \\log \\frac{p(X,Z|\\Theta)}{p(Z|X\\Theta)}\n",
    "     = \\mathbb{E}_{Z\\sim q} \\log \\frac{q(Z)}{p(Z|X,\\Theta)}\n",
    "     + \\mathbb{E}_{Z\\sim q} \\log \\frac{p(X,Z|\\Theta)}{q(Z)} \\\\ \n",
    "    &= KL\\bigl(q\\|p(\\cdot|X,\\Theta)\\bigr) + \\mathcal{L}\\bigl(q, \\Theta\\bigr)\n",
    "\\end{align*}  \n",
    "\n",
    "since the Bayes theorem posits that $p(X,Z|\\Theta) = p(Z|X,\\Theta) p(X|\\Theta)$. Call this equiation the **\"master equation\"**.\n",
    "\n",
    "Now note that since the Kullback-Leibler divergence is always non-negative, one has the following inequality:\n",
    "$$\\log p(X|\\Theta) \\geq \\mathcal{L}\\bigl(q, \\Theta\\bigr)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try to make the lower bound as large as possible by changing $\\Theta$ and varying $q$. But first note that the \n",
    "left-hand side of the **master equation** is independent of $q$, whence maximization of $\\mathcal{L}$ with respect to $q$ (with $\\Theta$ fixed) is equivalent to minimization of $KL\\bigl(q\\|p(\\cdot|X,\\Theta)\\bigr)$ with respect to $q$ taking $\\Theta$ fixed. Since $q$ is arbitrary, the optimal minimizer $q^*_\\Theta$ is $q^*(Z|\\Theta) = p(Z|X,\\Theta)$ for all $Z$.\n",
    "\n",
    "Now at the optimal distributuion $q^*_\\Theta$ the **master equation** becomes\n",
    "$$ \\log p(X|\\Theta)\n",
    "= \\mathcal{L}\\bigl(q^*_\\Theta, \\Theta\\bigr)\n",
    "= \\mathbb{E}_{Z\\sim q^*_\\Theta} \\log \\frac{p(X,Z|\\Theta)}{q^*(Z|\\Theta)}\n",
    "= \\mathbb{E}_{Z\\sim q^*_\\Theta} \\log p(X,Z|\\Theta) - \\mathbb{E}_{Z\\sim q^*_\\Theta} \\log q^*(Z|\\Theta)\n",
    "$$\n",
    "for any $\\Theta$. Thus the problem of log-likelihood maximization reduces to that of maximizing the sum of expectations on the right-hand side.\n",
    "\n",
    "This new problem does not seem to be tractable in general since the optimization paramters $\\Theta$ affect both the expected log-likelihood $\\log p(X,Z|\\Theta)$ under $Z\\sim q^*_\\Theta$ and the entropy of the optimal distribution of the latent variables $Z$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully using an iterative procedure which switches between the computation of $q^*_\\Theta$ and the maximization of $\\Theta$ might be effective. Consider the folowing :\n",
    "* **E**-step: considering $\\Theta_i$ as given and fixed find $q^*_{\\Theta_i} = \\mathop{\\text{argmin}}_q KL\\bigl(q\\|p(\\cdot|X,\\Theta_i)\\bigr)$ and set $q_{i+1} = q^*_{\\Theta_i}$;\n",
    "* **M**-step: considering $q_{i+1}$ as given, solve $\\mathcal{L}(q_{i+1},\\Theta) \\to \\mathop{\\text{max}}_\\Theta$, where \n",
    "$$ \\mathcal{L}(q,\\Theta) = \\mathbb{E}_{Z\\sim q} \\log p(X,Z|\\Theta) - \\mathbb{E}_{Z\\sim q} \\log q(Z) $$\n",
    "\n",
    "The fact that $q_i$ is considered fixed makes the optimization of $\\mathcal{L}(q_i,\\Theta)$ equivalent to maximization of the expected log-likelihood, since the entropy term is fixed. Therefore the **M**-step becomes:\n",
    "* given $q_{i+1}$ find $\\Theta^*_{i+1} = \\mathop{\\text{argmax}}_\\Theta \\mathbb{E}_{Z\\sim q_{i+1}} \\log p(X,Z|\\Theta)$ and put $\\Theta_{i+1} = \\Theta^*_{i+1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, if the latent variables are mutually independent, then the optimal $q$ must be factorizable into marginal densities and:\n",
    "\n",
    "\\begin{align*}\n",
    "    KL\\bigl(q\\|p(\\cdot|X,\\Theta)\\bigr)\n",
    "    &= \\mathbb{E}_{Z\\sim q} \\log q(Z) - \\sum_j \\mathbb{E}_{z_j\\sim q_j} \\log p(z_j|X,\\Theta)\\\\\n",
    "    &= \\sum_j \\mathbb{E}_{z_j\\sim q_j} \\log q_j(z_j) - \\sum_j \\mathbb{E}_{z_j\\sim q_j} \\log p(z_j|X,\\Theta)\n",
    "     = \\sum_j KL\\bigl(q_j\\|p_j(|X,\\Theta)\\bigr)\n",
    "\\end{align*}\n",
    "where $q_j$ is the marginal desity of $z_j$ in $q(Z)$ (the last term in the first line comes from the Fubini theorem).\n",
    "\n",
    "Therefore the **E**-step could be reduced to a set of minimization problems with respect to one-dimensional density functions:\n",
    "$$ q_j^* = \\mathop{\\text{argmin}}_{q_j} KL\\bigl(q_j\\|p_j(\\cdot|X,\\Theta)\\bigr) $$\n",
    "since the Kulback-Leibler divergence in this case in additively separable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correctness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that the **master equation** is an identity: for all densities $q$ on $Z$ and for all admissible parameters $\\Theta$\n",
    "$$ \\log p(X|\\Theta) = KL\\bigl(q\\|p(\\cdot|X,\\Theta)\\bigr) + \\mathcal{L}\\bigl(q, \\Theta\\bigr) $$\n",
    "\n",
    "Hence if after the **E**-step the Kulback-Leibler divergence is reduced:\n",
    "$$ KL\\bigl(q'\\|p(\\cdot|X,\\Theta)\\bigr) \\leq KL\\bigl(q\\|p(\\cdot|X,\\Theta)\\bigr) $$\n",
    "then for the same set of parameters $\\Theta$ one has\n",
    "$$ \\mathcal{L}(q,\\Theta) \\leq \\mathcal{L}(q',\\Theta) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just after the **E**-step one has $q_{i+1} = p(Z|X,\\Theta_i)$, whence $KL\\bigl(q_{i+1}\\|p(\\cdot|X,\\Theta_i)\\bigr) = 0$. In turn, this implies via the **master equation** that the following equality holds:\n",
    "$$ \\log p(X|\\Theta_i) = \\mathcal{L}(q_{i+1},\\Theta_i) $$\n",
    "\n",
    "After the **M**-step, since $\\Theta_{i+1}$ is a maximizer, or at least an \"_improver_\" of $\\mathcal{L}(q_{i+1},\\Theta)$ compared to its value at $(q_i,\\Theta_i)$, one has\n",
    "$$ \\mathcal{L}(q_{i+1},\\Theta_i) \\leq \\mathcal{L}(q_{i+1},\\Theta_{i+1}) $$\n",
    "\n",
    "Threfore the effect of a single complete round of **EM** on the log-likelihood itself is:\n",
    "$$ \\log p(X|\\Theta_i) = \\mathcal{L}(q_{i+1},\\Theta_i) \\leq \\mathcal{L}(q_{i+1},\\Theta_{i+1}) \\leq \\mathcal{L}(q_{i+2},\\Theta_{i+1}) = \\log p(X|\\Theta_{i+1}) $$\n",
    "where the equality is achieved between the **E** and the **M** step within one round. This implies that **EM** indeed iteratively improves the log-likihood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note, that in the general case, without attaining zero Kulback-Leibler divergence at the $E$-step, one cannot be sure that the real log-likelihood is improved by each iteration and one can just say that\n",
    "$$\\mathcal{L}(q_{i+1},\\Theta_i) \\leq \\log p(X|\\Theta_i)$$\n",
    "which does not uncover a relationship with $\\log p(X|\\Theta_{i+1})$. And without the guarantee that **EM** improves the log-likelihood to the maximum one cannot be sure about the consistency of the estimators. The key question is whether the lower bound $\\mathcal{L}(q,\\Theta)$ is any good."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application of the EM to NIST data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each image is a random element in a discrete probability space $\\Omega = \\{0,1\\}^{N\\times M}$ with product-measure\n",
    "$$\\mathbb{P}(\\omega) = \\prod_{i=1}^N\\prod_{j=1}^M \\theta_{ij}^{\\omega_{ij}} (1-\\theta_{ij})^{1-\\omega_{ij}}$$\n",
    "for any $\\omega\\in \\Omega$. In particular $M=N=28$. Basically each bit of the image is independent of any other bit and each one is a Bernoulli random variable with parameter $\\theta_{ij}$: $\\omega_{ij}\\sim \\text{Bern}(\\theta_{ij})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's apply the EM algorithm to this dataset. The proposed model is the following.\n",
    "\n",
    "Consider a mixture model of discrete probability spaces. Suppose there are $K$ componets in the mixture. Then each image is distributed according to the following law:\n",
    "$$p(\\omega|\\Theta)\n",
    "= \\sum_{k=1}^K \\pi_k p_k(\\omega|\\theta_k)\n",
    "= \\sum_{k=1}^K \\pi_k \\prod_{i=1}^N \\prod_{j=1}^M \\theta_{kij}^{\\omega_{ij}} (1-\\theta_{kij})^{1-\\omega_{ij}}$$\n",
    "where $\\theta_{kij}$ is the paramter of the probability distribution of the $(i,j)$-th random variable (pixel) in the $k$-th class, and $\\pi_k$ is the (prior) porbability of the $k$-th mixutre to generate a random element, $\\sum_{k=1}^K \\pi_k= 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose $X=(x_i)_{i=1}^n \\in \\Omega^n$ is the dataset. The log-likelihood is given by\n",
    "$$ \\log p(X|\\Theta) = \\sum_{s=1}^n \\log \\sum_{k=1}^K \\pi_k\n",
    "\\prod_{i=1}^N \\prod_{j=1}^M \\theta_{kij}^{x_{sij}} (1-\\theta_{kij})^{1-x_{sij}}$$\n",
    "where $x_{sij}\\in\\{0,1\\}$ -- is the value of the the $(i,j)$-th pixel at the $s$-th observation.\n",
    "\n",
    "If the source $Z=(z_i)_{i=1}^n$ components of the mixture at each datapoint were known, then the log-likelihood would have been\n",
    "$$ \\log p(X,Z|\\Theta) = \\sum_{s=1}^n \\log \\prod_{k=1}^K \\Bigl[ \\pi_k \n",
    "\\prod_{i=1}^N \\prod_{j=1}^M \\theta_{kij}^{x_{sij}} (1-\\theta_{kij})^{1-x_{sij}} \\Bigr]^{1_{z_s = k}}$$\n",
    "where $1_{z_s = k}$ is the indicator and take the value $1$ if $\\{z_s = k\\}$ and $0$ otherwise ($1_{\\{k\\}}(z_s)$ is another notation)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The log-likelihood simplifies to\n",
    "$$ \\log p(X,Z|\\Theta) = \\sum_{s=1}^n \\sum_{k=1}^K 1_{z_s = k} \\Bigl( \\log \\pi_k + \n",
    "\\sum_{i=1}^N \\sum_{j=1}^M \\bigl( x_{sij} \\log \\theta_{kij} + (1-x_{sij}) \\log (1-\\theta_{kij}) \\bigr) \\Bigr) $$\n",
    "and further into a more separable form\n",
    "$$ \\log p(X,Z|\\Theta)\n",
    "= \\sum_{s=1}^n \\sum_{k=1}^K 1_{z_s = k} \\log \\pi_k\n",
    "+ \\sum_{s=1}^n \\sum_{k=1}^K 1_{z_s = k} \\Bigl( \\sum_{i=1}^N \\sum_{j=1}^M x_{sij} \\log \\theta_{kij}\n",
    "+ \\sum_{i=1}^N \\sum_{j=1}^M (1-x_{sij}) \\log (1-\\theta_{kij}) \\Bigr)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The expected log-likelihood under $z_s\\sim q_s$ with $\\mathbb{P}(z_s=k|X) = q_{sk}$, is given by\n",
    "$$ \\mathbb{E}\\log p(X,Z|\\Theta)\n",
    "= \\sum_{s=1}^n \\sum_{k=1}^K q_{sk} \\log \\pi_k\n",
    "+ \\sum_{s=1}^n \\sum_{k=1}^K q_{sk} \\sum_{i=1}^N \\sum_{j=1}^M \\bigl( x_{sij} \\log \\theta_{kij} + (1-x_{sij}) \\log (1-\\theta_{kij}) \\bigr)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Analytic solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the **E**-step one must compute $q^*(Z) = \\mathbb{P}(z_s=k|X) = \\hat{q}_{sk}$ based on the value of $\\Theta = ((\\pi_k), (\\theta_{kij}))$.\n",
    "$$\\hat{q}_{sk}\n",
    "= \\frac{p(x_s|z_s=k,\\Theta) p(z_s=k)}{\\sum_{l=1}^K p(x_s|z_s=l,\\Theta) p(z_s=l)}\n",
    "\\propto \\pi_k \\prod_{i=1}^N \\prod_{j=1}^M \\theta_{kij}^{x_{sij}} (1-\\theta_{kij})^{1-x_{sij}}\n",
    "$$\n",
    "and\n",
    "$$ q^*(Z) = \\prod_{s=1}^n q_{s z_s} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to improve numerical stability and avoid numerical underflow it is better to use the following procedure for computation of the conditional probability:\n",
    "$$ l_{sk} = \\log \\pi_k + \\sum_{i=1}^N \\sum_{j=1}^M \\log \\bigl( \\theta_{kij}^{x_{sij}} (1-\\theta_{kij})^{1-x_{sij}} \\bigr)$$\n",
    "set $l^*_s = \\max_k l_{sk}$ and compute the log-sum\n",
    "$$ \\hat{l}_s = l^*_s + \\log \\sum_{k=1}^K \\text{exp}\\Bigl\\{ l_{sk} - l^*_s \\Bigr\\} $$\n",
    "This seemingly redunant subctration and addition of $l^*_s$ helps avoid underflow during the numerical exponentiation.\n",
    "However if the underflow still took place in the sum-exp, then those values of $\\hat{l}_s$ that evaluated to floating point _NAN_ can reliably be replaced with respective values of $l^*_s$.  After this sanitization the **E**-step's optimal distribution is calculated as:\n",
    "$$ \\hat{q}_{sk} = \\text{exp}\\Bigl\\{ l_{sk} - \\hat{l}_s \\Bigr\\} $$\n",
    "which would yield a numerically accurate distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $l^*_s >> l_{sk}$ for all $k$ but such that $l^*_s = l_{sk}$ (let it be $k^*$), then an underflow occurs at the sum-exp step, whence for some very small $\\epsilon > 0$ one has\n",
    "$$\\hat{l}_s = l^*_s + \\log (1+\\epsilon) $$\n",
    "whence\n",
    "$$\\hat{q}_{sk} = \\text{exp}\\bigl\\{l_{sk} - \\hat{l}_s\\bigr\\} = (1+\\epsilon)^{-1} \\cdot \\text{exp}\\bigl\\{l_{sk} - l^*_s \\bigr\\}$$\n",
    "For $k=k^*$ on has $\\hat{q}_{sk} = \\frac{1}{1+\\epsilon}\\approx 1$, and for $k\\neq k^*$ -- $\\hat{q}_{sk} = \\frac{\\eta}{1+\\epsilon} \\approx 0$ for some extremely small $\\eta>0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variables in the code have the following dimensions:\n",
    "* $\\theta \\in [0,1]^{K\\times (N\\times M)}$;\n",
    "* $\\pi \\in [0,1]^{1\\times K}$;\n",
    "* $x \\in \\{0,1\\}^{n\\times (N\\times M)}$;\n",
    "* $z \\in [0,1]^{n\\times K}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## A bunch of wrappers to match the task specifications\n",
    "def posterior( x, clusters ) :\n",
    "    pi = np.ones( clusters.shape[ 0 ], dtype = np.float ) / clusters.shape[ 0 ]\n",
    "    q, ll = __posterior( x, theta = clusters, pi = pi )\n",
    "    return q\n",
    "\n",
    "## The likelihood is a byproduct of the E-step's minimization of Kullback-Leibler\n",
    "def likelihood( x, clusters ) :\n",
    "    pi = np.ones( clusters.shape[ 0 ], dtype = np.float ) / clusters.shape[ 0 ]\n",
    "    q, ll = __posterior( x, theta = clusters, pi = None )\n",
    "    return ll\n",
    "\n",
    "## The core procedure for computing the conditional density of classes\n",
    "def __posterior( x, theta, pi ) :\n",
    "## Unfortunately sometimes there can be negative machine zeros, which\n",
    "##  spoil the log-likelihood computation by poisoning with NANs.\n",
    "## That is why the theta array is restricted to [0,1].\n",
    "    theta_clipped = np.clip( theta, 0.0, 1.0 )\n",
    "## Note that the power formulation is just a mathematically convenient\n",
    "##  way of writing \\theta if x=1 or (1-\\theta) otherwise. This comapct\n",
    "##  way of computing had to be abandoned, since it required huge amounts\n",
    "##  of memory.\n",
    "#     ll_skij = np.where( x[:, np.newaxis] > 0,\n",
    "#         np.log( theta_clipped ), np.log( 1 - theta_clipped ) )\n",
    "#     ll_sk = np.sum( ll_skij, axis = ( 2, ) ) + np.log( pi )\n",
    "#     del ll_skij\n",
    "## Iterate over classes\n",
    "    ll_sk = np.zeros( ( x.shape[ 0 ], theta.shape[ 0 ] ), dtype = np.float )\n",
    "    for k in xrange( theta.shape[ 0 ] ) :\n",
    "        ll_sij = np.where( x > 0, np.log( theta_clipped[ k ] ), np.log( 1 - theta_clipped[ k ] ) )\n",
    "        ll_sk[ :, k ] = np.sum( ll_sij, axis = ( 1, ) ) + np.log( pi[ k ] )\n",
    "## Find the largest unnormalized probability.\n",
    "    llstar_s = np.reshape( np.max( ll_sk, axis = ( 1, ) ), ( ll_sk.shape[ 0 ], 1 ) )\n",
    "## Compute the log-sum-exp of the individual log-likelihoods\n",
    "    ll_s = np.reshape( np.log( np.sum( np.exp( ll_sk - llstar_s ), axis = ( 1, ) ) ), ( ll_sk.shape[ 0 ], 1 ) )\n",
    "## The sum-exp could never be anything lower than 1, since at least one\n",
    "##  element of each row of ll_sk has to be lstar_s, whence the respective\n",
    "##  difference should be zero and the exponent -- 1. Thus even if the\n",
    "##  rest of the sum is close to machine zero, the logarithm would still\n",
    "##  return 0.\n",
    "    ll_s += llstar_s\n",
    "    del llstar_s\n",
    "## Normalise the likelihoods to get conditional probability, and compute\n",
    "##  the sum of the log-denominator, which is the log-likelihood.\n",
    "    return np.exp( ll_sk - ll_s ), np.sum( ll_s )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the **M**-step for some fixed $q(Z)$ one solves $\\mathbb{E}\\log p(X,Z|\\Theta)\\to \\max_\\Theta$ subject to $\\sum_{k=1}^K \\pi_k = 1$ which is a convex optimization problem with respect to $\\Theta$, since the log-likelihood as a linear combination of convex functions is convex. The first order condition is $\\sum_{s=1}^n \\frac{q_{sk}}{\\pi_k} - \\lambda = 0$ for all $k=1,\\ldots,K$, whence $ \\lambda = \\sum_{s=1}^n \\sum_{l=1}^K q_{sl} = n $ and finally\n",
    "$$\\hat{\\pi}_k = \\frac{\\sum_{s=1}^n q_{sk}}{n}$$\n",
    "For $\\theta_{kij}$, $i=1,\\ldots,N$, $j=1,\\ldots,M$ and $k=1,\\ldots,K$ the FOC is\n",
    "$$ \\sum_{s=1}^n q_{sk} \\frac{x_{sij}}{\\theta_{kij}} - \\sum_{s=1}^n q_{sk} \\frac{1-x_{sij}}{1-\\theta_{kij}} = 0 $$\n",
    "whence\n",
    "$$\\hat{\\theta}_{kij} =  \\frac{\\sum_{s=1}^n q_{sk} x_{sij}}{ \\sum_{s=1}^n q_{sk} } = \\frac{\\sum_{s=1}^n q_{sk} x_{sij}}{ n \\hat{\\pi}_k }$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## A wrapper to match the task specifications\n",
    "def learn_clusters( x, z ) :\n",
    "    theta, pi = __learn_clusters( x, z )\n",
    "    return theta\n",
    "\n",
    "## The E-step is simple: just compute the optimal parameters under\n",
    "##  the current conditional distribution of the latend variables.\n",
    "def __learn_clusters( x, z ) :\n",
    "## The prior class probabilities\n",
    "    pi = np.sum( z, axis = ( 0, ) )\n",
    "## Pixel probabilities conditional on the calss\n",
    "    theta = np.dot( z.T, x ) / np.reshape( pi, ( z.shape[ 1 ], 1 ) )\n",
    "## Return\n",
    "    return theta, np.reshape( pi, ( 1, pi.shape[ 0 ] ) ) / np.sum( pi )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A it has been mentioned eariler, the EM algorithm switches between **E** and **M** steps until convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## A wrapper for the core em algorithm above\n",
    "def em_algorithm( x, K, maxiter, verbose = True, rel_eps = 1e-4 ) :\n",
    "## Initialize the model parameters with uniform [0,1] random numbers\n",
    "    theta = rand.uniform( size = ( K, x.shape[ 1 ] ) )\n",
    "## Run the em algorithm\n",
    "    ll, theta, pi, status = __em_algorithm( x, theta_1 = theta,\n",
    "        pi_1 = None, niter = maxiter, rel_eps = rel_eps, verbose = verbose )\n",
    "## Return the history of theta and the final log liklihood\n",
    "    if verbose :\n",
    "        if status[ 'status' ] != 0 :\n",
    "            print \"Convergence not achieved. %d\" % ( status[ 'status' ], )\n",
    "    return theta, ll"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The procedure above actually invokes the true EM core, defined below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## The core of the EM algorithm\n",
    "def __em_algorithm( x, theta_1, pi_1 = None, niter = 1000, rel_eps = 1e-4, verbose = True ) :\n",
    "## If we were supplied with an initial estimate of the prior distribution,\n",
    "##  then assume the full model is needed.\n",
    "    full_model = pi_1 is not None\n",
    "## If the prior cluster probabilities are not supplied, assume uniform distribution.\n",
    "    pi_1 = pi_1 if pi_1 is not None else np.ones( theta_1.shape[ 0 ], dtype = np.float ) / theta_1.shape[ 0 ]\n",
    "## Allocate the necessary space for the history of model estimates\n",
    "    theta_hist = np.empty( ( 0, ) + theta_1.shape, dtype = np.float )\n",
    "    pi_hist = np.empty( ( 0, theta_1.shape[ 0 ] ), dtype = np.float )\n",
    "    ll_hist = np.empty( 0, np.float )\n",
    "## Set \"old\" estimated to zero. At this line the current estimates are in fact\n",
    "##  the initially provided ones.\n",
    "    theta_0, pi_0 = np.zeros_like( theta_1 ), np.zeros_like( pi_1 )\n",
    "## Initialize the loop\n",
    "    status, kiter, rel_theta, rel_pi, ll = -1, 0, np.nan, np.nan, -np.inf\n",
    "    while kiter < niter :\n",
    "## Dump the current estimators and other information.\n",
    "        if verbose :\n",
    "            print( \"Iteration %d: avg. log-lik: %.3f, $\\\\Theta$ div. %.3f, $\\\\Pi$ div. %.3f\" % (\n",
    "                kiter, ll / x.shape[ 0 ], rel_theta, rel_pi ) )\n",
    "            show_data( theta_1 - theta_0 if True else theta_1, n = theta_0.shape[ 0 ],\n",
    "                n_col = min( 10, theta_0.shape[ 0 ] ), cmap = plt.cm.hot, interpolation = 'nearest' )\n",
    "## The convergence criterion is the L^∞ norm of relative L^1 errors\n",
    "        if max( rel_pi, rel_theta ) < rel_eps :\n",
    "            status = 0\n",
    "            break ;\n",
    "## Overwrite the initial estimates\n",
    "        theta_0, pi_0 = theta_1, pi_1\n",
    "## E-step: call the core posterior function to get both the log-likelihood\n",
    "##  and the estimate of the conditional distribution.\n",
    "        z_1, ll = __posterior( x, theta_0, pi_0 )\n",
    "## M-step: compute the optimal parameters under the current estimate of the posterior\n",
    "        theta_1, pi_1 = __learn_clusters( x, z_1 )\n",
    "## Discard the computed estimate of pi if the model is discriminative (conditional likelihood).\n",
    "        if not full_model :\n",
    "            pi_1 = pi_0\n",
    "## Record the current estimates to the history\n",
    "        theta_hist = np.vstack( ( theta_hist, theta_1[np.newaxis] ) )\n",
    "        pi_hist = np.vstack( ( pi_hist, pi_1 ) )\n",
    "        ll_hist = np.append( ll_hist, ll )\n",
    "## Check for bad float numbers\n",
    "        if not ( np.all( np.isfinite( theta_1 ) ) and np.all( np.isfinite( pi_1 ) ) ) :\n",
    "            status= -2\n",
    "            break ;\n",
    "## Check convergence: L^1 relative error. If the relative margin is exactly\n",
    "##  zero, then return NaNs. This makes teh loop exhaust all iterations, since\n",
    "##  any comprison agains a NaN returns False.\n",
    "        rel_theta = np.sum( np.abs( theta_1 - theta_0 ) / ( np.abs( theta_0 ) + rel_eps ) ) if rel_eps > 0 else np.nan\n",
    "        rel_pi = np.sum( np.abs( pi_1 - pi_0 ) / ( np.abs( pi_0 ) + rel_eps ) ) if rel_eps > 0 else np.nan\n",
    "## Next iteration\n",
    "        kiter += 1\n",
    "    return ll_hist, theta_hist, pi_hist, { 'status': status, 'iter': kiter }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Miscellanea"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to be able to plot more flexibly, define another arranger."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## A more flexible image arrangement\n",
    "def arrange_flex( images, n_row = 10, n_col = 10, N = 28, M = 28, fill_value = 0 ) :\n",
    "## Create the final grid of images row-by-row\n",
    "    im_grid = np.full( ( n_row * N, n_col * M ), fill_value, dtype = images.dtype )\n",
    "    for k in range( min( images.shape[ 0 ], n_col * n_row ) ) :\n",
    "## Get the grid cell at which to place the image\n",
    "        i, j = ( k // n_col ) * N, ( k % n_col ) * M\n",
    "## Just put the image in the cell\n",
    "        im_grid[ i:i+N, j:j+M ] = np.reshape( images[ k ], ( N, M, ) )\n",
    "    return im_grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The procedure below loads the MNIST data from a comma-separated text file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def load_mnist_from_csv( filename ) :\n",
    "## Read the CSV file\n",
    "    data = np.loadtxt( open( filename, \"rb\" ), dtype = np.short, delimiter = \",\", skiprows = 0 )\n",
    "## Peel off the lables\n",
    "    return data[:,1:], data[:,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A procedure for pickling numpy arrays of floats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import base64, zlib\n",
    "\n",
    "def np_pack( data ) : return base64.b64encode( zlib.compress( data ) )\n",
    "def np_unpack( text ) : return np.frombuffer( zlib.decompress(\n",
    "    base64.decodestring( text ) ), dtype = np.float )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A handy procedure for one-line visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def show_data( data, n, n_col = 10, **kwargs ) :\n",
    "## Get the number of rows necessary to plot the needed number of images\n",
    "    n_row = ( n + n_col - 1 ) // n_col\n",
    "## Set the dimensions of the figure\n",
    "    plt.figure( figsize = ( n_col, n_row ) )\n",
    "## Plot!\n",
    "    plt.imshow( arrange_flex( data[:n], n_col = n_col, n_row = n_row ), **kwargs )\n",
    "    plt.show( )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a convenient procedure for running experiments. By setting relative error to zero the algorithm is forced to exhaust all the allocated iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def experiment( data, K, maxiter, verbose = True, until_convergence = False ) :\n",
    "## Run the EM\n",
    "    return em_algorithm( data, K, maxiter, rel_eps = 1.0e-4 if until_convergence else 0.0, verbose = verbose )\n",
    "\n",
    "def visualize( data, clusters, ll, n_col = 2, plot_ll = True ) :\n",
    "## Display the result\n",
    "    print \"Final conditional log-likelihood value per observation achieved %f in %d iteration(s)\" % ( ll[-1] / data.shape[ 0 ], len( ll ) )\n",
    "## Plot the first difference of average log-likelihood\n",
    "    if plot_ll :\n",
    "        plt.figure( figsize = ( 16, 9 ) )\n",
    "        ax = plt.subplot(111)\n",
    "        ax.set_title( r\"avg. log-likelihood change between successive iterations (log scale)\" )\n",
    "        ax.plot( np.diff( ll / data.shape[ 0 ]  ) )\n",
    "        ax.set_ylabel( r\"$\\Delta_i \\frac{1}{n} \\sum_{s=1}^n \\mathbb{E}_{z_s\\sim q_i} \\log p(x_s,z_s|\\Theta_i)$\" )\n",
    "        ax.set_yscale( 'log' )\n",
    "## Plot the final estimates\n",
    "    if n_col > 0 :\n",
    "        show_data( clusters[-1], n = clusters.shape[1], n_col = n_col, cmap = plt.cm.gray, interpolation = 'nearest' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function that produces and embed a video in HTML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Make simple animations of the EM estimatora\n",
    "## http://jakevdp.github.io/blog/2013/05/12/embedding-matplotlib-animations/\n",
    "## http://jakevdp.github.io/blog/2012/08/18/matplotlib-animation-tutorial/\n",
    "from matplotlib import animation\n",
    "from IPython.display import HTML\n",
    "from tempfile import NamedTemporaryFile\n",
    "VIDEO_TAG = \"\"\"<video controls><source src=\"data:video/x-m4v;base64,{0}\" type=\"video/mp4\">Your browser does not support the video tag.</video>\"\"\"\n",
    "def display_animation( anim ) :\n",
    "    plt.close( anim._fig )\n",
    "    if not hasattr( anim, '_encoded_video' ) :\n",
    "        ffmpeg_writer = animation.FFMpegWriter( )\n",
    "        with NamedTemporaryFile( suffix = '.mp4' ) as f:\n",
    "            anim.save( f.name, fps = 15, extra_args = [ '-vcodec', 'libx264' ] )# , writer = ffmpeg_writer )\n",
    "            video = open( f.name, \"rb\" ).read( )\n",
    "        anim._encoded_video = video.encode( \"base64\" )\n",
    "    return HTML( VIDEO_TAG.format( anim._encoded_video ) )\n",
    "\n",
    "def animate( theta, ll, n_col = 10, n_row = 10, diff = True, interval = 20, **kwargs ) :\n",
    "## First set up the figure, the axis, and the plot element we want to animate\n",
    "    fig = plt.figure( )\n",
    "    ax1, ax2 = fig.add_subplot( 211 ), fig.add_subplot( 212 )\n",
    "    bg = arrange_flex( np.zeros_like( theta[ 0 ] ), n_col = n_col, n_row = n_row )\n",
    "    if diff :\n",
    "        im = ax1.imshow( bg, vmin = -1.0, vmax = +1.0, **kwargs )\n",
    "    else :\n",
    "        im = ax1.imshow( bg, vmin = +0.0, vmax = +1.0, **kwargs )\n",
    "    line, = ax2.plot( [], linestyle = \"-\", color = 'blue' )\n",
    "    ax2.set_xlim( -0.1, theta.shape[ 0 ] + 0.1 )\n",
    "    ax2.set_yscale( 'log' )\n",
    "    ll_diff = np.diff( ll ) / data.shape[ 0 ]\n",
    "    ax2.set_ylim( np.min( ll_diff ) * 0.9, np.max( ll_diff ) * 1.1 )\n",
    "## Animation function.  This is called sequentially\n",
    "    def update( i ) :\n",
    "        if diff :\n",
    "            frame = theta[ i ] - theta[ i-1 ] if i > 0 else theta[ 0 ]\n",
    "            frame /= np.max( np.abs( frame ) )\n",
    "        else :\n",
    "            frame = theta[ i ]\n",
    "        im.set_data( arrange_flex( frame, n_col = n_col, n_row = n_row ) )\n",
    "        if i > 0 :\n",
    "            line.set_data( np.arange( i ), ll_diff[ :i ] )\n",
    "        return im, line,\n",
    "## call the animator.  blit=True means only re-draw the parts that have changed.\n",
    "    return animation.FuncAnimation( fig, update, frames = theta.shape[ 0 ], interval = interval, blit=True )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading and visualizing MNIST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Preamble"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obtain the MNIST data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if False :\n",
    "## Fetch MNIST dataset from SciKit and create a local copy.\n",
    "    from sklearn.datasets import fetch_mldata\n",
    "    mnist = fetch_mldata( \"MNIST original\", data_home = './data/' )\n",
    "    np.savez_compressed('./data/mnist/mnist_scikit.npz', data = mnist.data, labels = mnist.target )\n",
    "\n",
    "if False :\n",
    "## Fetch the data from the provided CSV (!) files and save as a compressed data blob\n",
    "    data, labels = load_mnist_from_csv( \"./data/mnist/mnist_train.csv\" )\n",
    "    np.savez_compressed( './data/mnist/mnist_train.npz', labels = labels, data = data )\n",
    "    data, labels = load_mnist_from_csv( \"./data/mnist/mnist_test.csv\" )\n",
    "    np.savez_compressed( './data/mnist/mnist_test.npz', labels = labels, data = data )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Study"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First of all load and binarize the training data using the values 127 as the threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "assert( os.path.exists( './data/mnist/mnist_train.npz' ) )\n",
    "with np.load( './data/mnist/mnist_train.npz', 'r' ) as npz :\n",
    "    mnist_labels, mnist_data = npz[ 'labels' ], np.array( npz[ 'data' ] > 127, np.int )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Case : $K=2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at some 6s and 9s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Mask\n",
    "inx_sixes, inx_nines = np.where( mnist_labels == 6 )[ 0 ], np.where( mnist_labels == 9 )[ 0 ]\n",
    "## Extract\n",
    "sixes = mnist_data[ rand.choice( inx_sixes, 90, replace = False ) ]\n",
    "nines = mnist_data[ rand.choice( inx_nines, 90, replace = False ) ]\n",
    "## Show\n",
    "show_data( sixes, n = 45, n_col = 15, cmap = plt.cm.gray, interpolation = 'nearest' )\n",
    "show_data( nines, n = 45, n_col = 15, cmap = plt.cm.gray, interpolation = 'nearest' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They do indeed look quite distinct. Now collect them into a single dataset and estimate the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data = mnist_data[ np.append( inx_sixes, inx_nines ) ]\n",
    "clusters, ll = experiment( data, 2, 30 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimate deltas show that the **EM** algorithm's **E**-step actually transfers the unlikely observations between classes, as is expected by constructon of the algorithm.\n",
    "\n",
    "Judging by the plot below, it turns out that 30 iterations is more than enough for the **EM** to get meaninful estimates the class ideals, represented by the probability porduct-measure on $\\Omega^{28\\times 28}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "visualize( data, clusters, ll )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see how well the **EM** algorithm performs on a model with more classes. But before that let's have a look at a random sample of the handwritten digits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "indices = np.arange( mnist_data.shape[ 0 ] )\n",
    "rand.shuffle( indices )\n",
    "plt.imshow( arrange_flex( mnist_data[ indices[:100] ] ), cmap = plt.cm.gray, interpolation = 'nearest' )\n",
    "plt.show( )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Case : $K=10, 15, 20$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The original size of the trainig sample is too large to fit in these RAM banks :) That is why I had to limit the sample to a random subset of 2000 observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = mnist_data[ indices[:2000]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the procedure that perfoems **EM** algorithm and return the history of the parameter estimates as well as the dynamics of the log-likelihood lower bonud."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "clusters, ll = experiment( data, 10, 50 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can clearly see, that $50$ iterations were not enough for the alogirithm to converge: though the changes are tiny, even on the log-scale, they are still unstable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "visualize( data, clusters, ll, n_col = 5 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Furthermore the returned templates of the mixture components and clearly suboptimal: the procedure seems to get stuck at individual examples. One way to remedy this is to use more iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "clusters_10, ll_10 = experiment( data, 10, 1000, verbose = True, until_convergence = True )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As one can see, increasing the number of iterations does not necessarily improve the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "visualize( data, clusters_10, ll_10, n_col = 5 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see if changing $K$ does the trick."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "clusters_15, ll_15 = experiment( data, 15, 1000, verbose = True, until_convergence = True )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "clusters_20, ll_20 = experiment( data, 20, 1000, verbose = True, until_convergence = True )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "visualize( data, clusters_10, ll_10, n_col = 10, plot_ll = False )\n",
    "visualize( data, clusters_15, ll_15, n_col = 10, plot_ll = False )\n",
    "visualize( data, clusters_20, ll_20, n_col = 10, plot_ll = False )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "display_animation( animate( clusters_20, ll_20, 10, 2, diff = False, cmap = plt.cm.gray, interpolation = 'nearest' ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "</hr>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Record good runs for the posterity\n",
    "clusters_b64_clob = \"\"\"<KEY>\"\"\"\n",
    "if not 'clusters_b64_clob' in locals( ) :\n",
    "    clusters_b64_clob = np_pack( clusters[-1:] )\n",
    "    print clusters_b64_clob\n",
    "elif False :\n",
    "    clusters = np_unpack( clusters_b64_clob ).reshape( ( 2, -1 ) )\n",
    "show_data( np_unpack( clusters_b64_clob ).reshape( ( 2, -1 ) ), 2, 2 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_train_samples = 1000\n",
    "indices = np.arange( X.shape[ 0 ] )\n",
    "np.random.shuffle( indices )\n",
    "\n",
    "train_indices = indices[ : n_train_samples ]\n",
    "X_train = X[ train_indices, : ]\n",
    "y_train = y[ train_indices ]\n",
    "\n",
    "# del X, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X, y = np.array( X_train > 0.5, np.int ), y_train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def theta_show( Theta, n_col = 10, **kwargs ) :\n",
    "    plt.figure( figsize = ( n_col, ( nclasses + n_col-1 ) // n_col ) )\n",
    "    plt.imshow( arrange_flex( Theta, n_col = n_col,\n",
    "        n_row = ( nclasses + n_col-1 ) // n_col ), **kwargs )\n",
    "    plt.show( )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nclasses = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "uPi = - np.log( np.random.uniform( size = (1, nclasses) ) )\n",
    "Pi = uPi / np.sum( uPi )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Theta = np.random.uniform( size = ( nclasses, X.shape[ 1 ] ) )\n",
    "Theta_init, Pi_init = Theta.copy(), Pi.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "theta_show( Theta, 5, cmap = plt.cm.gray, interpolation = 'nearest' )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Theta, Pi = Theta_init.copy( ), Pi_init.copy( )\n",
    "kiter, niter = 0, 1000\n",
    "while kiter < niter :\n",
    "## Run EM for a couple of iterations\n",
    "    Theta_old = Theta.copy( )\n",
    "    ll_hist, theta_hist, pi_hist, S = __em_algorithm( X, Theta, Pi, verbose = False, niter = 10 )\n",
    "    Theta, Pi = theta_hist[-1], pi_hist[-1]\n",
    "    kiter += S[ 'iter' ]\n",
    "    print \"Iteration %d: avg. log-lik: %.8e\" % ( kiter, ll_hist[-1] / X.shape[ 0 ], )\n",
    "    theta_show( Theta - Theta_old, 15, interpolation = 'nearest' )\n",
    "    if S['status'] != -1 :\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print Pi * X.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print Pi_init * X.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "theta_show( Theta, 5, cmap = plt.cm.gray, interpolation = 'nearest' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A random variable $X\\sim \\text{Beta}(\\alpha,\\beta)$ if the law of $X$ has density\n",
    "$$p(u) = \\frac{\\Gamma(\\alpha+\\beta)}{\\Gamma(\\alpha)\\Gamma(\\beta)} u^{\\alpha-1}(1-u)^{\\beta-1} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\log p(X,Z|\\Theta) = \\sum_{s=1}^n \\log \\prod_{k=1}^K \\Bigl[ \\pi_k \n",
    "\\prod_{i=1}^N \\prod_{j=1}^M \n",
    "\\frac{\\Gamma(\\alpha_{kij}+\\beta_{kij})}{\\Gamma(\\alpha_{kij})\\Gamma(\\beta_{kij})} x_{sij}^{\\alpha_{kij}-1}(1-x_{sij})^{\\beta_{kij}-1} \\Bigr]^{1_{z_s = k}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",<commit_msg>Commented the animator<commit_after>{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os, time, warnings\n",
    "warnings.filterwarnings( \"ignore\" )\n",
    "# from IPython.core.display import HTML\n",
    "from IPython.display import display, HTML\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "## This is the correct way to use the random number generator,\n",
    "##  since it allows finer control.\n",
    "rand = np.random.RandomState( 321 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# _EM_ and NIST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before proceeding, it seems pedagogically necessary (at least for myself) to revise the **EM**-slgorithm and show its \"correctness\", so to say.\n",
    "\n",
    "The $\\TeX$ markup used here uses the \"align*\" environment and thus should not be viewed though nbViewer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A brief description of the **EM** algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The EM algorithm seeks to maximize the likelihood by means of successive application of two steps: the E-step and the M-step.\n",
    "\n",
    "For any probability measure $Q$ on the space of latent variables $Z$ with density $q$ the following holds:  \n",
    "\\begin{align*}\n",
    "\\log p(X|\\Theta)\n",
    "    &= \\int q(Z) \\log p(X|\\Theta) dZ\n",
    "     = \\mathbb{E}_q \\log p(X|\\Theta) \\\\\n",
    "   %% &= \\Bigl[p(X,Z|\\Theta) = p(Z|X,\\Theta) p(X|\\Theta) \\Bigr] \\\\\n",
    "    &= \\mathbb{E}_{Z\\sim q} \\log \\frac{p(X,Z|\\Theta)}{p(Z|X\\Theta)}\n",
    "     = \\mathbb{E}_{Z\\sim q} \\log \\frac{q(Z)}{p(Z|X,\\Theta)}\n",
    "     + \\mathbb{E}_{Z\\sim q} \\log \\frac{p(X,Z|\\Theta)}{q(Z)} \\\\ \n",
    "    &= KL\\bigl(q\\|p(\\cdot|X,\\Theta)\\bigr) + \\mathcal{L}\\bigl(q, \\Theta\\bigr)\n",
    "\\end{align*}  \n",
    "\n",
    "since the Bayes theorem posits that $p(X,Z|\\Theta) = p(Z|X,\\Theta) p(X|\\Theta)$. Call this equiation the **\"master equation\"**.\n",
    "\n",
    "Now note that since the Kullback-Leibler divergence is always non-negative, one has the following inequality:\n",
    "$$\\log p(X|\\Theta) \\geq \\mathcal{L}\\bigl(q, \\Theta\\bigr)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try to make the lower bound as large as possible by changing $\\Theta$ and varying $q$. But first note that the \n",
    "left-hand side of the **master equation** is independent of $q$, whence maximization of $\\mathcal{L}$ with respect to $q$ (with $\\Theta$ fixed) is equivalent to minimization of $KL\\bigl(q\\|p(\\cdot|X,\\Theta)\\bigr)$ with respect to $q$ taking $\\Theta$ fixed. Since $q$ is arbitrary, the optimal minimizer $q^*_\\Theta$ is $q^*(Z|\\Theta) = p(Z|X,\\Theta)$ for all $Z$.\n",
    "\n",
    "Now at the optimal distributuion $q^*_\\Theta$ the **master equation** becomes\n",
    "$$ \\log p(X|\\Theta)\n",
    "= \\mathcal{L}\\bigl(q^*_\\Theta, \\Theta\\bigr)\n",
    "= \\mathbb{E}_{Z\\sim q^*_\\Theta} \\log \\frac{p(X,Z|\\Theta)}{q^*(Z|\\Theta)}\n",
    "= \\mathbb{E}_{Z\\sim q^*_\\Theta} \\log p(X,Z|\\Theta) - \\mathbb{E}_{Z\\sim q^*_\\Theta} \\log q^*(Z|\\Theta)\n",
    "$$\n",
    "for any $\\Theta$. Thus the problem of log-likelihood maximization reduces to that of maximizing the sum of expectations on the right-hand side.\n",
    "\n",
    "This new problem does not seem to be tractable in general since the optimization paramters $\\Theta$ affect both the expected log-likelihood $\\log p(X,Z|\\Theta)$ under $Z\\sim q^*_\\Theta$ and the entropy of the optimal distribution of the latent variables $Z$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully using an iterative procedure which switches between the computation of $q^*_\\Theta$ and the maximization of $\\Theta$ might be effective. Consider the folowing :\n",
    "* **E**-step: considering $\\Theta_i$ as given and fixed find $q^*_{\\Theta_i} = \\mathop{\\text{argmin}}_q KL\\bigl(q\\|p(\\cdot|X,\\Theta_i)\\bigr)$ and set $q_{i+1} = q^*_{\\Theta_i}$;\n",
    "* **M**-step: considering $q_{i+1}$ as given, solve $\\mathcal{L}(q_{i+1},\\Theta) \\to \\mathop{\\text{max}}_\\Theta$, where \n",
    "$$ \\mathcal{L}(q,\\Theta) = \\mathbb{E}_{Z\\sim q} \\log p(X,Z|\\Theta) - \\mathbb{E}_{Z\\sim q} \\log q(Z) $$\n",
    "\n",
    "The fact that $q_i$ is considered fixed makes the optimization of $\\mathcal{L}(q_i,\\Theta)$ equivalent to maximization of the expected log-likelihood, since the entropy term is fixed. Therefore the **M**-step becomes:\n",
    "* given $q_{i+1}$ find $\\Theta^*_{i+1} = \\mathop{\\text{argmax}}_\\Theta \\mathbb{E}_{Z\\sim q_{i+1}} \\log p(X,Z|\\Theta)$ and put $\\Theta_{i+1} = \\Theta^*_{i+1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, if the latent variables are mutually independent, then the optimal $q$ must be factorizable into marginal densities and:\n",
    "\n",
    "\\begin{align*}\n",
    "    KL\\bigl(q\\|p(\\cdot|X,\\Theta)\\bigr)\n",
    "    &= \\mathbb{E}_{Z\\sim q} \\log q(Z) - \\sum_j \\mathbb{E}_{z_j\\sim q_j} \\log p(z_j|X,\\Theta)\\\\\n",
    "    &= \\sum_j \\mathbb{E}_{z_j\\sim q_j} \\log q_j(z_j) - \\sum_j \\mathbb{E}_{z_j\\sim q_j} \\log p(z_j|X,\\Theta)\n",
    "     = \\sum_j KL\\bigl(q_j\\|p_j(|X,\\Theta)\\bigr)\n",
    "\\end{align*}\n",
    "where $q_j$ is the marginal desity of $z_j$ in $q(Z)$ (the last term in the first line comes from the Fubini theorem).\n",
    "\n",
    "Therefore the **E**-step could be reduced to a set of minimization problems with respect to one-dimensional density functions:\n",
    "$$ q_j^* = \\mathop{\\text{argmin}}_{q_j} KL\\bigl(q_j\\|p_j(\\cdot|X,\\Theta)\\bigr) $$\n",
    "since the Kulback-Leibler divergence in this case in additively separable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correctness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that the **master equation** is an identity: for all densities $q$ on $Z$ and for all admissible parameters $\\Theta$\n",
    "$$ \\log p(X|\\Theta) = KL\\bigl(q\\|p(\\cdot|X,\\Theta)\\bigr) + \\mathcal{L}\\bigl(q, \\Theta\\bigr) $$\n",
    "\n",
    "Hence if after the **E**-step the Kulback-Leibler divergence is reduced:\n",
    "$$ KL\\bigl(q'\\|p(\\cdot|X,\\Theta)\\bigr) \\leq KL\\bigl(q\\|p(\\cdot|X,\\Theta)\\bigr) $$\n",
    "then for the same set of parameters $\\Theta$ one has\n",
    "$$ \\mathcal{L}(q,\\Theta) \\leq \\mathcal{L}(q',\\Theta) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just after the **E**-step one has $q_{i+1} = p(Z|X,\\Theta_i)$, whence $KL\\bigl(q_{i+1}\\|p(\\cdot|X,\\Theta_i)\\bigr) = 0$. In turn, this implies via the **master equation** that the following equality holds:\n",
    "$$ \\log p(X|\\Theta_i) = \\mathcal{L}(q_{i+1},\\Theta_i) $$\n",
    "\n",
    "After the **M**-step, since $\\Theta_{i+1}$ is a maximizer, or at least an \"_improver_\" of $\\mathcal{L}(q_{i+1},\\Theta)$ compared to its value at $(q_i,\\Theta_i)$, one has\n",
    "$$ \\mathcal{L}(q_{i+1},\\Theta_i) \\leq \\mathcal{L}(q_{i+1},\\Theta_{i+1}) $$\n",
    "\n",
    "Threfore the effect of a single complete round of **EM** on the log-likelihood itself is:\n",
    "$$ \\log p(X|\\Theta_i) = \\mathcal{L}(q_{i+1},\\Theta_i) \\leq \\mathcal{L}(q_{i+1},\\Theta_{i+1}) \\leq \\mathcal{L}(q_{i+2},\\Theta_{i+1}) = \\log p(X|\\Theta_{i+1}) $$\n",
    "where the equality is achieved between the **E** and the **M** step within one round. This implies that **EM** indeed iteratively improves the log-likihood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note, that in the general case, without attaining zero Kulback-Leibler divergence at the $E$-step, one cannot be sure that the real log-likelihood is improved by each iteration and one can just say that\n",
    "$$\\mathcal{L}(q_{i+1},\\Theta_i) \\leq \\log p(X|\\Theta_i)$$\n",
    "which does not uncover a relationship with $\\log p(X|\\Theta_{i+1})$. And without the guarantee that **EM** improves the log-likelihood to the maximum one cannot be sure about the consistency of the estimators. The key question is whether the lower bound $\\mathcal{L}(q,\\Theta)$ is any good."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application of the EM to NIST data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each image is a random element in a discrete probability space $\\Omega = \\{0,1\\}^{N\\times M}$ with product-measure\n",
    "$$\\mathbb{P}(\\omega) = \\prod_{i=1}^N\\prod_{j=1}^M \\theta_{ij}^{\\omega_{ij}} (1-\\theta_{ij})^{1-\\omega_{ij}}$$\n",
    "for any $\\omega\\in \\Omega$. In particular $M=N=28$. Basically each bit of the image is independent of any other bit and each one is a Bernoulli random variable with parameter $\\theta_{ij}$: $\\omega_{ij}\\sim \\text{Bern}(\\theta_{ij})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's apply the EM algorithm to this dataset. The proposed model is the following.\n",
    "\n",
    "Consider a mixture model of discrete probability spaces. Suppose there are $K$ componets in the mixture. Then each image is distributed according to the following law:\n",
    "$$p(\\omega|\\Theta)\n",
    "= \\sum_{k=1}^K \\pi_k p_k(\\omega|\\theta_k)\n",
    "= \\sum_{k=1}^K \\pi_k \\prod_{i=1}^N \\prod_{j=1}^M \\theta_{kij}^{\\omega_{ij}} (1-\\theta_{kij})^{1-\\omega_{ij}}$$\n",
    "where $\\theta_{kij}$ is the paramter of the probability distribution of the $(i,j)$-th random variable (pixel) in the $k$-th class, and $\\pi_k$ is the (prior) porbability of the $k$-th mixutre to generate a random element, $\\sum_{k=1}^K \\pi_k= 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose $X=(x_i)_{i=1}^n \\in \\Omega^n$ is the dataset. The log-likelihood is given by\n",
    "$$ \\log p(X|\\Theta) = \\sum_{s=1}^n \\log \\sum_{k=1}^K \\pi_k\n",
    "\\prod_{i=1}^N \\prod_{j=1}^M \\theta_{kij}^{x_{sij}} (1-\\theta_{kij})^{1-x_{sij}}$$\n",
    "where $x_{sij}\\in\\{0,1\\}$ -- is the value of the the $(i,j)$-th pixel at the $s$-th observation.\n",
    "\n",
    "If the source $Z=(z_i)_{i=1}^n$ components of the mixture at each datapoint were known, then the log-likelihood would have been\n",
    "$$ \\log p(X,Z|\\Theta) = \\sum_{s=1}^n \\log \\prod_{k=1}^K \\Bigl[ \\pi_k \n",
    "\\prod_{i=1}^N \\prod_{j=1}^M \\theta_{kij}^{x_{sij}} (1-\\theta_{kij})^{1-x_{sij}} \\Bigr]^{1_{z_s = k}}$$\n",
    "where $1_{z_s = k}$ is the indicator and take the value $1$ if $\\{z_s = k\\}$ and $0$ otherwise ($1_{\\{k\\}}(z_s)$ is another notation)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The log-likelihood simplifies to\n",
    "$$ \\log p(X,Z|\\Theta) = \\sum_{s=1}^n \\sum_{k=1}^K 1_{z_s = k} \\Bigl( \\log \\pi_k + \n",
    "\\sum_{i=1}^N \\sum_{j=1}^M \\bigl( x_{sij} \\log \\theta_{kij} + (1-x_{sij}) \\log (1-\\theta_{kij}) \\bigr) \\Bigr) $$\n",
    "and further into a more separable form\n",
    "$$ \\log p(X,Z|\\Theta)\n",
    "= \\sum_{s=1}^n \\sum_{k=1}^K 1_{z_s = k} \\log \\pi_k\n",
    "+ \\sum_{s=1}^n \\sum_{k=1}^K 1_{z_s = k} \\Bigl( \\sum_{i=1}^N \\sum_{j=1}^M x_{sij} \\log \\theta_{kij}\n",
    "+ \\sum_{i=1}^N \\sum_{j=1}^M (1-x_{sij}) \\log (1-\\theta_{kij}) \\Bigr)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The expected log-likelihood under $z_s\\sim q_s$ with $\\mathbb{P}(z_s=k|X) = q_{sk}$, is given by\n",
    "$$ \\mathbb{E}\\log p(X,Z|\\Theta)\n",
    "= \\sum_{s=1}^n \\sum_{k=1}^K q_{sk} \\log \\pi_k\n",
    "+ \\sum_{s=1}^n \\sum_{k=1}^K q_{sk} \\sum_{i=1}^N \\sum_{j=1}^M \\bigl( x_{sij} \\log \\theta_{kij} + (1-x_{sij}) \\log (1-\\theta_{kij}) \\bigr)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Analytic solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the **E**-step one must compute $q^*(Z) = \\mathbb{P}(z_s=k|X) = \\hat{q}_{sk}$ based on the value of $\\Theta = ((\\pi_k), (\\theta_{kij}))$.\n",
    "$$\\hat{q}_{sk}\n",
    "= \\frac{p(x_s|z_s=k,\\Theta) p(z_s=k)}{\\sum_{l=1}^K p(x_s|z_s=l,\\Theta) p(z_s=l)}\n",
    "\\propto \\pi_k \\prod_{i=1}^N \\prod_{j=1}^M \\theta_{kij}^{x_{sij}} (1-\\theta_{kij})^{1-x_{sij}}\n",
    "$$\n",
    "and\n",
    "$$ q^*(Z) = \\prod_{s=1}^n q_{s z_s} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to improve numerical stability and avoid numerical underflow it is better to use the following procedure for computation of the conditional probability:\n",
    "$$ l_{sk} = \\log \\pi_k + \\sum_{i=1}^N \\sum_{j=1}^M \\log \\bigl( \\theta_{kij}^{x_{sij}} (1-\\theta_{kij})^{1-x_{sij}} \\bigr)$$\n",
    "set $l^*_s = \\max_k l_{sk}$ and compute the log-sum\n",
    "$$ \\hat{l}_s = l^*_s + \\log \\sum_{k=1}^K \\text{exp}\\Bigl\\{ l_{sk} - l^*_s \\Bigr\\} $$\n",
    "This seemingly redunant subctration and addition of $l^*_s$ helps avoid underflow during the numerical exponentiation.\n",
    "However if the underflow still took place in the sum-exp, then those values of $\\hat{l}_s$ that evaluated to floating point _NAN_ can reliably be replaced with respective values of $l^*_s$.  After this sanitization the **E**-step's optimal distribution is calculated as:\n",
    "$$ \\hat{q}_{sk} = \\text{exp}\\Bigl\\{ l_{sk} - \\hat{l}_s \\Bigr\\} $$\n",
    "which would yield a numerically accurate distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $l^*_s >> l_{sk}$ for all $k$ but such that $l^*_s = l_{sk}$ (let it be $k^*$), then an underflow occurs at the sum-exp step, whence for some very small $\\epsilon > 0$ one has\n",
    "$$\\hat{l}_s = l^*_s + \\log (1+\\epsilon) $$\n",
    "whence\n",
    "$$\\hat{q}_{sk} = \\text{exp}\\bigl\\{l_{sk} - \\hat{l}_s\\bigr\\} = (1+\\epsilon)^{-1} \\cdot \\text{exp}\\bigl\\{l_{sk} - l^*_s \\bigr\\}$$\n",
    "For $k=k^*$ on has $\\hat{q}_{sk} = \\frac{1}{1+\\epsilon}\\approx 1$, and for $k\\neq k^*$ -- $\\hat{q}_{sk} = \\frac{\\eta}{1+\\epsilon} \\approx 0$ for some extremely small $\\eta>0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variables in the code have the following dimensions:\n",
    "* $\\theta \\in [0,1]^{K\\times (N\\times M)}$;\n",
    "* $\\pi \\in [0,1]^{1\\times K}$;\n",
    "* $x \\in \\{0,1\\}^{n\\times (N\\times M)}$;\n",
    "* $z \\in [0,1]^{n\\times K}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## A bunch of wrappers to match the task specifications\n",
    "def posterior( x, clusters ) :\n",
    "    pi = np.ones( clusters.shape[ 0 ], dtype = np.float ) / clusters.shape[ 0 ]\n",
    "    q, ll = __posterior( x, theta = clusters, pi = pi )\n",
    "    return q\n",
    "\n",
    "## The likelihood is a byproduct of the E-step's minimization of Kullback-Leibler\n",
    "def likelihood( x, clusters ) :\n",
    "    pi = np.ones( clusters.shape[ 0 ], dtype = np.float ) / clusters.shape[ 0 ]\n",
    "    q, ll = __posterior( x, theta = clusters, pi = None )\n",
    "    return ll\n",
    "\n",
    "## The core procedure for computing the conditional density of classes\n",
    "def __posterior( x, theta, pi ) :\n",
    "## Unfortunately sometimes there can be negative machine zeros, which\n",
    "##  spoil the log-likelihood computation by poisoning with NANs.\n",
    "## That is why the theta array is restricted to [0,1].\n",
    "    theta_clipped = np.clip( theta, 0.0, 1.0 )\n",
    "## Note that the power formulation is just a mathematically convenient\n",
    "##  way of writing \\theta if x=1 or (1-\\theta) otherwise. This comapct\n",
    "##  way of computing had to be abandoned, since it required huge amounts\n",
    "##  of memory.\n",
    "#     ll_skij = np.where( x[:, np.newaxis] > 0,\n",
    "#         np.log( theta_clipped ), np.log( 1 - theta_clipped ) )\n",
    "#     ll_sk = np.sum( ll_skij, axis = ( 2, ) ) + np.log( pi )\n",
    "#     del ll_skij\n",
    "## Iterate over classes\n",
    "    ll_sk = np.zeros( ( x.shape[ 0 ], theta.shape[ 0 ] ), dtype = np.float )\n",
    "    for k in xrange( theta.shape[ 0 ] ) :\n",
    "        ll_sij = np.where( x > 0, np.log( theta_clipped[ k ] ), np.log( 1 - theta_clipped[ k ] ) )\n",
    "        ll_sk[ :, k ] = np.sum( ll_sij, axis = ( 1, ) ) + np.log( pi[ k ] )\n",
    "## Find the largest unnormalized probability.\n",
    "    llstar_s = np.reshape( np.max( ll_sk, axis = ( 1, ) ), ( ll_sk.shape[ 0 ], 1 ) )\n",
    "## Compute the log-sum-exp of the individual log-likelihoods\n",
    "    ll_s = np.reshape( np.log( np.sum( np.exp( ll_sk - llstar_s ), axis = ( 1, ) ) ), ( ll_sk.shape[ 0 ], 1 ) )\n",
    "## The sum-exp could never be anything lower than 1, since at least one\n",
    "##  element of each row of ll_sk has to be lstar_s, whence the respective\n",
    "##  difference should be zero and the exponent -- 1. Thus even if the\n",
    "##  rest of the sum is close to machine zero, the logarithm would still\n",
    "##  return 0.\n",
    "    ll_s += llstar_s\n",
    "    del llstar_s\n",
    "## Normalise the likelihoods to get conditional probability, and compute\n",
    "##  the sum of the log-denominator, which is the log-likelihood.\n",
    "    return np.exp( ll_sk - ll_s ), np.sum( ll_s )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the **M**-step for some fixed $q(Z)$ one solves $\\mathbb{E}\\log p(X,Z|\\Theta)\\to \\max_\\Theta$ subject to $\\sum_{k=1}^K \\pi_k = 1$ which is a convex optimization problem with respect to $\\Theta$, since the log-likelihood as a linear combination of convex functions is convex. The first order condition is $\\sum_{s=1}^n \\frac{q_{sk}}{\\pi_k} - \\lambda = 0$ for all $k=1,\\ldots,K$, whence $ \\lambda = \\sum_{s=1}^n \\sum_{l=1}^K q_{sl} = n $ and finally\n",
    "$$\\hat{\\pi}_k = \\frac{\\sum_{s=1}^n q_{sk}}{n}$$\n",
    "For $\\theta_{kij}$, $i=1,\\ldots,N$, $j=1,\\ldots,M$ and $k=1,\\ldots,K$ the FOC is\n",
    "$$ \\sum_{s=1}^n q_{sk} \\frac{x_{sij}}{\\theta_{kij}} - \\sum_{s=1}^n q_{sk} \\frac{1-x_{sij}}{1-\\theta_{kij}} = 0 $$\n",
    "whence\n",
    "$$\\hat{\\theta}_{kij} =  \\frac{\\sum_{s=1}^n q_{sk} x_{sij}}{ \\sum_{s=1}^n q_{sk} } = \\frac{\\sum_{s=1}^n q_{sk} x_{sij}}{ n \\hat{\\pi}_k }$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## A wrapper to match the task specifications\n",
    "def learn_clusters( x, z ) :\n",
    "    theta, pi = __learn_clusters( x, z )\n",
    "    return theta\n",
    "\n",
    "## The E-step is simple: just compute the optimal parameters under\n",
    "##  the current conditional distribution of the latend variables.\n",
    "def __learn_clusters( x, z ) :\n",
    "## The prior class probabilities\n",
    "    pi = np.sum( z, axis = ( 0, ) )\n",
    "## Pixel probabilities conditional on the calss\n",
    "    theta = np.dot( z.T, x ) / np.reshape( pi, ( z.shape[ 1 ], 1 ) )\n",
    "## Return\n",
    "    return theta, np.reshape( pi, ( 1, pi.shape[ 0 ] ) ) / np.sum( pi )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A it has been mentioned eariler, the EM algorithm switches between **E** and **M** steps until convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## A wrapper for the core em algorithm above\n",
    "def em_algorithm( x, K, maxiter, verbose = True, rel_eps = 1e-4 ) :\n",
    "## Initialize the model parameters with uniform [0,1] random numbers\n",
    "    theta = rand.uniform( size = ( K, x.shape[ 1 ] ) )\n",
    "## Run the em algorithm\n",
    "    ll, theta, pi, status = __em_algorithm( x, theta_1 = theta,\n",
    "        pi_1 = None, niter = maxiter, rel_eps = rel_eps, verbose = verbose )\n",
    "## Return the history of theta and the final log liklihood\n",
    "    if verbose :\n",
    "        if status[ 'status' ] != 0 :\n",
    "            print \"Convergence not achieved. %d\" % ( status[ 'status' ], )\n",
    "    return theta, ll"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The procedure above actually invokes the true EM core, defined below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## The core of the EM algorithm\n",
    "def __em_algorithm( x, theta_1, pi_1 = None, niter = 1000, rel_eps = 1e-4, verbose = True ) :\n",
    "## If we were supplied with an initial estimate of the prior distribution,\n",
    "##  then assume the full model is needed.\n",
    "    full_model = pi_1 is not None\n",
    "## If the prior cluster probabilities are not supplied, assume uniform distribution.\n",
    "    pi_1 = pi_1 if pi_1 is not None else np.ones( theta_1.shape[ 0 ], dtype = np.float ) / theta_1.shape[ 0 ]\n",
    "## Allocate the necessary space for the history of model estimates\n",
    "    theta_hist = np.empty( ( 0, ) + theta_1.shape, dtype = np.float )\n",
    "    pi_hist = np.empty( ( 0, theta_1.shape[ 0 ] ), dtype = np.float )\n",
    "    ll_hist = np.empty( 0, np.float )\n",
    "## Set \"old\" estimated to zero. At this line the current estimates are in fact\n",
    "##  the initially provided ones.\n",
    "    theta_0, pi_0 = np.zeros_like( theta_1 ), np.zeros_like( pi_1 )\n",
    "## Initialize the loop\n",
    "    status, kiter, rel_theta, rel_pi, ll = -1, 0, np.nan, np.nan, -np.inf\n",
    "    while kiter < niter :\n",
    "## Dump the current estimators and other information.\n",
    "        if verbose :\n",
    "            print( \"Iteration %d: avg. log-lik: %.3f, $\\\\Theta$ div. %.3f, $\\\\Pi$ div. %.3f\" % (\n",
    "                kiter, ll / x.shape[ 0 ], rel_theta, rel_pi ) )\n",
    "            show_data( theta_1 - theta_0 if True else theta_1, n = theta_0.shape[ 0 ],\n",
    "                n_col = min( 10, theta_0.shape[ 0 ] ), cmap = plt.cm.hot, interpolation = 'nearest' )\n",
    "## The convergence criterion is the L^∞ norm of relative L^1 errors\n",
    "        if max( rel_pi, rel_theta ) < rel_eps :\n",
    "            status = 0\n",
    "            break ;\n",
    "## Overwrite the initial estimates\n",
    "        theta_0, pi_0 = theta_1, pi_1\n",
    "## E-step: call the core posterior function to get both the log-likelihood\n",
    "##  and the estimate of the conditional distribution.\n",
    "        z_1, ll = __posterior( x, theta_0, pi_0 )\n",
    "## M-step: compute the optimal parameters under the current estimate of the posterior\n",
    "        theta_1, pi_1 = __learn_clusters( x, z_1 )\n",
    "## Discard the computed estimate of pi if the model is discriminative (conditional likelihood).\n",
    "        if not full_model :\n",
    "            pi_1 = pi_0\n",
    "## Record the current estimates to the history\n",
    "        theta_hist = np.vstack( ( theta_hist, theta_1[np.newaxis] ) )\n",
    "        pi_hist = np.vstack( ( pi_hist, pi_1 ) )\n",
    "        ll_hist = np.append( ll_hist, ll )\n",
    "## Check for bad float numbers\n",
    "        if not ( np.all( np.isfinite( theta_1 ) ) and np.all( np.isfinite( pi_1 ) ) ) :\n",
    "            status= -2\n",
    "            break ;\n",
    "## Check convergence: L^1 relative error. If the relative margin is exactly\n",
    "##  zero, then return NaNs. This makes teh loop exhaust all iterations, since\n",
    "##  any comprison agains a NaN returns False.\n",
    "        rel_theta = np.sum( np.abs( theta_1 - theta_0 ) / ( np.abs( theta_0 ) + rel_eps ) ) if rel_eps > 0 else np.nan\n",
    "        rel_pi = np.sum( np.abs( pi_1 - pi_0 ) / ( np.abs( pi_0 ) + rel_eps ) ) if rel_eps > 0 else np.nan\n",
    "## Next iteration\n",
    "        kiter += 1\n",
    "    return ll_hist, theta_hist, pi_hist, { 'status': status, 'iter': kiter }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Miscellanea"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to be able to plot more flexibly, define another arranger."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## A more flexible image arrangement\n",
    "def arrange_flex( images, n_row = 10, n_col = 10, N = 28, M = 28, fill_value = 0 ) :\n",
    "## Create the final grid of images row-by-row\n",
    "    im_grid = np.full( ( n_row * N, n_col * M ), fill_value, dtype = images.dtype )\n",
    "    for k in range( min( images.shape[ 0 ], n_col * n_row ) ) :\n",
    "## Get the grid cell at which to place the image\n",
    "        i, j = ( k // n_col ) * N, ( k % n_col ) * M\n",
    "## Just put the image in the cell\n",
    "        im_grid[ i:i+N, j:j+M ] = np.reshape( images[ k ], ( N, M, ) )\n",
    "    return im_grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The procedure below loads the MNIST data from a comma-separated text file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def load_mnist_from_csv( filename ) :\n",
    "## Read the CSV file\n",
    "    data = np.loadtxt( open( filename, \"rb\" ), dtype = np.short, delimiter = \",\", skiprows = 0 )\n",
    "## Peel off the lables\n",
    "    return data[:,1:], data[:,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A procedure for pickling numpy arrays of floats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import base64, zlib\n",
    "\n",
    "def np_pack( data ) : return base64.b64encode( zlib.compress( data ) )\n",
    "def np_unpack( text ) : return np.frombuffer( zlib.decompress(\n",
    "    base64.decodestring( text ) ), dtype = np.float )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A handy procedure for one-line visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def show_data( data, n, n_col = 10, **kwargs ) :\n",
    "## Get the number of rows necessary to plot the needed number of images\n",
    "    n_row = ( n + n_col - 1 ) // n_col\n",
    "## Set the dimensions of the figure\n",
    "    plt.figure( figsize = ( n_col, n_row ) )\n",
    "## Plot!\n",
    "    plt.imshow( arrange_flex( data[:n], n_col = n_col, n_row = n_row ), **kwargs )\n",
    "    plt.show( )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a convenient procedure for running experiments. By setting relative error to zero the algorithm is forced to exhaust all the allocated iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def experiment( data, K, maxiter, verbose = True, until_convergence = False ) :\n",
    "## Run the EM\n",
    "    return em_algorithm( data, K, maxiter, rel_eps = 1.0e-4 if until_convergence else 0.0, verbose = verbose )\n",
    "\n",
    "def visualize( data, clusters, ll, n_col = 2, plot_ll = True ) :\n",
    "## Display the result\n",
    "    print \"Final conditional log-likelihood value per observation achieved %f in %d iteration(s)\" % ( ll[-1] / data.shape[ 0 ], len( ll ) )\n",
    "## Plot the first difference of average log-likelihood\n",
    "    if plot_ll :\n",
    "        plt.figure( figsize = ( 16, 9 ) )\n",
    "        ax = plt.subplot(111)\n",
    "        ax.set_title( r\"avg. log-likelihood change between successive iterations (log scale)\" )\n",
    "        ax.plot( np.diff( ll / data.shape[ 0 ]  ) )\n",
    "        ax.set_ylabel( r\"$\\Delta_i \\frac{1}{n} \\sum_{s=1}^n \\mathbb{E}_{z_s\\sim q_i} \\log p(x_s,z_s|\\Theta_i)$\" )\n",
    "        ax.set_yscale( 'log' )\n",
    "## Plot the final estimates\n",
    "    if n_col > 0 :\n",
    "        show_data( clusters[-1], n = clusters.shape[1], n_col = n_col, cmap = plt.cm.gray, interpolation = 'nearest' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function that produces and embed a video in HTML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Make simple animations of the EM estimatora\n",
    "## http://jakevdp.github.io/blog/2013/05/12/embedding-matplotlib-animations/\n",
    "## http://jakevdp.github.io/blog/2012/08/18/matplotlib-animation-tutorial/\n",
    "from matplotlib import animation\n",
    "from IPython.display import HTML\n",
    "from tempfile import NamedTemporaryFile\n",
    "VIDEO_TAG = \"\"\"<video controls><source src=\"data:video/x-m4v;base64,{0}\" type=\"video/mp4\">Your browser does not support the video tag.</video>\"\"\"\n",
    "def display_animation( anim ) :\n",
    "    plt.close( anim._fig )\n",
    "    if not hasattr( anim, '_encoded_video' ) :\n",
    "        ffmpeg_writer = animation.FFMpegWriter( )\n",
    "        with NamedTemporaryFile( suffix = '.mp4' ) as f:\n",
    "            anim.save( f.name, fps = 15, extra_args = [ '-vcodec', 'libx264' ] )# , writer = ffmpeg_writer )\n",
    "            video = open( f.name, \"rb\" ).read( )\n",
    "        anim._encoded_video = video.encode( \"base64\" )\n",
    "    return HTML( VIDEO_TAG.format( anim._encoded_video ) )\n",
    "\n",
    "def animate( theta, ll, n_col = 10, n_row = 10, diff = True, interval = 20, **kwargs ) :\n",
    "## First set up the figure, the axis, and the plot element we want to animate\n",
    "    fig = plt.figure( )\n",
    "## Create two subplots\n",
    "    ax1, ax2 = fig.add_subplot( 211 ), fig.add_subplot( 212 )\n",
    "## Set the background\n",
    "    bg = arrange_flex( np.zeros_like( theta[ 0 ] ), n_col = n_col, n_row = n_row )\n",
    "## Initialize different ranges for the image aritst\n",
    "    if diff :\n",
    "        im = ax1.imshow( bg, vmin = -1.0, vmax = +1.0, **kwargs )\n",
    "    else :\n",
    "        im = ax1.imshow( bg, vmin = +0.0, vmax = +1.0, **kwargs )\n",
    "## Compute log-likelihood differences and initialize geomtery for the line artist\n",
    "    ll_diff = np.diff( ll ) / data.shape[ 0 ]\n",
    "    line, = ax2.plot( [], linestyle = \"-\", color = 'blue' )\n",
    "    ax2.set_xlim( -0.1, theta.shape[ 0 ] + 0.1 )\n",
    "    ax2.set_yscale( 'log' )\n",
    "    ax2.set_ylim( np.min( ll_diff ) * 0.9, np.max( ll_diff ) * 1.1 )\n",
    "## Animation function.  This is called sequentially\n",
    "    def update( i ) :\n",
    "## Compute the frame\n",
    "        if diff :\n",
    "            frame = theta[ i ] - theta[ i-1 ] if i > 0 else theta[ 0 ]\n",
    "            frame /= np.max( np.abs( frame ) )\n",
    "        else :\n",
    "            frame = theta[ i ]\n",
    "## Draw the frame on the image artist\n",
    "        im.set_data( arrange_flex( frame, n_col = n_col, n_row = n_row ) )\n",
    "        if i > 0 :\n",
    "## Show history on the line artist\n",
    "            line.set_data( np.arange( i ), ll_diff[ :i ] )\n",
    "## Return an iterator of artists in this frame\n",
    "        return im, line,\n",
    "## Call the animator.\n",
    "    return animation.FuncAnimation( fig, update, frames = theta.shape[ 0 ], interval = interval )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading and visualizing MNIST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Preamble"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obtain the MNIST data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if False :\n",
    "## Fetch MNIST dataset from SciKit and create a local copy.\n",
    "    from sklearn.datasets import fetch_mldata\n",
    "    mnist = fetch_mldata( \"MNIST original\", data_home = './data/' )\n",
    "    np.savez_compressed('./data/mnist/mnist_scikit.npz', data = mnist.data, labels = mnist.target )\n",
    "\n",
    "if False :\n",
    "## Fetch the data from the provided CSV (!) files and save as a compressed data blob\n",
    "    data, labels = load_mnist_from_csv( \"./data/mnist/mnist_train.csv\" )\n",
    "    np.savez_compressed( './data/mnist/mnist_train.npz', labels = labels, data = data )\n",
    "    data, labels = load_mnist_from_csv( \"./data/mnist/mnist_test.csv\" )\n",
    "    np.savez_compressed( './data/mnist/mnist_test.npz', labels = labels, data = data )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Study"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First of all load and binarize the training data using the values 127 as the threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "assert( os.path.exists( './data/mnist/mnist_train.npz' ) )\n",
    "with np.load( './data/mnist/mnist_train.npz', 'r' ) as npz :\n",
    "    mnist_labels, mnist_data = npz[ 'labels' ], np.array( npz[ 'data' ] > 127, np.int )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Case : $K=2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at some 6s and 9s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Mask\n",
    "inx_sixes, inx_nines = np.where( mnist_labels == 6 )[ 0 ], np.where( mnist_labels == 9 )[ 0 ]\n",
    "## Extract\n",
    "sixes = mnist_data[ rand.choice( inx_sixes, 90, replace = False ) ]\n",
    "nines = mnist_data[ rand.choice( inx_nines, 90, replace = False ) ]\n",
    "## Show\n",
    "show_data( sixes, n = 45, n_col = 15, cmap = plt.cm.gray, interpolation = 'nearest' )\n",
    "show_data( nines, n = 45, n_col = 15, cmap = plt.cm.gray, interpolation = 'nearest' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They do indeed look quite distinct. Now collect them into a single dataset and estimate the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data = mnist_data[ np.append( inx_sixes, inx_nines ) ]\n",
    "clusters, ll = experiment( data, 2, 30 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimate deltas show that the **EM** algorithm's **E**-step actually transfers the unlikely observations between classes, as is expected by constructon of the algorithm.\n",
    "\n",
    "Judging by the plot below, it turns out that 30 iterations is more than enough for the **EM** to get meaninful estimates the class ideals, represented by the probability porduct-measure on $\\Omega^{28\\times 28}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "visualize( data, clusters, ll )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see how well the **EM** algorithm performs on a model with more classes. But before that let's have a look at a random sample of the handwritten digits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "indices = np.arange( mnist_data.shape[ 0 ] )\n",
    "rand.shuffle( indices )\n",
    "plt.imshow( arrange_flex( mnist_data[ indices[:100] ] ), cmap = plt.cm.gray, interpolation = 'nearest' )\n",
    "plt.show( )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Case : $K=10, 15, 20$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The original size of the trainig sample is too large to fit in these RAM banks :) That is why I had to limit the sample to a random subset of 2000 observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = mnist_data[ indices[:2000]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the procedure that perfoems **EM** algorithm and return the history of the parameter estimates as well as the dynamics of the log-likelihood lower bonud."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "clusters, ll = experiment( data, 10, 50 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can clearly see, that $50$ iterations were not enough for the alogirithm to converge: though the changes are tiny, even on the log-scale, they are still unstable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "visualize( data, clusters, ll, n_col = 5 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Furthermore the returned templates of the mixture components and clearly suboptimal: the procedure seems to get stuck at individual examples. One way to remedy this is to use more iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "clusters_10, ll_10 = experiment( data, 10, 1000, verbose = True, until_convergence = True )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As one can see, increasing the number of iterations does not necessarily improve the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "visualize( data, clusters_10, ll_10, n_col = 5 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see if changing $K$ does the trick."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "clusters_15, ll_15 = experiment( data, 15, 1000, verbose = True, until_convergence = True )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "clusters_20, ll_20 = experiment( data, 20, 1000, verbose = True, until_convergence = True )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "visualize( data, clusters_10, ll_10, n_col = 10, plot_ll = False )\n",
    "visualize( data, clusters_15, ll_15, n_col = 10, plot_ll = False )\n",
    "visualize( data, clusters_20, ll_20, n_col = 10, plot_ll = False )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "display_animation( animate( clusters_20, ll_20, 10, 2, diff = False, cmap = plt.cm.gray, interpolation = 'nearest' ) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A random variable $X\\sim \\text{Beta}(\\alpha,\\beta)$ if the law of $X$ has density\n",
    "$$p(u) = \\frac{\\Gamma(\\alpha+\\beta)}{\\Gamma(\\alpha)\\Gamma(\\beta)} u^{\\alpha-1}(1-u)^{\\beta-1} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\log p(X,Z|\\Theta) = \\sum_{s=1}^n \\log \\prod_{k=1}^K \\Bigl[ \\pi_k \n",
    "\\prod_{i=1}^N \\prod_{j=1}^M \n",
    "\\frac{\\Gamma(\\alpha_{kij}+\\beta_{kij})}{\\Gamma(\\alpha_{kij})\\Gamma(\\beta_{kij})} x_{sij}^{\\alpha_{kij}-1}(1-x_{sij})^{\\beta_{kij}-1} \\Bigr]^{1_{z_s = k}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
{'fraction_non_alphanumeric': 0.2104779622866855, 'fraction_numerical': 0.012115574542364317, 'mean_word_length': 2.581434564448466, 'pattern_counts': {'":': 741, '<': 20, '<?xml version=': 0, '>': 40, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '256573', 'n_tokens_mistral': 33754, 'n_tokens_neox': 32169, 'n_words': 8549}