<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step1:  <p style = "font-size
Step2:  Preamble
Step3:  EM and MNIST
The $\TeX$ markup used here uses the "align*" environment and thus should not be viewed though nbViewer.
Before proceeding, it seems pedagogically necessary (at least for myself) to revise the EM-slgorithm and show its "correctness", so to say.
A brief description of the EM algorithm
The EM algorithm seeks to maximize the likelihood by means of successive application of two steps
Step4:  Classify using the maximum aposteriori rule.
Step5:  A procedure to compute the log-likelihood of each observaton with respect to each mixture component. Used in the posterior computation.
Step6:  The actual procedure for computing the E*-step
Step7:  Analytic solution
Step8:  A wrapper to match the assignment specifications.
Step9:  A it has been mentioned eariler, the EM algorithm switches between E and M steps until convergence.
Step10:  The procedure above actually invokes the true EM core, defined below.
Step11:  Define a convenient procedure for running experiments. By setting relative error to zero the algorithm is forced to exhaust all the allocated iterations.
Step12:  Miscellanea
Step13:  The folowing pair of procedures are used to plot the digits in a clear manner. The first one just creates a canvas for the image
Step14:  This procedure displays the images on a nice plot. Used for one-line visualization.
Step15:  MIscellanea
Step17:  Define a function that produces (using ffmpeg) and embeds a video in HTML into IPython
Step18:  Miscellanea
Step19:  Or obtain the data from the provided CSV files.
Step20:  Study
First of all load and binarize the training data using the value 127 as the threshold.
Step21:  Case 
Step22:  They do indeed look quite distinct. Now collect them into a single dataset and estimate the model.
Step23:  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.
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}$.
Step24:  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.
Step25:  Case 
Step26:  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.
Step27:  One can clearly see, that $50$ iterations were not enough for the alogirithm to converge
Step28:  Let's see if changing $K$ does the trick.
Step29:  For what values of $K$ was it possible to infer the templates of all digits?
Step30:  Obviously, the model with more mixture components is more likely to produce "templates" for all digits. For larger $K$ this is indeed the case.
Having run this algorithm for many times we are able to say that the digits $3$ and $8$, $4$ and $9$ and sometimes $5$ tend to be poorly separated. Furthermore due to there being many different handwritten variations of the same digit one should estimate a model with more classes. 
The returned templates of the mixture components are clearly suboptimal
Step31:  As one can see, increasing the number of iterations does not necessarily improve the results.
Step32:  Judging by the plot of the log-likelihood, the fact that the EM is guaranteed to converge to local maxima and does so extremely fast, there was no need for more than 120-130 iterations. The chages in the log-likelihood around that number of iterations are of the order $10^{-4}$. Since we are working in finite precision arithmetic (double), the smallest precision is $\approx 10^{-14}$.
Let's see the dynamics of the estimated of the EM iterations. You will have to ensure that ffMPEG is installed (Windows
Step33:  The parameter estimates of the EM stabilize pretty quickly. In fact most templates stabilize by iterations 100-120.
Choosing $K$
Among many methods of model selection, let's use simple training sample fittness score, givne by the value of the log-likelihood. Becasue the models are nested with respect to the number of mixture components, one should expect the likelihood to be a non decreasing function of $K$ (on average dut to randomization of the initial parameter values).
For large enough $K$ this method may lead to overfitting.
Step34:  Indeed the log-likelihood does not decrease with $K$ on average. Nevertheless the model with the highes likelihood turs out to have this many mixture components
Step35:  A nice, yet expected coincidence 
Step36:  ... and get the posterior mixture component probabilities.
Step37:  Use a simple majority rule to automatically assign lables to templates.
Step38:  Assign the labels $l$ to templates $t$ according to its score, based on the average of the top-$5$ log-likelihoods of observations with label $l$ and classfified with template $t$.
Step39:  Compare the label assignments. Here are the templates.
Step40:  These are the templates, which were assigned different labels by the majority and "trust" methods.
Step41:  Here are the pictures of templates ordered according to their label.
Step42:  Classification
Step43:  Let's see the best template for each test observation in some sub-sample.
Step44:  The digits are shown in pairs
Step45:  Not surprisingly, majority- and likelihood-based classification accuracies are close.
Let's see which test observations the model considers an artefact and for which it cannot reliably assign a template
Step46:  </hr>
Let's see how more pasrimonious models fare with respect to accuracy on thte test sample.
Accuracy of $K=30$
Step47:  Accuracy of $K=20$
Step48:  Accuracy of $K=15$
Step49:  Accuracy of $K=10$
Step50:  As one can see test sample accuracy of the model falls drammatically for less number of mixture components. This was expected, since due to various reasons, one being thet the data is handwritten, it is higly unlikely, that a single digit would have only one template.
Step51:  <br/><p style="font-size
Step52:  <hr/>
A random variable $X\sim \text{Beta}(\alpha,\beta)$ if the law of $X$ has density
$$p(u) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} u^{\alpha-1}(1-u)^{\beta-1} $$
$$ \log p(X,Z|\Theta) = \sum_{s=1}^n \log \prod_{k=1}^K \Bigl[ \pi_k 
\prod_{i=1}^N \prod_{j=1}^M 
\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}}$$
\begin{align}
\mathbb{E}q \log p(X,Z|\Theta)
&= \sum{k=1}^K \sum_{s=1}^n q_{sk} \log \pi_k \
&+ \sum_{k=1}^K \sum_{i=1}^N \sum_{j=1}^M \sum_{s=1}^n q_{sk} \bigl(
    \log \Gamma(\alpha_{kij}+\beta_{kij}) - \log \Gamma(\alpha_{kij}) - \log \Gamma(\beta_{kij}) \bigr) \
&+ \sum_{k=1}^K \sum_{i=1}^N \sum_{j=1}^M \sum_{s=1}^n q_{sk} \bigl(
    (\alpha_{kij}-1) \log x_{sij} + (\beta_{kij}-1) \log(1-x_{sij}) \bigr) \
\end{align}
Derivative of a Gamma function does not seem to yeild analytically tracktable solutions.
<hr/>
<p style="font-size
<ASSISTANT_TASK:>
Python Code:
## Add JS-based table of contents
from IPython.display import HTML as add_TOC
add_TOC( u<h1 id="tocheading">Table of Contents</h1></br><div id="toc"></div>
<script src="https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js"></script></br></hr></br> )
Explanation: <p style = "font-size: 24pt;text-align: center;"><strong>Expectation Minimization and MNIST</strong></p>
<p style = "font-size: 16pt;text-align: center;"><strong><i>Nazarov Ivan, 101мНОД (ИССА)</i></strong></p>
End of explanation
import os, time as tm, warnings
warnings.filterwarnings( "ignore" )
# from IPython.core.display import HTML
from IPython.display import display, HTML
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed( 569034853  )
## This is the correct way to use the random number generator,
##  since it allows finer control.
rand = np.random.RandomState( np.random.randint( 0x7FFFFFFF ) )
Explanation: Preamble
End of explanation
## A bunch of wrappers to match the task specifications
def posterior( x, clusters ) :
    pi = np.ones( clusters.shape[ 0 ], dtype = np.float ) / clusters.shape[ 0 ]
    q, ll = __posterior( x, theta = clusters, pi = pi )
    return q
## The likelihood is a byproduct of the E-step's minimization of Kullback-Leibler
def likelihood( x, clusters ) :
    pi = np.ones( clusters.shape[ 0 ], dtype = np.float ) / clusters.shape[ 0 ]
    q, ll = __posterior( x, theta = clusters, pi = pi )
    return np.sum( ll )
Explanation: EM and MNIST
The $\TeX$ markup used here uses the "align*" environment and thus should not be viewed though nbViewer.
Before proceeding, it seems pedagogically necessary (at least for myself) to revise the EM-slgorithm and show its "correctness", so to say.
A brief description of the EM algorithm
The EM algorithm seeks to maximize the likelihood by means of successive application of two steps: the E-step and the M-step.
For any probability measure $Q$ on the space of latent variables $Z$ with density $q$ the following holds:
\begin{align}
\log p(X|\Theta)
    &= \int q(Z) \log p(X|\Theta) dZ
     = \mathbb{E}q \log p(X|\Theta) \
   %% &= \Bigl[p(X,Z|\Theta) = p(Z|X,\Theta) p(X|\Theta) \Bigr] \
    &= \mathbb{E}{Z\sim q} \log \frac{p(X,Z|\Theta)}{p(Z|X\Theta)}
     = \mathbb{E}{Z\sim q} \log \frac{q(Z)}{p(Z|X,\Theta)}
     + \mathbb{E}{Z\sim q} \log \frac{p(X,Z|\Theta)}{q(Z)} \ 
    &= KL\bigl(q\|p(\cdot|X,\Theta)\bigr) + \mathcal{L}\bigl(q, \Theta\bigr)\,,
\end{align}  
since the Bayes theorem posits that $p(X,Z|\Theta) = p(Z|X,\Theta) p(X|\Theta)$. Call this equiation the "master equation".
Now note that since the Kullback-Leibler divergence is always non-negative, one has the following inequality:
$$\log p(X|\Theta) \geq \mathcal{L}\bigl(q, \Theta\bigr) \,.$$
Let's try to make the lower bound as large as possible by changing $\Theta$ and varying $q$. But first note that the 
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$.
Now at the optimal distributuion $q^_\Theta$ the master equation becomes
$$ \log p(X|\Theta)
= \mathcal{L}\bigl(q^\Theta, \Theta\bigr)
= \mathbb{E}{Z\sim q^_\Theta} \log \frac{p(X,Z|\Theta)}{q^(Z|\Theta)}
= \mathbb{E}{Z\sim q^\Theta} \log p(X,Z|\Theta) - \mathbb{E}{Z\sim q^\Theta} \log q^*(Z|\Theta) \,,
$$
for any $\Theta$. Thus the problem of log-likelihood maximization reduces to that of maximizing the sum of expectations on the right-hand side.
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$.
Hopefully using an iterative procedure which switches between the computation of $q^_\Theta$ and the maximization of $\Theta$ might be effective. Consider the folowing :
* 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}$;
* M*-step: considering $q_{i+1}$ as given, solve $\mathcal{L}(q_{i+1},\Theta) \to \mathop{\text{max}}\Theta$, where 
$$ \mathcal{L}(q,\Theta) = \mathbb{E}{Z\sim q} \log p(X,Z|\Theta) - \mathbb{E}_{Z\sim q} \log q(Z) \,.$$
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:
* 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}$.
Now, if the latent variables are mutually independent, then the optimal $q$ must be factorizable into marginal densities and:
\begin{align}
    KL\bigl(q\|p(\cdot|X,\Theta)\bigr)
    &= \mathbb{E}{Z\sim q} \log q(Z) - \sum_j \mathbb{E}{z_j\sim q_j} \log p(z_j|X,\Theta)\
    &= \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)
     = \sum_j KL\bigl(q_j\|p_j(|X,\Theta)\bigr) \,,
\end{align}
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).
Therefore the E-step could be reduced to a set of minimization problems with respect to one-dimensional density functions:
$$ q_j^* = \mathop{\text{argmin}}_{q_j}\,\, KL\bigl(q_j\|p_j(\cdot|X,\Theta)\bigr) \,, $$
since the Kulback-Leibler divergence in this case in additively separable.
Correctness
Recall that the master equation is an identity: for all densities $q$ on $Z$ and for all admissible parameters $\Theta$
$$ \log p(X|\Theta) = KL\bigl(q\|p(\cdot|X,\Theta)\bigr) + \mathcal{L}\bigl(q, \Theta\bigr) \,.$$
Hence if after the E-step the Kulback-Leibler divergence is reduced:
$$ KL\bigl(q'\|p(\cdot|X,\Theta)\bigr) \leq KL\bigl(q\|p(\cdot|X,\Theta)\bigr) \,,$$
then for the same set of parameters $\Theta$ one has
$$ \mathcal{L}(q,\Theta) \leq \mathcal{L}(q',\Theta) \,.$$
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:
$$ \log p(X|\Theta_i) = \mathcal{L}(q_{i+1},\Theta_i) \,.$$
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
$$ \mathcal{L}(q_{i+1},\Theta_i) \leq \mathcal{L}(q_{i+1},\Theta_{i+1}) \,.$$
Threfore the effect of a single complete round of EM on the log-likelihood itself is:
$$ \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}) \,,$$
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.
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
$$ \mathcal{L}(q_{i+1},\Theta_i) \leq \log p(X|\Theta_i) \,,$$
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.
Application of the EM to MNIST data
Each image is a random element in a discrete probability space $\Omega = {0,1}^{N\times M}$ with product-measure
$$ \mathbb{P}(\omega) = \prod_{i=1}^N\prod_{j=1}^M \theta_{ij}^{\omega_{ij}} (1-\theta_{ij})^{1-\omega_{ij}} \,,$$
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})$.
Let's apply the EM algorithm to this dataset. The proposed model is the following.
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:
$$p(\omega|\Theta)
= \sum_{k=1}^K \pi_k p_k(\omega|\theta_k)
= \sum_{k=1}^K \pi_k \prod_{i=1}^N \prod_{j=1}^M \theta_{kij}^{\omega_{ij}} (1-\theta_{kij})^{1-\omega_{ij}}$$
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$.
Suppose $X=(x_i){i=1}^n \in \Omega^n$ is the dataset. The log-likelihood is given by
$$ \log p(X|\Theta) = \sum{s=1}^n \log \sum_{k=1}^K \pi_k
\prod_{i=1}^N \prod_{j=1}^M \theta_{kij}^{x_{sij}} (1-\theta_{kij})^{1-x_{sij}} \,,$$
where $x_{sij}\in{0,1}$ -- is the value of the the $(i,j)$-th pixel at the $s$-th observation.
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
$$ \log p(X,Z|\Theta) = \sum{s=1}^n \log \prod_{k=1}^K \Bigl[ \pi_k 
\prod_{i=1}^N \prod_{j=1}^M \theta_{kij}^{x_{sij}} (1-\theta_{kij})^{1-x_{sij}} \Bigr]^{1_{z_s = k}} \,,$$
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).
The log-likelihood simplifies to
$$ \log p(X,Z|\Theta) = \sum_{s=1}^n \sum_{k=1}^K 1_{z_s = k} \Bigl( \log \pi_k + 
\sum_{i=1}^N \sum_{j=1}^M \bigl( x_{sij} \log \theta_{kij} + (1-x_{sij}) \log (1-\theta_{kij}) \bigr) \Bigr) \,,$$
and further into a more separable form
$$ \log p(X,Z|\Theta)
= \sum_{s=1}^n \sum_{k=1}^K 1_{z_s = k} \log \pi_k
+ \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}
+ \sum_{i=1}^N \sum_{j=1}^M (1-x_{sij}) \log (1-\theta_{kij}) \Bigr) \,.$$
The expected log-likelihood under $z_s\sim q_s$ with $\mathbb{P}(z_s=k|X) = q_{sk}$, is given by
$$ \mathbb{E}\log p(X,Z|\Theta)
= \sum_{s=1}^n \sum_{k=1}^K q_{sk} \log \pi_k
+ \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) \,.$$
Analytic solution: E-step
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}))$.
$$\hat{q}{sk}
= \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)}
\propto \pi_k \prod_{i=1}^N \prod_{j=1}^M \theta_{kij}^{x_{sij}} (1-\theta_{kij})^{1-x_{sij}}
$$
and
$$ q^(Z) = \prod_{s=1}^n q_{s z_s} $$
Note that the denominator is actually the log-likelihood of the data.
In order to improve numerical stability and avoid numerical underflow it is better to use the following procedure for computation of the conditional probability:
$$ l_{sk} = \sum_{i=1}^N \sum_{j=1}^M \log \bigl( \theta_{kij}^{x_{sij}} (1-\theta_{kij})^{1-x_{sij}} \bigr) \,,$$
set $l^s = \max_k l{sk}$ and compute the log-sum
$$ \hat{l}s = \log \sum{k=1}^K \text{exp}\Bigl{ ( l_{sk} - l^s ) + \log \pi_k \Bigr} \,,$$
and then compute the consitional distribution:
$$ \hat{q}{sk} = \text{exp}\Bigl{ l_{sk} + \log \pi_k - ( \hat{l}_s + l^*_s ) \Bigr} \,.$$
This seemingly redunant subctration and addition of $l^_s$ helps avoid underflow during the numerical exponentiation. After this sanitization the E*-step's optimal distribution would be numerically accurate.
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
$$ \hat{l}_s = l^s + \log (1+\epsilon) \,,$$
whence
$$ \hat{q}{sk} = \text{exp}\bigl{l_{sk} - \hat{l}s\bigr} = (1+\epsilon)^{-1} \cdot \text{exp}\bigl{l{sk} - l^_s \bigr} \,.$$
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$.
The variables in the code have the following dimensions:
* $\theta \in [0,1]^{K\times (N\times M)}$;
* $\pi \in [0,1]^{1\times K}$;
* $x \in {0,1}^{n\times (N\times M)}$;
* $z \in [0,1]^{n\times K}$.
Wrappers required for the assignment.
End of explanation
## Classifier
def classify( x, theta, pi = None ) :
    pi = pi if pi is not None else np.ones( theta.shape[ 0 ], dtype = np.float ) / theta.shape[ 0 ]
## Compute the posterior probabilities of the data
    q_sk, ll_s = __posterior( x, theta = theta, pi = pi )
## Classify according to max pasterior:
    c_s = np.argmax( q_sk, axis = 1 )
    return c_s, q_sk, ll_s
Explanation: Classify using the maximum aposteriori rule.
End of explanation
def __component_likelihood( x, theta ) :
## Unfortunately sometimes there can be negative machine zeros, which
##  spoil the log-likelihood computation by poisoning with NANs.
## That is why the theta array is restricted to [0,1].
    theta_clipped = np.clip( theta, 0.0, 1.0 )
## Iterate over classes
    ll_sk = np.zeros( ( x.shape[ 0 ], theta.shape[ 0 ] ), dtype = np.float )
## Make a binary mask of the data
    mask = x > 0
    for k in xrange( theta.shape[ 0 ] ) :
## Note that the power formulation is just a mathematically convenient way of
##  writing \theta if x=1 or (1-\theta) otherwise.
        ll_sk[ :, k ] = np.sum( np.where( mask,
            np.log( theta_clipped[ k ] ), np.log( 1 - theta_clipped[ k ] ) ), axis = ( 1, ) )
    return ll_sk
Explanation: A procedure to compute the log-likelihood of each observaton with respect to each mixture component. Used in the posterior computation.
End of explanation
## The core procedure for computing the conditional density of classes
def __posterior( x, theta, pi ) :
## Get the log-likelihoods of each observation in each mixture component.
    ll_sk = __component_likelihood( x, theta )
## Find the largest unnormalized probability.
    llstar_s = np.reshape( np.max( ll_sk, axis = ( 1, ) ), ( ll_sk.shape[ 0 ], 1 ) )
## Subtract the largest exponent
    ll_sk -= llstar_s
## In the rare case when the largest exponent is -Inf, force the differences
##  to zero. This effective treaks such observations as having unfiorm likelihood
##  across classes. This way the priors don't get masked by really small numbers. 
##  I could've used ``np.nan_to_num( ll_sk - llstar_s )'' but it actually copies
##  the ll_sk array.
    ll_sk[ np.isnan( ll_sk ) ] = 0.0
## Don't forget to add the log-prior probability (Numpy broadcasting applies!).
##  Adding priors before dealing with infinities would mask then and yield
##  incorrect estimates of the log-likelihoods!
    ll_sk += np.log( np.reshape( pi, ( 1, ll_sk.shape[ 1 ] ) ) )
## Compute the log-sum-exp of the individual log-likelihoods. Negative infinities
##  resolve to 0.0 while the largest exponent resolves to a one. This step cannot
##  produce NaNs
    ll_s = np.reshape( np.log( np.sum( np.exp( ll_sk ), axis = ( 1, ) ) ), ( ll_sk.shape[ 0 ], 1 ) )
## The sum-exp could never be anything lower than 1, since at least one
##  element of each row of ll_sk has to be lstar_s, whence the respective
##  difference should be zero and the exponent -- 1. Thus even if the
##  rest of the sum is close to machine zero, the logarithm would still
##  return 0.
## Normalise the likelihoods to get conditional probability, and compute
##  the sum of the log-denominator, which is the log-likelihood.
    return np.exp( ll_sk - ll_s ), ll_s + llstar_s
Explanation: The actual procedure for computing the E*-step: the conditional distribution and the log-likelihood scores.
End of explanation
## The E-step is simple: just compute the optimal parameters under
##  the current conditional distribution of the latent variables.
def __learn_clusters( x, z ) :
## The prior class probabilities
    pi = z.sum( axis = ( 0, ) )
## Pixel probabilities conditional on the calss
    theta = np.tensordot( z, x, ( 0, 0 ) ) / pi.reshape( ( pi.shape[ 0 ], 1 ) )
## Return: regularization should be done at **E**-step!
#     return np.clip( theta, 1.0/784, 1.0 - 1.0/784 ), pi / x.shape[ 0 ]
    return theta, pi / x.shape[ 0 ]
Explanation: Analytic solution: M-step
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
$$ \hat{\pi}k = \frac{\sum{s=1}^n q_{sk}}{n} \,.$$
For $\theta_{kij}$, $i=1,\ldots,N$, $j=1,\ldots,M$ and $k=1,\ldots,K$ the FOC is
$$ \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 \,,$$
whence
$$ \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 } \,.$$
This M-step procedure is implemented below.
End of explanation
## A wrapper for the above function
def learn_clusters( x, z ) :
    theta, pi = __learn_clusters( x, z )
## Just return theta: in the condtional model the pi are fixed.
    return theta
Explanation: A wrapper to match the assignment specifications.
End of explanation
## A wrapper for the core em algorithm below
def em_algorithm( x, K, maxiter, verbose = True, rel_eps = 1e-4, full = False ) :
## Initialize the model parameters with uniform [0.25,0.75] random numbers
    theta_1 = rand.uniform( size = ( K, x.shape[ 1 ] ) ) * 0.5 + 0.25
    pi_1 = None if not full else np.ones( K, dtype = np.float ) / K
## Run the em algorithm
    tick = tm.time( )
    ll, theta, pi, status = __em_algorithm( x, theta_1 = theta_1,
        pi_1 = pi_1, niter = maxiter, rel_eps = rel_eps, verbose = verbose )
    tock = tm.time( )
    print( "total %.3f, %.3f/iter" % ( ( tock - tick ), ( tock - tick ) / len( ll ), ) )
## Return the history of theta and the final log liklihood
    if verbose :
        if status[ 'status' ] != 0 :
            print "Convergence not achieved. %d" % ( status[ 'status' ], )
    if full :
        return ( theta, pi ), ll
    return theta, ll
Explanation: A it has been mentioned eariler, the EM algorithm switches between E and M steps until convergence.
End of explanation
## The core of the EM algorithm
def __em_algorithm( x, theta_1, pi_1 = None, niter = 1000, rel_eps = 1e-4, verbose = True ) :
## If we were supplied with an initial estimate of the prior distribution,
##  then assume the full model is needed.
    full_model = pi_1 is not None
## If the prior cluster probabilities are not supplied, assume uniform distribution.
    pi_1 = pi_1 if full_model else np.ones( theta_1.shape[ 0 ], dtype = np.float ) / theta_1.shape[ 0 ]
## Allocate the necessary space for the history of model estimates
    theta_hist, pi_hist = theta_1[ np.newaxis ].copy( ), pi_1[ np.newaxis ].copy( )
    ll_hist = np.asarray( [ -np.inf ], dtype = np.float )
## Set "old" estimated to zero. At this line the current estimates are in fact
##  the initially provided ones.
    theta_0, pi_0 = np.zeros_like( theta_1 ), np.zeros_like( pi_1 )
## Initialize the loop
    status, kiter, rel_theta, rel_pi, ll = -1, 0, np.nan, np.nan, -np.inf
    while kiter < niter :
## Dump the current estimators and other information.
        if verbose :
            print( "Iteration %d: avg. log-lik: %.3f, $\\Theta$ div. %.3f, $\\Pi$ div. %.3f" % (
                kiter, ll / x.shape[ 0 ], rel_theta, rel_pi ) )
            show_data( theta_1 - theta_0 if True else theta_1, n = theta_0.shape[ 0 ],
                n_col = min( 10, theta_0.shape[ 0 ] ), cmap = plt.cm.hot, interpolation = 'nearest' )
## The convergence criterion is the L^∞ norm of relative L^1 errors
        if max( rel_pi, rel_theta ) < rel_eps :
            status = 0
            break ;
## Overwrite the initial estimates
        theta_0, pi_0 = theta_1, pi_1
## E-step: call the core posterior function to get both the log-likelihood
##  and the estimate of the conditional distribution.
        z_1, ll_s = __posterior( x, theta_0, pi_0 )
## Sum the individual log-likelihoods of observations
        ll = ll_s.sum()
## M-step: compute the optimal parameters under the current estimate of the posterior
        theta_1, pi_1 = __learn_clusters( x, z_1 )
## Discard the computed estimate of pi if the model is discriminative (conditional likelihood).
        if not full_model :
            pi_1 = pi_0
## Record the current estimates to the history
        theta_hist = np.vstack( ( theta_hist, theta_1[np.newaxis] ) )
        pi_hist = np.vstack( ( pi_hist, pi_1[np.newaxis] ) )
        ll_hist = np.append( ll_hist, ll )
## Check for bad float numbers
        if not ( np.all( np.isfinite( theta_1 ) ) and np.all( np.isfinite( pi_1 ) ) ) :
            status= -2
            break ;
## Check convergence: L^1 relative error. If the relative margin is exactly
##  zero, then return NaNs. This makes teh loop exhaust all iterations, since
##  any comprison agains a NaN returns False.
        rel_theta = np.sum( np.abs( theta_1 - theta_0 ) / ( np.abs( theta_0 ) + rel_eps ) ) if rel_eps > 0 else np.nan
        rel_pi = np.sum( np.abs( pi_1 - pi_0 ) / ( np.abs( pi_0 ) + rel_eps ) ) if rel_eps > 0 else np.nan
## Next iteration
        kiter += 1
    return ll_hist, theta_hist, pi_hist, { 'status': status, 'iter': kiter }
Explanation: The procedure above actually invokes the true EM core, defined below.
End of explanation
def experiment( data, K, maxiter, verbose = True, until_convergence = False, full = False ) :
## Run the EM
    return em_algorithm( data, K, maxiter, rel_eps = 1.0e-4 if until_convergence else 0.0, verbose = verbose, full = full )
Explanation: Define a convenient procedure for running experiments. By setting relative error to zero the algorithm is forced to exhaust all the allocated iterations.
End of explanation
## A more flexible image arrangement
def arrange_flex( images, n_row = 10, n_col = 10, N = 28, M = 28, fill_value = 0 ) :
## Create the final grid of images row-by-row
    im_grid = np.full( ( n_row * N, n_col * M ), fill_value, dtype = images.dtype )
    for k in range( min( images.shape[ 0 ], n_col * n_row ) ) :
## Get the grid cell at which to place the image
        i, j = ( k // n_col ) * N, ( k % n_col ) * M
## Just put the image in the cell
        im_grid[ i:i+N, j:j+M ] = np.reshape( images[ k ], ( N, M, ) )
    return im_grid
Explanation: Miscellanea: visualization
In order to be able to plot more flexibly, define another arranger.
End of explanation
def setup_canvas( axis, n_row, n_col, N = 28, M = 28 ) :
## Setup major tick marks to the seam between images and disable their labels
    axis.set_yticks( np.arange( 1, n_row + 1 ) * N, minor = False )
    axis.set_xticks( np.arange( 1, n_col + 1 ) * M, minor = False )
    axis.set_yticklabels( [ ], minor = False ) ; axis.set_xticklabels( [ ], minor = False )
## Set minor ticks so that they are exactly between the major ones
    axis.set_yticks( ( np.arange( n_row + 1 ) + 0.5 ) * N, minor = True )
    axis.set_xticks( ( np.arange( n_col + 1 ) + 0.5 ) * M, minor = True )
## Make their labels into cell x-y coordinates
    axis.set_yticklabels( [ "%d" % (i,) for i in 1+np.arange( n_row + 1 ) ], minor = True )
    axis.set_xticklabels( [ "%d" % (i,) for i in 1+np.arange( n_col + 1 ) ], minor = True )
## Tick marks should be oriented outward
    axis.tick_params( axis = 'both', which = 'both', direction = 'out' )
## Return nothing!
    axis.grid( color = 'white', linestyle = '--' )
Explanation: The folowing pair of procedures are used to plot the digits in a clear manner. The first one just creates a canvas for the image: it sets up both axes properly and adds labels to them.
End of explanation
def show_data( data, n, n_col = 10, transpose = False, **kwargs ) :
## Get the number of rows necessary to plot the needed number of images
    n_row = ( n + n_col - 1 ) // n_col
## Transpose if necessary
    if transpose :
        n_col, n_row = n_row, n_col
## Set the dimensions of the figure
    fig = plt.figure( figsize = ( n_col, n_row ) )
    axis = fig.add_subplot( 111 )
## Plot!
    setup_canvas( axis, n_row, n_col )
    axis.imshow( arrange_flex( data[:n], n_col = n_col, n_row = n_row ), **kwargs )
## Plot
    plt.show( )
def visualize( data, clusters, ll, n_col = 2, plot_ll = True ) :
## Display the result
    print "Final conditional log-likelihood value per observation achieved %f in %d iteration(s)" % (
        ll[-1] / data.shape[ 0 ], len( ll ) )
## Plot the first difference of average log-likelihood
    if plot_ll :
        plt.figure( figsize = ( 12, 7 ) )
        ax = plt.subplot(111)
        ax.set_title( r"avg. log-likelihood change between successive iterations (log scale)" )
        ax.plot( np.diff( ll / data.shape[ 0 ]  ) )
#         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)$" )
        ax.set_ylabel( r"$\Delta_i \frac{1}{n} \sum_{s=1}^n \log p(x_s|\Theta_i)$" )
        ax.set_yscale( 'log' )
## Plot the final estimates
    if n_col > 0 :
        show_data( clusters[-1], n = clusters.shape[1], n_col = n_col, cmap = plt.cm.spectral, interpolation = 'nearest' )
Explanation: This procedure displays the images on a nice plot. Used for one-line visualization.
End of explanation
def animate( theta, ll, pi = None, n_col = 10, n_row = 10, interval = 1, **kwargs ) :
## Create a background
    bg = arrange_flex( np.zeros_like( theta[ 0 ] ), n_col = n_col, n_row = n_row )
## Compute log-likelihood differences and sanitize them.
    ll_diff = np.maximum( np.diff( ll ), np.finfo(np.float).eps )
    ll_diff[ ~np.isfinite( ll_diff ) ] = np.nan
## Set up the figure, the axis, and the plot elements we want to animate
    fig = plt.figure( figsize = ( 12, 12 ) )
## Create three subplots  and sposition them specifically
    if pi is None :
        ax1, ax3, ax2 = fig.add_subplot( 311 ), fig.add_subplot( 312 ), fig.add_subplot( 313 )
    else :
        ax1, ax4 = fig.add_subplot( 411 ), fig.add_subplot( 412 )
        ax3, ax2 = fig.add_subplot( 413 ), fig.add_subplot( 414 )
## Initialize different ranges for the image aritsts
    setup_canvas( ax1, n_row = n_row, n_col = n_col )
    ax1.set_title( r"Current estimate of the mixture components" )
    setup_canvas( ax2, n_row = n_row, n_col = n_col )
    ax2.set_title( r"Change between successive iterations" )
## Initialize geomtery for the delta log-likelihood plot.
    ax3.set_xlim( -0.1, ll.shape[ 0 ] + 0.1 )
    ax3.set_yscale( 'log' ) #; ax3.set_yticklabels( [ ] )
    ax3.set_title( r"Change between successive iterations of EM (log scale)" )
    ax3.set_ylabel( r"$\Delta_i \sum_{s=1}^n \log p(x_s|\Theta_i)$" )
    ax3.set_ylim( np.nanmin( ll_diff ) * 0.9, np.nanmax( ll_diff ) * 1.1 )
    ax3.grid( )
## Setup a plot for prior probabilites
    if pi is not None :
        classes = 1 + np.arange( len( pi[ 0 ] ) )
        ax4.set_xticks( classes )
        ax4.set_ylim( 0.0, 1.0 )
        ax4.set_title( r"Current estimate of the mixture weights" )
        ba1 = ax4.bar( classes, pi[ 0 ], align = "center" )
## Setup the artists
    im1 = ax1.imshow( bg, vmin = +0.0, vmax = +1.0, **kwargs )
    im2 = ax2.imshow( bg, vmin = -1.0, vmax = +1.0, **kwargs )
    line1, = ax3.plot( [ ], linestyle = "-", color = 'blue' )
## Animation function.  This is called sequentially
    def update( i ) :
## Compute the frame
        frame = theta[ i ] - theta[ i-1 ] if i > 0 else theta[ 0 ]
        frame /= np.max( np.abs( frame ) )
## Draw frames on the image artists
        im1.set_data( arrange_flex( theta[ i ], n_col = n_col, n_row = n_row ) )
        im2.set_data( arrange_flex( frame,      n_col = n_col, n_row = n_row ) )
        if i > 0 :
## Show history on the line artist
            line1.set_data( np.arange( i ), ll_diff[ :i ] )
        if pi is not None :
            [ b.set_height( h ) for b, h in zip( ba1, pi[ i ] ) ]
            if i > 0 :
                [ b.set_color( 'green' if h > p else 'red' ) for b, h, p in zip( ba1, pi[ i ], pi[ i-1 ] ) ]
## Return an iterator of artists in this frame
            return ( im1, im2, line1, ) + ba1
        return im1, im2, line1,
## Call the animator.
    return animation.FuncAnimation( fig, update, frames = theta.shape[ 0 ], interval = interval, blit = True )
Explanation: MIscellanea: animating the EM
This function creates an animation of successive iterations of a run of the EM.
End of explanation
## Make simple animations of the EM estimatora
## http://jakevdp.github.io/blog/2013/05/12/embedding-matplotlib-animations/
## http://jakevdp.github.io/blog/2012/08/18/matplotlib-animation-tutorial/
from matplotlib import animation
from IPython.display import HTML
from tempfile import NamedTemporaryFile
def embed_video( anim ) :
    VIDEO_TAG = <video controls autoplay muted loop><source src="data:video/x-m4v;base64,{0}"
        type="video/mp4">Your browser does not support the video tag.</video>
    plt.close( anim._fig )
    if not hasattr( anim, '_encoded_video' ) :
        ffmpeg_writer = animation.FFMpegWriter( )
        with NamedTemporaryFile( suffix = '.mp4' ) as f:
            anim.save( 'myanim.mp4', fps = 12, extra_args = [ '-vcodec', 'libx264' ] )# , writer = ffmpeg_writer )
            video = open( 'myanim.mp4', "rb" ).read( )
        anim._encoded_video = video.encode( "base64" )
    return HTML( VIDEO_TAG.format( anim._encoded_video ) )
Explanation: Define a function that produces (using ffmpeg) and embeds a video in HTML into IPython
End of explanation
if False :
## Fetch MNIST dataset from SciKit and create a local copy.
    from sklearn.datasets import fetch_mldata
    mnist = fetch_mldata( "MNIST original", data_home = './data/' )
    np.savez_compressed('./data/mnist/mnist_scikit.npz', data = mnist.data, labels = mnist.target )
Explanation: Miscellanea: obtaining the data
Try to download the MNIST data from the SciKit repository.
End of explanation
if False :
## The procedure below loads the MNIST data from a comma-separated text file.
    def load_mnist_from_csv( filename ) :
## Read the CSV file
        data = np.loadtxt( open( filename, "rb" ), dtype = np.short, delimiter = ",", skiprows = 0 )
## Peel off the lables
        return data[:,1:], data[:,0]
## Fetch the data from the provided CSV (!) files and save as a compressed data blob
    data, labels = load_mnist_from_csv( "./data/mnist/mnist_train.csv" )
    np.savez_compressed( './data/mnist/mnist_train.npz', labels = labels, data = data )
    data, labels = load_mnist_from_csv( "./data/mnist/mnist_test.csv" )
    np.savez_compressed( './data/mnist/mnist_test.npz', labels = labels, data = data )
Explanation: Or obtain the data from the provided CSV files.
End of explanation
assert( os.path.exists( './data/mnist/mnist_train.npz' ) )
with np.load( './data/mnist/mnist_train.npz', 'r' ) as npz :
    mnist_labels, mnist_data = npz[ 'labels' ], np.array( npz[ 'data' ] > 127, np.int )
    
assert( os.path.exists( './data/mnist/mnist_test.npz' ) )
with np.load( './data/mnist/mnist_test.npz', 'r' ) as npz :
    test_labels, test_data = npz[ 'labels' ], np.array( npz[ 'data' ] > 127, np.int )
Explanation: Study
First of all load and binarize the training data using the value 127 as the threshold.
End of explanation
## Mask
inx_sixes, inx_nines = np.where( mnist_labels == 6 )[ 0 ], np.where( mnist_labels == 9 )[ 0 ]
## Extract
sixes = mnist_data[ rand.choice( inx_sixes, 90, replace = False ) ]
nines = mnist_data[ rand.choice( inx_nines, 90, replace = False ) ]
## Show
show_data( sixes, n = 45, n_col = 15, cmap = plt.cm.gray, interpolation = 'nearest' )
show_data( nines, n = 45, n_col = 15, cmap = plt.cm.gray, interpolation = 'nearest' )
Explanation: Case : $K=2$
Let's have a look at some 6s and 9s.
End of explanation
data = mnist_data[ np.append( inx_sixes, inx_nines ) ]
clusters, ll = experiment( data, 2, 30 )
Explanation: They do indeed look quite distinct. Now collect them into a single dataset and estimate the model.
End of explanation
visualize( data, clusters, ll )
Explanation: 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.
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}$.
End of explanation
indices = np.arange( mnist_data.shape[ 0 ] )
rand.shuffle( indices )
show_data( mnist_data[ indices[:100] ] , n = 100, n_col = 10, cmap = plt.cm.gray, interpolation = 'nearest' )
Explanation: 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.
End of explanation
sub_sample = np.concatenate( tuple( [ rand.choice( np.where( mnist_labels == i )[ 0 ], size = 200 ) for i in range( 10 ) ] ) )
train_data, train_labels = mnist_data[ sub_sample ], mnist_labels[ sub_sample ]
# train_data, train_labels = mnist_data, mnist_labels
Explanation: Case : $K=10, 15, 20, 30, 60$ and $90$
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.
End of explanation
clusters_10, ll_10 = experiment( train_data, 10, 50 )
Explanation: 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.
End of explanation
visualize( train_data, clusters_10, ll_10, n_col = 10, plot_ll = True )
Explanation: 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.
End of explanation
clusters_15, ll_15 = experiment( train_data, 15, 50, verbose = False, until_convergence = False )
clusters_20, ll_20 = experiment( train_data, 20, 50, verbose = False, until_convergence = False )
clusters_30, ll_30 = experiment( train_data, 30, 50, verbose = False, until_convergence = False )
Explanation: Let's see if changing $K$ does the trick.
End of explanation
visualize( train_data, clusters_15, ll_15, n_col = 10, plot_ll = False )
visualize( train_data, clusters_20, ll_20, n_col = 10, plot_ll = False )
visualize( train_data, clusters_30, ll_30, n_col = 10, plot_ll = False )
Explanation: For what values of $K$ was it possible to infer the templates of all digits?
End of explanation
clusters_60, ll_60 = experiment( train_data, 60, 500, verbose = False, until_convergence = True )
Explanation: Obviously, the model with more mixture components is more likely to produce "templates" for all digits. For larger $K$ this is indeed the case.
Having run this algorithm for many times we are able to say that the digits $3$ and $8$, $4$ and $9$ and sometimes $5$ tend to be poorly separated. Furthermore due to there being many different handwritten variations of the same digit one should estimate a model with more classes. 
The returned templates of the mixture components are clearly suboptimal: the procedure seems to get stuck at individual examples. This may happen for any $K$ and allowing for more iterations does not remedy this.
Some possibilities do exist: add regularizers to the E step, that tie neighbouring pixel distributions together.
End of explanation
visualize( train_data, clusters_60, ll_60, n_col = 15 )
Explanation: As one can see, increasing the number of iterations does not necessarily improve the results.
End of explanation
## If you want to see this animation ensure that ffmpeg is installed and uncomment the following lines.
anim_60 = animate( clusters_60, ll_60, n_col = 15, n_row = 4,
    interval = 1, cmap = plt.cm.hot, interpolation = 'nearest' )
embed_video( anim_60 )
Explanation: Judging by the plot of the log-likelihood, the fact that the EM is guaranteed to converge to local maxima and does so extremely fast, there was no need for more than 120-130 iterations. The chages in the log-likelihood around that number of iterations are of the order $10^{-4}$. Since we are working in finite precision arithmetic (double), the smallest precision is $\approx 10^{-14}$.
Let's see the dynamics of the estimated of the EM iterations. You will have to ensure that ffMPEG is installed (Windows: and is on the PATH environment variable).
End of explanation
## Test model with K from 12 up to 42 with a step of 4
classes = 12 + np.arange( 11, dtype = np.int ) * 3
ll_hist = np.full( len( classes ), -np.inf, dtype = np.float )
## Store parameters
parameter_hist = list( )
for i, K in enumerate( classes ) :
## Run the experiment
    c, l = experiment( train_data, K, 50, verbose = False, until_convergence = False )
    ll_hist[ i ] = l[ -1 ]
    parameter_hist.append( c[ -1 ] )
## Visualize the final parameters
    show_data( c[-1], n = K, n_col = 13, cmap = plt.cm.hot, interpolation = 'nearest' )
Explanation: The parameter estimates of the EM stabilize pretty quickly. In fact most templates stabilize by iterations 100-120.
Choosing $K$
Among many methods of model selection, let's use simple training sample fittness score, givne by the value of the log-likelihood. Becasue the models are nested with respect to the number of mixture components, one should expect the likelihood to be a non decreasing function of $K$ (on average dut to randomization of the initial parameter values).
For large enough $K$ this method may lead to overfitting.
End of explanation
print classes[ np.argmax( ll_hist ) ]
Explanation: Indeed the log-likelihood does not decrease with $K$ on average. Nevertheless the model with the highes likelihood turs out to have this many mixture components:
End of explanation
# clusters = parameter_hist[ np.argmax( ll_hist ) ] * 0.999 + 0.0005
clusters = clusters_60[-1]
Explanation: A nice, yet expected coincidence :)
Classification: label assignment
Select a model ...
End of explanation
## Compute the posterior component probabilities, and use max-aposteriori
##  for the best class selection. 
c_s, q_sk, ll_s = classify( train_data, clusters )
Explanation: ... and get the posterior mixture component probabilities.
End of explanation
template_x_label_maj_60 = np.full( clusters.shape[ 0 ], -1, np.int )
for t in range( clusters.shape[ 0 ] ) :
    l, f = np.unique( train_labels[ c_s == t ], return_counts = True)
    if len( l ) > 0 :
## This is too blunt an approach: it does not guarantee surjectivity of the mapping.
        template_x_label_maj_60[ t ] = l[ np.argmax( f ) ]
Explanation: Use a simple majority rule to automatically assign lables to templates.
End of explanation
## there are 10 labels and K templates
label_cluster_score = np.full( ( clusters.shape[ 0 ], 10 ), -np.inf, np.float )
## Loop over each template
for t in range( clusters.shape[ 0 ] ) :
## The selected templates are chosen according to max-aposteriori rule.
    inx = np.where( c_s == t )[ 0 ]
## Get the assigned lables and their frequencies
    actual_labels = train_labels[ inx ]
    l, f = np.unique( actual_labels, return_counts = True )
## For each template and each associated label in the training set the
##  score is average of the top-5 highest log-likelihoods.
    label_cluster_score[ t, l ] = [ np.average( sorted(
        ll_s[ inx[ actual_labels == a ] ].flatten( ), reverse = True )[ : 5 ] ) for a in l ]
## For each template choose the label with the highes likelihood.
template_x_label_lik_60 = np.argmax( label_cluster_score, axis = 1 )
Explanation: Assign the labels $l$ to templates $t$ according to its score, based on the average of the top-$5$ log-likelihoods of observations with label $l$ and classfified with template $t$.
End of explanation
show_data( clusters, clusters.shape[ 0 ], 10,
          cmap = plt.cm.spectral, interpolation = 'nearest' )
Explanation: Compare the label assignments. Here are the templates.
End of explanation
mask = np.asarray( template_x_label_maj_60 != template_x_label_lik_60, dtype = np.float ).reshape( (-1,1) )
show_data( clusters * mask, clusters.shape[ 0 ], 10,
          cmap = plt.cm.spectral, interpolation = 'nearest' )
print "\nLikelihood based: ", template_x_label_lik_60[ mask[:,0] > 0 ]
print "Majority bassed: ", template_x_label_maj_60[ mask[:,0] > 0 ]
Explanation: These are the templates, which were assigned different labels by the majority and "trust" methods.
End of explanation
show_data( clusters[ np.argsort( template_x_label_lik_60 ) ], clusters.shape[ 0 ],
          10, cmap = plt.cm.spectral, interpolation = 'nearest' ) 
show_data( clusters[ np.argsort( template_x_label_maj_60 ) ], clusters.shape[ 0 ],
          10, cmap = plt.cm.spectral, interpolation = 'nearest' ) 
Explanation: Here are the pictures of templates ordered according to their label.
End of explanation
## Run the classifier on the test data
c_s_60, q_sk, ll_s = classify( test_data, clusters )
Explanation: Classification: test sample
Shall we try running the classifier on the test data?
End of explanation
## Show a sample of images and their templates
sample = np.random.permutation( test_data.shape[ 0 ] )[:64]
## Stack image and tis best template atop one another
display_stack = np.empty( ( 2 * len( sample ), test_data.shape[ 1 ] ), dtype = np.float )
display_stack[0::2] = test_data[ sample ] * q_sk[ sample, c_s_60[ sample ], np.newaxis ]
display_stack[1::2] = clusters[ c_s_60[ sample ] ]
## Display
show_data( display_stack, n = display_stack.shape[ 0 ], n_col = 16,
    transpose = False, cmap = plt.cm.spectral, interpolation = 'nearest' )
Explanation: Let's see the best template for each test observation in some sub-sample.
End of explanation
print "Accuracy of likelihood based labelling: %.2f" % (
    100 * np.average( template_x_label_lik_60[ c_s_60 ] == test_labels ), )
print "Accuracy of simple majority labelling: %.2f" % (
    100 * np.average( template_x_label_maj_60[ c_s_60 ] == test_labels ), )
Explanation: The digits are shown in pairs: each first digit is the test observation (colour is determined by the confidence of the classifier -- the whiter the higher), and each second -- is the best template.
Let's see how accurate the classification was. Recall that the component assignment was done using 
$$ \hat{t}_s = \mathop{\text{argmax}}_k p\bigl( C_s = k \, \big\vert\, X = x_s\bigr) \,, $$
i.e. maximum aposteriori rule.
Then, with the classes assigned, labels were deduced based on either :
* simple majority;
* observations with top largest likelihoods in each class.
By accuracy I understand the following socre:
$$ \alpha = 1 - \frac{1}{|\text{TEST}|} \sum_{s\,\in\,\text{TEST}} \mathbf{1}{ l_s \neq L(\hat{t}_s)}\,, $$
where $l_s$ -- is the actual lablesof an observation $s$, $\hat{t}_s$ -- is the inferred mixture component (class) of that observation, and $k\mapsto L(k)$ is the component to label mapping.
End of explanation
## Now display the test observations, which the model cold nor classify at all.
bad_tests = np.where( np.isinf( ll_s ) )[ 0 ]
show_data( test_data[ bad_tests ], n = max( len( bad_tests ), 10 ), n_col = 15, cmap = plt.cm.gray, interpolation = 'nearest' )
# print q_sk[ bad_tests ]
Explanation: Not surprisingly, majority- and likelihood-based classification accuracies are close.
Let's see which test observations the model considers an artefact and for which it cannot reliably assign a template: i.e. the posterior class probablity for these cases coincides with the prior. This happens when the likelihood of an observation is identical within each class.
End of explanation
clusters = clusters_30[-1]
c_s, q_sk, ll_s = classify( train_data, clusters )
template_x_label_maj_30 = np.full( clusters.shape[ 0 ], -1, np.int )
for t in range( clusters.shape[ 0 ] ) :
    l, f = np.unique( train_labels[ c_s == t ], return_counts = True)
    if len( l ) > 0 :
        template_x_label_maj_30[ t ] = l[ np.argmax( f ) ]
label_cluster_score_30 = np.full( ( clusters.shape[ 0 ], 10 ), -np.inf, np.float )
for t in range( clusters.shape[ 0 ] ) :
    inx = np.where( c_s == t )[ 0 ]
    actual_labels = train_labels[ inx ]
    l, f = np.unique( actual_labels, return_counts = True )
    label_cluster_score_30[ t, l ] = [ np.average( sorted(
        ll_s[ inx[ actual_labels == a ] ].flatten( ), reverse = True )[ : 5 ] ) for a in l ]
template_x_label_lik_30 = np.argmax( label_cluster_score_30, axis = 1 )
show_data( clusters[ np.argsort( template_x_label_lik_30 ) ], clusters.shape[ 0 ],
          10, cmap = plt.cm.spectral, interpolation = 'nearest' )
print template_x_label_lik_30[ np.argsort( template_x_label_lik_30 ) ].reshape((3,-1))
c_s_30, q_sk, ll_s = classify( test_data, clusters )
print "Accuracy of likelihood based labelling: %.2f" % (
    100 * np.average( template_x_label_lik_30[ c_s_30 ] == test_labels ), )
print "Accuracy of simple majority labelling: %.2f" % (
    100 * np.average( template_x_label_maj_30[ c_s_30 ] == test_labels ), )
Explanation: </hr>
Let's see how more pasrimonious models fare with respect to accuracy on thte test sample.
Accuracy of $K=30$
End of explanation
clusters = clusters_20[-1]
c_s, q_sk, ll_s = classify( train_data, clusters )
template_x_label_maj_20 = np.full( clusters.shape[ 0 ], -1, np.int )
for t in range( clusters.shape[ 0 ] ) :
    l, f = np.unique( train_labels[ c_s == t ], return_counts = True)
    if len( l ) > 0 :
        template_x_label_maj_20[ t ] = l[ np.argmax( f ) ]
label_cluster_score_20 = np.full( ( clusters.shape[ 0 ], 10 ), -np.inf, np.float )
for t in range( clusters.shape[ 0 ] ) :
    inx = np.where( c_s == t )[ 0 ]
    actual_labels = train_labels[ inx ]
    l, f = np.unique( actual_labels, return_counts = True )
    label_cluster_score_20[ t, l ] = [ np.average( sorted(
        ll_s[ inx[ actual_labels == a ] ].flatten( ), reverse = True )[ : 5 ] ) for a in l ]
template_x_label_lik_20 = np.argmax( label_cluster_score_20, axis = 1 )
show_data( clusters[ np.argsort( template_x_label_lik_20 ) ], clusters.shape[ 0 ],
          10, cmap = plt.cm.spectral, interpolation = 'nearest' )
print template_x_label_lik_20[ np.argsort( template_x_label_lik_20 ) ].reshape((2,-1))
c_s_20, q_sk, ll_s = classify( test_data, clusters )
print "Accuracy of likelihood based labelling: %.2f" % (
    100 * np.average( template_x_label_lik_20[ c_s_20 ] == test_labels ), )
print "Accuracy of simple majority labelling: %.2f" % (
    100 * np.average( template_x_label_maj_20[ c_s_20 ] == test_labels ), )
Explanation: Accuracy of $K=20$
End of explanation
clusters = clusters_15[-1]
c_s, q_sk, ll_s = classify( train_data, clusters )
template_x_label_maj_15 = np.full( clusters.shape[ 0 ], -1, np.int )
for t in range( clusters.shape[ 0 ] ) :
    l, f = np.unique( train_labels[ c_s == t ], return_counts = True)
    if len( l ) > 0 :
        template_x_label_maj_15[ t ] = l[ np.argmax( f ) ]
label_cluster_score_15 = np.full( ( clusters.shape[ 0 ], 10 ), -np.inf, np.float )
for t in range( clusters.shape[ 0 ] ) :
    inx = np.where( c_s == t )[ 0 ]
    actual_labels = train_labels[ inx ]
    l, f = np.unique( actual_labels, return_counts = True )
    label_cluster_score_15[ t, l ] = [ np.average( sorted(
        ll_s[ inx[ actual_labels == a ] ].flatten( ), reverse = True )[ : 5 ] ) for a in l ]
template_x_label_lik_15 = np.argmax( label_cluster_score_15, axis = 1 )
show_data( clusters[ np.argsort( template_x_label_lik_15 ) ], clusters.shape[ 0 ],
          15, cmap = plt.cm.spectral, interpolation = 'nearest' )
print template_x_label_lik_15[ np.argsort( template_x_label_lik_15 ) ].reshape((1,-1))
c_s_15, q_sk, ll_s = classify( test_data, clusters )
print "Accuracy of likelihood based labelling: %.2f" % (
    100 * np.average( template_x_label_lik_15[ c_s_15 ] == test_labels ), )
print "Accuracy of simple majority labelling: %.2f" % (
    100 * np.average( template_x_label_maj_15[ c_s_15 ] == test_labels ), )
Explanation: Accuracy of $K=15$
End of explanation
clusters = clusters_10[-1]
c_s, q_sk, ll_s = classify( train_data, clusters )
template_x_label_maj_10 = np.full( clusters.shape[ 0 ], -1, np.int )
for t in range( clusters.shape[ 0 ] ) :
    l, f = np.unique( train_labels[ c_s == t ], return_counts = True)
    if len( l ) > 0 :
        template_x_label_maj_10[ t ] = l[ np.argmax( f ) ]
label_cluster_score_10 = np.full( ( clusters.shape[ 0 ], 10 ), -np.inf, np.float )
for t in range( clusters.shape[ 0 ] ) :
    inx = np.where( c_s == t )[ 0 ]
    actual_labels = train_labels[ inx ]
    l, f = np.unique( actual_labels, return_counts = True )
    label_cluster_score_10[ t, l ] = [ np.average( sorted(
        ll_s[ inx[ actual_labels == a ] ].flatten( ), reverse = True )[ : 5 ] ) for a in l ]
template_x_label_lik_10 = np.argmax( label_cluster_score_10, axis = 1 )
show_data( clusters[ np.argsort( template_x_label_lik_10 ) ], clusters.shape[ 0 ],
          10, cmap = plt.cm.spectral, interpolation = 'nearest' )
print template_x_label_lik_10[ np.argsort( template_x_label_lik_10 ) ].reshape((1,-1))
c_s_10, q_sk, ll_s = classify( test_data, clusters )
print "Accuracy of likelihood based labelling: %.2f" % (
    100 * np.average( template_x_label_lik_10[ c_s_10 ] == test_labels ), )
print "Accuracy of simple majority labelling: %.2f" % (
    100 * np.average( template_x_label_maj_10[ c_s_10 ] == test_labels ), )
Explanation: Accuracy of $K=10$
End of explanation
print "Model with K = 10: %.2f" % ( 100 * np.average( template_x_label_lik_10[ c_s_10 ] == test_labels ), )
print "Model with K = 15: %.2f" % ( 100 * np.average( template_x_label_lik_15[ c_s_15 ] == test_labels ), )
print "Model with K = 20: %.2f" % ( 100 * np.average( template_x_label_lik_20[ c_s_20 ] == test_labels ), )
print "Model with K = 30: %.2f" % ( 100 * np.average( template_x_label_lik_30[ c_s_30 ] == test_labels ), )
print "Model with K = 60: %.2f" % ( 100 * np.average( template_x_label_lik_60[ c_s_60 ] == test_labels ), )
Explanation: As one can see test sample accuracy of the model falls drammatically for less number of mixture components. This was expected, since due to various reasons, one being thet the data is handwritten, it is higly unlikely, that a single digit would have only one template.
End of explanation
( clusters_full, pi_full ), ll_full = experiment( data, 30, 1000, False, True, True )
anim_full = animate( clusters_full, ll_full, pi = pi_full, n_col = 15, n_row = 2, interval = 1, cmap = plt.cm.hot, interpolation = 'nearest' )
embed_video( anim_full )
Explanation: <br/><p style="font-size: 20pt;font-weight: bold; text-align: center;font-family: Courier New"> Ignore everything below </p><br/>
Let digress for a moment and consider the full model and create a video to see how the estimates are refined.
End of explanation
from sklearn.neighbors import KernelDensity
from sklearn.decomposition import PCA
from sklearn.grid_search import GridSearchCV
pca = PCA( n_components = 50 )
X_train_pca = pca.fit_transform( X_train )
params = { 'bandwidth' : np.logspace( -1, 1, 20 ) }
grid = GridSearchCV( KernelDensity( ), params )
grid.fit( X_train_pca )
print("best bandwidth: {0}".format( grid.best_estimator_.bandwidth ) )
params
kde = grid.best_estimator_
new_data = kde.sample( 100 )
new_data = pca.inverse_transform( new_data )
print new_data.shape
plt.figure( figsize = ( 9, 9 ) )
plt.imshow( arrange_flex( new_data ), cmap = plt.cm.gray, interpolation = 'nearest' )
plt.show( )
Explanation: <hr/>
A random variable $X\sim \text{Beta}(\alpha,\beta)$ if the law of $X$ has density
$$p(u) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} u^{\alpha-1}(1-u)^{\beta-1} $$
$$ \log p(X,Z|\Theta) = \sum_{s=1}^n \log \prod_{k=1}^K \Bigl[ \pi_k 
\prod_{i=1}^N \prod_{j=1}^M 
\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}}$$
\begin{align}
\mathbb{E}q \log p(X,Z|\Theta)
&= \sum{k=1}^K \sum_{s=1}^n q_{sk} \log \pi_k \
&+ \sum_{k=1}^K \sum_{i=1}^N \sum_{j=1}^M \sum_{s=1}^n q_{sk} \bigl(
    \log \Gamma(\alpha_{kij}+\beta_{kij}) - \log \Gamma(\alpha_{kij}) - \log \Gamma(\beta_{kij}) \bigr) \
&+ \sum_{k=1}^K \sum_{i=1}^N \sum_{j=1}^M \sum_{s=1}^n q_{sk} \bigl(
    (\alpha_{kij}-1) \log x_{sij} + (\beta_{kij}-1) \log(1-x_{sij}) \bigr) \
\end{align}
Derivative of a Gamma function does not seem to yeild analytically tracktable solutions.
<hr/>
<p style="font-size: 20pt; text-align: center;font-family: Courier New">Non-parametric approach</p>
End of explanation
<END_TASK>