{
 "metadata": {
  "name": "",
  "signature": "sha256:8daf891c3a257420768717dc0e3c4dd8d262a1adf5bbcff1fbf559176a325f69"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Gibbs sampling for the beta-bernoulli mixture model"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Bayesian beta-bernoulli (fixed size) mixture model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Given a fixed $K$ and $D$, the beta-bernoulli mixture model is a generative model describing $D$-dimensional binary vectors ($y_i \\in \\{0,1\\}^{D}$) drawn from a $K$-component mixture.\n",
      "\n",
      "We can describe the probabilistic model as:\n",
      "\n",
      "\\begin{align*}\n",
      "  \\pi|\\alpha &\\sim \\text{Dirichlet}(\\{ \\frac{\\alpha}{K} \\}_{k=1}^{K}) \\\\\n",
      "  c_i|\\pi &\\sim \\text{Discrete}(\\pi) \\\\\n",
      "  p_d^{k} | \\beta,\\gamma &\\sim \\text{Beta}(\\beta, \\gamma) \\\\\n",
      "  y_i^{d} | c_i{=}k,p_d^{k} &\\sim \\text{Bernoulli}(p_d^{k})\n",
      "\\end{align*}\n",
      "\n",
      "Let us clear up some of this notation. $\\pi$ is a $K$-dimensional vector living in the $(K-1)$-dimensional probability simplex. For each $k \\in [K]$, $\\{p_d^k\\}_{d=1}^{D} \\in [0,1]^{D}$. Given $N$ data points, the assignment vector $\\{c_i\\}_{i=1}^{N} \\in \\{0, ..., K-1\\}^{N}$. The hyperparameters of this model are $\\mathcal{H} = (\\alpha, \\beta, \\gamma)$.\n",
      "\n",
      "As we will see later on, this model has nice analytical properties for Gibbs sampling, since the beta distribution is a nice conjugate prior for the bernoulli distribution.\n",
      "\n",
      "The inference problem we will consider is the following. Given a dataset $\\mathcal{Y} = \\{y_i\\}_{i=1}^{N}$, we want to learn the posterior distribution on the assignment vector $\\mathcal{C} = \\{c_i\\}_{i=1}^{N}$. That is, we want to be able to estimate and draw samples from $p(\\mathcal{C} | \\mathcal{Y}; \\mathcal{H})$."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Gibbs sampling for estimating the posterior"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This problem can be solved efficiently with [Gibbs sampling](http://en.wikipedia.org/wiki/Gibbs_sampling), which is another [Markov-chain monte carlo](http://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo) (MCMC) variant. Let the notation $\\mathcal{C}_{\\neg i} = \\{ c_j \\in \\mathcal{C} : j \\neq i \\}$. The Gibbs sampler assumes we can efficiently sample from the distribution $p(c_i | \\mathcal{C}_{\\neg i}, \\mathcal{Y}; \\mathcal{H})$. Every iteration of the sampling then works by\n",
      "sampling for each $i \\in [N]$,\n",
      "\\begin{align*}\n",
      "    c^{(t)}_{i} \\gets p(c_i | \\{ c_j \\in \\mathcal{C}^{(t)} : j < i \\}, \\{ c_j \\in \\mathcal{C}^{(t-1)} : j > i \\}, \\mathcal{Y}; \\mathcal{H})\n",
      "\\end{align*}\n",
      "\n",
      "Let the notation $\\mathcal{Y}^{k}_{\\neg i} = \\{ y_j \\in \\mathcal{Y} : c_j = k, j \\neq i \\}$.\n",
      "It turns out, for the beta-bernoulli model, we can indeed [derive](http://homepage.tudelft.nl/19j49/Publications_files/TR_1.pdf) that\n",
      "\\begin{align*}\n",
      "  p(c_i{=}k|\\mathcal{C}_{\\neg i}, \\mathcal{Y};\\mathcal{H}) \\propto\n",
      "    \\frac{ |\\mathcal{Y}^{k}_{\\neg i}| + \\frac{\\alpha}{K} }{ N - 1 + \\alpha } \\prod_{d=1}^{D} \\frac{(\\beta+\\sum_{y_k\\in \\mathcal{Y}^{k}_{\\neg i}} y_k^{d})^{y_i^d} (  \\gamma + |\\mathcal{Y}^{k}_{\\neg i}| - \\sum_{y_k\\in \\mathcal{Y}^{k}_{\\neg i}} y_k^{d})^{(1-y_i^{d})} }{ \\beta + \\gamma + |\\mathcal{Y}^{k}_{\\neg i}| }\n",
      "\\end{align*}\n",
      "We can construct $p(c_i{=}k|\\mathcal{C}_{\\neg i}, \\mathcal{Y};\\mathcal{H})$ exactly by then enumerating through all $K$ clusters and normalizing."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Implementation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We'll first implement a sampler from a discrete distribution. Note while `scipy.stats.rv_discrete` [(Docs)](http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_discrete.html) implements this functionality, it is quite heavy weight to call within a loop, so we roll our own."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# workspace setup\n",
      "import numpy as np\n",
      "import scipy as sp\n",
      "import scipy.io\n",
      "import scipy.misc\n",
      "import scipy.special\n",
      "\n",
      "%matplotlib inline\n",
      "import matplotlib.pylab as plt\n",
      "\n",
      "import itertools as it\n",
      "import math"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def discrete_sample(pmf):\n",
      "    coin = np.random.random()\n",
      "    acc = 0.0\n",
      "    n, = pmf.shape\n",
      "    for i in xrange(n):\n",
      "        acc0 = acc + pmf[i]\n",
      "        if coin <= acc0:\n",
      "            return i\n",
      "        acc = acc0\n",
      "    return n-1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Ok enough playing around! Now the actual Gibbs sampler implementation. Note we make use of [numpy broadcasting](http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html) for performance:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def gibbs_beta_bernoulli(Y, K, alpha, beta, gamma, niters):\n",
      "    N, D = Y.shape\n",
      "    alpha, beta, gamma = map(float, [alpha, beta, gamma])\n",
      "\n",
      "    # start with random assignment\n",
      "    assignments = np.random.randint(0, K, size=N)\n",
      "\n",
      "    # initialize the sufficient statistics (cluster sums) accordingly\n",
      "    sums = np.zeros((D, K), dtype=np.int64)\n",
      "    cnts = np.zeros(K, dtype=np.int64)\n",
      "    for yi, ci in zip(Y, assignments):\n",
      "        sums[:,ci] += yi\n",
      "        cnts[ci] += 1\n",
      "\n",
      "    history = np.zeros((niters, N), dtype=np.int64)\n",
      "\n",
      "    # precomputations\n",
      "    nplog = np.log\n",
      "    npexp = np.exp\n",
      "    nparray = np.array\n",
      "    logsumexp = sp.misc.logsumexp\n",
      "    lg_denom = nplog(N - 1 + alpha)\n",
      "    alpha_over_K = alpha/K\n",
      "    beta_plus_gamma = beta + gamma\n",
      "\n",
      "    for t in xrange(niters):\n",
      "        for i, (yi, ci) in enumerate(zip(Y, assignments)):\n",
      "            # remove from SS\n",
      "            sums[:,ci] -= yi\n",
      "            cnts[ci] -= 1\n",
      "\n",
      "            lg_term1 = nplog(cnts + alpha_over_K) - lg_denom - D*nplog(beta_plus_gamma + cnts)\n",
      "            lg_term2 = nplog(beta + sums)\n",
      "            lg_term3 = nplog(gamma + cnts - sums)\n",
      "\n",
      "            lg_dist = lg_term1 + (lg_term2*yi[:,np.newaxis] + lg_term3*(1-yi[:,np.newaxis])).sum(axis=0)\n",
      "            lg_dist -= logsumexp(lg_dist) # normalize\n",
      "\n",
      "            # reassign\n",
      "            ci = discrete_sample(npexp(lg_dist))\n",
      "            assignments[i] = ci\n",
      "            sums[:,ci] += yi\n",
      "            cnts[ci] += 1\n",
      "        history[t] = assignments\n",
      "\n",
      "    return history"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Testing"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us first generate a small toy dataset (using the beta-bernoulli model as the generative process). Then we will test our Gibbs sampler on this small dataset."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "alpha, beta, gamma = 2., 1., 1.\n",
      "K = 2\n",
      "D = 3\n",
      "N = 4\n",
      "\n",
      "pis = np.random.dirichlet(alpha/K*np.ones(K))\n",
      "cis = np.array([discrete_sample(pis) for _ in xrange(N)])\n",
      "aks = np.random.beta(beta, gamma, size=(K, D))\n",
      "\n",
      "print 'Pi:', pis\n",
      "print 'C :', cis\n",
      "\n",
      "def bernoulli(p):\n",
      "    return 1 if np.random.random() <= p else 0\n",
      "\n",
      "Y = np.zeros((N, D), dtype=np.int64)\n",
      "for i in xrange(N):\n",
      "    Y[i] = np.array([bernoulli(aks[cis[i], d]) for d in xrange(D)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Pi: [ 0.17704708  0.82295292]\n",
        "C : [1 1 0 1]\n"
       ]
      }
     ],
     "prompt_number": 56
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "niters = 50000\n",
      "chain = gibbs_beta_bernoulli(Y, K, alpha, beta, gamma, niters)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 57
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Posterior distribution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now the natural question is, how do we verify that the Gibbs sampler is indeed drawing samples from $p(\\mathcal{C} | \\mathcal{Y})$. This is a bit trickier in the Bayesian case than the non-Bayesian case, since what we are really after is a *distribution* instead of, e.g. a point estimate.\n",
      "\n",
      "Luckily, for small problems (small $K$ and $N$), we can actually calculate $p(\\mathcal{C} | \\mathcal{Y}; \\mathcal{H})$ exactly by brute force enumeration. This is because $p(\\mathcal{C} | \\mathcal{Y} ; \\mathcal{H})$ is a discrete distribution of size $K^{N}$, and we can analytically calculate the joint distribution $p(\\mathcal{C}, \\mathcal{Y}; \\mathcal{H})$; from the joint, the posterior follows by $p(\\mathcal{C} | \\mathcal{Y} ) = \\frac{p(\\mathcal{C}, \\mathcal{Y})}{\\sum_{c} p(c ,\\mathcal{Y})}$, where the summation in the denominator is over all possible $K^N$ assignment vectors (from here on, we drop the hyperparameter dependence to ease notation).\n",
      "\n",
      "For an arbitrary $\\mathcal{C}, \\mathcal{Y}$, we have $p(\\mathcal{C}, \\mathcal{Y}) = p(\\mathcal{C}) p(\\mathcal{Y} | \\mathcal{C})$. From [page 2](http://homepage.tudelft.nl/19j49/Publications_files/TR_1.pdf), we have \n",
      "\\begin{align*}\n",
      "  p(\\mathcal{C}) = \\frac{\\Gamma(\\alpha)}{\\Gamma(|\\mathcal{Y}|+\\alpha)}\n",
      "      \\prod_{k=1}^{K} \\frac{\\Gamma( |\\mathcal{Y}^{k}| + \\frac{\\alpha}{K})}{ \\Gamma(\\frac{\\alpha}{K})}\n",
      "\\end{align*}\n",
      "where $\\Gamma(\\cdot)$ is the [Gamma](http://en.wikipedia.org/wiki/Gamma_function) function.\n",
      "We can derive $p(\\mathcal{Y}|\\mathcal{C})$ as follows.\n",
      "\\begin{align*}\n",
      "  p(\\mathcal{Y} | \\mathcal{C}) &= \n",
      "    \\prod_{k=1}^{K} p(\\mathcal{Y}^{k} | \\mathcal{C}^{k}) \\\\\n",
      "    &= \\prod_{k=1}^{K} \\int_{\\Theta_k} [\\prod_{y_i \\in \\mathcal{Y}^k} p(y_i | \\Theta_k) ] p(\\Theta_k ; \\mathcal{H}) \\; d\\Theta_k \n",
      "\\end{align*}\n",
      "Now for each $k$, we evaluate the inner integral as:\n",
      "\\begin{align*}\n",
      "  \\int_{\\theta_1,...,\\theta_D} \\prod_{d=1}^{D} \\left(\\prod_{y_i \\in \\mathcal{Y}^k} \\theta_d^{y_i^d} (1-\\theta_d)^{1-y_i^d}\\right) \\frac{1}{B(\\beta,\\gamma)} \\theta_d^{\\beta-1}(1-\\theta_d)^{\\gamma-1} \\; d\\theta_1...\\theta_{D}\n",
      "\\end{align*}\n",
      "where $B(\\beta,\\gamma)$ is the [Beta function](http://en.wikipedia.org/wiki/Beta_function). Noting that $\\int_{0}^{1} x^{(m-1)} (1-x)^{(n-1)} \\; dx = \\frac{\\Gamma(m)\\Gamma(n)}{\\Gamma(m+n)}$, we can simplify the above integral to:\n",
      "\\begin{align*}\n",
      "\\frac{1}{B(\\beta,\\gamma)^{D}} \\prod_{d=1}^{D} \\frac{\\Gamma( \\sum_{y_i \\in \\mathcal{Y}^k} y_i^d + \\beta) \\Gamma( |\\mathcal{Y}^{k}| - \\sum_{y_i \\in \\mathcal{Y}^k} y_i^d + \\gamma)}{ \\Gamma( |\\mathcal{Y}^{k}| + \\beta + \\gamma )}\n",
      "\\end{align*}\n",
      "And therefore,\n",
      "\\begin{align*}\n",
      "  p(\\mathcal{Y} | \\mathcal{C}) = \\frac{1}{B(\\beta,\\gamma)^{KD}} \\prod_{k=1}^{K} \\prod_{d=1}^{D} \\frac{\\Gamma( \\sum_{y_i \\in \\mathcal{Y}^k} y_i^d + \\beta) \\Gamma( |\\mathcal{Y}^{k}| - \\sum_{y_i \\in \\mathcal{Y}^k} y_i^d + \\gamma)}{ \\Gamma( |\\mathcal{Y}^{k}| + \\beta + \\gamma )}\n",
      "\\end{align*}\n",
      "\n",
      "Therefore, we can compute the posterior distribution of the data as follows:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def all_assignment_vectors():\n",
      "    return it.product(range(K), repeat=N)\n",
      "\n",
      "def lg_pr_joint(C, Y, K, alpha, beta, gamma):\n",
      "    N, D = Y.shape\n",
      "    nks = np.bincount(C, minlength=K)\n",
      "    assert nks.shape[0] == K\n",
      "    assert C.shape[0] == N\n",
      "\n",
      "    # log P(C)\n",
      "    gammaln = sp.special.gammaln\n",
      "    betaln = sp.special.betaln\n",
      "    term1 = gammaln(alpha) - gammaln(N + alpha) - K*gammaln(alpha/K)\n",
      "    term2 = sum(gammaln(nk + alpha/K) for nk in nks)\n",
      "    lg_pC = term1 + term2\n",
      "\n",
      "    # log P(Y|C)\n",
      "    term1 = K*D*betaln(beta, gamma)\n",
      "    term2 = D*sum(gammaln(nk + beta + gamma) for nk in nks)\n",
      "    sums = np.zeros((K, D))\n",
      "    for yi, ci in zip(Y, C):\n",
      "        sums[ci] += yi\n",
      "    def fn1(nk, sum_yid):\n",
      "        assert nk >= sum_yid\n",
      "        return gammaln(sum_yid + beta) + gammaln(nk - sum_yid + gamma)\n",
      "    term3 = sum(sum(fn1(nk, yid) for yid in row) for nk, row in zip(nks, sums))\n",
      "    lg_pYgC = -term1 - term2 + term3\n",
      "    \n",
      "    return lg_pC + lg_pYgC\n",
      "\n",
      "def brute_force_posterior(Y, K, alpha, beta, gamma):\n",
      "    N, _ = Y.shape\n",
      "\n",
      "    # enumerate K^N cluster assignments\n",
      "    lg_pis = np.array([lg_pr_joint(np.array(C), Y, K, alpha, beta, gamma) for C in all_assignment_vectors()])\n",
      "    lg_pis -= sp.misc.logsumexp(lg_pis)\n",
      "    \n",
      "    return np.exp(lg_pis)\n",
      "\n",
      "# generate an ID for each K^N element\n",
      "idmap = { C : i for i, C in enumerate(all_assignment_vectors()) }\n",
      "revidmap = { i : C for i, C in enumerate(all_assignment_vectors()) }\n",
      "\n",
      "actual_posterior = brute_force_posterior(Y, K, alpha, beta, gamma)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 58
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's pause for a second and take a look at the [MAP](http://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation) estimator:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'P(C=actual|Y):', actual_posterior[idmap[tuple(cis)]]\n",
      "print 'max_C P(C|Y):', actual_posterior.max()\n",
      "print 'argmax_C P(C|Y):', revidmap[actual_posterior.argmax()]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "P(C=actual|Y): 0.0864616067769\n",
        "max_C P(C|Y): 0.132805028009\n",
        "argmax_C P(C|Y): (0, 0, 0, 0)\n"
       ]
      }
     ],
     "prompt_number": 59
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice how the ground truth cluster assignment is *not* the MAP estimator!\n",
      "\n",
      "To measure the distance between the actual posterior distribution and that produced by our Gibbs sampler, we compare the KL-divergence of the actual and the empirical. Recall for discrete distributions the KL-divergence (relative entropy) is defined as\n",
      "\\begin{align*}\n",
      "  D(P||Q) = \\sum_{x} P(x) \\log\\frac{P(x)}{Q(x)}\n",
      "\\end{align*}"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "smoothing = 1e-5\n",
      "skip = 100\n",
      "skipped_chain = chain[::skip]\n",
      "window = 10000\n",
      "\n",
      "def kl(a, b):\n",
      "    return np.sum([p*np.log(p/q) for p, q in zip(a, b)])\n",
      "\n",
      "def histify(history, K):\n",
      "    _, N = history.shape\n",
      "    hist = np.zeros(K**N, dtype=np.float)\n",
      "    for h in history:\n",
      "        hist[idmap[tuple(h)]] += 1.0\n",
      "    return hist\n",
      "\n",
      "def fn(i):\n",
      "    hist = histify(skipped_chain[max(0, i-window):(i+1)], K) + smoothing\n",
      "    hist /= hist.sum()\n",
      "    return kl(actual_posterior, hist)\n",
      "\n",
      "kls = map(fn, xrange(skipped_chain.shape[0]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 60
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's take a peek at the posterior distribution of the *last* window compared with our brute force posterior distribution:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "H = histify(skipped_chain[skipped_chain.shape[0]-1-window:], K) + smoothing\n",
      "H /= H.sum()\n",
      "print \"C\\t\\tP_g(C|Y)\\tP_a(C|Y)\\t|Diff|\"\n",
      "for c, (a, b) in zip(all_assignment_vectors(), zip(H, actual_posterior)):\n",
      "    def f(x): return '%.7f' % (x)\n",
      "    print \"\\t\".join(map(str, [c,f(a),f(b),f(math.fabs(a-b))]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "C\t\tP_g(C|Y)\tP_a(C|Y)\t|Diff|\n",
        "(0, 0, 0, 0)\t0.1240000\t0.1328050\t0.0088050\n",
        "(0, 0, 0, 1)\t0.0900000\t0.0864616\t0.0035384\n",
        "(0, 0, 1, 0)\t0.0860000\t0.0864616\t0.0004616\n",
        "(0, 0, 1, 1)\t0.0200000\t0.0227718\t0.0027718\n",
        "(0, 1, 0, 0)\t0.0280000\t0.0288205\t0.0008205\n",
        "(0, 1, 0, 1)\t0.1200000\t0.0910871\t0.0289129\n",
        "(0, 1, 1, 0)\t0.0280000\t0.0227718\t0.0052282\n",
        "(0, 1, 1, 1)\t0.0340000\t0.0288205\t0.0051795\n",
        "(1, 0, 0, 0)\t0.0260000\t0.0288205\t0.0028205\n",
        "(1, 0, 0, 1)\t0.0180000\t0.0227718\t0.0047718\n",
        "(1, 0, 1, 0)\t0.0620000\t0.0910871\t0.0290871\n",
        "(1, 0, 1, 1)\t0.0380000\t0.0288205\t0.0091795\n",
        "(1, 1, 0, 0)\t0.0180000\t0.0227718\t0.0047718\n",
        "(1, 1, 0, 1)\t0.0980000\t0.0864616\t0.0115384\n",
        "(1, 1, 1, 0)\t0.0660000\t0.0864616\t0.0204616\n",
        "(1, 1, 1, 1)\t0.1440000\t0.1328050\t0.0111949\n"
       ]
      }
     ],
     "prompt_number": 61
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This looks reasonable. We can also get a sense for the trend as we increase the number of iterations:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(np.arange(0, chain.shape[0], skip) + 1, kls)\n",
      "plt.xlabel('Iterations')\n",
      "plt.ylabel('Posterior KL-divergence')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEPCAYAAABCyrPIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG2lJREFUeJzt3X2YXFWd4PFvdTovnYQESUIgISYRUPGFBVYEBaVAcGCH\nAWdlGFjXGVGMLig8MI+jgK49+zysDI47g/gyMDKoLIgi4MAiIIjF8CYoJBASggSEgfCSBEJeSeik\na//43aKqO93Vt5u6dStd38/z3KfuvXXr3lMnnfurc84954AkSZIkSZIkSZIkSZIkSZKUu3OAJcBi\n4CpgfL7JkSTlYR7wFNUg8FPgr3NLjSSprs4Mz70O6AEmAtuS1xUZXk+S1MIWAOuBlcAVOadFkpST\nPYGlwDSiJHI98IlcUyRJGlSWVUbvA+4FXk62rwM+CFxZOWDGjD3Lq1Y9mWESJGlUehLYq9En7Wj0\nCWssAw4GuoACcCRRYnjDqlVP0tNTplx2+frXv557GlplMS/MC/Oi/kLUwDRclgHhYeDHwO+BR5J9\nl/Y/aMuWDFMgSUotyyojgAuTZVBbtsCkSRmnQpI0pCxLCKlYQgjFYjHvJLQM86LKvKgyL7JXyPn6\n5aeeKjN/fs6pkKQdSKFQgAzu35YQJEmAAUGSlMg9IGzenHcKJEnQAgHBEoIktQYDgiQJaIGAYJWR\nJLWG3AOCJQRJag0GBEkS0AIBwSojSWoNuQcESwiS1BoMCJIkoAUCglVGktQacg8IlhAkqTXkHhAs\nIUhSa8g9IFhCkKTWkHtA2LABYopQSVKecg8IP/gBfOc7eadCkpR1QHgHsLBmWQuckfE1JUkj0Jnx\n+R8H9k/WO4AVwPX9DxozJuNUSJKG1MwqoyOBJ4Fna3eedRa89loTUyFJGlAzA8JJwFX9d06Y4KOn\nktQKmhUQxgF/BlzT/42uLksIktQKsm5DqDgGeBBY1f+Nu+7qZv166OyEYrFIsVhsUpIkacdQKpUo\nlUqZX6eQ+RXC1cDNwI/67S9ffHGZxx6D7363SSmRpB1coVCADO7fzagymkQ0KF830JtWGUlSa2hG\nldFGYPpgbxoQJKk15N5T2aeMJKk15B4QLCFIUmswIEiSgBYICFYZSVJryD0gWEKQpNZgQJAkAS0Q\nEKwykqTWkHtAsIQgSa3BgCBJAlogIFSqjJxXWZLylXtA6OyEjg7o6ck7JZLU3nIPCADjx8OWLXmn\nQpLaW0sEhAkTDAiSlLeWCAjjx/voqSTlrSUCgiUEScpfSwQESwiSlL+WCAiWECQpfy0RECwhSFL+\nWiYgWEKQpHxlHRB2Bn4OPAYsBQ4e6CCrjCQpf50Zn/8i4JfACcm1Jg10kFVGkpS/LAPCVOBDwF8n\n21uBtQMdaAlBkvKXpspoN+Ay4JZk+13AZ1J8bj6wCrgceAj4F2DiQAdaQpCk/KUpIfyQuKmfl2w/\nAfyMCBJDnfsA4AvA74B/Ar4C/M/ag7q7u1myBNasgblzixSLxdSJl6R2UCqVKJVKmV+nkOKY3wPv\nAxYC+yf7FgH7DfG53YD7iJICwKFEQDi25phyuVzmi1+EvfeGM85InW5JaluFQgHS3b+HJU2V0QZg\nWs32wQzSFtDPi8CzwNuT7SOBJQMdaBuCJOUvTZXR3wA3Am8D7gVmEE8NpfFF4EpgHPAkcMpAB9mG\nIEn5SxMQHgQOA95BFFGWAWmns3kYOHCogyqzpkmS8pOmyugLwGTgUWBxsn5aIxNhCUGS8pcmIHwW\nWFOzvQZY0MhE2IYgSflLExA6+h03BhjbyEQ4lpEk5S9NG8KtwNXAJUQbwueodlJrCNsQJCl/aQLC\nl4kqov+RbN8G/KCRibCEIEn5SxMQtgHfT5ZMzJoFy5dndXZJUhpperodCnwdmEc1gJSJfglvVrlc\nLtPTA7vtBo88ArNnN+CskjSK5dlT+TLg/xCB4cBkeX8jEzF2LBx2GNx9dyPPKkkajjRVRq8CN2ed\nkKlTYePGrK8iSRpMmoDwG+CbwHVAbdPvQ41MiH0RJClfaQLCwUSbwfv67T+8kQnxSSNJyleagFDM\nOhHg8BWSlLcsZ0wbFquMJClfaQLCD4FfAbOS7SeAsxqdEKuMJClfaQLCdOCnRAc1iKGvtzY6IVYZ\nSVK+spwxbVgsIUhSvrKeMS012xAkKV/DnTEN4HHSz5iWmiUEScpXmoDwcaIfQsXbiSqjxcDKRiXE\nNgRJyleagPBp4ANEj2WIfgkPAfOB/wX8eIjPPw2sIxqlexhkHCRLCJKUrzQBYSywD/BSsj0TuAI4\nCPh3hg4IZSKIvFLvINsQJClfaZ4ymkM1GEBUE80BXgZeT3mdIYdptcpIkvKVdnC7m4CfETf2jwMl\nYBIxEupQysDtRJXRJcC/DHSQVUaSlK80AeF0Iggckmz/CLiWuNGnGeDuEOAF4nHV24BlwF2VN7u7\nuwF44QVYvbpIk4ZOkqQdRqlUolQqZX6doapyOoFHgXc26HpfJzq6fSvZLpfL8QDT0qVwwgnxKkka\nXF4zpm0l+h3MHeH5JwI7JeuTgI8Sj6tuxzYEScpXmiqjXYAlwANAZU6zMnBcis/OBK6vudaVxEB5\n27ENQZLylSYgfG2AfeUB9g3kj8B+aQ40IEhSvtIEhBIwD9iLeFpoYsrPDcuECVYZSVKe0vRDWABc\nQzwyCrAH1WqghrGEIEn5ShMQTgcOJYafAPgDsGujEzJ2LGzbFoskqfnSBIQtyVLRSfo2hNQKhQgK\nWxs+9Y4kKY00AeFO4Dyi7eAoovroxiwS09kJPQ0fWFuSlEaajg0dwKlEHwKAW4Ef0JhSwhsd0wB2\n3hmefjpeJUkDy6pjWpqnhT5GDFdxaaMv3p8lBEnKT5oqo+OAJ4ghr48lg0dOK2xDkKT8pAkInyL6\nIPwcOBl4Crgsi8RYQpCk/KT9tf86cDPQSzQufwz4TMMT02kJQZLykqaE8F+AHxLVRicQ8xnMzCIx\nVhlJUn7SlBD+Crga+DyQ6eASVhlJUn7SBISTMk9FwhKCJOWnXpXRPcnrBmB9v2XdYB96MywhSFJ+\n6pUQKlNmTm5GQsASgiTlqV5A2GWIz77SyISAJQRJylO9gPAQMTxFAXgrsCbZ/xbgGWB+wxPjY6eS\nlJt6bQjziJv+bUQP5WnJ8qfJvoazykiS8pOmH8IHgF/WbN8MfDCLxFhlJEn5SRMQnge+SrXEcB6w\nYhjXGAMsJMWQ2ZYQJCk/aQLCycQMadcD1yXrJw/jGmcCS0kxXLYlBEnKT5qOaS8DZyTruwMvDOP8\nexBDX5wPnD1kYmxUlqTcpCkh1LppmMf/I/AlYlC8IVllJEn5Ge7cBsOZoedYYCXRflAc7KDu7u43\n1levLtLTM+ihktSWSqUSpVIp8+sMdwq204DvAbMZumH5fwOfBLYCE4ApwLXEYHkVfabQPPVUOPjg\neJUkDSyrKTSHW2X0veT1vhTHngvMIZ5MOgm4g77BYDs2KktSfoYbECpGEpmGfMrINgRJyk9m8yP3\nc2ey1GUJQZLyUy8gXFznvZ0bnRDwsVNJylO9gPAgA1fzFIDfZ5EYq4wkKT/1AsKiZBnIaRmkxSoj\nScpRvUbl64D3DbD/74DPZpEYSwiSlJ96AeEvgJ9RHdm0A/hn4LBkaThLCJKUn3oB4UHgY8AVwNHA\nNcAM4E/IaE5lSwiSlJ96AWEX4DngU8CVQA/wOWASQ0+vOSKWECQpP2mm0ARYDxwE/C7ZLgNva3hi\nfOxUknJTLyDMa1YiKqwykqT8jHToikxYZSRJ+WmpgGAJQZLy01IBwRKCJOVnqIDQCTzejISAjcqS\nlKehAsJWYBkwtwlpYexYeP31ZlxJktRfmuGvdwGWAA8AG5N9ZeC4Ridm+nRYvTrW77wTOjrgQx9q\n9FUkSQNJExC+lrxW+iQUSDHZzUjMmQPPPhvrn/gErFgBmzZBV1cWV5Mk1UrTqFwiqo2mADsBS0kx\n2c1IzJoFL74IjzwSpYMpU2Dz5iyuJEnqL01AOBG4nxjs7kSi6ugvskjMuHEwbRqcf36UECZMsE1B\nkpolzdzIjwBHAiuT7RnAr4F9G3D9crnct/Zp771h+XJ46in48IfhnnvgrW9twJUkaZQoFAowsrnt\n60pTQigAq2q2Xx5GQiYQpYtFRFXTN4b6wMsvw777wvz5MH68JQRJapY0jcq3ALcCVxGB4C+Bm1Oe\nfzNwOLApudbdwKHJ64AeewymTo31ceMMCJLULGkCwt8C/5W4kZeBS4Drh3GNTcnrOGAM8Eq9g2fO\nrK4bECSpedIEhDJwbbKMRAcxlPaewPeJqqNUxo2DLVtGeFVJ0rDUCwj3AIcAG9i+30GZeAw1jV5g\nP2AqUfVUJB5lBaC7u/uNA4vFIsVi8Y1t2xAkCUqlEqVSKfPrNLyVeghfA14D/iHZ3u4po1of+Qic\ne268SpJCXk8ZdRKd0kZqOrBzst4FHAUsTPthq4wkqXmGakPYSox2Ohd4ZgTn3x34ERF4OoAriD4M\nqVhlJEnNk/XgdouBA0aWNJ8ykqRmGs7gdrUyGdyuPwOCJDVPmoBQAuYBewG3AxNTfu5Nsw1Bkpon\nzdAVC4BriA5pAHswvI5pI2YbgiQ1T5qAcDrRS3ldsv0HYNfMUlTDKiNJap40AWFLslR00sQ2BKuM\nJKk50gSEO4HziLaDo4jqoxuzTFSFVUaS1DxpAsKXieGvFwOfA34JfDXLRFVYZSRJzZPmaaEvAhcB\nl9bsOzPZl6lx42JOZUlS9tKUED41wL5TGpyOAdmGIEnNU6+EcDLw34D59G0z2ImYNS1ztiFIUvPU\nCwj3Ai8Qcyj/A9WR9dYR8yxnzjYESWqeegHhmWQ5khiyehvwjmRZnH3SrDKSpGZK+9jpeGA2McHN\nJ4EfZpimN9RWGf3iFzHf8h13wG9+04yrS1J7SfOUUQcxL/JngO8BFwIPZ5moinHjYM0a+OY34fzz\nYY89YOXKePLomWdg2rRmpEKS2kPaQeo+AHyCCAqQrmTxps2YAffdB729sHAh3HADTJ8eAeI//sOA\nIEmNlGYKtsOAvyHmWP57YE+iH8IZDbh+3Sk04wAo9EvlscfCggVwXJoZGSRplMlqCs00JYQ7k2Un\nYDLwJI0JBqn0DwYQVUfPPdesFEhSe0hT9fNeYh7kJcBS4EHgPVkmaihz5kRAuOoqOOIIuPPOPFMj\nSaNDmhLCpcDZQOXZnmKy74MZpWlIs2fDKafAPvvAO98JpRIcdlheqZGk0SFNCWEi1WAAMYPapJTn\nn5N8dgnwKA2qatp3X5g1Cx54AI46ClasaMRZJam9pQkIfyTmVZ5HDGPxVeCplOfvAc4C3g0cTEy2\ns8+wU9nPAQdEEJg8OdoTDAiS9OalCQinEDOkXQdcSwxl8emU538RWJSsbwAeA2YNM411zZ5tQJCk\nRqjXhtAFfB7Yixi76GziF/9IzQP2B+5/E+fYjgFBkhqjXkD4EfA6cDdwDPAuov/BSEwGfp58fkPt\nG93d3W+sF4tFisXisE48YwasXg3f+hacffbAj6lK0o6sVCpRKpUyv0692+di4pFTiMDxO+IX/nCN\nBf4fcDPwT/3eG7JjWhrf/S5cfjm8//1w/PGw334wc+abPq0ktaSsOqbVa0PYOsj6cBSAy4j+C/2D\nQcOcfjr8+tfw6qtw9NFw2WVZXUmSRq96VUb7AutrtrtqtsvAlBTnPwT470QbxMJk3znALcNL5tCm\nTo2OagcdBMuXN/rskjT61QsIYxpw/rtp0kB4FbNn23NZkkaiqTfrZpg1C55/Pu9USNKOx4AgSQJG\nYUDYffcICOecY9WRJA1H2glydhjjx8O2bXDBBdDT46B3kpTWqCshQEyx+e1vw+bNeadEknYcozIg\ndHXFVJurV+edEknacYzKgAAGBEkarlEdEFatyjsVkrTjGNUBwRKCJKWX99igDRncbiCvvQY77xwN\ny46AKmk0yWpwu7xvlZkFBIhAcMcdMeXmhg1wyy0wfz4sWxb7hjnStiS1BAPCCCxYEEHg1Vdhp53g\nwANh5crozbx8OSxaNPQ5JKnVGBBGfIF4ra02WroUTjghXiVpR5NVQBh1PZX7G6j9YPx42LKl+WmR\npFY2ap8yqmfcOHj99bxTIUmtxYAgSQLaNCBYZSRJ22vLgGAJQZK2Z0CQJAHZB4R/BV4CFmd8nWEZ\nMwZ6e2PeBElSyDogXA4cnfE1hq1QiHYESwmSVJV1QLgLWJPxNUbEaiNJ6qst2xAgAoJPGklSVe49\nlbu7u99YLxaLFJs04pxVRpJ2FKVSiVKplPl1mjGW0TzgRuC9A7yX+VhGg3nb2+D22+NVknYkWY1l\nZJWRJAnIPiD8BLgXeDvwLHBKxtdLzSojSeor6zaEkzM+/4j5lJEk9WWVkSQJaOOAYJWRJPXVtgHB\nKiNJ6qutA4JVRpJU1bYBwSojSeqrbQOCVUaS1FdbBwSrjCSpqm0DglVGktRX2wYEq4wkqa/cRzvN\ny9ixcNttsHQpvPwyXHQRzJqVd6okKT9tGxAOOABuvRUOPTSm0jz1VDjtNJgzB7Zuhddeg7lzY3vt\nWrjuOnjwQVi1Cl56CSZOjGk4p06FnXeO12OOgcMPhw0b4hqFQhxXaMaYspL0JuV9q8pt+OtamzbB\nuefCY4/Bc89F+0JXFyxbFlVLmzbBEUdAsQi77gozZ8a+jo4IFmvXwiuvRClj82bo6YkgsG1bBIrj\nj49jAd7znljfsAE2boS994a99oJ16+LV4bglDSWr4a8NCHWUy/D88zBpUpQChrJhQ5Qupk6tlgoW\nL4a77qoGiIcfjiAzeTJMmBBVVs88EyWSd78brr462+8kacdnQBjlfvUruPDCmLRHkupxgpxRbtq0\naNyWpLwYEFqEAUFS3gwILcKAIClvBoQWMXlyPJ20eXPeKZHUrrIOCEcDy4AngC9nfK0dWqFgKUFS\nvrIMCGOA7xBB4V3E/Mr7ZHi9HVqpVDIgJEqlUt5JaBnmRZV5kb0seyq/H1gOPJ1sXw0cDzyW4TV3\nWKVSicmTi1x8MZx4Yoyz9Pzz0Xnt8MNhl11gypSoWhozpvq5nh548UV49VWYMSM6zkH0a9i0KY5d\nt65aFdXbG/0ldtkl1js64C1vib4T27bF+bq6+l6j2UqlEsViMb8EtBDzosq8yF6WAWE28GzN9nPA\nQRleb4f36U/DHXfAN74RndZ23z32X3IJrF8fN/aNG6s36xkzYiiN6dPjpr5yJaxZEzf6rq5Ytm6N\nQDJxYnymUIigsmpVjOfU2xufWbs2zjt2bASPrq743JQpESzGj4/PVpaOjur6TjvB7NkRxCpLT08c\nM3bs9ktnv7+6/ue7++4ITrX7KusTJsR37R+wBhoeJM2+Vjum8j07OmJ5/HG46aZqPlQ+09EReTBm\nTKyXy/WX3t6hj2nGUvmOY8bE30HHIHUUnZ3xb137vZctgxtuqH7vwYaEqde1aaD3Bju+Nr8HWx/q\n/Xrrg71fe+5661nIMiDY42yYFiyIpZ7e3rjZbtsWN/VZs+ImW9HTE/+Z3swfTm9vBJ5162JZuzbm\njhjsJrNmTZRSxo+PXtjjxlWDTU/P9svWrdX0DXS+zs74fGVf5TzlcgwRsmRJ3//Eaf+T99/XisfU\n5kNvb9wE162L9f7HbNsWS2/v9jeZgQJ3KyxQTfvWrYPfjCsPWFTeL5cjLzZurH7veur9/acJzLXX\nHWx9qPfrrQ/2fu25B1vPsi9vlvHmYKCbaEMAOAfoBf6+5pjlwJ4ZpkGSRqMngb3yTsRwdBKJngeM\nAxZho7Ikta1jgMeJksA5OadFkiRJUisbjZ3W/hV4CVhcs28X4DbgD8CvgNqBtM8hvv8y4KM1+/9z\nco4ngItq9o8Hfprs/y0wt7HJb6g5wG+AJcCjwBnJ/nbMjwnA/US16VLgG8n+dsyLijHAQuDGZLtd\n8+Jp4BEiLx5I9rVdXowhqpHmAWMZPe0LHwL2p29AuBD422T9y8AFyfq7iO89lsiH5VQb+R8g+nEA\n/JJqw/xpwPeS9b8k+na0qt2A/ZL1yUTV4T60b34kD/7SSfzHPJT2zQuAs4ErgRuS7XbNiz8SAaBW\n2+XFB4Bbara/kiyjwTz6BoRlwMxkfbdkGyLS15aMbiGezNqdvp33TgL+ueaYSl+OTmBVoxLdBL8A\njsT8mAj8Dng37ZsXewC3A4dTLSG0a178EZjWb19ueZHX4HYDdVqbnVNasjaTqEYiea38Q88ivndF\nJQ/6719BNW9q820rsJbtf120onlEyel+2jc/Oohfdy9RrUpr17z4R+BLxGPoFe2aF2UiOP4e+Gyy\nL7e8yLJjWj3t2mmtTPt998nAtcCZwPp+77VTfvQSVWhTgVuJX8e12iUvjgVWEnXmxUGOaZe8ADgE\neAGYQbQbLOv3flPzIq8Swgqi0bFiDn0j3GjyElHsgyjarUzW++fBHkQerEjW+++vfOatyXoncXN5\npfFJbpixRDC4gqgygvbOD4hfaDcRjYDtmBcfBI4jqkp+AhxB/H20Y15ABAOIqpzriXaAtsuL0dxp\nbR7bNypX6v2+wvYNROOA+UR+VBqI7ifq/Qps30D0/WT9JFq0gShRAH5MVA/Uasf8mE71SZEu4N+B\nj9CeeVHrMKptCO2YFxOBnZL1ScA9xJND7ZgXo7LT2k+A54HXiXq7U4j6utsZ+BGyc4nvvwz4k5r9\nlUfIlgPfrtk/HvgZ1UfI5mXwHRrlUKKaZBFRPbCQ+CNtx/x4L/AQkRePEPXn0J55Ueswqk8ZtWNe\nzCf+JhYRj2ZX7oPtmBeSJEmSJEmSJEmSJEmSJEmS1Cwbkte5wMkNPve5/bbvafD5JUkNVBkfqUi1\n12taQ43l1X/sJUlSC6vctH8LvEr0iD6TGKfrm8Q48Q8DC5LjisBdwL9RHUjsF8Rok49SHXHyAmJ0\nyIXEWDtQLY0UknMvJnoen1hz7hJwDTEc8f+tSecFxAinDyeflSQ1WCUg1I6LAxEAzkvWxxNzEMwj\nbtob6Dtz1FuS1y7iJl/Z7l9CqGx/nBhSoADsCjxDDEZWJILSrOS9e4nRLKfRdxTLKWm/nNRMeY12\nKjVaod/2R4G/In7h/5YYH2av5L0HiJt4xZnEeDL3EaNJ7j3EtQ4FriKGJV4J3AkcmGw/QIxnVU7O\nOZcIEpuBy4A/B14b7peTmsGAoNHsC8TEPPsDexIDhgFsrDmmSIw8ejAxX8FCYg7kespsH4AqY9Zv\nqdm3jRgCfBsxrPHPifkAbkFqQQYEjRbrqQ4lDDEJzWlUG47fTnVe41pTgDXEL/h3EoGhooeBG57v\nIuan7SAmNvkwUTLoHyQqJhEjVt5MzCX8n4b8NlIO8poxTWqUyi/zh4lf4ouAy4khgOcRw04XiKqd\nP2f7GahuAT4PLCWGY7+v5r1LiUbjB4FP1nzuemJe8IeTfV9Kzr8P289uVSYC1b8RJY8CcNaIv60k\nSZIkSZIkSZIkSZIkSZIkSZIkSVKr+P+DJR/nhH+rSgAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x114d38690>"
       ]
      }
     ],
     "prompt_number": 64
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Posterior predictive distribution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Another distribution we can look at is the posterior *predictive* distribution $p(y|\\mathcal{C},\\mathcal{Y};\\mathcal{H})$, that is, the distribution over the *next* data point. To derive this, we first note that $p(y|\\mathcal{C},\\mathcal{Y}) = \\sum_{k=1}^{K} p(y,c{=}k|\\mathcal{C},\\mathcal{Y}) = \\sum_{k=1}^{K} p(c{=}k|\\mathcal{C}) p(y|\\mathcal{Y}_{k})$.\n",
      "\n",
      "From [Eq. 9, pg. 3](http://homepage.tudelft.nl/19j49/Publications_files/TR_1.pdf), we have $p(c{=}k|\\mathcal{C}) = \\frac{|\\mathcal{C}_{k}| + \\frac{\\alpha}{K}}{|\\mathcal{C}| + \\alpha}$. Next, $p(y|\\mathcal{Y}_{k})$ is simply the posterior predictive of the beta-bernoulli likelihood model. That is, $p(y_d{=}1|\\mathcal{Y}_{k})= \\frac{\\beta + \\sum_{y_i \\in \\mathcal{Y}_{k}} y_i^d}{ \\beta + \\gamma + |\\mathcal{Y}_{k}|  } $. Putting it together, we have\n",
      "\\begin{align*}\n",
      "  p(y|\\mathcal{C},\\mathcal{Y}) = \\sum_{k=1}^{K} \\frac{|\\mathcal{C}_{k}| + \\frac{\\alpha}{K}}{|\\mathcal{C}| + \\alpha} \\prod_{d=1}^{D}  (\\hat{p}_d^k)^{y_d} (1-\\hat{p}_d^k)^{1-y_d}\n",
      "\\end{align*}\n",
      "where $\\hat{p}_d^k = \\frac{\\beta + \\sum_{y_i \\in \\mathcal{Y}_{k}} y_i^d}{ \\beta + \\gamma + |\\mathcal{Y}_{k}| } $.\n",
      "\n",
      "We can also look at the non-Bayesian case, where we know the exact values for $\\pi$ and $\\{ p_d^k \\}$.\n",
      "\\begin{align*}\n",
      "  p(y | \\pi, \\{p_d^k\\}) = \\sum_{k=1}^{K} \\pi_k \\prod_{d=1}^{D} (p_d^k)^{y_d} (1-p_d^k)^{1-y_d}\n",
      "\\end{align*}\n",
      "\n",
      "This makes explicit how our prior information affects our estimates of $\\hat{\\pi}$ and $\\hat{p_d^k}$. Now for the reference distributions, we know both the ground truth cluster assignment and the exact model parameters, so we can compute $ p(y|\\mathcal{C},\\mathcal{Y}) $ and $ p(y | \\pi, \\{p_d^k\\})$ exactly. For our Gibbs sampler, we estimate $p(y|\\mathcal{C},\\mathcal{Y})$ by averaging over our samples of $\\mathcal{C}$. Details below:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def posterior_predictive1(C, Y, K, alpha, beta, gamma):\n",
      "    N, D = Y.shape\n",
      "    assert C.shape[0] == N\n",
      "    \n",
      "    nks = np.bincount(C, minlength=K)\n",
      "    sums = np.zeros((D, K), dtype=np.int64) \n",
      "    for yi, ci in zip(Y, C): \n",
      "        sums[:,ci] += yi\n",
      "        \n",
      "    lg_term1 = np.log(nks + alpha/K) - np.log(N + alpha)\n",
      "\n",
      "def posterior_predictive(C, Y, K, alpha, beta, gamma):\n",
      "    N, D = Y.shape\n",
      "    assert C.shape[0] == N\n",
      "\n",
      "    nks = np.bincount(C, minlength=K)\n",
      "    sums = np.zeros((K, D), dtype=np.int64) \n",
      "    for yi, ci in zip(Y, C): \n",
      "        sums[ci] += yi\n",
      "\n",
      "    def fn(yvalue):\n",
      "        def fn1(nk, sum_yid, yd):\n",
      "            assert nk >= sum_yid\n",
      "            theta = (beta + sum_yid) / (beta + gamma + nk) \n",
      "            assert theta >= 0.0 and theta <= 1.0 \n",
      "            return np.log(theta) if yd else np.log(1.-theta)\n",
      "        def fn2(nk, row):\n",
      "            assert len(yvalue) == row.shape[0]\n",
      "            term1 = np.log(nk + alpha/K) - np.log(N + alpha)\n",
      "            term2 = sum(fn1(nk, sum_yid, yd) for sum_yid, yd in zip(row, yvalue))\n",
      "            return term1 + term2\n",
      "        return sp.misc.logsumexp([fn2(nk, row) for nk, row in zip(nks, sums)])\n",
      "\n",
      "    yvalues = it.product([0, 1], repeat=D)\n",
      "    lg_pr_yvalue = map(fn, yvalues)\n",
      "    return np.exp(lg_pr_yvalue)\n",
      "\n",
      "def posterior_predictive_nonbayes(pis, aks):\n",
      "    K, = pis.shape\n",
      "    assert aks.shape[0] == K\n",
      "    _, D = aks.shape\n",
      "    \n",
      "    def fn(yvalue):\n",
      "        def fn2(yd, theta_kd):\n",
      "            return np.log(theta_kd) if yd else np.log(1.-theta_kd)\n",
      "        def fn1(pi_k, theta_k):\n",
      "            return np.log(pi_k) + sum(fn2(yd, theta_kd) for yd, theta_kd in zip(yvalue, theta_k))\n",
      "        return sp.misc.logsumexp([fn1(pi_k, theta_k) for pi_k, theta_k in zip(pis, aks)])\n",
      "    \n",
      "    yvalues = it.product([0, 1], repeat=D)\n",
      "    lg_pr_yvalue = map(fn, yvalues)\n",
      "    return np.exp(lg_pr_yvalue)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 65
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# reference distributions\n",
      "actual_posterior_predictive = posterior_predictive(cis, Y, K, alpha, beta, gamma)\n",
      "nonbayes_posterior_predictive = posterior_predictive_nonbayes(pis, aks)\n",
      "\n",
      "posterior_predictives = np.array([\n",
      "    posterior_predictive(assignment, Y, K, alpha, beta, gamma) for assignment in skipped_chain])\n",
      "    \n",
      "def fn1(i):\n",
      "    posteriors = posterior_predictives[min(0, i-window):(i+1)].mean(axis=0)\n",
      "    return kl(actual_posterior_predictive, posteriors)\n",
      "\n",
      "posterior_predictive_kls = map(fn1, xrange(posterior_predictives.shape[0]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 66
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "PH = posterior_predictives[posterior_predictives.shape[0]-1-window:].mean(axis=0)\n",
      "print \"y\\t\\tP(y|C) [gibbs]\\tP(y|C) [actual]\\tP(y|params)\\t|gibbs-actual|\"\n",
      "for y, ((a, b), c) in zip(it.product([0,1],repeat=D), zip(zip(PH, actual_posterior_predictive), nonbayes_posterior_predictive)):\n",
      "    print \"\\t\".join(map(str, [y, a, b, c, math.fabs(a-b)]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "y\t\tP(y|C) [gibbs]\tP(y|C) [actual]\tP(y|params)\t|gibbs-actual|\n",
        "(0, 0, 0)\t0.129465256173\t0.130765432099\t0.145796887467\t0.00130017592593\n",
        "(0, 0, 1)\t0.0661453919753\t0.0707160493827\t0.0602427431159\t0.00457065740741\n",
        "(0, 1, 0)\t0.107173910494\t0.0973827160494\t0.271569844867\t0.00979119444444\n",
        "(0, 1, 1)\t0.0656987746914\t0.0566913580247\t0.135681228984\t0.00900741666667\n",
        "(1, 0, 0)\t0.196573169753\t0.177382716049\t0.0713598709559\t0.0191904537037\n",
        "(1, 0, 1)\t0.107571737654\t0.110024691358\t0.0312461999697\t0.0024529537037\n",
        "(1, 1, 0)\t0.197682108025\t0.216691358025\t0.189236544177\t0.01900925\n",
        "(1, 1, 1)\t0.129689651235\t0.140345679012\t0.0948666804639\t0.0106560277778\n"
       ]
      }
     ],
     "prompt_number": 67
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's take a look at the posterior predictive distributions of the *last* window compared to the ground truth distributions:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice how our Gibbs estimate of the distribution does not quite match the non-Bayesian one, which makes sense. Since $N$ is small in our case, our prior distribution is still weighing in non-negligibly. Finally, let's take a look at the KL-divergence (w.r.t. the *Bayesian* distribution with ground truth clusters) as the number of iterations increases:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(np.arange(0, niters, skip) + 1, posterior_predictive_kls)\n",
      "plt.xlabel('Iterations')\n",
      "plt.ylabel('Posterior predictive KL-divergence')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 68,
       "text": [
        "<matplotlib.text.Text at 0x11546be10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEPCAYAAAB7rQKTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYFOW59/HvLKKiyKogi44sKiCiIAiIMu6Eo7gkqKhx\nwX0/UeP6qpyYRaPGSEwMx4gSI+ISNxJEzIkDRgWUHdmRwQEURAWRZZyl3j/uartpZrqfxq6pXn6f\n6+pruqqrup8pmLr72e4HRERERERERERERERERERERERERERERPLCYGAxsAy4vZ5jRvmvzwWOjNk/\nBlgHzK/nvFuAWqBFWkoqIiJZpwhYDpQAuwFzgK5xxwwBJvrPjwamxbx2LBZ46go0HYBJwEoUaERE\nMlphgO/dFws05UAVMB44I+6YocBY//l0oBnQxt9+F/i6nvf+HXBbGssqIiIBCTLQtAMqYrZX+/tS\nPSbeGf5x835oAUVEJHjFAb6353hcQQrnNQbuAk5OcL6IiGSQIAPNGqwvJaIDVhNJdEx7f199OmF9\nPnNjjp+JNdOt3+HATp28FStWpFxoEZE8twLoHHYhXBVjBS4BGpF8MEA/dhwMgH9ufaPOIPFgAE/M\nfffdF3YRMoauRZSuRZSuRRTurVHOguyjqQauB94CFgIvAIuAq/wHWJD5BBs0MBq4Nub854H3gYOx\nfpxL6/iMtF8QERFJryCbzgDe9B+xRsdtX1/PucMd3r9jyiUSEZEGFWSNRjJEaWlp2EXIGLoWUboW\nUboWwcrlEVt+c6OIiLgqKCiANMcG1WhERCRQCjQiIhIoBRoREQmUAo2IiARKgUZERAKlQCMiIoFS\noBERkUAp0IiISKAUaEREJFAKNCIiEigFGhERCZQCjYiIBEqBRkREAqVAIyIigVKgERGRQCnQiIhI\noBRoREQkUAo0IiISKAUaEREJlAKNiIgESoFGREQC5RJo9gLuAZ70t7sApwVWIhERySkugeZp4Dtg\ngL+9FvhVYCUSEZGc4hJoOgEPYsEGYEtwxRERkVzjEmgqgT1jtjv5+0RERJIqdjhmJDAJaA+MA44B\nLgmuSCIikktcajSTgR8Dl2KBpjfwTgqfMRhYDCwDbq/nmFH+63OBI2P2jwHWAfPjjn8IWOQf/wrQ\nNIXyiIhIA3IJNGcD1cA//Ec1cKbj+xcBj2PBphswHOgad8wQoDM2mu1K4ImY1572z403GegO9ASW\nAnc6lkdERBqYS6C5D9gYs70Ra05z0RdYDpQDVcB44Iy4Y4YCY/3n04FmQBt/+13g6zre922gNuac\n9o7lERGRBuYSaArq2Ffk+P7tgIqY7dX+vlSPSWQEMDGF40VEpAG5BJqZwO+w0WadgUf9fS48x+Pi\ng5nreXdjw67HOR4vIiINzGXU2Q1YZoAX/O23gesc338N0CFmuwNWY0l0THt/XzKXYP07J9Z3wMiR\nI79/XlpaSmlpqcPbiojkj7KyMsrKygL9jLqaxdKpGFiCBYO1wAxsQMCimGOGANf7P/sBv/d/RpQA\nE4AeMfsGA48Ag4AN9Xy253muFSMREQEoKCiANMcGlxrNIcCt2A0/crwHnOBwbjUWRN7C+nWewoLM\nVf7ro7H+lSHYoIEt2DDqiOexYNIS68e5FxuJ9gegEVa7AvgAuNahPCIi0sBcotY8bMjxLKDG3+fh\n3k8TFtVoRERSFESNxuXNZmKTNLONAo2ISIqCCDQuo84mYJ3/+wMtYh4iIiJJuUStcuoebnxQeouS\ndqrRiIikKKyms2ylQCMikqKwms60wqaIiOwyrbApIiKB0gqbIiISKK2wKSIigdIKmyIiEijXkQWt\niOYfm0b9+cUyiUadiYikKKzhzb2JzqMp8J9vAlZhucwylQKNiEiKwgo007BgM8/f7gF8DDQFrsES\nZmYiBRoRkRSFNY9mLXAEFmx6+88/AU4GfpvOwoiISO5xCTSHYDWYiIXAocAK3FfCFBGRPOUy6uxj\nbJmA8Vh16hws2OwOVAVXNBERyQUu7XB7Ytmbj/G33wP+BGzH0tNsDqZoP5j6aEREUhTGYIBibBXL\n49P5oQ1EgUZEJEVhDAaoBmqBZun8UBERyR8ufTRbgPlYzSaS58wDbgyqUCIikjtcAs0r/iN+0qaI\niEhSru1wjYEDgMUBliXd1EcjIpKisCZsDgVmY4k1AY4E3khnIUREJHe5BJqRwNHA1/72bKBjUAUS\nEZHc4hJoqoCNcftqAyiLiIjkIJdA8zFwATZwoAvwB+D9IAslIiK5wyXQ3AB0x1bVfB74BvjvIAsl\nIiK5w2VkQS9gVtAFCYBGnYmIpCis9WjKgDbAS8ALwIJ0FiBACjQiIikKa3hzKZbrbAMwGssScE86\nCyEiIrnLJdAAfAY8BlwNzAXudTxvMDbJcxlwez3HjPJfn4vN0YkYA6zDAlusFlg6nKXAZJSHTUQk\no7kEmm7YXJoFwOPYiLN2DucV+ccP9t9jONA17pghQGdsNNuV2Lo3EU/758a7Aws0BwP/52+LiEiG\ncgk0Y7B5NKcCg7C1aNY7nNcXWA6UY3NxxgNnxB0zFBjrP5+O1U7a+NvvEp0kWt85Y4EzHcoiIiIh\ncUmq2W8X37sdUBGzvRrLMJDsmHbA5wnetzXWpIb/s/Uulk9ERBpAokDzEjCMnftIwLI3H57kvV2H\nfMWPbkhlqJiX6PiRI0d+/7y0tJTS0tIU3lpEJPeVlZVRVlYW6GckGsLWFlgLlNTzenmS9+6H9e1E\n+lnuxFLXPBhzzJ+x4dPj/e3FWPNcpMZSAkwAesScsxgbCfc5sD/wDnBoHZ+v4c0iIilq6OHNa/2f\n5fU8kvkI6+QvARoB57Jz1uc3gIv85/2wvqB1JPYGcLH//GLgNYeyiIhISBJFrW+pv1nKA/ZxeP8f\nAb/HRqA9BfwGuMp/bbT/MzIybQtwKdEsBM9jtZuW2OCDe7GRaC2AF7H1ccqBc9g56SeoRiMikrKw\nMgP8Eqvd/M3fvgBrVsv0SZsKNCIiKQor0Mxj547/uvZlGgUaEZEUhZWCZgtwIdb8VYTVaL5NZyFE\nRCR3uQSa87F+kHX+4xx/n4iISFJprR5lGDWdiYikKKyms1jZuC6NiIiEKNVAk8s1IBERCUCqgeaf\ngZRCRERy1q7WUD7FJkxmMvXRiIikKBP6aCLUhCYiIk52NdCIiIg4SbRMwC0JXts73QUREZHclCjQ\nNKH+pJq/D6AsIiKSgxIFmjFYp39dTg+gLCIikoMS9dG8DRxUx/4RwGPBFEdERHJNokDzM2AycHDM\nvjuBm4HjgiyUiIjkjkRNZxOBSuBN4AzgcqAvcCzwdfBFExGRXOAyH+Y44FXgPSxz8/ZAS5Q+mrAp\nIpKihl74LHYp5z2A74Baf9t1KecwKdCIiKQorBU2s5UCjYhIijIpBY2IiIgTBRoREQmUAo2IiATK\nNdCUACf5zxuT+QMBREQkQ7gEmiuBl4DR/nZ7bLiziIhIUi6B5jpgIPCNv70U2C+wEqVRTU3YJRAR\nEZdAU+k/IoqpP6tzRtmeLVNLRURymEugmQLcjfXNnIw1o00IslDpsm1b2CUQERGXSTmFWJ6zU/zt\nt4C/kPm1Gq+iwqN9+7CLISKSPcKasHkmMBb4if94EvcgMxhYDCwDbq/nmFH+63OBIx3O7QvMAGYD\nHwJ96vvw1atByQFERMLlEmiGYjf7Z4HTSJzxOVYR8DgWMLoBw4GucccMAToDXbDRbU84nPtb4B4s\nKN3rb9epf3945RXH0oqISCBcAs0lWDB4GbvhfwI85XBeX2A5UA5UAeOx5QZiDcVqSwDTgWZAmyTn\nfgY09Z83A9YkKsQ33yR6VUREguZaO/kOW5emFhsUcCZwWZJz2gEVMdurgaMdjmkHtE1w7h3Af4CH\nsUDZP1EhGjdOUkoREQmUS6AZgq1DczxQhvXRDHM4z7V3JNVOp6eAG7FJo8OAMdhouDqM5LXXYNEi\nKC0tpbS0NMWPEhHJbWVlZZSVlQX6GS43+fH+YxKpLXrWDxiJ9bOALQNdCzwYc8yfseA13t9eDAwC\nDkpw7jdEU+AUABuJNqXF8sBjwgQ47bQUSi0iksfCGnV2HvAaqa+s+RHWyV8CNALOBd6IO+YN4CL/\neT8saKxLcu5yLBgBnIBlKqhXVVWKpRYRkbRK1HT2HnAMO660GeGywmY1cD0276YIa/JaBFzlvz4a\nmIg1zS0HtgCXJjkXbHTaH4HdgW3+dr2++y5JKUVEJFA5vcImeDz7LFx4YdhFERHJDmE1nT3ruC8j\nqelMRCRcLoHmsLjtYqB3AGUJhAKNiEi4EgWau4DNQA//Z+Sxnp079TOWAo2ISLgSBZpfA02Ah/yf\nkUcLbNJkxmvZUoFGRCRsLk1nH2KpXiKaYZkBMt4llyjQiIiEzSXQ3IfNb4nYiE2mzHiNGinQiIiE\nzSXQ1DXMrSjdBQnCbrsp0IiIhM0l0MwEfgd0wrI4P+rvy3gKNCIi4XMJNDdgqfpfwHKSbQeuC7JQ\n6bLbbsoMICISNpfszd9S/+qYGU01GhGR8CUKNI8BNwET6njNwxYty2gKNCIi4UsUaCJpZh5piIIE\nQYFGRCR8iQLNR/7PsgYoRyAUaEREwpco0MxP8JoHHJ7msqSdAo2ISPgSBZrT/Z/X+j+fxebUXBBo\nidJIEzZFRMKXKNCU+z9PAY6I2T8PmE0WjERTjUZEJHyumQEGxmwfQ5YsmKZAIyISPpd5NCOAp4Gm\n/vZGoksuZzQFGhGR8LkEmplYx39TrCazMfHhmUOBRkQkfC5NZ22Ap7AUNBuBbsBlQRYqXZSCRkQk\nfC6B5hlgMtDW314G/CyoAqWTajQiIuFzCTStsNpMjb9dBVQHVqI0UqAREQmfS6D5FmgZs90P2BRM\ncdJL82hERMLnMhjgFiyxZkfgfWBf4CdBFipdVKMREQlfskBTBBznPw7FRp0tAbKii12BRkQkfMma\nzmqA87E+mQVY/rOsCDKgUWciIpnAZYb/o8Bu2ICALf45HjArwHKlg7d1q0eLFvDNNxZ0REQksYKC\nAkhz9heXNyvDAku849NZkAB4nudx8MHw2mvQrVvYxRERyXxBBBqXUWelWFCJf7gYDCzG5t7Ul4Rz\nlP/6XOBIx3NvABZhzXkPJirAYYfBggWOpRURkbRznUfzByxj8yxsieeWCc8wRcDjWMDoBgwHusYd\nMwToDHQBrgSecDj3eGwZ6cOBw4CHExVCgUZEJFwugWY8sB44GxvW/AXWX5NMX2A5ttxAlf8+Z8Qd\nMxQY6z+fDjTDUt4kOvca4Df+fvzy1KtlS9iYNdnZRERyj2uus/uBlcAnwC+B1g7ntQMqYrZX+/tc\njmmb4Nwu2HDraVj/0VGJCqFJmyIi4XKZsDkZa7qK1GKG+fuSqWsAQV1S7XQqBppjGQr6AC9ik0l3\nMnLkSGbNgooKKCsrpbS0NMWPEhHJbWVlZZSVlQX6GS43+W+BxkCtv12IDXMGCyb71HNeP2Ak1s8C\ncKf/HrGd93/GaiXj/e3FwCDgoATnvgk8AEzxX1sOHA18Gff5nud5/PWv8Pbb8OyzSX9PEZG8F9ao\ns73944r9RyHQxH/UF2QAPsKauUqARsC5wBtxx7wBXOQ/74ctQ7AuybmvASf4zw/2X48PMt9T05mI\nSLhcms52VTVwPfAWNorsKWxI8lX+66OBidjIs+VYLenSJOcCjPEfkSwFkUBVJ2UHEBEJV1qrRxnG\n8zyPCRNg9Gj4xz/CLo6ISOYLq+ksq6npTEQkXMkCTTGWrTlrqelMRCRcyQJNNTYS7MAGKEsgGjVS\noBERCZPLYIAWwMfADHYc1jw0qEKlk5rORETC5RJo7vF/RiZgFuA+GTN0ajoTEQmXS6Apw9LQ9MEC\nzAws91lWUNOZiEi4XEadnYMlvBzmP5/hP88KajoTEQmXS43m/2G1mUgtZl/g/4CXgipUOqnpTEQk\nXC41mgJ2TMX/JVk00VNNZyIi4XKp0UzCUsGMwwLMuVhiy6ygpjMRkXC51EwKsEXPBmKDAd4FXg2y\nUGnieZ7Hpk3QoQN8803YxRERyXxBpKDJmiawXeB5nse2bdCiBWzbFnZxREQyX0PnOnvP//ktsDnu\nkTX1A/XRiIiEK+drNACFhbB8udVquncPuVQiIhksjOzNxVius6zWqBE89RT84Q9hl0REJP+4JNVc\nQhYn1QQLNOXl8PXXYZdERCT/5HxSTbBJm+XlsNdeYZdERCT/pJJUM1bWJNUEq9GsWgX77x92SURE\n8o9rUs0SoDPwL6Cx43kZo1Ej+PRT2GOP4D+rosLm7YiIiHFJQXMlltdstL/dnuyYsPm9Qv+3DLqP\nZvVq6NhRfUEiIrFcAs11WFaAyNyZpcB+gZUoAF/4mdo2bYLa2uA+58UXoboa5s1zP6e21solIpKr\nXAJNpf+IKCbL+mi2bIHOnWHPPWHz5mA+Y8oUeOABOPpomDPH9i1YAF9+ac+rq+3n9u2wdm30+XXX\n2TlBBkARkTC5BJopwN1Y38zJWDPahCALFYTjj4fmzdPbrFVZCZ5nedTOOguefx5GjIBZs2z/8OHw\nzDMwfrylwZk6FQYMgDPPhK++gtJSGw3XqBFMmmTvuXkzXH55dFtEJNu5BJrbsWUC5gNXAROxNWqy\nxksvWW0jnYHm2WdtcMG4cfDmm9CvH5x4Ipx8sm2XlVmN5oMP4LbboHdvC3bnngsLF8Khh0KfPjBx\nItx6Kzz2mDWhDR9uk0ufeCI95RQRyQY3Oe7LNF68k07yvDff3Gl3yjZs8Lx99vG8n//c87p29bzu\n3T3vySejr/fr53l77+15DzzgeeB5vXt73qpVnvf00/b6oEGed+GF0eO3b/e8/ff3vJ49bf/q1XZe\nr16et2SJ573++g8vs4iIC0LqGpldx745DV6K1O10AUeM8LzRo3/4P8TLL3vekCGet26dBYQ//tHz\nqqujr1dUeN7nn3teTY3nnX++5y1evOP5M2d63po1O+57+GHPa9vW87Zts+2XXvK8wkLPKy62z6is\n/OHlFhFJhgACTaLEacOB84FjsTVoIpoANcCJ6S5MmvnXLOoXv7BMzq1aWRNV69apv2llJVxwAfTv\nD7fcAlu3QuPGP7ywNTWwbh20bRvdd8klsM8+MG2aDRZ46y1o2fKHf5aISH0aej2aA4GDgAewfprI\nsd8A87A8aJlsp0DzzDNw//3wySfw9NN2I09FTQ2ceioUFVkfzX4BD/KurLT0Oe+/Dw8+CO3awe9/\nb/1NzzwDr7wChxyi1Doikj4Nnb15FZYV4CTgP/7zz7AJm66FGIxlf16GBau6jPJfnwscmcK5twC1\nWC42JwceaEGmRQsb7ZWqt9+20WITJwYfZAB2390mmw4caIHt9ddtmYP58+HYY21Y9K23Bl8OEZEf\nwiWVzBSs+aw58BbwIXAucEGS84qAx7FAtcY/7w1gUcwxQ7DUNl2Ao4EngH4O53bAhlqvcij/9w70\nc1Dfd5+NBkvF5s1WG7rqKqvRNLRmzeBvf7M1dU47DZYssXk4L7wADz+sWo2IZC6X4c2FwFbgbOBP\nwDDgMIfz+gLLgXKgChgPnBF3zFBgrP98OtAMaONw7u+A2xzKsIOOHWHFChuKvGSJ+3meBz/+sTVT\nXXZZqp+aPieeaEEGrCwvvWTzcl5+ObwyiYgk4xJoAPpjNZh/pnBeO6AiZnu1v8/lmLYJzj3D304h\n0UtUx47QtSssXWo1AhfLl8OiRfDkk1CcYelEr7gCRo2KZh6oz2efWYA999xoFoLt222i6ZxsGEMo\nIlnLJWD8N3AnlkjzY6AT8I7Dea5D5FLpdNoTuAu4bxfPB6BJE+jWDaZPdzt+1iybXBlGk1kyp59u\nfU6jRtV/THm51ciOPdYCzqBB8POf2/bs2TB4sDULiogEwbWPZgo2rHlvYAVwo8N5a7C+lIgOWE0k\n0THt/WN2q+fcTtiSBXNjjp+JNbWtjy/AyJEjv39eWlpKaWlpzLbN3h80KPkvMnMm9OqV/LgwFBbC\nb35jNZMXX7Sg4XkWfG66yWpiJ58M11wDt99u2QfGjIFXX4VLL7WUOePH21DtYcMsY4GI5I+ysjLK\nysrCLgY9sEmbn/qPmbj10RRjQakEaIRN8uwad8wQLKUN2CCAaSmcC7CS+kedJZyUNGaM5118cfLJ\nS5s2eV7Hjp43ZcoPnwgVpFtv9bxHHvG8Jk08b7/9PG/ffT3v0089r0MHzxs7Nvn5Tz/teQce6Hlb\ntgRdUhHJZAQwYdOlRvO/wM1Em8tK/X0DkpxXDVyPjVQrAp7CRo1d5b8+GgsyQ7CO/y3ApUnOjbfL\nF6RVK9iwIflxf/2r5Sk77rhd/aSG8dBD9rNNGxsOff/91jx4xRVw0UXJz7/kEsvR9sgjcE9da6pm\nsOrqzOs7E5Eol/6NuUBPh32Zxg/OdXv/fWsuSjbM+fLLLdBcc02aSxewykrLKHDccVDg2Iu1ciUc\ndZRlIDj4YMtKkMlmzLCEpe+9ZwF1yRK44QY45hjYd9+6z9m2zbJoH3SQnXf22dC0afLPWrUK3nnH\nVlC98kprsiwutmUgCgttkMmuWLHCVn8dNCi6QJ9ImBo6M0DEa1hz2bP+8RcAvYGz0lmQACQMNEuX\nwn/9FyxblvhN+vSxzMoDktXfcsRf/mKTQBs3tn6pyLIGTZtaZoKFC6FLF8tYkC61tdavVFRk2Rf+\n8x/bjulS+96HH1ptrUsXm8T68MNw2GG2RMMRR9ggh3XrLFAOGmTD0V991WpskybB//wPlJRYUO3e\n3VIS3XWXpSO6+24bOHHooXDzzfa7//vf8NxzMHmyDS9v2tTmLnmeTagFK3OTJlYTvPzyxL/r0qUW\nWL79Fv70J/joI0s7tGED7L+/vdfFF8PVV9tnt2kD69fbJOH16+GXv7Sh7Yls2BD90hD/JaOqyj6j\nIZY1l+wUVqBpDvwCOMbffhcYCWT6gsUJA82XX9rN6quvdty/dKktyXzCCVYraNECPv/cbiT5ZPp0\nq+29+qp98+/c2W7Q//633YT79IEePeyGvmCB3TSnTrXsCfvvv+N7LVtmNYKTTtpx/+rVlgpozBi7\n8R11lN1QO3Sw0XFPPGHX/ogjbGG5ceNs7Z8rrrAJtCNG7HzTXbnS/k0/+sjKPm2azT36+9/t8x98\n0JoUwYLFPffYe5eXwx13WC65xYvt9ykosKAzdKh9ZqSGt2mTlffjj61mtGWLlffiiy0w/PrXFsT2\n3tvSBI0aZaP7DjjAFr3r2NE+Z9gwuPBCC1iffGK1pYICuPdeK3f//hZc2re350VFVq5evWzS8bvv\n2mqup59uZZ41yx4bN9qSGFu32rlNmtiE3ooKC3K7726pjIYPd6/tSv5o6ECzJ3A1NnN/HjAGmzyZ\nLRIGmpoa+4Pbvn3H9v1zzrEULwsX2h/z/Pl2w8pXkewD8+bZDXLoUGuO277drmG/fnZTv/JKu/lP\nnmxB+sEH7Wb8zDNw550WtMeOhX/9y2pGU6da8+V551mtY/JkqyWdd57dHKdOhVNOsfeaNAnOP9++\n5Q8Y4N7E5HlWxuJi+5nK8PQFC6zMvXu7n1NVZbWfn/3MbvJHHGE1l1/9Co480pYUP/xwW+guWbk3\nbKi7+a+qyvrRnnvO3qtHDwuUhx9uZe3Vy/6damosqKxda2XYts2+ALRvbzW+Sy+FTp1s2PuAARb4\n5syJDuXfc08rv/q+ctP27Vabr6iwL3PNmlltvWNHKCpq2EDzIvAdludsMJbuJRvWoYlIGGjABgQs\nWhT9g/70U/vj2mcf++Z44YV2g23fvgFKm0U2bbIb6QsvWG3nuefs5ul5FpSffNKCSU2N3QTHj7ea\nxnnnwY03Wo3j5JPhJz9JnPk6khm7osJqOdli61ar4SxaZMEyWWAJQ2Ul/Pa39v976lQrc/fu9v9/\n+nRbILB5c/jRj+CnP7VJzmC1pQUL7HdbuNBqS6edZrVR9TGlbssWu9F//rnVYmfPtvvRGWdYk+oB\nB9h1njPHvnxs3WpfCvr3ty/K1dW2fcAB0fdcu9a+TGzdav9W06bZ32njxtYqsWiRtTAUFNjfaYcO\n9je5caN9iRk2rGEDzXxsaDPY6LQP2THpZaZLGmgOOcQSVUbmjlxzjf1xtW5tzSiXXQaPP94AJc0x\nmzbZN+uBA+HRR635Eax5aFc7zSU4tbV204ltRvM8m5f14YdWK23e3JotKyvt76VbNws+X3wBEybY\nTeqUUyzwFBZak+qpp1rztOfZQIpx4yzorllj/yeuvdaaGr/4wj6zvNwGcqSyfEdtrQ0KmTPHasgX\nXWR9afFNguvWWRmmTLGAeuaZ1kz+8cf2u/Xvb5+/cKFlAikshCFD7AvSRx9Z8PU8++Y/YoTV3ufP\nt1rBoYdazXL9emtiLiqyoLBli73n11/b//tVq+y9iouttv/++/Y+rVvbdTj8cCvHypV2X1q3zt6z\nSxerEbdpY9e3vNy+IGzfbsFmzRr73E2b7PeOLDdSXGy13J49rRWistJ+5x497BoUFe18nTwPCgsb\nNtDMZsfAEr+d6ZIGmmOOsSaegQNtu2dP6zPo2tWGC992mzpNd9XWrdb8oj6A7PfFF3bD22sv63dq\n1WrnY+bNs9GKGzdaTXb1ahsu36qV3WhbtrSmz9pa+xa9YIEFnshNsaDAWg5mz7Zm2CFDrFbYu7c1\n8UVs22Y36MJCC3ATJtjNtk8fu3n+5S/2jb51a/u8jRvtpl1RYYNDjj3Wgs0HH9jfeffuFhQ++MDK\ncdRRdtP+6iur6XXqZPsimUE+/dSagNu0sRt2VZV9gVqwwM5futR+x1at7P9/x452D6mosPft2TPa\nXN+rl9UwEv2NVFUlH3gTuSb77WfBpGfPHzZYp6H7aGqwZJoRewLb/OcekOGDX5MHmrPOsmaBs8+2\nSN60qX3raN68gUooksNqauxbf6tWFlzquqFWVlotJ/JaRYUtUDh3rvUpzZhh53fpYjfod96xG3xh\noTXrDR4MfftGz/c8CxKLFllNrE0bq1kNGJCZKaQyUVijzrJV0kBz9dVWXb32WvvP2bGjfQMSkcxQ\nXW19DIte40tVAAAHDklEQVQWWW1q6FAtiRG0IAJNXo8pad3aqu5g7aIlJaEWR0TiFBdb03akeVuy\nU16PE2nTJhpoyssVaEREgpDXgaZ1axtWCDYR8chsGuogIpIl1HS2zvpnnn/eOiBFRCS98r5GU1Fh\n+al++tPsmhQoIpIt8nrUWVUVXHCBDXscM8bGvYuI5DMNb05N0kAjIiI7CiLQ5HXTmYiIBE+BRkRE\nAqVAIyIigVKgERGRQCnQiIhIoBRoREQkUAo0IiISKAUaEREJlAKNiIgESoFGREQCpUAjIiKBUqAR\nEZFAKdCIiEigGiLQDAYWA8uA2+s5ZpT/+lwgdp3L+s59CFjkH/8K0DS9RRYRkXQJOtAUAY9jAaMb\nMBzoGnfMEKAz0AW4EnjC4dzJQHegJ7AUuDOw3yAHlJWVhV2EjKFrEaVrEaVrEaygA01fYDlQDlQB\n44Ez4o4ZCoz1n08HmgFtkpz7NlAbc077IAqfK/RHFKVrEaVrEaVrEaygA007oCJme7W/z+WYtg7n\nAowAJv7gkoqISCCCDjSuS1zu6mpudwPfAeN28XwREcly/YBJMdt3svOAgD8D58VsLwZaO5x7CfAe\nsEc9n70cC3R66KGHHnq4P5aTZYqBFUAJ0AiYQ92DASJNX/2AaQ7nDgY+BloFU2wREckmPwKWYFEy\nMjrsKv8R8bj/+lygV5JzwYY7rwJm+48/BVFwERERERGRBucyQTQbjQHWAfNj9rXAhnovxeYWNYt5\n7U7sGiwGTonZ39t/j2XAYzH7dwde8PdPAw5Mb/HTpgPwDtZ0ugC40d+fj9diD2x4/xxgIfAbf38+\nXouIIqyVY4K/na/XohyYh12LGf6+fL0WaVeENbOVALtRd59QtjoWy5oQG2h+C9zmP78deMB/3g37\n3XfDrsVyoiP7ZmBzlMD6xgb7z68l2gR5LjZvKRO1AY7wn++NNa12JT+vBUBj/2cx9gc/kPy9FgA3\nA88Bb/jb+XotVmKBJVa+Xou068+OI9Xu8B+5ooQdA01khB7YDXix/zx+hN4kbKDF/ljqnojzsFF/\nkWOO9p8XA1+kq9ABew04CV2LxsCHWMaMfL0W7YF/AccTrdHk67VYCbSM2xfatci1pJouE0RzSWus\nOQ3/Z+Q/UVvsd4+InQQbu38N0esTe+2qgU3s/I0o05Rgtbzp5O+1KMS+ja4j2qSYr9fiUeDnRLOG\nQP5eCw8Luh8BV/j7QrsWxbvyG2QwL+wChCgyBj5f7A38HbgJ2Bz3Wj5di1qsKbEp8Bb2bT5WvlyL\n04D1WJ9EaT3H5Mu1ADgG+AzYF+uXWRz3eoNei1yr0azBOosjOrBjRM4167AqMFg1d73/PP46tMeu\nwxp2zAsX2R855wD/eTF24/oq/UVOi92wIPMs1nQG+XstIjYB/8Q6b/PxWgzA8iauBJ4HTsD+f+Tj\ntQALMmBNWq9i/Sz5ei3SzmWCaDYrYefBAJG21TvYuXOvEXAQdk0inXvTsbbVAnbu3Itkzj6PzO3c\nKwD+ijWTxMrHa9GK6MihPYGpwInk57WINYhoH00+XovGQBP/+V5YBpVTyM9rEZj6Jnlmu+eBtVhu\ntwrgUqxN9F/UPVzxLuwaLAZOjdkfGa64HFsHKGJ34EWiwxVLAvgd0mEg1lw0h+iE3cHk57XoAczC\nrsU8rH8C8vNaxBpEdNRZPl6Lg7D/E3OwKQCR+2A+XgsREREREREREREREREREREREREREREREUnu\nW//ngcDwNL/3XXHb76X5/UVEJAtE8qeVEp1l7ipZDsH43GwiIpKHIsFgGrARy0BwE5Yf8CFsnY65\nwJX+caXAu8DrRBMYvoZlz11ANIPuA1i229lYLi6I1p4K/Peej830PyfmvcuAl7C07X+LKecDWMbm\nuf65IiKSJSKBJjZvFlhgudt/vju2BkwJFgy+ZceVBpv7P/fEgkdkO75GE9n+MZYapADYD1iFJUEs\nxYJdW/+197HsvC3ZMSvvPq6/nEhDyrXszSLpVhC3fQpwEVYjmYblj+rsvzYDCw4RN2H5pj7AsuN2\nSfJZA4FxWPr29cAUoI+/PQPLdef573kgFny2A08BZwHbUv3lRBqCAo1I6q7HFlw7EuiEJSoE2BJz\nTCmWSbkftl7MbGCPJO/rsXNgi6wZUhmzrwZbKqEGS//+MrYeyyREMpACjUhim4mmXAdbXOxaoh3+\nB2Np2ePtA3yN1TgOxQJORBV1Dxh4F1t/vRBbsOo4rCYTH3wi9sIy8L4J3Az0TPrbiIQg11bYFEmX\nSE1iLlZzmAM8jaVKL8HS8xdgTVxnsfOKhZOAq4GF2LIVH8S89r9YZ/9M4Kcx570K9Pc/08PS/q/H\n1lSKXw3RwwLg61hNqQD42S7/tiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIikrn+P5cM1q5B\nVgayAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x114ae5350>"
       ]
      }
     ],
     "prompt_number": 68
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that this distribution seems to converge more slowly than the posterior in the previous section. "
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "USPS digit dataset example"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this section, we use our Gibbs sampler to learn clusters for digit recognition. Our digit dataset is the [USPS digit dataset](http://www.gaussianprocess.org/gpml/data/). We'll just focus on the digit 2 for this example."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "usps_digits = sp.io.loadmat('data/usps_resampled/usps_resampled.mat')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "usps_class_2 = usps_digits['train_labels'][2,:] == 1\n",
      "usps_Y0 = usps_digits['train_patterns'][:,usps_class_2]\n",
      "usps_Y = np.zeros(usps_Y0.shape)\n",
      "usps_Y[usps_Y0>0.1] = 1.0 # make the image binary\n",
      "plt.imshow(usps_Y[:,3].reshape((16,16)), cmap=plt.cm.binary)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "<matplotlib.image.AxesImage at 0x112abe950>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAAD7CAYAAABOrvnfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnW2MLNlZ3389/VZVPTsz966FDXiVJZaQSJQ4kAQRB7Sb\n2JYcxzGKEiVBJDHeJF8CiQnEeA2K2f0QyQsiJi/iA0mMbIEdgokckAjBlrIkUYgBYxsHMGDz5rXl\nBfnembm3q/p98qHmqTl9uqqn+3RVdVfV85NKXf0yU9Wnz7/OOc85/3pAURRFURRFURRFURRFURRF\nURRFUZSq8sQTT1wBuummW8nbtfZSaWW9kQNXV1fpx33mmWd45plnCjy0Hm8TwjDk3r173Lt3jy9+\n8YvJvr3Z70VR5HSeQRBw9+5d7t69y6OPPprs25v93nPPPVeJ8jyUY7VaLcjQ9tEO//d1wKeA3wbe\ntsP/URSlJFwF3wb+HbHo/wTwTcBX5XVSiqIUg6vgvxb4NPB7wBT4T8A3bvrHTz75pONh3dDjVZs6\nl2fZ381V8F8OfNZ4/sL1axtR5x+wCccrmzqXZ9nfreP4d5lRQGU3JNB5dXW1tH/bey6Mx2MmkwnT\n6ZTZbMZisVj6f61Wi3a7TafTodvt0uv16Pf7zsfs9/v0ej263S6dToejoyMJMCXfZT6fM5vNmE6n\nTCYTxuMxrVYr+Zy9bz/arynLuAr+c8BjxvPHiFv5Jczo45NPPln7lmhXrq6uWCwWyWY/z3rPVYCj\n0YjLy0uGwyGj0YjxeMxsNuPq6oqjoyM6nQ69Xg/f95lOp8zncyC+ULjgeR4nJycMBgN836ff79Pp\ndGi1WiwWC2azGZPJhCiK6Ha7tNvt5O+Ojo6WtlartfKa/XpTRP/888/z/PPPb/RZ1xLpAL8JvBr4\nPPCLxIG73zA+kzktp6SzWCyYz+fJNpvNlp5nvb5YLJyON5lMePjwYbINh8Ol52mvDYdDJpOJ0/F6\nvR7Hx8fJNhgMlp6nvTYYDPA8j3a7nfQ2ZN/c0l4/OtplEqq6rJuWc23hZ8C3Af+dOGL/H1kWu+KI\niF662Zs8ugp+Op0SRVGymS18q9VaauEBjo6O6Ha7zGYzp+N1u1183082s4W/urpaauGlLKbTafI5\nGQrIo7m/WCzodrvAzVBEWcVV8AD/7XpTckK66dKKTyaTZJMxbdpr0tXelvl8nozj5dHu0vf7feBG\n7J7nOR+v3W7T7/eTsXyv11vp0stwQcQ+Go2Sz8om8YRer8d8Pqfb7SbDGunSa+8ynV0Er+SMCN4O\nWo3H42SMnba5trhyLHtbLBZJCw/LLfsuQwi5iNibCFS+h4hduukSLOz3+3iel+zb5yNib7fbKvgM\nVPAHhNnCm4KPoojRaLTyKNt0OnU+phkAlOi/2cK32+2kBZXP7YIZWJOoutnCy3eX182ehe/7TCYT\nfN9PPmu27DKWV7Fno4I/MOwuvbTuYRgShiFRFK3su7bwwJKwTAFKC5/23i6YFw7zAiMtfNp7MvY3\npxBF7HJOcmGSFl9Fn44K/oDI6tKbgpeouewPh0PnFl66v9K1tvez3nMVvXkxE9Ga+1nvdbtdxuNx\nMjVot+zSG5Gx/a69kDqjgj8g0rr00oUPwzCZGnvw4MHSVJnrNJk5PpZHuOl2S9DO/Eyv13OOgM/n\n86UAoRmgkxbeDiJOJhPa7XYSULTH7ObCoH6/ry38LajgD4x1XXpp3R88eMDl5WXy6Cp4c2wcBAFw\nE6Czp+XMTYJ52yJRd7mAyfedTqcr03Lm1m63N2rZpQeggs9GBX9ArAvaRVGUCP7y8pKLi4vk0XXl\nW6/XYzAYJFF5Ebs9Lef7/tKiGJnv3pbJZMJwOEyi8iL2tJV25qKfdrudrCo0x+xyQTKj9jLLoKSj\ngj8g0sbwZos4HA6TVv3i4oLz83POz88ZjUZOxzPn1E2Bm9NyvV6PIAg4Pj7m5OSEk5OTpOu/LaPR\nKBkOmPPu5rTceDxOhi+Xl5dcXl5ydHSUtNqybNYWvOd52sJvgAr+gLAX3kyn05UuvS34+/fvO9+B\nxlxB1+v1kgtAVgt/cnLCnTt38DzP6XjmCjq5mLXb7ZWFN9LCX15ecv/+/WSJrLTsEpGXc5OAnjnG\nV8Gno4K/hbLda+Ym01Dmlra23rULa563aTqRaLwE6aQFNcfxrsjaAXu1nSzAMU0vcn5pXgK7XNIW\nJKnLbhUV/Brq7l7zfZ+zszNOT0955JFHknG67/vJijZZyipmlF0EYa6EM+22chEZDAZLFzazHNVl\nlw8q+DWYXeyy3WthGCZTUeaYut/vEwRBcoxOp+Mcpfc8j9PTU05PTzk5OVkRvQhe5uGLErzneQwG\ng8ylsldXV4l7zvf95JzSxv7S/Z/NZlu77OSYKvgGU2f3Wr/fT4R+cnKS2sJ3u92lFn4XTMGbATdZ\nRZc1z75YLNRllxMq+DXU3b3W7/eTllNa9yAICIKgkBZezlsEJ4E36cabwyFznn0+n6vLLidU8Guo\nu3tN5uHtLQiClTF82i2ptsXu0puOtzQjjAh0Pp+ryy4nVPBrqLt7rdvtJq25tOzyPKuF3wW7S28v\nhwWWZgpE8LPZTF12OaGCv4U6u9fEhSaiMB/TovR5BLTSjC72mN1s3WVBjbrs8kEFv4a6u9c6nc5S\nV9fu+krQrqgWvtfrpY7Z5T3P85JYhrrs8kEFv4a6u9fMyr5uK2oeXsQn3Xhzqs4OiKrLLh92Efxj\nwHuBLyG+T/0PA/8mj5M6JOrsXjMrvT11Zb5mBsny6NKL4GG5Gz+bzRLXmznVKbERddntzi6CnwL/\nDPg4cAx8FPgQNbp7bd3da2b3Om0hir1YJa8uvVygROzS/ba76LI/Go3UZZcTuwj+C9cbwENioX8Z\nNRR8Xd1rWctM1215dOllXwTZ7XYzlywvFouktQZ12e1KXmP4x4GvBj6S0/87COruXgNWov1p+2nv\nuR7LnF6DVeORvX91dZV8P3XZ7U4egj8GPgC8hbilTzikVFMurjdzbl268ea+TMWZU3PymisyJWYG\nn8wup70azRSR+Zrs24/7doW5XjSiKFoJIApZgdVWq5V022W6MQxDPM8jiiJ8309+P9/3k9+6ai67\nbVJN7Sr4LvCTwI8CH7TfNAW/T1xdb2EYcn5+zsXFRRKJHw6HifBNC2seLURWzEAq8HA4XFrbLpVs\nNBqpKyyFppSn3Zg+++yzmZ/dRfAt4hRTvw784A7/p3BcXW9hGHJxcZEE5GzRm+vdixS8tG72fLh8\n3vd9dYWloOW5yi6C/4vA3wV+FfjY9WtvB35215MqAhfXWxRFS9NtaS286d/OI8qbNg24roLO53OC\nIFBXWAZansvsIvj/DVQiPaer6y2KoqRVl0dznF5UC2/e0WU8HifdTvtuMPJZ2xmmrrAYLc9VGrHS\nztX1Zi6XNbcwDFfG8BJUy7NLPx6Pl7qOZuU0u6lidFFX2DJanqs0SvDbut7SDDJmZN5u4Xf90bOm\nAc1upwQU7bUB6gpbRctzlUYIHtxcb/ZUnH1RyDtKL+dpVj6JAMNNS2RPQYmlVV1hq2h5LtMIwbu6\n3szgnLmZr9kLN/I4T/PClDbGlPfk7q/muagr7AYtz1UaJfhtXW9hGKYG9dKCfEVMy0nllG6nebEy\ng0pRFKkrLAUtz1UaIXhwc71FUbQU4bWn7WTfXAmXR5deKigsdzvlDrX2dFEQBOoKy0DLc5lGCN7V\n9RZF0cqCnKz9PLv0csccqZydTmfpXm32wo/xeKyusBS0PFdplOC3db1FUbTWxZW1NHfX85T9+XzO\n0dER0+l0rYNtNBqpKywFLc9VGiV4F9eb6dqS/5Xm6LLf2+U8zYq6iXtN/PfqCltGy3OVRggeVkVv\njsPNWyiZLf8urrddzxXiCmdXJHNqSPbtHoZZ0fO6GO0b06Ri3s/PXPJq58Jbd6EGllp/+z2Zhzdd\ni3bcpmpihwYJvgrYjqx1N6gw3wuCoNQccWVj3wtv05x0nuc5uSTrjAr+gDAr9qZOrXa7TRAEpeaI\nK5sswd+Wk06Crps6JGG37L9VQAV/YJgVO82xlfbo+36pOeLKxhT8NjnpJEeAPaVq7ksQD24Ce3VG\nBX9AmHO9UrGzXFvmc7nHXZk54srGnAPfNCed53mZzsh2u720GEe69FUrl21RwR8QWV3XtAQR5hYE\nQak54srGLpdNc9KJwcleDp1144v5fF6pcnFBBX9ArAtOZaWD8jxvqTUvI0dc2dhd+k1z0nmel8y6\nyAVAysC2xs5ms9qLHVTwB0fWIo80Ucu+CLvMHHFlk2ZYscfsdk46mVo1ezfy/WUqTlZfVnWosy0q\n+AMibb223YpL8gmzG28G58rIEVc2aUG7tDG7nZPO/O5pN71Is87WHRX8AbGuS58WmJMtCIJSc8SV\njV0uInbpxmflpPN9f+2969JuilGlcnFhV8G3gV8GXgD+2u6no6zr0kvrLtNv8ij55crMEVc2puBh\ns5x0YRhu1LJX9ULowq6CfwvxbaofyeFcGs8mK8pE5LLI5vT0NLmtclk54srG7LbL801y0vX7/aQX\nYI7ZTYu0RO3lQlh3dhH8y4HXA/8S+I58TqfZZI3hpUsvLbwI/ezsjLOzM3zfX+v+ylqaWxWkXGR/\n05x0nucl31VeS7sngrbwm/Eu4K3ASU7n0niypp/MLr3Zwp+dnXHnzh183890fpn7ae9VASkXU/ib\nOBglyaZ97zq5F0JaMLNK5eKCq+DfAPwhcQKKJ7M+dEi55VzJEtBt77lgzq3bU232FJw93153XC9S\nMg8vU5MSvDTjGlUXexm55V4FvJG4S+8Rt/LvBf6++aFDyS3niqt7zbXi+L5fa9ebUgx2Y1pEbrnv\nvt4AngD+OZbY64Cre801+ON5Xq1db8r+yWsevraeQhf3mqvg+/1+rV1vyv7JQ/A/f73VDlf3mmti\nwX6/X3vXm7JfdKXdGlzdazJfvC29Xq/Wrjdl/6jg1+DqXpPVYNvS7XZr7XpT9o8K/hZc3GuuLbzc\nvabOrjdlv6jg1+DqXnNt4eUuLXV1vSn7RwW/Blf3Wq/XczqefbvlurnelP2jgr8FF/eaq+DNJbV1\ndb0p+0UFvwZX95qs4d4W8+JSR9ebsn9U8Gtwda95nrfT8erqelP2jwp+Dbu413Y5Zl1db8r+aYzg\nXXKTtVotda8ptaIRgnfNTdZqtdS9ptSKRgv+ttxkrVZL3WtKrWic4LfJTQaoe02pFY0QPLjlJmu1\nWupeU2pFIwTvmpus1Wqpe02pFY0T/Da5yVqtlrrXlFrRCMGDW24yQN1rSq1ohOBdc5MB6l5TakXj\nBL9NbrJWq6XuNaVW7CL4M+A/AH+S+CaWTwH/N4+TKgKX3GSAuteUWrGL4P818DPA37z+P4NczqgA\nXHOTAWsda+peU6qGq+BPgW8A3nT9fAZc5HJGBeCamwxuuv3qXlPqgKvgvwL4I+BHgFcCHyXOJBvm\ndF654pqbTP5W3WtKXXAVfAf4GuDbgF8CfhB4GniH+aFDyi2nYqw2WRfnde+NRiPG43EShJXYjLnJ\nWgxztWXV2Ca3nKsCXgb8AnFLD/D1xIJ/g/GZq6oWoHJYXF1dLQ217OdZ74VhyPn5Oefn59y/fz/Z\nl+fymv3e0dERd+/e5e7duzz66KPJvr3Z7x2KLfq6YUvVtmsL/wXgs8BXAr8FvAb4Ncf/pShrERGn\nBVbNzX49DEMuLi64uLjg8vKSBw8e8PDhQ4bDIVEUJa2/2dLXvZHaJUr/T4AfA3rAZ4A353JGipKC\nCN6eOl33GEVRIvTLy8slsUt337wPggRq68wugv8E8OfzOhFFycJs4Wez2dLiqOl0uvTcfC2KoqRV\nl8cwDAnDUFt4RTlURPDSck8mE8bjMePxOGmp07YwDBkOhytbGIZLAT3pFajgFeUAMFt4U/DSNbcf\nZTNbc9mX52ktfN3FDip4pSLYXXpp3dNELfsi7KyLgtnCa5deUQ6ErC69KfjhcJiM0aXrbgbnzM18\nTbrzGrRTlAMhrUsvrXUYhjx8+HApMCdbGIapQb20IJ+28IpyQKzr0kvrLtNv8hhFUdIrSJu2s1fe\nqeAV5QBYF7SLoigR/OXlZbLI5uLigiiKVhbkZO1rl15RDoS0MbzZpR8Oh0mrfnFxkSyRjaJorRsy\na2lunVHBKwePvfBmOp2udOltwd+/f58oilbcj+uckeZ+XVHBK5XAFr05Dpfgm1wEzGk5ZRm9RYui\nNAgVvKI0CBW8ojQIFbyiNAgVvKI0CBW8ojQIFbyiNAgVvKI0CBW8ojSIXQT/duI71X4SeB/Qz+WM\nFEUpDFfBPw78I+JkFH8KaAN/J6dzUhSlIFzX0l8CUyAA5tePn8vrpBRFKQbXFv4e8APAHwCfB86B\nD+d1UoqiFINrC/8K4NuJu/YXwE8A30ycmCLhkHLLVQWXHGp1Z5cccVlJQNe953ke/X6fXq9Hr9ej\n2+3S6XSWNkkPfggZg8vILfe3gdcC//D6+d8Dvg74VuMzmltuS1xzqNW9nKMocsoRNx6PV1J5r0vz\nLftBEHB2dsbZ2Rl37txJ9uW5vGa/53nevosKKCa33KeAfwH4wIg4t9wvOv4v5RrXHGp1vzXTaDRy\nyhEnIm6327Tb7aRltjf79SAIOD095fT0lJOTEx555BGOj48ZDAb4vp+0/mZLv+9WflNcBf8J4L3A\nLwML4FeAH87rpJqMSw61ugt+PB4754gTwUu3fJNH3/cToZ+cnCyJXbr73W6XbrebCL4q7HLHm++7\n3pSccM2hNp/P933qhTKZTJxyxEkLL+NuGZObY/O050EQcHx8zPHxcdK6B0FAEASNbeGVAnDNoTab\nzfZ96oUymUyccsSZXXoRc7/fp9/vJy112hYEAYPBYGULgmApoCe9AhW84oRrDrXpdLrvUy+U6XTq\nlCMuS/DSNbcfZTNbc9mX51ktfFVQwR8YLjnU6t7CS673bXPEmYKXLr207mmiln0RdtZFwWzh2+32\n0rTeoaOCPyBcc6jVvYWfzWZOOeLMMby08HYrPhgMkjG6dN3N4FzWEEC689rCK8645lCbTCb7PvVC\nMS+A2+SIW9elTwvMyRYEwUpAL20zo/TawitOuORQq7vgzQQU2+aIW9ell9Zdpt/k0ff9pFeQNm1n\nr7xTwStOuOZQG4/H+z71QjEvgtvkiFvXwvu+z2AwSEQui2xOT0/xfX9lQU7WvnbpFWdcc6iNRqN9\nn3qhrFtivG7pcdYYXrr00sKL0GWJrO/7mUtw1y3NrQIq+ANilxxqdcclR5wdpe92uytderOFl/Xx\nvu8nkfcsg03We4dOYwRfBReaHYWWQJS5pa2tN8esh/z99oEtenMcLsE3uQiY03J1pRGCr4oLbTQa\ncXl5yXA4TKacZrMZV1dXSde01+vh+34ifvm7Knw/Zf80SvCH7kKbTCbJVFsYhsnqscViQavVotPp\nJEs/5dw6nQ6j0Wir7yZlooJvHo0QPFTDhSYryswlo9LCi+ClhYd4yqnb7Safs7+DuX90dJQs0Lm6\nuqq94UZJpxGCr4oLbT6fJ2N3c3242aXv9+ObA4vYPc/LPP/JZEK73WY6nSZBJSmLqgSZlHxplOAP\n3YUm52hvZpcebsQu3fSs72Mv/TQvfCr4ZtIowVfBhWYG1szpJmnhZRGJvL5YLJhOp8n593o9oihK\nPmu27NLDUbE3l0YIHqrjQpM5XVnMYW6dTif1PRn7mx5tEbtcGORiV7W130q+NELwVXGh2fPF9n7W\ne9PpNHFwyefsMbv53au0FFTJl0YJ/tBdaO12O1kIIo9A0iJL0M78TK/XYz6fp96QwezGy6o9beGb\nzW2CfzfwV4E/JE4pBXAX+HHgjwG/B/wt4kQUB00VXGgSdZcVX3AToLOn5cxtPp9v1LJX0c6p5Mtt\ngv8R4N8S36FWeBr4EPENLN92/fzpQs4uJ6riQuv1egwGgyQqL2K3p+V830/upHp8fMx8Pk9MHOaY\n3by4SdRe7JxKM7lN8P+LOLuMyRuBJ6733wM8T0UEf+guNM/zkrl/U+DmtJx5V9WTkxNOTk5YLBZJ\nqy3LZtN6M9rCKy5j+JcCL17vv3j9/KCpigvNXEEnVk65KWNaC39ycsKdO3eSFYHSsktPRnoxabdk\nUsE3k12DdlfXWypF5JZzcYWZc+vSjTf3zTug2ndFdclNtst3s++4mnZfddPd5fs+V1dXSbddvlsQ\nBIxGI3zfJ4qi5LNRFCX3ddtlLX1V3HmmWUgu+OZmrkqURmCdLdZ+zON335Vtcsu5CP5F4GXAF4Av\nJQ7opWIKPg9cXW9hGHJ+fs7FxcVKmqJN7nq6bW4y1x/f933Ozs44PT1dSW+UdrdU81i33dnFtNeK\n+FxtoFVxH2bFbuRiPxwOlzLISFmORqO1v2/ev/uu2I3ps88+m/lZF8H/FPAm4Lnrxw86/A8nXF1v\nYRiWmpvMNSjmeZ5zTrM0wXuex2AwyLybq2tQctvfQH67QxG8LFLKWnost7ja9DeH6twE4zbBv584\nQPcS4LPAO4B3Av8Z+AfcTMuVhovrLYqiUnOTuQq+3+875TTLuv+67/srFzL5fLvddha8GQs5dHde\nWgBzneDn8zlBEKT+tub+YrGg2+0CN+VZBW4T/DdlvP6avE9kE1xdb1EUlZqbzPXH7/f7TjnNzKCe\nXBD6/X5yETO70+Z3cl1nIIadQ3fn2bMzEsiUMjJ7SOZnR6NR5u88n8+TqVK4Kc+q3FugUivtXF1v\n5nLZMnKTiattW2Qe3iWnmXme/X5/qUttVk7zu7guHZ7NZpVw59ldevMc04xFUqfk4mr/zllDo3a7\nrYIvAlfXW5pBpsjcZNLV25Zut+uU0yzrZo1m5TQj/vJdXM1BVXHnZU3HmmUnAUV7jYb8rnIBSKsb\nUt5VETtUTPDg5nqzp+KKzk3m2sJLbvJtc5rZXfper5dUZLipnOaFy/M85xa+Su68NOOQPWa3PRby\nW5q9PhG7lLmUpX0T0UOnUoJ3db2ZwbkycpO5tvCdTscpp1nahSltzC7veZ630w0+JpNJJdx5aTGf\ntPOU90ajEf1+f6lOpA2L0i6uVaGSgt/W9RaGYam5yXq9ntP3MyvRtjnNzPOUyimtmXnu5vd1jZxP\nJpNKuPPs+iLnIL0fs9EwyzeKorWzG2nDJm3hC8LF9RZFUam5yVwFb1ambXKaSZdeBA/L3fjZbJYE\n6czv7ip4EXIV3Hmm4OU8pdWWmQq7rIMg2KhltxuIKlApwbu63qIoWlkMUmRuMvGxb4t5cdkmp5lZ\nGeW5CN289539XV27ouZKtEN255kXIXkuQp9Op5llOx6Pk16AOWY3L/oStZcGoipUUvDbut6iKMpc\n+rluSahrbjLP85y+37qlu+uW9AJL4m+328nCkNu+owtRFFXCnSf1RfbFRjydTteWqST2kJ5TVi9P\nW/iC2cX1Zt4QUv6XvW8/z5ru2iQ3mSvr8pate89c1531/dK+qwthGCb/65DdeVJfTOFvUp6yAtGe\n3ZDeXVqQVwVfELbozXG47Xoyp+Vucz+l/fjm3Lp04819c47cnjcvG9e13C6uN7uHIRearIvRPnG9\nuMnvLFO2Zp4/qXdVEztUUPAuuLregiBwdq9VgbLdh8r+aZTgt3W9BUHg7F6rAmW7D5X90wjBg5vr\nzfd9J/dalSjbfajsl0YI3tX1lra4ZhP3WlUo232o7J9GCX5b15u5XHZb91oVKNt9qOyfRgv+Ntdb\nmkFmE/daVSjbfajsn0YIHtYvkc1yvdlTcZu616pEme5DZf80QvCurjczOLeNe60qlO0+VPZPowS/\nrestCAJn91oVKNt9qOyfTQSfll/u+4E3ABPgM8CbgYsiTjAvXFxvvu87udeqRNnuQ2W/bCL4tPxy\nP0ecV25BfBfbt3PA6aZcXW9yu+Jt3WtVoWz3obJ/NhF8Wn65Dxn7HwH+Rl4nVASurjff953da1Wg\nbPehsn/yGMM/RXz/+oNlF9ebq3utCpTtPlT2z66C/x7icfz70t4sIrfcrqQJ1HZ9pbXYWS478//u\nExfX2y4591zZpDzT3nPFpVzk0fb6H2pOuqJzywnfArweeHXWB/LOLedK1li1LjnGquJ6c3Utupan\na7lUrb7YjWneueUAXge8lTgNVblJ1B1YF5yqQ46xqrjeXF2LrsHQbctEyrLO9WUTwdv55b6XOCrf\n4yZ49wvAPy7iBPMibfqpTjnGquJ6c3EtugrejE1smwOvrvVlE8Gn5Zd7d94nUiR2NLpuOcaq4npz\ndS26imI+d8uBV+f60oiVdnYXrW45xqrienN1Lbpm8nHNgVfn+tI4wdcxx1hVXG+urkXXTD6uOfDq\nXF8aIXiof46xqrjeXFyLri38Ljnw6lpfGiH4tDFu2hisqjnGquJ6c3Uturbwrjnw6lxfGif4OuYY\nq4rrzdW16Jq6yzUHnrngpm71pRGCh/rnGKuK683Ftegq+F1y4NW1vjRC8OaPLM/lh6tDjrGquN5c\nXYuuufpcc+DVub40SvCyP5/XK8dYVVxvrq5F11x9rjnw6lxfGiV484fcxPVWlRxjVXG97eJadME1\nB16d60vlBG+aDqTimJu5ukkq07oKDSxdze33pBsmc8NhGOJ5HlEU4ft+0m32fT+Zzip7Lf0urrfb\n3F1plduVfeTqk267lEkQBIxGo+S3k+NFUZSc2y715dBz0lVK8JuMAc213+aPUKZrahf7qAtRFDm5\n3sp2r/m+X2quvrLrSxWoheA9z2MwGGQuYZTgVFmuKdcgkyuj0cjJ9Va2e83zvFJz9ZVdX6pAZQVv\nBkJkdVPW/GcYhqW6psoW/Hg8dna9lele6/f7pebq20d9OXQqJXhgKcorARGpzOa4yIwIe55XqmvK\ndd7Ylclk4uR6K9u91u/3S8/VV2Z9qQKVErzdRTOdSGkGBamgUunLck25LgV1ZTKZOLneynav9Xq9\nUnP1lV1fqkBlBZ+2TBFIAkv2XG+ZrilXQbgynU6dXG9lu9e63W6pufrKri9VoFKCB1Z+HDNCas9/\nyo8nlb9wXGdxAAAGEUlEQVQs11TZLby4wrZ1vWWNcYtyr3W73dJz9ZVdXw6dSgk+rYKmjcHkPc/z\nGI/HS3ngynBNld3Cy91WtnW9le1ek/FxWbn6yq4vVeC2mpmWZkr4TuKUUy8B7uV/aqvYXVD58aTV\nNbulZtDF9/1SXVNl36fMvm3ypq63st1r9sKoonP1lV1fqsBtgk9LMwXwGPBa4PeLOKl1mD8gLHfL\nJEJuT6eEYViqa6pswZuVcFvXW5nuNXMsXVauvrLry6Fzm+DT0kwB/Cvgu4D/mvcJrcPshslz+eG6\n3W6m06vf75fqmiq7i2cOO7ZxvZXtXjMvLmnlltcCH6Hs+lIFXAab3wi8APxqzudyK/IDyn673Wax\niG/7u87h5Xleqa6psq/265Z+rlsSWrZ7bd3S3SLKsuz6UgW2FXwAfDdxd14o7ZvKD2j+kJs4vaRF\nKss1tQ9cXG9lu9dgfZnlXZ5l15cqsK3gX0Hcxf/E9fOXAx8FvpY4sLfEMwXklnOtBK6uqW3Y1T6a\nB9u42vbhXiubMutLq9VKphbNKT1zkwtEnheJInPLfRJ4qfH8d4E/S0aU3hT8PnF1TVWpYrtQtnut\nKrjWl1artZfytBvTXXLLSZqpR4nTTL2DOHIvVMIi5Oqakhsa1JWy3WtVwbW+tFqtgy/P2wSflmbK\n5I/ndSJF4uqaqrvgy3avVQXX+gIcfHlWaqXdLri4piaTyZ7Pulj24V6rCi71pdVqHXx5NkLwrq4p\nWWxTV8p2r1UF1/rSarUOvjwbJ/htXFOy4KaulO1eqwqu9aXVah18eTZC8ODmmqp7C78P91pVcKkv\nwMGXZyME7+qaqnsLX7Z7rSq41hfg4MuzcYLfxjVVlfuUuVK2e60quNaXVqt18OXZCMGDm2uq7oLf\nh3utKrjUF+Dgy7MRgnd1TVXlXuOulO1eqwqu9QU4+PJslOBlf1PX1L7XxRdN2e61quBaX4CDL88i\nj3Z1SILZxE2W9rm6U6Z7rUq41pdDKM/r/5n6jxsjeEVpCusE34xBmaIogApeURrFXgS/qVlfj6fH\nq/vxyv5uKng9nh5vj8drhOAVRdkPKnhFaRBFTss9T3x7LEVRyuXngSf3fRKKoiiKoiiKotSC1wGf\nAn4beFvBx3oM+B/ArwH/D/inBR8PoA18DPjpEo51BnwA+A3g14GvK/h4bycuy08C7wPcksxl827g\nxev/L9wFPgT8FvBzxN+5yON9P3F5fgL4L8BpwccTvhNYEH/f2tAGPk2cvaYLfBz4qgKP9zLgz1zv\nHwO/WfDxAL4D+DHgpwo+DsB7gKeu9zvkWzltHgd+hxuR/zjwppyP8Q3AV7MsiO8jTlwKcQPxzoKP\n91puZq/eWcLxIG6YfpY4sUutBP8XiL+Y8PT1VhYfBF5d4P9/OfBh4C9RfAt/SizAsrhLfMG8Q3xx\n+WngNQUc53GWBfEpbrIdvez6eZHHM/nrwI+WcLyfAP40JQi+7Hn4LyfOYCO8cP1aGTxOfHX9SIHH\neBfwVuKuWdF8BfBHxJmAfgX498TJPoviHvADwB8AnwfOiS9uRfNS4m4w148vXfPZvHkK+JmCj1Fq\nNuayBb8vv+wx8Vj3LcDDgo7xBuKEmh+jnIy6HeBrgB+6fhxSbG/pFcC3E184v4y4TL+5wOOlcUV5\ndeh7gAlxrKIoJBvz9xqvFVp3yhb854jHK8JjxFe3IukCP0ncNftggcd5FfBG4m7Z+4G/DLy3wOO9\ncL390vXzDxALvyj+HPB/gC8CM+KA1qsKPJ7wInFXHuBLSclSXADfArye4i9oZjbm3+UmG/OXFHzc\n0ugAnyH+kj2KD9q1iEX3rgKPkcYTlBOl/5/AV17vPwM8V+CxXkk80+ETl+t7gG8t4DiPsxq0k9mc\np8k3iJZ2vNcRz0S8JOfjZB3PpHZBO4C/Qhz8+TTxNE+RfD3xePrjxF3tjxH/oEXzBOVE6V9J3MIX\nMYWUxndxMy33HuLeU568nzg+MCGO9byZWAAfpphpOft4TxFPF/8+N/Xlhwo43pib72fyO9RQ8Iqi\nKIqiKIqiKIqiKIqiKIqiKIqiKIqiKIpSef4/hHklgi1HoiIAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x111640450>"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "usps_Y = np.array(usps_Y.T, dtype=np.int64) # convert dataset to our dimension"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 70
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We are ready to run our Gibbs sampler on this dataset. First, we use the hyper-parameters used by [Maaten](http://homepage.tudelft.nl/19j49/Publications_files/TR_1.pdf)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# USPS priors\n",
      "usps_alpha, usps_beta, usps_gamma = 50., .5, .5\n",
      "usps_K = 40\n",
      "usps_N, usps_D = usps_Y.shape\n",
      "usps_niters = 1000\n",
      "print usps_N, usps_D"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "475 256\n"
       ]
      }
     ],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "usps_chain = gibbs_beta_bernoulli(usps_Y, usps_K, usps_alpha, usps_beta, usps_gamma, usps_niters)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's plot the posterior predictive of each cluster separately, to see what each cluster looks like."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def posterior_thetas(C, Y, K, beta, gamma):\n",
      "    N, D = Y.shape\n",
      "    assert C.shape[0] == N\n",
      "    thetas = np.zeros((K, D))\n",
      "    for k in xrange(K):\n",
      "        thetas[k] = (beta + Y[C==k].sum(axis=0)) / (beta + gamma + N)\n",
      "    return thetas\n",
      "\n",
      "usps_chain = usps_chain[::skip]\n",
      "pts = np.array([posterior_thetas(c, usps_Y, usps_K, usps_beta, usps_gamma) for c in usps_chain])\n",
      "pts = pts.mean(axis=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 39
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "digits_per_row = 5\n",
      "assert not (pt.shape[0] % digits_per_row)\n",
      "img = np.vstack([\n",
      "    np.hstack([\n",
      "        p.reshape((16, 16)) for p in pt[i:(i+digits_per_row)]]) \n",
      "    for i in xrange(0, pt.shape[0], digits_per_row)]) \n",
      "plt.rcParams['figure.figsize'] = (10.0, 8.0)\n",
      "plt.imshow(img, cmap=plt.cm.binary)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 42,
       "text": [
        "<matplotlib.image.AxesImage at 0x1140aa510>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAToAAAHeCAYAAAALl8fIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvVmsLd16HTR2s5q9z7n//R0hbJNc2XQRIYIIRCIZpPzn\nwQ8oQPIUR5HyYswTCiQgkK9frOMHUHwlRPcWYYIThWDTCAkpUjCgcwEhQxKSgJREJEgWcSJfg3P/\n++9mdXvvxcPeo/aosb5ZzerPWd+QSmutWlWzZlXNOebXzW8CiUQikUgkEolEIpFIJBKJRCKRSCQS\niUQikUgkEonEXvHPAvjrAP4GgJ88cF0SiURi67gA8DcB/DCAAYC/DOC3HbJCiUQicb7l8n4Xnonu\nVwAsAPxnAH7flq+RSCQSvbBtovvNAP6W/P7Vl30VvvjiiyWA3HLLLbetbi/cEuKy9MeaKF6I+Pa3\nv43lcon379/j/fv3W778ejh0XZbLZbX9zM/8DH76p3+6eOzZ2Vnjfv2/dGxXdHkuy+XrK9/0epvW\nJYI+W60r8Fxf3XZdl21jk3ro84i+6/b09ITHx0c8PT1hOp1iMpnUPqfTKX7u534Ov//3/37M53PM\n53M8PDzg8vISl5eXuLi4wPX1Nd6+fYs3b97g+voaw+Gw2ngcjz0/P8fFxQUuLi4AdG9XTcdtm+j+\nNoBvyO9v4Fmqq+H9+/f48OED3r9/j3fv3uHdu3dbrsbHAzYk3RaLRfWfoqlj+n9nZ2dYLpc7JR/t\nIPu4XmJ74HvSd/b09FRrjyS3xWKBh4cHPDw8YDKZ4Pb2Fnd3d7i/v8dsNsNsNsNkMsF3v/tdPDw8\nYLFY4OnpCZeXlxgMBri8vKyR5sPDA8bjMcbjMR4fHzEcDmvtiG3p6ekpHMCJDx8+4MOHD53ud9tE\n9xcA/MN4dkb8HQB/AMAf9IM4Eh3DqHgMWC6XeHx8xMPDAx4fH7FYLEIJ5Pz8HOfn5zg7O8P5ed3q\nwH26f18kp4008fFA35e+S5IRN0ppi8UCt7e3+Oqrr3Bzc4Pb21vMZjPM53NMJhN8+eWXeHx8xOPj\nI5bLJQaDQUV0JE22b7Z1kivrc3FxEQ7YXl8AK0LSz/zMzxTvddtE9wDgDwP4c3j2wP4cgL8WHXhM\nUtw+6+LkxZfPRvUjP/IjmM/nobpFIosIjfso+pekvz5k1OW5+P3sCqfaXpqwjXookRDaJufzeSW1\nzWYz3N3dVWR3e3uLxWKBxWKBb3zjG5hOpyvEBTy3EW2zlCJJbGy3JET+d35+jqenp5VBfa373LiE\n/ljuq3McGnqfJRsIVQMd5TjSceOx2jBUsuMnbRw+KpYIsisicnRpLjr2GCS8XdjoPkWo2urkNp1O\na0RHsptMJivPVj+1/V1eXmI4HGIwGGA0GuHq6gpXV1e4vr6u1NjRaITRaFSRHtuytvcmvPwfHrRt\niS7xgiZDr5IYRXm1gyjhqcjvUpt+Z+NgA1ES1P/6oMnWp/sjQk+b3ccFvie2GZpT2C7ZRp+engAA\nl5eXK6Tk71rtfOpg0PYL1O3UlOhUsmP9NmlLSXQ7ROTBYgPixkakjSkivYeHh4q0nNz4SXvI5eXl\nymgKrN9Y2jy9vFf9bDovcXzwQW25XK60ycViURHQxcUFRqNRTRJTkwkALBaLyrbnEp5LaBHZ0Rmx\njXaURLdlRNKbf1cC84akBOeSndo86HrXawKoNRQ2Jo7ChKtrmzYk99ppPZuOTxwPIseS2s2GwyGA\nZ4mP6udyuazUzvF4XHuvT09PmM1m1aDsJhgAVfvnoE2i24ZNzpFEt0WUVFRXV0liSnAc+dxGF9mW\nouvS06Uqrqohemxk99gm+bj6GiHJ7vig7U3takp62o7Pzs4qaW40GlVlqOfV4/A46LIfqO2ZYSjE\nNu2nSXRbRmSH801JTl338/l8RQKMRkK/lh4L1L2w5+fnoVdN1dtdhIc0EVyGoxwftB0BqIiOAyNJ\nTwfe8/NzDIfDiug05k5DpEiKj4+PAFBr1w8PD5WUqB5bIInuo0BEQupYUNucEl/kGVTbhzsjtCFo\nI6QKEElsg8GgIsBtOAxKUmdKcx8HIs/02dlZZduNQj/YFnWGA2dDzOdzAKjauQcMA3XpUePstE1u\n0xueRLcDROTjUp2PXCQwLyc6RkNIIrd7dL2Hh4daWWy8PL503SaUQky6ntf3euuU30X1P3VEYTbu\n7WSbAV49sxqz6eeUBmE353hbLUl0m7aRJLotw196Sarz0VNVBR81leQ0vojnRtdUVZjX5fUuLi4w\nGAyqYyNbWlvDisJnumDd6/WFk3CSXTdEpML2x3enZhGSnZ8flaP2ukjj8baUEt2RovRSopFLPaE6\nWvK3S3AaRnJ5eVk1rmg0XS6XmEwmeHx8rNn9+B/tIcS6ISGu7vQ5b53r9UVKc/3gUpzuZ9uJSC46\n1snOB3vu88gEHfzTRvcRQEe+i4uLqpFo8CRxfn5eeUuV5CjlKclpHJ1fj6CzgZsagT1sxc+Nfu8K\nEQlt69pNJLcLT/PHitI7UCnO/4tm5kRbpJZG5fnAvm2SA5LodgZVRyNnAY9hg2Ackdo+VJqjFOdR\n5YSqaGqDox1OVWaNzQMQlrdP7FqNjcouSS+nBG2PJUJrIromkuMz5eDqE/gBVG1O7c46kG/zHSXR\n7Qg+2rnX6ezsrJLuVHx3h4NKc03ZS1Q9BlCR3GAwqJWvCQQ424L1OUSHV5VabXe7JrtTl+oikovU\nxlKWHCc1t+GpRKexoaVr+BTGKKpgEyTR7Qj6gqiyUqWkmuq2CUp2rqbyd5NIr0GcSpi83sPDQ81G\notPK1OZ3CGzTZtfFLnfKUpyj9LxKJAcgJDo/hmU3EZ2bZ9w0k6rrRwidc0qbnYvyauz1ift9bBY+\nSjIoU0dZxu65TVClq03Q16kRfV8H6XzYHCWJrklFjaILdEDV38vlssosPBwOq/myzF8X2aA3RRLd\nnqBeUldjiUgdiBpVX7JzmweJbj6fV+WSfLdJcm3EHMW76fNJHBYuqZVIDlj1qvqmTjAAFdGR5JjC\nSVXXbSKJbg9QiQnAijTXpgJEv7tck+SlEqFKdEqGj4+PtQj2TQmvqwTqhJYEdxxoIzlty2rjUy+r\nkxyJjmWT6JTsaEbZtnkhiW5P8JemDaWPWuroEkLhZXtD3GSOYXRM1/uIjtsV0Z2686GEtufRJMWp\n9L1cvqZ10mSdXChHA9bVrKLp1ndhmyOS6A4AtWts8lJ9dkIfla8peNnV5CYyK0l+fe7Lj9s22fVV\n+08BTe/PzQf+7FQ7UDPMw8NDtVDO3d1djezoKHOS82SxuyA5IInuYNhUNSXWITn1hnUJ5Cxh3f/8\nOA8xYR23gYjkkuye0fb+SqE+kW1O15iYTCa4v7/HZDLBbDarVrXzuFAlu12mDQOS6A6KTV6mE1wp\nJio6Ljqvzdi8i3uIytil6poEt4p1ByuX5kh0VFu55itVVw0jYbYT9bJG08m2iSS6PaFvB26TbKL9\nbHiaDirKh+eJAjwi3cvfJ0F0VWP7SIxJcNuDE5ymGtPFdEhyTCbLRBLj8Rhv3rzB1dVVFVJSWnNi\nm0ii2yNc2oqg9g8nu6bzIo9X5OrXCf4ap+cZUVTS20bISReUrtXF4dKl7MTm8DASzaWoJDeZTGor\n252dnWEwGOD6+hpv376tVv9S1TWJ7hNCV1ta1NnbzmsiO903GAyq8BOfiRHZ+w4l0XV9Tkli+4WS\nHImOmbKV7Dx27vLyEldXVzWJzkNKdoUkuh0hsoW54T+aEaBqpWYB9uMdTUGaOqna6xbVNboXYl+k\n0sUpkSS3P/A9uDRHL+tkMqnWep3NZlWiV9rluI4r4+Z2GTMXIYluR3DpKkpXo/FrHmPUJTo8CtT0\nZRRVfYiIrRRmouWXPKOJ00CpnS0WC0ynU9ze3uL29hZ3d3eYTCZVOIl6VSnF6XQvj5tLie4jhMYW\nlZwDLuUtl8sV9zsQ26KctFxlVZWCBmGXICPC82tFYSeJ04O2FR1Mp9Mp7u7u8NVXX+H+/r4WN8c1\nJaiu0i5Hac4TVewSSXQ7AhsECadkM3PpThcKIfGUGkKkFutoqwvwONHp+UqSUdxZSnGnCx8QtU3P\nZjPc399XRKezIIDXKV4kOUp0w+GwmKF4V0ii2yK0QTB4cjqdVhKVk1ubbcxVyaZrukRXUlvV0+oT\nqFVKVILdV2NMHA8ijWG5fF6LhAHBVFkpydG76nY5kpyGk2i7irSHDBg+Ymij0ODJ2WzWmE6a5wL1\nBJpRTFt0TmSj0/Q4nkBAPa1KYlo3J7gkutNBk1mERHdzc4OvvvoKNzc3uL+/x3Q6rUJIBoMBRqNR\nRXIeN6fTvUrXBrbb5pLotgRvFCrRTSaT2v9N6qhmHPb5hH28rp5xmFJaFDvHevDYaPpPEt1poWTD\npQPi5uYG3/ve93B/f18R3WAwqGY8qDSnEh21CAArA6y3uW0iiW7LUDuZhniU1MKSTYxEF6m7ei3C\n7X8lciypzVEGk+j7vtEUV3eIsJdPGSWvPM0gi8UCd3d3lcp6e3tb2eWA1+Syw+EQo9GolkwzipVz\njYTXbOsf6yCJbouIwjWisA1H9GIjtaHJtqcSnEpxferoK5EdQ5xa0yyJ9AhvDxHJKdFx7urNzQ1u\nbm4q2xzbDqd4UW2NMgaXPPneDt1+t43Fm5LotoiInEphG0T0MiPbm4eklIjPzy/VM7LpRSsvHVqa\n4wivnl8d+fWT5yTWQ0nSZ0aSu7u7iuQYN+er1XHCvs56cJtcJJmrNkJNhvZkYPPg8CS6LSKyabgj\nQOGkojMhVJwvSWCRitqnMTRJiH3L2hVSdd0forbHRdDpgCDRMThYQ0UowZHsnOTUZqzQ2E864Vxl\nbbJPd0ES3RahhBVlBNEX5wvg+IhVkg5LpOYSY8m+cX5+vpLwMFojtosDZNdwNbWpLn3qmaTYDBLc\n4+NjFSt3c3OD7373u/jqq6+qpJoPDw8YDoe4uLioSXGqrrIsAFUWE5fsouB5nwK5qbMiiW6LUKlM\nF/l1ldTJkPt85IokxJJdKpIO+en7PQ+YqxXHQHJEF3Vc1demOh+Lg+XYobGYJLqvvvoKX375ZRVO\nwvmsy+WyIjpd4IaxciQ6viOuVdJkavEsOpTyNpHgk+i2CLe36aYvMyI6nl+yz0UzGyKS86QALjG6\nROeZXXntYyO7tnp0rWcSXDtITowFpUT35Zdf4u7urso7x2Od6DRWbrlcriyUzvM4cOtAGw3M2hbX\nfX9JdFuCvhgA1YTm4XC4on6WgnQBhDMotPySVOLERtXBc/HTwButiK6IRty+2JQkt022UTlJfM/Q\nZ61p0bmVZtrwXODVe8opYpwp4aYbIhqggVeJrjRzIlXXA8NtdHS1A/WAXuA1MJiNRlVUJTsiCvvQ\nRqDEpkRGFUJRCiFpc070wTYIalfSZEkyPlVE7S/KfhM51ngO40Xn83mtzZVMKzrIkghV81GtZBvv\nJ4luy+DLoqtdbRMaPBzZ3zyExKU5V0sjKS1aeIRwKc3L9+PWJattSINa1jah9jw1FSTZrQa6M/tN\nKSkEn5mT42w2W2ln2lbVfDIYDGplqvRXctStgyS6LUJHLRIPX/JisahenEp3TnhOLlHgpAZg8tPJ\nTb+r8TciVa1/k3rcF8di43NsYtT+VFGS6iKSU3VUj18sFjVvqw7a3obpsaWNTx1j2i59YNJ9fZBE\ntyOQ7FwUv7i4WJmX6iTn0pafr2Er0dxV34jl8tkbHM3YcNUCQI1Q2+A2m30T3DqEtQ1J4WOGS12R\nJqGZbmhvdpsZY+20rXjuRZbnZdLEQ2cGyY9wW7ISaZ93l0S3I+gL9d8cJZXogFjl05fZRGYu4ZVC\nW5TsIjLVz75qg9d9X2QXkfQ6558SIpKLzBUcsOlN1YFZg4qpsRA+mLMsfmpbZaaT8Xi8Ij1yelk0\nw6IPkuh2BCcdSnJ8cZpGqSsiR0MUC8frl6Ajd+SUiO6lKw6lrrbdQ2IVkckkUlFVotPkEQQlOi2T\nCWTZzkvXBoDxeFzNitBrU41VO56qzX2QRLdDRB1PpSXa66L/I7hndZOFf1383yZJqJF/l3AJNLE+\noogBEiAzZAPAaDSqLWeoNmG1C2tbdYcDwTbCJAAah+eDuceE9kUS3Y6gdgwfgThSqWrr50Uv01/2\nJh4pl/62QRZexi7JrhS2kOgPJTkAFTFp+2RQMPPPjcdjzOfzFc3CnWXRptfk9ZTsRqNRtYiOTlXc\nZFnEJLodwsmEL4nkF6l5bapk1GDWqZdLc5tCy3QP2a6w6XNIvELtuvxOJwHVVl2ykDMkogSuLE8l\nPSc+vZabY5T0XDpc10SRRLdjqKraJTati81sm/XaVZn7ILiU4rYHHzQYMcDpXCQ6VTFns1nNG6o2\nv7Ozs5o31Um0aYqiTieLAovXQRLdntHWMT+VzvuxkXQb2oi7VKdtEv6277vJ+aQhIbrsJqX1i4sL\nzOfzFU0FeL1nl+ja5mHzP8+ss40+kUS3R/RRsz5mstMO8THfhyMKm3F1ret5fcDnuE6gbFOZpXek\nIUg6U4HEQ1U28qbyfODVxleS4FySc7V2k3ASRxLdnvApdfgu+FTvt2RXXee8vtiXREeoZ56qqYaZ\nMGKgNOXPHQ9R/sXSvOtt216T6BKJAKWAWv0PiL2+bqvcRWzhrgcStbfxt+7nb01CEd2j2qiZrt9n\n7TQR3baQRJdIFBCly+piq9NtF7NEtq3GOpTkPMmErzPC7Ca60p3Wk5+qturaEpqxZJf3l0SXSATQ\njl6SVvRYnx/Kz115n/ch0WnCV19MSdMyMV+drwcRhZMwhISe2SiN2C7uL4kukQigUk3TND1XbUlu\nfUlu0469TWLwe9JF0TWjCVMyzWYzzOfzXkRHkmNwsh6zi3tKokskAmgHBVZtUyWV1GPJ+l5TP7VM\nt1+1/e6CrkTM+9Gcip7Cqel5RPGVzHbCuaw67WvdDCVNSKJLJApQ6UKN8m6Yd2j+tT5wu57u907v\n9i8lZf2/hKZ7UDWc5XiOOs887GV6ObwOy5vNZhV56gpiw+FwJ2FJSXSJRAAnEgA11azkTdVOvg7R\nRdH/ShgOHr+OBNR2D0q8TnQq0TU9FyVp/aRKPJ/PMR6Pa2EspeewCZLoEokGUL1y9dCxDQlEyTVS\nBaOwFiUWDfDtWq8S0ZU2r4cGFEfTHKOyVQJkfkbd3Evb5T7akESXSARQInGjemR016wa68bNOZE2\nxaa5dKRxbV3Jrkmia4r/8/mtJYJTAi4Rp3pzqRLzeWqcXdN9dEESXSJRgEtYkVTnXkVg/SDhyGhf\nKiMKSiZ0Sc2S2ttEZk2SXOSkieqpJOZrE1NCVrIj0SmJ6op5h5LovgHgTwL4ewEsAfxxAP8BgN8E\n4BcA/BCAXwHwYwC+3KiGicQBEZGcSzA6p3NdkoscEE3/+3Gu1jbZ9dquXSI6J7kmSYs2OCUxV1t1\n4zG6pgoJ2x0Z62BdolsA+NcA/GUAbwH8RQC/BODHXz6/BeAnAXzzZUskPlq0STGquqqk0hdarjoZ\nSqnOI7tc11kFStzRsdECSu748PL0c7lcVs4KnfwfpfCP1N1tY12i+7WXDQBuAfw1AL8ZwO8F8MXL\n/p8H8AFJdIlPAOzATnb8z9XZdYnOpSeSgs/QiFRV7u8j+UQeTqqMek23Q/Jcraur+mpzU8ktijFU\nG+MuyG4bNrofBvBPAPhfAXw/gO+87P/Oy+9E4qNHKcZN/1dsQnSuqmksW6lsV/O6SnVtdXGi43lN\nqdL5nEhyrJ/aDqPr70qaAzYnurcA/ksAfwTAjf23fNkSiU8GXaWlvvYkt6lFjoltB9GWyoxyx3mI\nTRei44p3Uar1Uh2anCybYBOiG+CZ5P4UgP/6Zd93APwAntXaHwTw69GJ79+/r76/e/cO796926Aa\nicTHDw9n8f+2TXJt19PwDpUWAayQnEtpJEYnS8KlVSdOHtNms/vw4QM+fPjQ7X47PpfovJ8H8Bt4\ndkoQ33rZ97N4ts19jlUb3XJX4mki8bEj6hsaYxYtj+khLqU05V2vp1lK1GbWRarTcpnZZD6fYzKZ\nVCuIcfUweqo515WbTgnzFcCa8HL98EbXlej+GQB/CMD/AeAvvez7KQB/DMAvAvgJvIaXJBKJjlhX\ncnP1bxMV2+1xGkaj/zd5bktxh/zPnS5tXuBNsS7R/c8ASvT6o2uWmUgkXuCzB9zLS/j6C20E1AXR\n+eqIaSo7IrVSeIyeoymcNl2cPULOjEgkjhBKclHyTyWgpkVn9Ni+cEJrcpj4fg2DaUte6kSnK4BF\nNr51kESXSBwZnOSaJKGI6LahAvJctYt18f5G5Nh2H7yOSnLbXOoQSKJLJI4CPjvA54m6XatJVd0W\nOZQktiYw1k8dKPP5vJolofXn+hGRuurZSzZFEl0icSTwxWhKJFfysm6b4PrOstA6a5p1plrnpP6z\ns7OamuqZhSPVO1XXROITgJKbZvlQklNPJlW9bdrkFOuWQZKjJEeie3h4qO5HM5QwhTrDR/qExvRB\nEl0isWdEdioSG+PXdO4nEIdhbNMmF6FreVrXaHUwXziHEp3GzrkTYtv3lESXSOwJTdH+Jc+kxq5x\n2xUZrAMGFVOSI8HNZrNqXQlXu6m2DgYDjEajlcVxtqmGE0l0icSe4La30n8R2bkDYhfqal+oFEop\nTm1ynMmhUqmGkSjRaQzdLsg7iS6R2CNIDBHRkTgArJBYySZ3SGkOeFVXXZpzotM609vKVb9IdNtc\nI8KRRJdI7ACeIBNAbQ1Ul9o0tZJPnyoFA++T5KLkmFzFSzclOTogWG+qqKPRqJLkPDi4yz2tM1c+\niS6R2DLcU+oT40vBvzy+RHC7sF11vR+9FzpMVF1V54OSHM+9uLjAcDjEeDyuNpXk+jo++pJdEl0i\nsUUoyXlMXFt+tX1M61oXKsVRVfVYORIcN70vqqpXV1e4urqqspSQ6Prc0zp56pLoEoktwyfkt831\nJCK73K5DSLTOpX0qiWroiDsfOPvB7XJqkxuNRhiPxyse1q51bMtRV0ISXSKxIbwDRh1SO3O0r0uc\n3C7q3aX+GkKiJEepzlVVhpDo7AeqqgwMVlVdn0Gferadp0iiSyQ2hNvfmjqjOhx0f5OzYRfqajS3\nVgnNZ2pQWlPvKr2qmqBTSZrqKW1zTKQZEZ06YZrqGQ0gXaaqJdElEhvAVdS2MBGeo4imdenx+6g7\nCUuJi/vV9qYSnQYDk2zUtkgP63g8rmUMXkeiKxFd03mKJLpEoiOikBB++qZSWyRxuPS2b6+q2w+V\n5HQ9Vic6VVU9rXtpDqtO2tf78xCc0n1Hc349BCclukRiC/CwhqjjAeXMH5EEcsiwEZXkaF9TwlNC\n803n4irU8eAkFy1+w88mx0KTA8JtnE1IokskOqKkPjV1QD8mCgRuWut0VyDRkbxUutOZDk52aot0\n+DxWT78U2d+anqMe40kB3M6ZRJdIbIguXlWHdkSWoWTmhLdv0JvqEt1yuaypqB420iR9KXkrwfl1\nI6dN0/NUYi2pq0l0icSaaAu/aDtP4fY4/exTly7oYtSPbHOqtqqk10bouvFYEijrwuuRAH0QiKQ2\nr7M6PFyqa0MSXSIRoAvJlVSt0n+lIOC+dSmhTcKJvKxqp3M1NVrbNbqmkzbLVLK6uLioEZ3WMSIz\nJ071Yi+Xr8k7254JkUSXSBTQV13184D22Q7r1KUJbdKcr0Wh6itVVQ8xaZPm3M5IiY7fSXBKdHqu\nq7P+fPQ/Pef8/HwlC0wJSXSJRAO6klwkyel37bT6uytKUpWSqEtDTQ6OrgSu50a2RL8WSU6J9eHh\nARcXF3h4eKjFCipp6fOLiJPPTW2AbevFKpLoEokC+trkSiomvzNlkTomNq1HRGxOeF3KBbCiDmo9\nS/Ywvy7vkZ/8L0pQULJR+gpgWkd6cnlc13nESXSJRICuNjk/3o8jWUTE0bcukZqmqqNeO7KdtdWb\n4SGRzax0/1621rFEygBCO2UUduPna31dqmtCEl0i0YIuBNdFrW0ra51rKBnod6qATlYlKMlQYnIS\nbKtjVNcmgi5tpUScKgk2zSuOkESXSAQoqYT+2VXiW7cO0ffomEhi7OrZ1Vg+l7QuLy9rxzmcxKMk\no022NK13qc6RxOcLXbfdYxJdItEA73Rqb4s67zbIzju5OxpKUFtbV4LzjdchmbgqHF1Tid6JzvPx\nlSQxvT9Kol5HX/dVZ1y0BV0n0SUSBZQM8S6xbFuS02uXfjtcxewrzak0qCEfbWVFaiq/a5weyY6O\nimig4DUYhuL1VLJz9TYlukRiA6hNyPcBdU/qOoRXIrSuRNV0XiQZqh3s8vKyGIdWmnYVwVV5nuvp\nnjx+TwlPSVaTdroEx1XDNIlnlyl0SXSJRAC1fTnZqQTUJ8ShdA0nqei/tnKUKLTzl7yaqpZGRNGF\n6LyObpcrbS7psQ5aP62nbpr+SYOP25BEl0gU4B2IpKaqrNuY1iG7SArrI8mVttI9aB1L9q2uRKfE\n5AQXEZ+THKeKRSEly+UyTPnkGVE+OaLr4u3yl6jfm0bOxP4RkYWii51q1++vRBQuxaxbjyaSKtUh\nKsO9pU3PyvuF28OiOLpSuT5TwYmu5JhQotMylDCXy+WKFKfH9sHREl1EXsvlciXLaenBAljR76n3\n62dX0TexHbiHTudbup3LJYbS730OWkpwTfFim9alj42urz1Pz3EJD1gVKJrK91x6LuFx8zmtru6r\nY8HPUwmuT8YXxVESnT9o/b5YLDCZTDCZTDCdTmvZFvzh6SpEXKRDt9Fo1HtkSGwONVbrSlLuxXTv\nWuk7UJ/ruSsoyRFRqMQ26qKk0bVefa/XlP2jTbr2epbKbNLCIo+rmwVUrV2H0ImjJDqgPOl4Pp/j\n/v4eNzc3uLu7q7KgMusCJQQAlYeGqxBx8VyOvBSLU6LbH1yVWSwWmE6nmM1mK1Kd2mNcEqda42rY\nruDXicIpSmaTda6lHbxP3fpcg5/reIuja3ctc12njX/vg4MTXaSiRisS8fv9/T2++uqrkOhcjVVX\nNBfZ5YL7nnE5AAAgAElEQVS7V1dXlSThRk4fIZMIN4PHVqkUriu9u0Sn5Pb09ITLy8sauQCxfYr7\nt/Hemtqn52zzxWKaOqir3pHq5ufuitR30b6byjxEfzoo0bkoywbsOet1u7+/x+3tLW5vb3F/f1+z\n2XlHIYldXl5WJDcajTCZTHB1dVV9+pJs0WiaZLce3Ban64Lq+6XpgefouSXJSW1MSmyb2HK87qX2\nqYkqozrynCYnmBMd1XGdWK8pjVy9yzbZHQcjOm0QbqCkOsNNR/3JZIK7uzvc399jMpmsLNyhI6FK\nabPZrJLwRqMRrq6uKnX2+voab9++rRqP2u18OkqiP5iTbLFYYDab1d6rOyMIN0r7poNRNCBtSgZN\n7VPXN10sFo3EHJFUtGliShKcGuRz4N0MB1ddPYr68fGxkr7odKDjQT8nkwlms1lNfSDYCNR9TWni\n8vKypsbSoQG8Gr+pInFCcxc7SWIVasdSKZ3PfTabhfYtoD4dyD/1vUbzMSOj9aaEpyERrnEotH5K\nyG1EF0mDvrYCgJrzJdEdR6G6uk1O1dHFYrGSu/78/LxGRjp6liLDo4bFRqq2uaenJ4zHY4zH45WO\nko2sG5S81CbHd6nvEFh9nnyXumyexllpvNVisVjJYqEqsEbZb/LeVJpT6ZNlOtm6Sq3HleyLOlir\nyhoNBIl+OKjq6gGE0WrgulCHThdxbylHQI2g1g7njgo2Kn6yQfE67Ghq99mWV+1Th9vl1J6lz5sE\npoOIPmN6zPmu3SjPa5A01dvH96gEuck70wHZ1VI9hvUrqZnuSCEZO0jYUZxboj8OLtGpjcZJTsnO\niY5kpiOmxstppHZUtgcZA6/eNJY1Ho/3Gqf1KcEHr0ii87mMQD1IlURHB5FLirwG24FOtHdj/qbm\nB1XBtS0qufJ3tI/10vqp5Kcqq/6XZLcdHFSi46dLXNyIyG7G81VtUQmAKi0bpXr72OFUDaEay2sw\noFivuSv3/qeGkhSt74zPzxc61g6tcxyV6NTgr5I4UJfQea1tqH4kHtbF03irpKb34U4NXR2L5bpa\nG/UHtVPy2GyD3XFwic5frNs/dA1I/e4ZDkh0JDttiE503OiUcAlBZ08sl0uMRqOa/S8bWDv0nZZU\nOiBOEum2tSi+zMMwlER01sS27FtKchx0Wddo02eg9WO9dCCPbHlA3bFBG55LgdkWu+EoiM7d96qi\nemyRTwHSvFUkOpXoSkQ3mUwqtWY+n1eq1dnZWW2KmLr/2cBTjWiGD2CR9BLZ3PS7kly0hoGGYVAy\np0mC738wGGzlXenAGk0ZLJGdD6CqiuogQEnNyT4iOVWZE91x8PASwkdDNnLNWe8OB58WpF46Jbrl\nclnFPjGWS+05JEJ2mMlkUgUPs6NRutO6ZGMrw98nUA/TaQr/cO+5qrV87mwL/K32M5ahEt+miOru\nXlEldm13KslpuJMO6iXvrDpcVI1N9MNBiU4lNRqVdQRnA9aG6kSnJKdTubyDKHFGS7oxzo6SwXQ6\nrcqhYyKluTr8WThZqWmh6ZxIuvOQEScUlehIHloOiU9NGH1RIpSSB9fVStavqfzo3pvCUTb1Hp8q\nDkZ02iiUnHRxWickAKHkxt/aMdRTpdegauskN5/Pq2s/PDxgMplguVxWJKcTzpPsmnPIuUSmcIlH\n1ULd1DxRkp74jvWd8z2pN38dolMvexPhlb53aScRcfm1SkSXZNcPByc6wr1vGhXuI6Ta4fQzsvkQ\ny+USg8EAAGrTxXTKGW12/P/h4QGDwQBv3ryp9iVe4Z1QEUlztDVFKpuHCinRqQPIBy+V1l2iA7A2\n0UX31OU4r6fu9/Ii2xxRIsokufVwcBsdR2W3qbnE4BKdq6xt3lD9j8Q4Go0qaW4ymWAwGFQdk9KA\nhqOovSTxij6Si6uiKsn7cdE7dUnLj1Hj/tPTU61dHRoR+TXZep3s2qS/RBkHleiAV/sZDf4kPm+c\narvRzrFuuIc6GGazWZW9ZDqd1gJd6ZFVyeAYOs0xwiVwIM7SAaD2GamoKuWVUDLg6/vR4N59oeSB\ndQdDH/izSZLrh4NKdKpqeu4xJTsduaMO0ffF83okWaZvoqd1Op1WKq16Y9dpoKcC9YZGg0H0riJb\nnB7fFLcYmSkiUolStO8Sro24NzaavN8VSXLr4+CqK1+a2nI8ot7jqbYBxt1dXFxUJKcTyGmj8/m2\n6sXLBldHk4MCQDhIKdE5VFIrkZ13fCcUD1reF5pIrmt9InW/yQ6dKOPgRBdBX7B61XZRPlCOzmej\n1Hmy7g1meYlXlKRtta1FUjrhJBD99k33ezCuS1Hbfl9aP5ceef2S+aN036q667zuxHo4SqID6nMU\ndyGuq5ThMVsEVR/NjEuvr3rMUrqLoYMJf7eZHVQF5m/91OP8dzRf2leE17psE6X6ahtqqkvkdNNg\n+JLUm+iGoyQ6Fc93RSJKUpFqANQzVpDoOOE/JbpmuL2tieBU4os+u6Bk8I/i9naFyBaobUilTKAc\nM8dnw+iAbaSZOnUcJdEpdv1yI2LzjqJe2EPZfI4Vbe/HCS3qsCXvZJM9qs3o79d2EtpGu4pUZt1f\nIltuLqG5k6bJI53oh6Mnun2hSfXxLUmu2VGg3leHklCJrJTomtS1aEDyQGQ9nzGQ2yI5rXtEZtGg\n6P+zfk2mlCZVP9ENSXQop4tyslOJLtEsFTXZ2UodP5J4IiKNjlXCizyVkfNiXbSRnJNvicwB1Miu\nD8kl2fXDpkR3AeAvAPhVAP8CgN8E4BcA/BCAXwHwYwC+3PAaFXYlSUXSXEnCy6DhOtrIDsCKVAO0\nB9X6MU3qq78XleZKEtw27HUlktPUTF1NHVpndT6kNLcdbEp0fwTAXwXwtZff3wTwSwC+BeAnX35/\nc8NrhFPDgFhC8BGv1Di0LF+ImJ0gyoYShQNk42tHRGK636Xkts7tYRtatmeq2QVBlAg6Mnm4ZMf7\nU+lNU8pHba5r/SNCzfa5GdH9FgC/B8C/BeBff9n3ewF88fL95wF8wIZEVxq1m0Z9HR0jO483zNK6\nFJF731WhbETdUVLdXBLrIr0sl6thGzrwOGHs4j21ER1NHT4NLZLeIpJzkl6H7HYZufAxYROi+3cB\n/JsAPpN93w/gOy/fv/Pye2O4StCk7rjaEtl5XC31JRVZTml0ZRltqlUiDvsovT+3rQHlzq2DVER0\nKtFtm+ja1O7IS69p0N2r6gNqpLauU0ci2+b6RPfPA/h1AH8JwLvCMcuXrRElVVQ3l7i0szRJBKWp\nY9oYHx8fcX9/j+l0isViURGdp4PSxqdEeeqNKLK9RcdEpgciCqtosq251BR5W9ukw3VVwZIUpyTn\nHno1h6gk58li+5BcyeandUqb3ivWJbp/Gs9q6u8BMMazVPen8CzF/QCAXwPwg3gmwxW8f/8ewPNL\n+eKLL/Du3buVBqMSV7S+a3S8k5x2HMVyuawtqXhzc4PJZILFYrFin9NVxXQaTjolVr2P+qnQd+nn\nlUguCitxqSkiushzuW4cWun+IlNKKSQpsstFWbIjb2tTvUrPWvsFf3+qZPfhwwd8+PCh07HbeAJf\nAPg38Ox1/RaA3wDws3i2zX2OVRvdsmSn8dCOaF3QyJvlHUcbSyQZLJfL2kpg9/f3uLu7w93dHabT\nae28t2/f4vPPP8fnn3+Or33ta1UCAE0CQBI8JTgJlCQ2l3y4j4iIjgNKVD6Aal2PyWRSSeJsI8vl\nsraKm7+vkt22z/0pmbhHnnOimd5L75EzHXwRJ82nqG21SZorDbKu4TQNHJ8aXp5X+NC2FUfHJ/7H\nAPwigJ/Aa3hJ+SQbDXUk1Mar0lcp1EOJLpLsFE9PT5jNZtU2nU5rqqumjBoOh7URN7LRnTJKdipF\nm/TbpGpGBKPtRL2uWlYpU0pfdc41hyayc7VVTRuR2qqSXV8ps21Q+ZQluXWwDaL79ssGAH8XwI92\nOckbrUpuuiyh5oSL0pn7i/WG5Z0GeO54Xr5exwM4o5F2E0Pxp4K2e9d3VSI6JSh9l+fn5yuqIclD\nV3RTkvNsH5HHvA9KJB4RnJpWdCUyJdiSLa5rXbROLrklmnGwmRHaeFXc5/oN3GazWY0EAaw0kpIq\npJtLjxGxkuh0gZ6S2J8k154GXNVP7itJdSQ3SjwuRan05glReS7JpGT/Wtd7WSI7bUteHz4jtxk6\nAXepU0RyavM89bbYBQcnOjYSVSUnkwnu7+9xf3+PyWRSGzHPz89rzoEmONGVDMYq2UWG7ZKkyN+n\nii4SXZtaWzrGJTo3Y2iKe9ajtBTmuipr0xZJc/P5vEb+KtFFMXN96lW6Pq9zCja4TXAQolM7HCW5\n6XRaERuNzNPpFPP5PBy9St5U7UzqfXInh9pRVELQRhQZm90WxO+nCJU0CLet8bs+v+h56Xk8ruSU\ncg98W867dVTEaL/Wi9fXuuj9lkiuLU7QpbeSuszvpXt2U4C221PEwYhOR2h6Pm9vb3F7e7siYVEl\nof2FHitKdEpmpU3tJ/rdN5X+1LtbIrpTbTglw7erqy5F85iShMyyXSVUYvFpVRE5rfN+SiTjhOPS\npQ7E0f1FRNf2TJvacom82E9UE/F3dKo2vYMSHaU5Ljd4e3uL733ve7UkhUA926q65n2Ba+1MOvpq\nZ3GJQDuLqkyRiptEV4dKGlGQa5OqFUldbpPT9qHTqfydObm6JLeJyur3qtoI7ccunen11C7XFgit\nROsZc3TQ5TEOtl16cVOie8VBiM69rYxno21Op8voBG0G72oMkttxgLohXNd7cCN25CGLHBce31dS\nb04JkbQTqWH+HAHUVE0/zw38bBs62Chh+rW2eX9+n9oO3EvflDxASZflaFtzUnKp0TULLdevA2Bl\n0EiiOxDRuceJL0EbDEVwEpwH6VKi0wbA8vhdG6WHqkRiv39qY4zCYE5ZqnNJzKU4JYUmAtJ35zY5\nVQ1dwvG5zLx2pNr1JUBtRyUzh5dZIlrWmfdXstH5Peg9u9oaSc+n2Ab74CBEd3l5icViURvJ1Jah\nE7IHg0G1uPR4PF6JJtfOAdQTGVKi85g5SoylxqIjI8uJGj2J9tRmRRAqORAl+1KJbPR/DwR2E4NO\nlOczL6nJkfrZBS5dcqAsmTFKkhPLctsk66rPK/qM7HEuELB9nrLtrSsOJtFFcXAlp8N4PK42n2Tv\n0fE6gvJ3RHRNmS0iW5OP6Iy3o53w1ODOB5egSsbzyOngXky3TTnRKXFECRuayKfLPRGqQqpkFZkw\nStdhGZEa31bHSGpk3yCU7BJlHITo1GbRNJFZGwcbOoCawVVHfp3ONZvNKiO2djIlUpUMtV7cp43R\n67SuxPApwqUq7mvzruo5SgKeCNWJQc9r8khGx2t92+ASopOclu0OEe4v3aOXHT3T6Jm57Vr7Tylh\n5ymbWIiDrhnhHil3BFDFZCCmSmIkLLW7+fxVGrHVrsHPUuS82gnVxuQu+6jzJV6hHdKfT6nDacdv\nIjkl1rbQCx6ng906al50Lf8/KldJ0u+xBH12astWu2g044KDtGdDWcf7/KnhYESnpOKL87pdYz6f\nr6iLfHmqklKK46c6Hlg+G4Gqxt4o2Gh1loSSJJCT+pugEoRO53LJpKSulbzdPK90bIlE/F2t0+FL\nEhn/0+P8Wfj9Nh1PqJ265HxQKS7atM2eutPiKIjOp+too1LngY9gAGoeVf1UldU9uZpUk95bnWmh\ndiKv76k2lL5wj6x7RyOicHXOBxKWo98jacXJcB2UyFKJNFLDo3JKpF6SspTANLVUdHyUvCAix1Nv\ntwclOr5IelSvr69rROWBw9pp2OAj1z9Q90y55KjS3HA4XOkYGmDMNE38j3V3Z8Um0sKnCLeTlexn\nqqI6GTSpvqU09yrNdFWZtc7Rd56rqqPeo393RHZKHTQjSU3bbCSZ8ViV/KJjsz0+42BEp5LV4+Mj\nxuNxJb1xzivV0Eh9cXJy6YDZgbl5DjCS3GAwAICamjQYDCqy08zCJUKLJI9EbMz3zybbWsmoroOV\neiD13SsZAd0zfLSplVE5TeEsXnclZrcbR4QY3bse44N4pKZmezwCiW65XGI0GtVyiz09PVWT+TVN\nk6Y6J7zBqAQ3HA5xdXWF6+vrWlgKg5C5eSQ6CW6xWFSdSiWEaLRcx8D9qSKS3NzGVXI0EG7aiN4t\nBytKhjrJ3yW6PtJNE2GpOq5E4kHrfl5kT3OtozSANxGfEyevmajj4Da6y8vLiuyULPj/YDCoeVPd\n21VaXGQwGFQkp0SnG9VXJTqdMsZ4PJ+lES2YkyjnAixJV9xIgBoETNtqtLau2ljZfnTGgl/brx/V\nOdrP69OZpfeiUBL3MrweTnAlZ0PkQClJeEqUiRgHJzrgmaxGo1GtcQ2HQ1xfX9fy0nEeLPAqQfmU\nMCWlq6urWqCxr/HALQpMJdEBqNl/RqNRkexOvaHx/tWh5MTj9jq10elgAmBl8HKC4MZj3VHQhei0\n3vpdSffi4qKS8tXJ1VQmz3Vtw8nJJTEluuh+vNxsd91wiCe0jIy4GiSqqdRnsxlubm5we3uLm5ub\nWuaKs7OzFTucEp0vjKJEF3U+DVj1HGMAqutRZXK1g8ecKvisPCwkiiXjZyn9uKp6ruK5bYtlNTkf\novfiEqjXTSU1n7GhRNREOlEsm0tk0XldbX7Z7l7x8gzCB3FQiU4/VQoYjUY1yUolKB9Rlcx8tFdS\nUuN107ShpjmNAIorN2VDWyUOlVRICAr1orvNyYktkoC6wG2EpTpH96BlnJ+fV/ObncCb1Ee9F79O\nU9sp1bvtvESMg86MiBDZQa6vrwGg8tDqi1fV1TtDKaV26bpAPYUQO6her+ThSrzCpZU2oolCI0qf\n6zzvpnfetbyIuF2iK5FZ1O76XHed8w6FJqn6kDgqonOy0UZCm1tEPBw1I7Wgyaul11VHiHdQoJ5K\nquTGT7xCJXSgHEzL4+ix9Oe6qQSj77ZtkGsqw0k3chREZW3STj629lV618dQ/4PZ6DocVH1vm1Pa\n1Mh8f5eH3qV+x/DyPgZ0eZZ+XJM97ZDoei+OY6j7PtDkKd4Hmmx0R0t0iUTieBA5aaL9ug+IZ310\nkXI3MFEclzMikUh8HFAC6xrw7aDKH4XZRMeqKWkbSKJLJBKt8KgET0aqwdLusWbwdxT/6Mc37dsE\nSXSJRKIGl848wLu0noerparmKsnpVEx3IgIoSnmbIIkukUiEtjeV4jRwWjdPuAGgRnSl2Soe+qWz\nlaIckJtOcUuiSyQSIalFNjlfuMjTyysZNWWn0fhXpmpjggZ1VHjo2LpIokskThxKctGmK+2pfS6a\nQQTUA6p9phED/nVK5ng8DjNIA6gFkm+CJLpE4oQRqaq+6VzfErG5xKVqKqUyLvbN4HCSHrPRcKkE\nLS+aCreOZJdEl0icKJrIDUCNbLisQbRQkYeJRDFykWToqjD3+7xhrde6oSdJdInECaLN8eDHkpCi\n6V2eZsqnYUaqLFfq8zV7mf+P5Xn9dH8fJNElEieMaGYD4fFw6nBQEouI7eysnolGzydhMmO4S3a6\n1nLX+rUhiS6ROFEoUajqqdKXqowM8i3l2IsyzniCDDodNKaO9Xh8fMRisajqxvNdomu6jxKS6BKJ\nE0SUsYfQBKnqDY2WDXWiKyUZVfIEUIuh4/+U6lQK5HkkvLb7KCGJLpE4UTg5qPQF1InO09hHsxrU\nw+oOCY+pU4mO+3SpUtrqIidJ231ESKJLJE4cavTnbyUv7lMi0/mqJfXVr+HHNSWw7eN86JINKYku\nkUjU4DMSIltdlHkkiqcrla8OC1+vOSoj2u8Byk1IokskEjW4JEd7GX/73FMlwa6e0FL8XbTQlJNc\nRHZtSKJLJBIA6iqsez31GCUd/6/v9dqkQie7CKm6JhKJXiD5eECxo0QubbF4JaeCEpo6PKL1WzJg\nOJFIrA0lEPe+RqQXzaCI7GWlwGMnrpLDw4OHvcwuSKJLJBIrJAc0rxNRksxKZNQWKlIKU3Fb4DrS\nHJBEl0icLLqopqVJ/55oMwoP8YDhJrKLAo65vw0ZR5dIJEI0SVclYov2KzmVrsMQFQ0YVpTscnp+\nW9hJG5LoEokThRKYwrMLNxGeBw57+dxfyjSsZMmySmErTZ7YnAKWSCRWEBEX92vuuJIk52pvm2Tl\n2YqdYF111bLbyu8i1SXRJRInAveYKtl5unPd/FxCZzYwa3AUY6fSoeef0/Lb5rKuq7YCSXSJxEnB\nyYefDw8PmM/n1TKGTkTRLAWqmyQxn+fq11CJTsvXaWYlRMHF+l8bkugSiRND5Ggg0U2nUywWixrh\nRaEfalNTolOvaZeFdpg6vUma6zqHtglJdInEiUAlLJIYSW02m9XSm7s66llGWNbZ2VlN6nMbGzdf\nJ0IdDU6eSqKl7CZ9kUSXSJwAIgmOxDafz2sbCUiXGnSVlJ+apLPJC6q2OSdQz3GnSTk9rm5dJNEl\nEicCdTwsFgtMp1Pc39+vqKvn5+cYDAZVAkyVuCidqTdW7WuRJ5VSny56rRKbTvdSkmOa9U1JDkii\nSyROCiQntcnd39/XbGaDwaCS6Eh43FyC0zJ91oSfp+ory49sftH0r02RRJdInAiikBJ1OFCKGgwG\nGA6H1ebOABJRqXyd2O8zIlTK4/+lYOTSPTjS65pIJFagHlBdYnAwGODy8rIiuNFohMvLOkVwNkRp\nfqySVmnal0/x6kJwfh2CpJoBw4lEAkCzREe7GAmOZOeL1LCcJonOpbnFYrHifHC1V+Pt1iG7NiTR\nJRInhEiNpANApbrIPkZyoVPC/1cJjsdrmIon9NTwFDpINIW7enojVbdPhpMkukTiBBA5C87Ozipy\n083XWtXzeV6b+hrNkdVj3CNLqY9QYvR0Tx6gvGvV9XMA/xGA3w5gCeDHAfwNAL8A4IcA/AqAHwPw\n5QbXSCQSW4KqrMvlsrZ+6nA4rIjOA375nWgK+XBpziU4Pd9JjOfxt3phSWguNXb1ysaTx7rh3wfw\nZwH8NgD/OIC/DuCbAH4JwG8F8N+//E4kEkcAJyt6WZXkNFCX50SxcZ7qXOehRiTXVi9KdQxaZlyf\nTkXT8BTPrNKGdYNUvg7gLwH4B2z/XwfwBYDvAPgBAB8A/COr99VtibJEIrEdLJfLikBIJtPpFLPZ\nDI+PjxiNRtXmjgJ+j8pU0nFy8qlkPsPC57HyPw1t8dkRapdjnB/J+eW/kNPWVV3/fgD/L4A/AeB3\nAPiLAP4ogO/HM8nh5fP71yw/kUhsGZTEBoMBnp6eqk/OhFDV1QnO7W0qten0Lg8K1oDgKBsxVVid\nE7tYLKqgZp2GFkmSXSW6dYnuEsA/CeAPA/jzAP49rKqpy5ctkUgcATTQVyWxiOiIiNw8No7SnKqY\nLJdEpTMe3MbnKilJTiU5fqeaDaDKnNIF6xLdr75sf/7l938B4KcA/BqeVdZfA/CDAH49Ovn9+/fV\n93fv3uHdu3drViORSHSBh2tcXl5WBKP7ozg2j3Pz7yxfiYzODrfjeQwdUJ8qpplNHh4ewrrznG9/\n+9v45V/+5ZpNsXj/Gzy7/xHAvwTg/wLwHsD1y/7fAPCzeJbwPkcg6aWNLpHYP5SwaEubz+eV9BTN\nL42CeSO1lRIZzwGwomZ6mIiqq1H24YeHh5UZFePxuGZP5MawGGzZRgcA/wqAPw1gCOD/xnN4yQWA\nXwTwE3gNL0kkEkcC2rWo9ukUL5XQiJK66iEj7jBQaS7KQqLlafJOxu2pCszreQIASqVdZlNsQnR/\nBcDvDPb/6AZlJhKJHcE9lm5ni0ijRHRapk72d6eBeURXyqSq6hJitBqYq8x95sjmzIhE4kThzgna\nyIB4YRsP0o1CQ/RTHRAO2gaZcl1nOIxGo8pBokQIYCXer2s6pyS6ROIEoVOsAFQ54zy0xCUpnhtl\nCPY4N7fLESQ0SnS8tqqng8EA4/F4RXrzdOtRyEqEJLpE4gShXkydXtUUL+dTtnR2haunep0IJYlu\nuVxWNj3gddaE2w95XNc1JZLoEokThhKMe0O5qZQHYIXo1skGXAoa9vLpMIlscj5drQlJdInECcIn\n0XOfq638rp/ugOgiUUXXV6nSiU6PUbLVuqjK3IYkukTiROFpkYDVOa1N3sw++eAiuKTI63l5pfm2\nkf2vWNe1argZMmA4kThyeIiJS1OlPlzyyOp/u8IuJvUnEolPCCVi01kUOpeVyxYqPH7O58/2kcC2\njSS6RCIBoD7bgd7Wx8dHTKdTTCaTKq2T5oxTaG670WiEq6srXF1dVf/5XNd9IokukUhU0JASSm7T\n6RS3t7e4vb3F/f09JpMJJpMJ5vN57dzhcIjxeIzxeIyrq6tqGUWuEwugk4d0F0iiSyROGE5smnKJ\nk/5vbm6qjURH6U7tb8PhENPpFKPRCNPptJY4gJPvmVBTA433gSS6ROJEoXY4EtvDwwNmsxmm02ml\nst7e3uLm5ga3t7eYTCaYzWaVCqu2N+4fDofVcTz/zZs3uL6+xps3b3B1dVUtkr0vVTaJLpE4Qajn\n9OnpqSK62WyG+/t73N3d4f7+vlJZuU2n08pGp7nsNKMIHRGTyaRKq/TZZ5/hs88+q+aujsfjanGe\nJLpEIrFTaE44qpnT6RT39/eVFHd7e4u7uzvc3d3VFq7RJQeZEMAXuZlOpxgMBpWX1oOAmS141yEo\nSXSJxImilOLIc8oNh8OKpBaLxcq6EJ65JPq+WCxwd3cHAJjP55jNZlVizdFoVEmCu7LZJdElEicM\nn7RPKNFxZTAmutRz9Hh+6iwLXSHs7u6uUo1JcpQEx+NxNXd2F0iiSyROENHkfSJKv0TC88BiLYvn\ncp+vzzqbzapJ+iQ5ZjhWNXYXSKJLJE4UJWmOBDQcDmtSWZsEp3h6eqq8sLpxdsXt7W1FcpQSNe36\ntr2xSXSJxAkiyjcH1BeGptSlyS896aZnEWY5nFHBeLu7u7ualHd/f1+z+TEObzQaVWVtk+yS6BKJ\nE7ZreoQAACAASURBVEUkrZHoGOPGXHDRotQ6p5WSmHpxOYNiMplUJHd2doaHhwfc399jsVjg9vYW\ny+USw+EQb968wZs3b6oMJtucRZFEl0icKFSiI9wuFxEdpTiqt5zf6kSnk/op4dEOxxCUp6cnDAYD\nvH37Fl//+tcxm81655rrgiS6ROJEofNPNTW5xsRpqnOu7+AZShgaokTHMulkYFAyva0aYqKxe1Rx\nKVVuC0l0icQJQ1Oi0x6ndjhdvIYkyHg3Ehw3hZIVVxbT+Lvb29sqiHixWNSI7vz8vHKEbAtJdInE\niUKdCiQVXf9BZzHoYtO6II4G+rqaqVlLAFQSHXPaTSYTPD4+VhLd3d1dNQ/26uoqiS6RSGwGleTa\njvOFcnRNVV1y0IlOyZIe1aurK1xfX+P29raSFJlIYDKZ4P7+HtfX15WKG6VWXwdJdInECUK9mr4o\njsIXtFYbnZNcRHQEbXWj0Qjj8bjmvKCdjkTHFE/uDd4ESXSJxImCBFWaIQHUCUa9oeoVbVpAWgOL\nLy8vq7x0TK9OG51KdHRS0Bu8DbJLokskThCeLUSJziUxPy+S5NpISMNRKNGxbCe72WxWBRKrc2QT\nJNElEokK0Xqr3M/PddZyVaKjRMfQFJ1NoZmNZ7NZbT7sJhP+k+gSiQSAVTIDXtdZ9WUL+xKdTi1j\nSvXBYICLi4tKRdUEoCQ6nQ62CZLoEonEirrqE/id2PqqkiqZcTYFQ1MeHh4AxETnZLgukugSiUQI\nEl5EalEmk+g/3de08LU7RBhczC1KBdUHSXSJRKIGdVC41zU6Vo8rkRmTB9AOR0dDdLzOwdVjeJ11\nYuuS6BKJxAoigitJYxHJReTlUpoTmR6rqaGiDCt9kUSXSCRqKElkUXZht9+VYvI49YvZhinVNamz\nkUSn6EN4SXSJRCIkkyg5Z0Q8SnKl85ToPFbOc+F1IbAmG2GEJLpEIhFmG3ZbmW5UPfX8qCwep7Y5\nEh1/a7qo0rxZryvR1WaXRJdInDiajP8RwZGgVPWMbHI8Tp0PTMDJlOqPj4+1WLmmubMs15ESXSKR\nCNEktSkplTbORXVy9KUO/RxdS2KxWKyQZIncvPy+ActJdInEiUJJSkkrchxoHjn1mOp3ggSkZWt5\n0+kUk8mkSsXUFZpBpS/hJdElEicIJTdflJr54Tg7gQTFNR6UwFRiA8rEw0n7nLivEl2XukYSnX62\nIYkukThB+AwEJS+SGgmJ9jQSnTsbKOkR0dxYZhFWouPUL69TkxNEc+f1QRJdInGiiMhKs/2qLY1E\n57a9aJpWFDis81dJnKrush68Fo/nXFcu3LNc1hfe6ZpgIIkukThBuOpKYlEbGiUw2umoapYcDv5d\nJTO10XHT0BKf0M+NRKeL9+hi2V3CUYAkukTiZOFENJvNKilON0p6VDVVglJnBImKx0YODg01cYlO\nyZBEp1lOLi4uatJinxROSXSJxAmBJOHkpI4Cfp/P5yueVk2dzk2DjCnV8VoaWqJEB7zOhPA1KNQh\nwvViKbU9PT1hOBzW0rd3sdkl0SUSJwK1m6ltjhIUVVaqsCqVaWiHruWqKiqnb6kHV9ek0LI0k7GW\nRynt4eEB0+m0IjiV2rheRdM8WUcSXSJxQog8ppToSHS0n7mtjQTFhakHg0FVJvCav07P5f8er6cL\n7JDkmFodQCXR+UwJnsMFrpPoEonECpS0SEi0i5HsaIvzifu+7sNoNFrxsD48PFQe0uiaBAmLa0cw\n6zDVUqCe8eTy8hKLxQKDwWAlk0oXJNElEieIKDxEw0Q0k4jGw41GI1xfX+Pq6grD4bBmw1MPqk/9\nood0OBzi7Oxshdy4kfgo3enyiJT+mhbNLiGJLpE4Ufj8VpISgIpQlFTOz88xHo9xfX2N6+trDIfD\nWiiI2uN8fiuAqiwSWER0JDuqsirxaZ36kl0SXSJxgnBbnXpGqTrSFkZiuby8xNXVFa6vr/HmzRsM\nh0NMJpNKZdV1Wr1cEpOqvdyc6JTklNycfPtM7E+iSyROFD7dyr2jKuExtIP7OEuCsXaTyaQKTaHX\nliowiWs8HmM8HuPq6qqR7FxqU6Kj/S+zlyQSiVZ4qIn/p0sQelgI95+fn+P+/r7abm9vcX9/X00b\no3dVbXvcSHDj8biS8qi2Kokp0elUsFwcJ5FIdILPTdWgX376LAad+H92doa7uzvc39/j7u6ukupI\ndLShqcr79u1bvHnzpiI5l+wuLy9rUmUpqLgvkugSiROC53HTmQ4kEjoTKNmpVHdxcYHZbFYdR2lO\nZ1TwHKqsw+GwsuuR6KjGqupKotMEAxo7p599kUSXSJwYSB6DwQCj0WglqaZPAZvP5yteThImPa46\nSZ+2NNrj6Kl98+YN3rx5g+vr64roxuNx5XSgs0HJl/XtmqWkhCS6ROKEQKJgsO5oNFrJ+QagmgI2\nnU4BlOeTqiMDQEWgFxcXFbGR3K6urqpNic4dDCrFebqndckuiS6ROBEoQZDo1B5HwqJ9jes7lHLO\nsRydykWb3Hg8xps3b/DZZ5/ha1/7Wo3cRqNR7TenjgEoeoDdGdIXSXSJxAlC54wCrw4IZiPhdn5+\nXssh55KfzlFVB8PV1RXevn1bbVG8nM5tVaLTwGOV8jT/XJ+km0ASXSJxkqA3U0FpTOe0vnnzprZ+\nhC91SMIaDoeVWupq6tXVVc0OR3Lk9TWRpzpLuhJdF2xCdD8F4A8BeALwfwL4cQBvAPwCgB8C8CsA\nfgzAlxtcI5FI7AAkOvVkUiojcY3H4youTjOb0HlBoqPnVMNHOEVMg4D1Wj59SwlLiU/Jjv85yXUh\nu/UUXuCHAfwPAH4bgBmeye3PAvjtAP4/AN8C8JMAvg/AN+3c5TqLW5wq3Bjbhuilr2vX2DX8nrye\n0T0f6718zHBbGHPT+awH5qrTtOoaA8fwkbdv3+L6+roW5Luu59Rtgk1lvOwLC19XovsKwALANYDH\nl8+/g2cp74uXY34ewAesEl2iBW6YjRYK1hcPvK512eaWPwRRKGG1EbfWzwNHPQnjp0Z6fQSAbd67\nejoBVGorpT5O37q+vq5N1F8ul7WJ97TNMR6OEtymoSHbuNd1ie7vAvh3APw/ACYA/hyAXwLw/QC+\n83LMd15+J3qglAXWc+2z87s73ic9e9LCY7ivpoSJvh6BpgDivUVqzscOfz4l7OredcCkBEbC0wBe\nH3Dd66qT8bcx0KqtbhOsS3T/IIA/imcV9nsA/nM82+sUy5ct0RNKBprqOsrfz02NvGq09e+HJAYP\nHYhAQmZ0Pu9fDdfrzHX8GNA2CBC7uHc3+OtKW20kHM222CYhH1Ki+6cA/C8AfuPl938F4EcA/BqA\nH3j5/EEAvx6d/P79++r7u3fv8O7duzWr8WnA00xrehtdFUnVBlfnPB21p7XhRuxaKoo6h3eS6D+d\nb+mJGynVRUGk+7gnxzZtzV3L6nLcJvd+SBNHX3z48AEfPnzodOy6d/M7APxpAL8TwBTAfwLgf8Oz\nt/U3APwsnm1znyOdEY1QyYUbp9ToepoqzUVER1uKEp6qEdEUGx+Jt3lPJekkmkiuRO/3yP20AdEr\n6HFX+7ZHdpG8dlVm0/19TES1bezCGfFXAPxJAH8Bz+El/zuAPw7gawB+EcBP4DW8JNGCKG//dDpd\nkeKUACKiU5KjW//x8bFGDABqxnyNSt822Xk0e1S+kp2SPNMEEZRaXVL0Jfj0+F2hSUrdtNyux1HN\njHCKJNeGQzyRlOhQV9FIblxXU5MY6gro7oF1olOy04wQ+snFR0oOi20gmrajNiAlQb0/Ej3JTs+j\nR49eveg5RBHzu1LLu9rT9o1d3/sxYxcSXaInPMSC2+PjYxV1zlxe6njwdTWjEBItl/n5gbqdS89V\nW90unBRRR1Nid3ukepXVFud1UvWW5flqU0p8Xqe+KJHYMRKcQ+t3aoQXIYluj4g6O9W16XSK+/v7\ncPFgzw4BvLrdvXwlR5atCQw5t9CxCwnA66jSnNZR75dePz+Xv0l0HjLD//QZ+bPqc39dwj2Oleyi\n53bqZJdEt2doZ9d1NSnR6eLBni1COzcnX7tRn9dQKY7xTrroMP/blYqjdrnIrubS3Hw+r47XcBi3\nhz0+Pq5MJYpUyYgo++JY1dMmNIV/nDKS6HYM7+Te2T10gmTlNi3uK3Vwt+HpfxqLR4O+EyelPl5n\nG1ASjeLn/Flo6m3NokHS5n1G3uJIxQcQqvhN2IWTYddE01bXVGOT6HaKkndOJToSGyUu5gFT25oT\nHe1XLuloiIoSi67cxDg0JVaSjHrzthGL5d7BNiO+hsnwOXAgAFCRsaquqtIq0Sm62iH9XW1DmttX\nuEdUT5duT1mNTaLbESKS0093FFxcXGA0GuHi4qIWTqLQ3P7eWOmlZNl+LUp07qV1FRjYTqeMHAlt\nqqCHyXDiuKvkLo1GToqmujShNDhtgl16QJW4mshO63KKSKLbMaKO7SqbehlJZiWiK62E5EZ5JQgN\n39AwFGaY3ZUdyqWZEonofessDlVd3VsbEV0k0amK3wfbeB77CvVI1bUdSXR7QjSy+opLy+WyZqPy\nxIhNS745majtykkw8uZqGevek6MUIByp7fpJ26X/r3nQXPrUrLdNwbRdEIXFdD2+y/7E/pFEt0dE\nUgzJjiqkqpbeyXwRX0Wktj49PdXUXHVauJTTV+posmOp7Sw6z4nOA321jm5LpFdan6GW6/v7wkmu\nK2nuyxaXWA9JdHsGO43HeLlDIOowPlFfoR5atwGyPJfm3Eu5jkTXZnOLzlGiUynVVVBVZynhcaqc\nDhTunV5XFS/Vt8/5SXTHiSS6HSNq+Epo0VzTkl1Hc3850ZE0ojTVTU4GV2f1+pvC7XFui+Oz8Clo\nqsaX7FxKmLpfn2HfBVQidDlvH3a4xGZIotsRPMSC+6JQCwC1zl9SvSISI3yBYcLDCSLCoEqrRLRJ\np/WwBr2uk5yThOZBi+qtdSPReTn6nNa9l+j9tR2bOF4k0e0QkaQReWDZ+dVmVZIQmtRanajv6pye\n2yYZ8dxtkF1kl+R33m9UP62/297ULkdHjqefiib590Wfc5LsjhtJdDuGd4AmVdalnz6dx1VVJZqS\nOkcooWxiyHdENjwnIz2Gjgf+r3n3eK47IEjQpXCVTVXKT4XATl21TqI7ImziuWOnjzKBROrcNtIa\n+bFKNqWgW79H9QB7nbo6FZQs2yTYTwH+nJue0SZt6lNCEt2RoIvk1wRKZD4FrKTORQTXhxhK6m0X\n76dKjCxH57qybirNdX0GpZCZT62j+0ASPWeXhE8ZSXRHgG00Qs0GwlizJonOCXBTiS6yPTbF2Lma\nXXJW+OyQEiIP76fewbtIvKfwHLogiW6P6BqT1bVRaliIrjGhwbiRvapJle3bMbqoTb5PN53HGpXp\n9S/Z3bjPA6SbQmu2iW2/26ZrtMUvJlaRRLcnNKkYEdo6hNvkNAW5rrfgnliWrdLTuh5KV1O97pFq\nG6nJTnJKfiUPqm8aX6j2P63Prshu2++2dI2+10m8Ioluj+iqanSRQkh08/m8ylCsKdgj1VSv4ZKc\nEmLfe2kiu+j+SjbBUnklsvP70RAbJcxN4wK7YJvvtqnsJLv1kES3B7hXMLJnEU1OAT3u6empWmtC\nVw2jfS5a8LmL6rruvUUoqa5+bikpJ4+Pnh+lPT0ukuj0/F2QXdN7jbCp6poEtx6S6HYIt6loqiHd\np8dEap1vdDyQ4LiCGNVXX1NBJR7Old3G9Kgu9064Osn7KKVi8vN98SBdq3YfUk6p7L7XXKeOfQjO\nJd1Td0IQSXQ7QkRyvpU6N4mJ5/M/X/hGSY5qK0NLtAxgNZfdNqZ6dbl3BYlb710/fe1aBdXz6XSK\nxWKB4XBY3Vd0rV1gU0JVKXRdwusCt0km2SXR7RSRJKad2pc09En1Z2dntf/1PBIdVVeGk6hdyolz\n3xJdZI9yJ4oSHfdpKiaCavpsNqucLefn5zWJbh/YVH3cV11TmqsjiW6HaCI3Xd/B14jQkdglQC2H\n6dE1SWXJ00lS4LYp2bmUyn1O6m5/1OUNNbDZyU4TbC6XyxoJqpofTX3blYfymG1k7lFPkqsjiW5H\n0M7rHVg7uv92oipJhD4LgogauhrqB4NBLafduuqr319JLXd7pN+zIlJho2Uf6XCI4ur8eXn4yqcI\nDx1KrCKJbkdwKYUBvbrpfnpMo07JfW7Pc29lqbGrRKdEp57WdSW6SCqNtuhZeBiMk7kee3Z2ViNo\n9xwT7mk9ZilsG4hILsluFUl0O4B7R1XNJKHpxv3z+byxUzapZJHq4qEk7ozYRMWhOqnqcyS9Ogkq\nwTMMRtfA4D1qjCCXaRwOhxgOh7i8fG22pXAVN/pvalc7drJMdbUZSXRbhEpYLq3RmM7v+r9mG+nS\nodzuxU+Sl64toeopCYRzSF2a69NRHh8fMZvNcHd3h+l0WlNbXR2PHCr8XQoa5nclLn2+kaOHqr+q\ntno+77ErKewjbCWxHyTRbQluR3MpTuPcdKrWYrFYsWV1RdR5XXpTaUkz8p6fn6/Mh+0Den1vb29x\nd3dXewZO4q5mu7qtHmKWQygh+nPSZ67n8VidIcFn1VVVV5JLovv4kUS3RUTeQ1dNGQemdqpNJIcm\ngvMVw9xTSi8vCWFdie6rr76qSUkPDw81yVWfT6RWk5Si50kp1G14bUHXLDeSVruS+inY+E4FSXQ7\nAqU6VVlJbFHkfxsig3MkyXniSvf8toWUtKmxpXpHEhYdIF5f/R05E0huGjjsUhxteK6e+9Q3z333\nqdqyInNG4hVJdFuG2oZofGeQq9vhVG3zmQBann6P0i1FZAegSt+k//m50fVKEo9Lnufn5xgMBhgO\nhytrsZKAlsvlyrWjcvm81JOrcXZ6XR1ESG6evUQTFWzLAXOM0LbEZ99XOj8FJNFtES4Nkeg4/zTK\nEwfUQyKayi3ljivlaosM/kp6HsfGazR1FFUPqRoPBgOcnZ1V4TH8TwlHbYVaThSKos4MnyGhdkA+\nYw07GQwGtVkgKlHuI4vJIVDywCdekUS3Zaj3UJ0Sal/SUdg3V0GaQkQiCc+vTzIpqa0lG1Sb5EMy\no0QH1J0d6vX1RbdVxdXAYbfFAa/2tJJzxZ0dZ2dnK57d6P70HXwKSNW1GUl0W0akml5eXtZSm0dk\nxu8KVQEj0vCJ+W7DUocDSUUJmMHDg8Gg+q5SkJYdqdCXl5cYDod4eHhYIWUPb4nSQFGtdefJYDCo\n1FZuvNZwOMRgMFh5fryOPx8AlT1PpdvSAOHP3+u7TfQlpHSKrI8kui0iCp1gp2uyw/FYVxl5vgbK\n6qYE4gZ82rgIl+4uLi4wn89rNrbBYFCzv7EObtDn/yQfTQulJOj54SJiIbnx+g8PDxiNRjVJ7enp\nqaonn4OH46gzQp8N71ljCN1uF6m0JVvXrsiuifTUXJBktx6S6LYM73zsVAyhKC32EkkVSibj8biS\nvCh9aYdVNZAOCNbHwzFUonFHAuuiC9RoWAhBQhkOhzWHA8tVooukXCUmnRtbmt+qEh2J0W14rCOv\nrbZSSp2s93K5rN1j6Z3oe90VujpIkuTWRxLdFsFOzE6kEo8SQTRX1ckAeCaT8XhcbapeukdRbWMX\nFxdYLBbV9TwlFKUbhZIN68HraZ30XnmPVAlJzCrJ8vosX4kOeCVSJRuW40RHktd7jzKk6DMCyimj\n+rxXRxvpdLlOk9rcVmaSXj8k0W0JqrYBqDrmcDisbEyURKKgV3dQsCOPRiOMx2OMRqMVKU43tXFR\nNeVvkp7PVFBjvs/SAF7VXVdfCY2To2MiIjNXo6MO6yqz74/skiUic1ujmg02SU/VVY3tWm5fgvPw\nkSS77kii2yK0g9PmRMnJM3noRkTG+dFohNFohOFwWJvD6UZ1SpFKqkpy2lk0SaerdUoaGipCUlGw\nrmdnZzX7Hp8F6zWfz6vyPO4r2u9EwXvzwOCS9ziS6Nzm2DemzglGr+vHdbG76fFdr69IkuuHJLot\nQhs5JTqdpxmRnMbVAViRzuhlHA6HKx1ISUElw6enp0avLPCqqqp0pyTATQ3/0f2qwyEiDp3Z4MSu\n98LvJRJy6VU/tXw+ew9n4Webp7UJfVTHdcrf5vUTdSTR7QgkAQbTug1MO6YSl3ZyDdFQ9bFEeJF3\nziUMkheJOMpuTJU2clLo/SlxRv8DdYeM2t2AejZkPb4ktZbITu87uueoXvrftux20fUSx4Ekuh1B\nbUh0QETzXL1DO1xla1KNlHw0xMLtfoyh82y/OkPD00fpNfx6hBMiP5Xkogn4rJtfJ7JFKuHxf5bL\n+kRkqYNBdC99oOprE+EljgdJdDuCehU9hAKIZz3wHHVUuETTZgNSo76rsJTmPI07Y+poq+O1fa4p\n4bax6NPvj6TvEm1J/Vbp04neJTs1D6hkx3ooSfo1NiG7xMeDJLodIlKh2ohOO2CTTamLJKGhLkqA\nqkrqnFCd4QCgCmmJnCBt9fA6eVyexhW6tOWe1kiiVduj/o6eWUlF/ZjJKiXJfkii2xFU4lG1rKsd\nSUmvdF7XOqhExbIpXZLgPN7v7Ow5mJhhLV2u6/+7dKf2NEqcZ2dntTmu/gxc/XbSKtkidV8TwW0i\n1R0SSXT9cIinsjwlj1Gk2pU6XqQObiKFaDnu6fV9moRTJSmfN9p0fbfXld6zqq/qidZ7jMiqjeSa\nJE7dtw073aHR1Ic+1nvaFC/3Hd78UUh0TS+tFNbQ9PuY0GRPU0TE0GSD63tNtfOpnUo3oK4u9n2u\nTkC6P1JP3QPrx/q9sP48XstTCbBrPT9mfAr3sE8cjOiiYFHf12QIj2xcHxNK91zq/PzsqmqVCEP/\nj67VN5C26foKjd9zSZJBy36u2t4iG6bXNzt/ooSDSnTRXM9IzVK47UZVqo+loUfk5vNf1YMYEXvb\nvfpziex9TnJ6/CbQ67k6GxGcLt7t9+Axc6WQk4/l3ScOg4NKdNrwtdNrJ3Ci89ACZpT92BARezQP\nNpJcacxvgx7TZqtSbIM0ouvo/UbZhD31ks4QUc8xy+9K+onEQSU6b/RRqp7ISO3ZK3hsU6DooVFy\nDET33kZ0pVAL/e2xasB+bZtd7tcXuPYYN58HzDJIgBonGKnqx/DePwVEZhb/L0JkZ43+3wcOKtHp\npHJdAtClGmB1SpEGwGqyRZ8HqcRwaGiKpGg1e8+q6+pmpLa7WhdNvm+y1e0CLq36mhARwaskp6q1\nemb1/lXai973sbzzTwUlM1NEfsBqe9V9/B5FFewKR0N0ushzqeHrg+N3TUmkwa+METuWxq4dn/fL\nhXN0loLa6BxunPeMw0wxzkbEsg7R8V0ld+ktUtVZd62nS3U8JkoxrxlL9t2RPmVEZiY3OfE4wjWv\nSMM6CYkOKJNd5ICIVDUGvqp0xxfAyfQ8TsvaN6JOz2UQSXiRnarJ48z7ZTJKJQq143k4yaEkOs93\nF0nsvkVOqkgl5jza6Fnp70R/RDZk/SxFRwD1LDKamMKlvH3gKOLogHracB8piChkQaUfRthH0Ae9\nr4frnZJSjZKbTp53O1VUnpKAT6HiMVwCkBJtpGrs8hmUPMVs3Eq+Xhd/PyV1KbJJcqDQBAKHeO8f\nO7zdRnZVldCV9Ag+c9W0fBEmzym4SxwF0amEAqA4UrjNivsA1CQa3a/n7vPBRo1FyU3JrmSriohe\nO7mTHK9zeXmJ0WhUqXhaJ2+Mu0TkLFHVsm3TZ9Am2ZHo+Bwo1e/7vX8qUBVVCc7br+6PiI4b1z25\nurrCaDSqjolsyrvAURAd8CrR8eZdvWmSctj5lfQionM1dtfQeqt67uqqL01Y8mS5mqcSLM+n3ZLP\n0UlzX+qrOo/chOBOA7e7qtSq9+62okjKiAaLfb/3jx1udlDTA81L+sn90cDM93l1dVXLeajt85OX\n6CKvDCd7u1HeG71+dxIksUR2mn3aB5TkIi+j1omqaFt5/qnk6KQf2eb2obo6IidSaePxJTMF8Erq\nLvkpydFUoNc/lCH8Y0Fk/6Qm4gO0SnP8L+qrfOb6H8t3CV8/t42DEx1vVjt8pLa63YafkTSkBAO8\nqrUkUSWDXT1Ylz6isBG+ZHbaklHXy/SNqn6bbetQndtVlBLp8dionkpoTcfqcSQ6t2WmdLeKEsmp\nJDedTlckZ5euS9qXhlbxuPPz8ypSQNvALtBGdP8xgH8OwK8D+Mde9v0mAL8A4IcA/AqAHwPw5ct/\nPwXgXwTwCOBfBfDflgp2Q7IbqEv6vndeV1vUK6QjBzubLuKyD4nOSU6NtiplRKmKmsoj/Ls2lmMg\nu4hYIjW1y4geOS08+QCfpaaH9zi7Q5L+McPNAx4hMJ1OV+znelwULqSDsTrkSHJXV1etg9w20EZ0\nfwLAfwjgT8q+bwL4JQDfAvCTL7+/CeAfBfAHXj5/M4D/DsBvBVBcnr7r6KoN1V+GRsZ70K2SJm2A\nfOC7Vl+9nj7KeciLSpul8kjawGoAp6qrPF4/9VqHIjslsxLRdSlLP4nIdOEeQ3/mSXYxfIBWO2hT\n6BfP1fO9LBLicDjE1dVVtUodA793hTai+58A/LDt+70Avnj5/vMAPuCZ6H4fgD8DYIFnSe9vAvhd\nAH45KrhNjWjqFKqKqhODqkqkNqo94ezsrNG93SRVNdU1uj+VRnXtBD3H7ZFaJ20wnuVDz/H7jgh2\nn527TTLbZriHq8YqaQDPz5A2JIY3HDq28lih7Zbv6eLiAsPhMLSD8hiuQeIOItrzgFcTDR0Zk8kE\n9/f3WC6XGA6HYR22hXVsdN8P4Dsv37/z8hsA/j7USe1X8SzZFdGV5PymKf0AWIm+LqmLuj6CSnMl\n93ZJIvJ6ttn6VG3WVbCiUZD7ldS98Sj4DFTSK5Gd1mXXaHsm7hjYxvX0nVLSp1fPk3vyfbAeqco+\nQ/udviMdnNUO9/j4WJmDIpKbTqcAUAv7IelRFZ5MJtX1dxlMvKkzYvmyNf0fwh+q/1ey3ah987Gs\nyAAAIABJREFUK2qgSm48nuEWupgzG3tTI4+Mqk33UtqvEp3fY3SOkrGOiqqeqg3Ey3K7XJe6bhtt\nxL8LouMAqKqqkpyG3ejgkCRXfwYu0WnYl9vhSrY6LlpOe5w7Kkh09/f3NUGAcZ/bVmPXIbrvAPgB\nAL8G4Afx7KgAgL8N4Bty3G952beC9+/fVx3wd//u340vvnjWhL3hdekMLulpJDYbPAmCL0LVSA80\nBuKQleiabeqZkngpzZCX4USn9XFnQ/SMXFVdxw62TUT3t8vrKOlxMNRO5rGZSXKviNq7S8sa4+ra\ng9pD2bdUhdXwqsiO7tdvez8fPnzAhw8fOt1bl7f8wwD+G7x6Xb8F4DcA/CyebXOf49UZ8Z/i2S5H\nZ8Q/hFWpbtnkDWwjuSZpxe0C0+m0ShSgBEJvDzdXpZrU1nXqqZ5gDQWJzvMAWo1Z4v1Mp9PK5sRN\nz+P9XV9f4/r6ulrJq5TdZNfw57lL0qV6xNgunaq0XC4xHo8xHo8xGo1qnthTJzx9R56EwZ1c+j1y\nOJDgbm5ucHNzg9vb26rNUtL7/PPP8X3f9334/PPP8ebNm6q9jsfjlT7QFS/Hhie0SXR/Bs+Oh78H\nwN8C8NMA/hiAXwTwE3gNLwGAv/qy/68CeADwL6NZrQ0ruok6oaL2crmspLoo/ZGuVK/XdrJbt559\npFKFq7n8rtJcm8Spz+LQEh3rEg1mu7qWStPeWaJBNiW7Z5SECDcv6bPyfsXN58SSDJnYQcvddQwd\n0E50f7Cw/0cL+//tl60XtmW30RdCMqNdzkcfHbFczexzva7/q/jfpBK7uulEFRGxk2pbZz8E9n19\nn3Ghz73LQHGKKNl0S20RKBPd+fk5RqNRTbLWvhi9n122kYPPjNg3+CI5tYVpoTZR7boQXpNNLjoO\nWF1IxoOF9TyVBEuLPu/zeR/KCcJrlTbWjRJykl0MDsxRG4raLB0Oaj7SuDsNU2GYT6mN7qK9Hozo\nXJ3ZNVws9xFG89c11Tna1+c+Stfw/W7j81HT6+Sxet6I9omSjXPf9WgjOo0DS7xCpTe1J0cmG0pn\n2l7pdVUb39nZWWUuorYVJeTcVXs9eYlOX9y6toI+JFfa7x3RSa4p5Ti/K8m5of0QZBfVcx/oovar\nEydRB58V31/XNkR7d5tEByCU6Hz2zDZxNGmamuASQmRjiaQfJwl3idNWpyqf2tL8YXP0UmxKcv6f\nk5sGYjrZaQf2dPKbkFwkMfY5LzL0R9JA27PYNnwAafKuJ2KzChH1SQ3MZziJkiXj8c7OzqoZKiUz\ny7bxURAdgJWGGT3oiOBK6ZEid7jaJUrSnUsBfdXWJmgoSpRrze9DR0IlOm88vHaXekadv+95kffO\nva77JDuvj9rnkuT6Q5+b9yFKckp0ah5gvxoOh41kt218FESnDZS//T+XhJwodKYExXI9/vz8vJqq\nok6J6KFHoR5tL6fJludk3ZSyOprS5fY5tdGt03j0nvqonxHJRSpsaeTue70+aKpXojsikmNf81x1\nOiCr4MBpY0p03k5PTnV1EvPGqWqIqqmRFBSdx3O84/HllCQbemr7oCvJuQqgSQ79OTRJc31Jromc\n+pxXkuhYXmQz439tNrNNOoDXJ8muDn8fipKZSPuZtlNVW90GR6Ir2ZN3gaMhulKDc5XUj/N0PPrQ\n5/N57YF7uZEk2OQB0hemZax7r24vYsOJ0q37Ajqsj9sXu5Bb07OOCKBUVpONywcFD41xg7eP/D6y\nc986nUGfU2R37UvsnzIiO5z+dvuxryPB4OBIWue7VSlO2+0ucTREB5Q7jhOC7icBePwObQSe8cPL\nZYgBXx6wOgIp1NawiTSg5KYNyHPze3S5qpWRfa5E0tG966fu79LxnawVrqYvl6+G6qgTsIwogNTv\ns68a7mVE77QrsX/qcJLzgQjAyvzVaFBmn1JE4U/7cEIQR0V0wGqj81FEwX18wP7Qo/L8OqoyEnwp\nnqyRYAPYRKJzSVVHSV3MW4mc1/bOv65NLqp/pHJG5BQd689Ij9GpePzfJUjtDN4BNN1+X+h5kdrc\ndI+ngjazChENxtHC85FE55EB+5xrfPAFrPnZtGlnjyS6JhtWGxGVrqmS1jYMpFoftSW6yq1krU4I\nbTyR1Oajccn2GD330jPpQyqR2qvvz50oXmcnn+i5dyFwkmrkaY/KOFVyK70rYDUiQf/TpQ59sRyf\nBha119JAtmscjOiayMWlnYjoXFKIkiu6+N107dL/RIlk+rwk3l8UOuL2Ds8c4eqXN5bouUQNjsfo\nZ9s7akMknfp71HvQuiiaCFn/d8eSn+NE59dXs4R/PxU0vTOP39T/neA8TpW/l8tlTYLj8z9UEPvB\nJTp/yCqNuLTjdio930dv79Sl65ekG5atUk0kgXR9WX5vTmzuYaWNUa8XTYQuXSOS6rrWVcmzC9Hx\n2qVMMdHzin63SfZ6j231UMdNdLwT3Dak9o8JbB/RojbeLl3zUI1Dn5v21eVyWYWQaB86OaKLRnwn\nOSW4aGaAd2pt1F0eYmRHYN20jtFx65JcyWulWR6UMEqGeZfmXN2L1LVSnXVf1Pm7DBY+OEUje0Qu\nui+SMPwaXp+obnr96D4jqfxUCI6gSUG1B42HY3tk5EIU9sT2SQ1D2zTfy/l5edEnfa+fpNc16vhO\nAhGxES5ZlTqiNmJXg1wF7DrSROpOdHykSkb3VwoM1vt0aU6vqYSiDpTomUWSi+7zZ6L3Ed2f348T\nlAaIRoGhXh/WW1V8N0fwnkt188FT79+fw6kRnAoIularJihVZxiJzt+Jvhdtm3qcX08Hd2/DLpzs\n4t0chOhKxOZG+sh+QNKKHoo3/OiB8hiP4+lCXooux6i06Sq2S6oldYv1UomI9VQph40uuteo3pFE\nxRgnHtcmyfl9lFRVn61RIjpXVSkN+DOJVFp/7l4Xl8xPjeSA1UBfrsQ1mUxqqioJcDabVdELQPyO\nSv3Fp0lqW/E6sWxth9y3LRyNRKcSjquy3pi9wUbQTkxjqBKDurkjKaeErhJB6R7bJDonuUhdVftc\nF7WgVH+NVNegTaqRbfawSPpW6UvLZzR818GEqhVtjCVHh9evRIIlaeFUyM7bIyW3yWSCu7u7mqrK\ndP3T6bRKZVYaoIBVzcOXnlSi8+dNsnRNYtuq7MGIzg3mkVPCG2rp5t2W5qqYqkV8CZQylOj0Ov5y\ntZy+YRcu+bhbXtXV0r3x2UTSSemciKic/HWU1vO73JdeS8vmd3+WTqZtJKMquF/P20t073qd0sZj\nTgHaFrmmyv39/QrRefAvEb1fb5v6H9u8Hsd9l5eX1afG5Pn6JiVCBeoL27f1yaOR6EqSCbB6g/qQ\ntSyfrdAUaR8RXURq2iH6Bjdq/XyqjBOdk5ffW3T/us/P8+8uyfnz0fta5/6aRvuIWLpcQ+vqarpu\n0X3z0+/L3+upkBxQJrrb29tie4ycSNpeeKyaTpTo+HwfHx9rpOTBw0zKSROHmjqAeMAaDAYYDocY\nDofHS3SuzjXZXCK1w4kuMsBHL4cPz4nOt8gxUerMTffpjYt2D3XRl+7XB4SmukSSj5epkm1Jcu3T\n8Z3olFQiabzPNaJ3ET2TqO2UVNZowDtVoqMERaJTD6z3F2BVS1JPuRMd2xi93ho94OeXiI+kF/VL\nHj8ajQC8OruacHCJztWQEqJGqiN2JI040flo4QbyJjvEJveobnlfLMTj3HhfXg5RIn0eW1KB/V58\nZO7qdS6B56k9tPQ+u6jGWm6TtKvtKLqmjvSnTnLeHul0mEwmNdMK7al8n0CZ6HwgjqICIolPy9Tz\nVAhxm64T4nL5nJKdhNeEo5rr2mRb8mP4vUSQpYcTpYYpqaptHcE7S5OqqvYP9a76vemnw6VSr2fT\nc9NGxgzE3K9ldp0zq+c1PQvg1WnRl2Dcjuu2QJUkI5W9rV2cGuHps/QwIJIG371qPN4/tO2RFBnW\n1NQ+WQd++mDF4zQExdsmg5B5zTYBiTgqogNiW0v0P7eIMAiOSG4LiDpzqXN3VbO0bqqq6sRndz64\n5Ob36WADpME2qlvTs+CnOwV0MOgj2bUZgVWtjYiuqfySisrrqrmCZSmxNj0LHn8qBEdEz1Ml8fPz\n53Ae7SsqGJTaHOGqKa/JT+8fLqhEAoxeezAYYDQa1dpUV+3gKIjOyaup8q7qNoVBaCf2l+Yj+qaq\nm3dMVQ80+FIX8G2T4KL7IdENBoPedWQZ+j0iui5kpOeRZEobR2sPFm0jOn2WSpolFVSfY5MKn+pr\nvc/ogNHkJGCb03cbObSU5KIBS0nO+zFDX6hG6+A+GAyqtnRxcdHJ5EUchOgiIzhQj/In+D8buT8g\nfWEutajI66qqS3PbaviRLaQ0c0Dvse3z7OysetnMt78N+ADg1yxBnx0Qr+nRtHV5zt4RvFOwHoR2\nWL2G1+tUoQP65eUlhsMhRqMRFotF7Rk50SnR6KAGINQE1H5HONH5VEcnOnWMaH2Gw2FNne7aX4+G\n6JpUUe1UbdKflh2lFu+irq4LF70jiaTpmbiE6Z9qnyDRda13iVz9GfWBj95nZ/VAY4/982fUBg8j\niZ5lSfp0Qu2iLXzqUDvXcDjEeDzG9fX1ipSnfcclOtcI2voYj/M2UZrm+fT0VEu7ptfiYK+Dc9e+\ne3CiU6Ml2TwiL3YiQo/RRt+F5ErX3xbZReJ6W5xgVC+XSrXR6VStLvX2Th5dMyKkEvyZeTxbkye0\nC+FE6q92jEgSV1LzcJymwfEUoLY4Skbj8biS5tw76rYxtjvvLxrBELVfHwz56f3CJTrPhKKDJst1\np0cTDkp02tGA1U5C6A1G4I0qOUQ2OR9tdmmncbKLVGytvxO03kf06YGcXeqiUk7pGbjKEdXX97Pu\nni2kS12i/0vHRqqvS708T9XYqOxTBNsUgIroqB5GE/XZDlWic42oKZJB27TDB7FoQIvmgut5KtW1\n4aDOCH1gl5eX1cN2vb3JpqPkFYWRlLysuyA57Wz6AtjBos7HY11d0PpH3q8ukmikLjaRi9vnuo6W\nPEdDTUh89OKVpCuvR6mO+mwje5DW1TvWLgayjxH6DIfDYTUw0civEp0+Vw89cm2gzQbepU7UCFg/\n3V/KYkQVvIvJ5WBEFz2oy8vL6oZ9cn9EePoinCyc7HYtxZVsa8Cqah09h0hFVTW1aV0Iv5eI2Hyf\nP1O3een76dpYgbp0R3ODlunqi9fb1RSvh3aELu/01KU4wtvnYDDAcvnseBgOhyvtpPR8o/YRaWhd\n+pceo2Sn9mKX+PQcqtWRxOg4uETn6qY2dqDdnqN2hyayi0hh22SnEpur5Hq/Uf1d5eZopXa5iERL\naFIR1bCv9WN99Bpdn5FLmNx4DSX9tnqxYfMZ8nuJ0LTOvr8p/OjUoO2T7W44HHYK0SiZPbxN+qDW\nJKGX6tfU7qKB7ahVV1b04uIifBjaUZpGFPfCuG2gZLDclTqj6hMDHJsM83wGag9xac69on0lUlcb\nS8SgZZf+73Id71BR+W2gtMH32Ca5+X5KCC65RrGCpwR9R/zdRHJd7btOPOs+267vuW/5B1ddAax4\nEFW60XATHeXZYHWmgKt1SnrR9XcJjpZv3ryp7oObXt8l0JJ9cZ1GVGooTfBntk6j0ndbIjqVtoC6\nw0ntNVFwd+ma+kyfnuL03H3jBT9l8Hl3kebcnFCSqvx5dnm2Kv15WdGx6+AQb3jptqPlcrniXfE4\nG36S3CKDfZcHr9hlA/ekmprjSzv12dnZitRW8l71rXOpAZckzEiaW+cZNTk/+K5LQd96vBNdCdF/\nfl96f10cOaeA6JmXjnMiKg08/ruPyaCkBTha+nP450ElOn66DUund7ir2dXVKKr/GKD1fHp6wnw+\nx8XFRRW3pB3ObYneCdUb1tZw2kbByKGj39dVDfx6JVsLByutR8mG06Qe9anLNu7pU0SJ7PV9qCBC\nSdufZdOgEqm0hxA8jmKuayRFqOpJIouCGY+5AesLps0OiMnGvUvqkSzdY1ND6iO5bEO6iTpH5PFV\n6TxaH0SPVyl33UHtmNvHMcFtceoYYlKK+Xy+Ym5QzYvtnH2WUxXpVPP3t893cxREB6wasvnpIQc8\nVonwGBuzSjV88brfDbxKAjoHkOdE5SuZuXs/imdqkva28Qz5fjzavdQxojmPSvJqs6RjZzQa1Z5l\n6fmU7jMRQ/uZvp/FYoG7u7tq0xhXzdDDtSV0PvbV1RWur69xfX2N8Xhcvb9NTSPr4KiIrqSaubSg\nn8fcgLVuHMVULfepYTqhmVvT/UX2PA9VYT20cUX13NZg4R3GB6toUSBdW1Sfidpgx+Nx1ZGiuifW\nh0txmih2Npvh5uYG3/ve9/Dll1/W1nllOnYupnNxcVGR2Xg8xmeffYavf/3rtalcOuj39epvgoMT\nXWQM9//UlqNk2Gbn0XMPAVfZohHTU9b4SuiRMZ8oOS50fiKn7fRxaqzzzPQeI5WU/zuRK9Ep+ev9\nXVxc1AiSUp3PFCl52BOr8Lap74GqKpdD/O53v1ttXAeWW0R04/EYV1dXFdF9/etfx2effYavfe1r\n+Oyzz/D27dtqrYf/v70ri5HsOstfTberu5ZeZjxSArYlR2GRE0VhUWJhYmzjEDmI5cEPIQ8oEIkX\nkEBCCiHw4DdAIARPeUEkQggsQQgoERJKIhhrFJE4AWwCSUgiPFES4nGv1VVd3VPunuKh+rv93a//\nc6u6p7pvt+d+0lVV3brLWb/zb+ccXVvxNAetUonOPTlOXJGNx/9LhRFEzzxLpOwdbFR6KCl45y/K\neyoyXUfW/f39cB2+lCfzJGUW5TXlNdct9Sgd+FaPSnQ8dDMXdhI9lNQrTAatKyWsnZ0d9Pv9bIew\n9fV1bGxsYH19Hf1+P9sKkQdJkcuaz83NodFoZMS2uLiI5eVlXLlyBZcvX8by8jJarVZ2nIVDsTSi\ni4gsRXJFIQjAeLIrC6mO71KbTl6O9nnls9xWGdnnqB77yg+cdeKBx9OKMYyk1WhStnamW7du5e7R\nOta0uANqdnYWjUYjOxibOMlUoAojuBZx69atzA7X6/WyY2trC+vr69nR7XYzAuz3+7kVtJXo5ufn\nsbCwkB333nsver1etlfsYDDIzBPEaQ5S50J1jdQz79hOdgQ7a6S6njT+bBpI2aVo21AjrpOCftdn\neXmkXPja4Wu1Wk4d5POA9CIDeu+d5FvtPTx0c2QSndsqU+YIJT1fkp6e2TI8ehcRrCcdeHX7w62t\nLXQ6HXQ6Hayvr2NtbQ3r6+vY2tpCr9dDt9vF9vZ2ri2T6Lgyyvb2Nra2ttBqtbCzs3PE4USb63A4\nxNzcXM5+N22UTnSEq6PeASJPnkooWkAkgeOu9jHt/EQkp5KWkjTnsvqcX+20RcSvxEJyYwNmnqMZ\nDz7pPpLyJvVqan5UnR4OhzkpVsOEeA/P6ftS0jyft7u7m53jJ1V137O3wlHogKSDTNQedP41Vybe\n39/P5srevn07F1Kii8MOBgPs7Oyg2+1m9lV9N+122iamPVCdO6JjY40M2nqousKRXJ9Vqx0uPa4d\n6CxtOKoauKSmDYtErFKX/8f8pYz+auTn83WFVnVU8NksC7XJ8ZzPhRzX6Px5Gtai6pHumcG8Oanz\nOVrXarvTvKlUx+upOp23IPLzBLcda3v0UCUN8VGi80HSdw5jvQ4GA/T7fWxtbWXtz4UAANlzT0ML\nOxdE54XuUlDKzjUcDpMbd6h3VjvvWTkniiQ69ywynSrh8Ijy504NhgGomqoES8LkLI0iKYrfVW2c\nhORUQnWboRKdzvWNJEclY5ZXrXa4CornTYNYVfLjKroV0ojsvpE050uH1ev1bIaPkptqJbVaLTNP\nvPrqq+j3+7l61UG6VquhXq+j3W7n+sU0UaozIuqUasNyQ71LeUp0atRkp1OijK47bTtOJEWx40c2\nR1X5fNECplsb5quvvprzVqnUw/c5Ug6daACYpGxSTqNILXKboJKsdirmJxVY7IMfbUQ6OFDiKMtG\nexHgM488REmXDKOUDIwWvKQTSKU8leKHwyH6/X7mweVzB4MBer1ebpOndrudeeH39vYyk8M0hZLS\nic7VrmjDZ99Ji/erHc6J7tKlS7nOQSMpO6FKHqcBtW0AyAUKu62Nadb7VGXloXnzWDMAR4iFUE/r\nJOl250aqjDz9KYmc6eSIr4g6l6ryrHNvIzz4TkbmM7RGJeeztM1eFGj7pH3Ng82V5NgPdf26Wq2W\n9SvOeGB72Nvby0JU6J2l4LKzs4P5+Xns7Oyg0WjkJD+qse5Uu1OUTnRswBpbRc8c3dbqtQPSc2EV\nJDqNVVM7EJ9zWqqsqmbaACKpR8+rXS4K9FUvI/POZ3iYhqvEx8nnOINwkfMoMjHQhhORsA5WujeB\nkifVcy6OwPeq15plrt7s0x7QLiq0fVJj8CBstcvNz89naqbW0/z8fDa9C0BOytapY91uF91uN+vb\nJDqGGbFuORtIBZdpoFQbnasfbMwejKhBtEB+jwX3JhKUCDy0wq9129S04EZdRaTmKdFp3iIJiwSg\nXk6qEdF+G5N09sjbNqnqGqmrLoFTenDVVSd7+4KjmgdX0zlIsjxUcnT19jQHtIsKbZ+RrU3tceo8\n43Q8Ho1GI/ukk4gedsbizc/Po1ar4datWxgOh7mZF05ye3t72eD9mlBd3SajB2NrKN3opGEgHSjr\nnY4Eqh2d17ES9XlOKNPIY+p8ZEtziVPhKmLkwHGj8rh08D+XkE9qu4zsjnyHx/bxUEki2h/D8877\nWLckTy1DLZNIiqwwAsuMgxBJTQep4XCY1U29Xs9N8+IMCB7scwwMB5DVE+1xar6gRO6zgXTvmGmh\nVImOHYuNm2xOsZX/qf0gUj/VyK8kwEJW8D8aVjXuy6WvOyE7lSL8OWpnm4SQtKOqlBORXJGxP5VO\nVY+nQXKaDiVR7VjRQOdEFxEc79UgVapUmg4dBCKp+m6Htk/WA23YkcOvXq9nZgMNMWk2mxnRAcgR\nHYlsMBjk5rWyLtw+z+8UdKaJcyHRcdRQ4qP3zO0sZPpxpMDRRSUgnlcvoHpipx1+Mu45TtwpW57m\nLaUiqpqo755EFVXnxyTpTiGS5jSf/O4OCCU8DVPwZ1AyrNVqmS3XCdHthep5rZCHthGVjNmu1LFA\nB4+aGOr1erYMU6PRyCRtEh1DinZ2dnJr0qlEFy3y4NrJNFB6HJ2O1FGwoXsXdVqU3l+r1XJxZRzl\nVa2lTUefp0TrkfSnbdOJyIh503d7hUeE51KU3qdlHHlzj2OTi/Kgz1cJKvKauQTps1d0wPM8a13p\n9S7J+wCgZFghD20j7IP0tJLsdCCNnBQ8VBoEkFt0wXfjI8nRG+sD9bRROtERbLS6mS6ArONowKs2\nbO0YLDh+B/JkRVHaja+aBlezTgNFNjn/Tz95PnIAuMrq90dSk5bdcfMb5UEldL+O0DKOiC6VHldd\n9V6fHaEeX3VqVMjDtQYKGEp2ahK6fft26KyYm5vL1E0KGD4VzyU5V1tVUzuNvncuiI6ZYyd0m5Qa\nQ3WkZqXwf3VYRBKR3uvzIbXzpVSwaec59Tuy6fEzIrjIJqb3KTEoyWne7zQPautUSc4ls4isJpmM\n79Ki1pdKCao+uQe6wiG0zaSkOtanXu9Ep/2I5U/tSq9VaY6E6BId03AaZHcuiA5AKA0Ah3YX2us8\nOl5nDvB+JwA34FON5WhCycA7xDTtdRGKnq1SnCIytqdUV353aU4lqDvNn6rf7sghOfn1LgGqCj3J\nu/Q5moeisqkQw9uMalYUPhRKXpTS3KzE63xgZZtmH2bd6XJkUQTENFAK0Y3zzvG3Xk9SYwGR7HTG\ngBu8tZGrpKEV6SqOjkplwSU45iOaCuezRniPSjvq2leCi/I4jhSie6K64+/oeU52k6Ql5Wgokr7L\nrMOLBh2wVODwMvTZKzrQaL26LRY4NB35ijocDKm1Rc6oO0UpRBeNsi4R8JxLNt5R/RrvdG6jUolG\n1Rwli1QHPQtEJBcZ2FMLBXiALhuPEt04IneVBhgfl6cSAYk2IqFIKpskLSmSS4XV+FFhPLT/qYDg\n14wbMJ3oPMh7MBjk/tcZFx5rNy2ULtEBR2N6eC46VMXhfd7QVQ3WDkgDqxqoI6mIdsGyENninOyi\nSe68Vm0tNBTzYL7Gdf5xJOWIBqVJMWlaPG6wSC11kqvIrhjRwBCZHpy8IpLT56UkOv2f72CMHonu\nwquuUQYoAWgQrRakj+yRLS6SziKbX8r2xU6ko0wU5sDPaUoN+lz3IEbSXMo2xwboU6qOM981KqOU\nKuOYRuM8jlTrda2dy2dYVDiKSO0E8usR+vVF5gZ9Fm3hvkgHg4L5Hg8WP42oh9KITqUx4JDAUiSn\nTgRdvin1H8NInAh1MUpXfZVceN5tdh6IOqkhvQhOcqlO7cSnRO5Sr3vETjtkZlpwb6A7kliH3g54\nHfOZmk5WYQTtWynJLDLhTEJw/J8kxxVMdCmm4XCYG4x9YYFpt9Vxtf8RADcBfEnO/RGArwB4EcDH\nASzJfx8C8HUAXwXwrtRDI2Olx1A5yUUN3NerU5JTItTn8H+VAIgiiWGSd9wJXFX192kgtDseVJKJ\nJmZPy8N62piE5KJycdJXz+FFIvmzRkq9P07/TD0TQDZDguvScatEEt2lS5dyu7jdaWxnEcYR3UcB\nPGXnPgXgzQDeCuBrGJEbALwJwHsOPp8C8OHU8yexnUQqqRNMRAS6np0u0RSRlHak1PN84U+9ZhpT\nVSL7U/S+FMkBebuJxzhNM5TkrBCZJVg+0WKsOtiwg7p3sEKMlOklRXbjzB86UO3v72f7RXA5JrZj\nANkCAamogGlinOp6HcCDdu7T8v3zAJ4++P7zAJ4F8CqAGwC+AeDtAD4XPdhtUVHIQGR7S6ktLFRu\npaerlrDh63eGlPA7lzhylc9HGCeXk3YibxDaYZU8VZpx5wpwOFBEE+MvArnpIKGk5t+SDOgDAAAT\n+UlEQVRT0rRG1LspYZo21Arj4dqX9kcnudnZ0ZaV7XYbrVYrcxKeVj3dqY3u/RiRGwB8L/Kk9m0A\n96VudHUkIj6/Rkd0JzudQKxEF82j5DkGCzM+j0HJlIoYpxdNT7qTwMaIxJgX90TrNTzvMwOi0JnT\ncNGfForUVZV0U+qqSh8axFoR3NlC26oupqsqq8bMzc/Po91uo9lsnmui+10AAwB/XXBNUqdTm1mK\n6NS7qEQXje7cyZ0jCJfwSU0a108u9cPnkuRo3/LneGc6CSInCQMp3SXP/+lo4DXRbAf1sF4EiY6I\n1NVIpddVLtzT6nVdSXNnC7eD+2rhKtFxdWKV6E7DCUGclOh+CcBPA3hSzn0HwAPy+/6Dc0fwzDPP\nZI33kUcewcMPP1wYK+ZE53Y0uq11dWIuKeNHRHZUW0lumhaSIO+nWxxIT9NSpMJbIpWMHmElOnZk\njnaRWu1S51l28FQYyrjreK7ocGmdXrtofiTrJlopo0IaJ7EvqwSu7Zn9b2dnB51OB2tra1hbW8Pq\n6io6nQ52d3cxHB4uGjA/P49ms4n5+flsUc/j4Nq1a7h27dpE107SEh4E8EkAbzn4/RSAPwbwGIBV\nue5NGEl3b8dIZf0MgO/DUaluyJFZG69Lah4Qq4fuEKaNX0lPycmJwM/rTkY0jOpCgUooep1u1qur\noGQZNS8iP53UNe1AfpaINiZ1NvhiBKrK+vnjwG2ECn9mKn/HeWZkkwTyq89SWqdZYnd3NyfJUw3i\nireLi4tYXFzE0tJSrpwq4stD6++490Se8H6/n+0Nsb6+jhs3buCll17CjRs3sLGxkS2tPhwO8cY3\nvjE7HnjgAdx33324//77cfXq1RPn56B+w0oeJ9E9ixGhXQXwLQDPYORlrePQKfGvAH4VwJcB/M3B\n597BubAUVZIh4fnie/oZeSD1Xh4ec1Ykxam6pzsbeedVNVBnFShhFkkr+j1qIFFcoHqkPc3qoeKz\nee1ZduSi/EXXkdij/6J7SXQcyHwfEbVbajiJDlAV0ogGqUnvc7sp+22328XGxgY2NjawsrKCl19+\nGTdv3sTNmzexvb2d2aC5tBP3muDsndO0KY8juvcG5z5ScP3vHRyFUF1ejZUqkanUplKdb3dHktvd\n3c3Z+AAcITc33JPsItsQn6P2HrXNUaUdF17iqmqkmrvqqp5VVa+9I6fI5TRIL/XMlG3N/1di8rKJ\nOt3t27czG4/aedg+9B06CNCofVEcMWVj0vbL765d6c59m5ubWF1dxerqakZwKysrWF1dza0SpG1Z\np3xp31UUxexNitLmuu7t7WF3dxfb29vo9XrhUsrufY1GkpQzQ981HOa3B/SKpT1PQ0j4Dl9mJkWK\nKllF+dV0R3FyzDPtTUrGGvKiKq1jXGziSVD0zGhg0PL3evP4v5StR1V6t8u5A6JWqx2Jrnf1etpl\n8lqBlrmf1zqL4lSjHfs2Nzexvr6O9fV1rK6uYmVlJbPNAYcrn3Cz6/39/cwkwdkTvhF2FMxc9H8K\npa1eQqLr9XrodDo5yQY4ujKwq31RmIFWkHceJTlXa9wmxhAT7l7Eg9JfKt5Pn5fKs3oMI7JjR1XD\nehQXF0lx/JxWp57kmV7uUYiQDlzjJD8tD50jqfUN5IOkfRFIDkp8R0VyR+F1pvB2SWlaoxqcoPr9\nPjqdDjY3NzP1dWNjIyM6XYKpXq+jVqtlz46ILgpSTgUwTyK9ly7R9Xo9bG1t5Rq5exJTRBdtpuEd\niiSnkl2RqseOSVGbjYDqUZFEV0QGJLrBYJCzIzrhsaOqGsaRkM/Szygf05boomemVP2i0KAiqQ44\n3BpPpQceWmYsayc5l+gmrZ+7FVpnCu1fg8EgtxF1r9c78p2fW1tb6HQ66HQ62Nrayv7b2dnJNs9x\niU5jX7e3tzOBwsOE+Kl9AZh8Q6tSiM6jpvv9fq4DaKxYSkUE8ks/U62cnT1cij0qoKKloKPgYC7t\nxKkqGtQYkY124Mgux0akBMC0+qqsqUnORTaV4yJSH1PSot/n+UodTnbRO4A80fnhKms02nsZ0RRQ\nEdwhUuYCrSuSD/snvaU8nPCU+Oh13d7ezqRAesZV9d3d3UW/3882uKbtuVar5SIftJ7ZNyKeGIdS\niU7ZXEHJS+cpeufg/xy5yfKuy6sdx0NDVCVNdR4lHRKkOgL4rqiTqY0jJcWpWn3p0qWcuuqBr9Mk\nOEVK0mJ5RCOmSqlROFDRnNSUVFoUP+l2uSh0JpIsIwn+bkU0CKsGQzOBEla32z1Ccjw4WV9VTz3v\ngxaJTwmOHleNYOCuYo1GIxQ8WLckxknqtxSi86jpnZ2dI6OySmYAQmLRwFqeizysJBDdmo2FPDc3\nlxsR3LjptoBIQnAiikbK6KC7XSvRJU1Ny2kgsq953lL3ReE+HsztBOh1qlBiU1KMJLkoni+lRldk\nd4hUOXHyfb/fzzkVNjc3jxAdSY1zWNUhoXNb9T3s8zs7O7m4R11dh/XUarUyW6xGOnikgWp041AK\n0XlMHLco5JEyXHsGi8JGWEC8hiMEDyU8t+ekEKVDiTFSB6IQGbc3AYf7WKRWcDhpRx13n47qbqsp\nGi1VooviH+mZc2lvXPmmDOQewxiVj0rLbsD2gfJuB+uc0hzt5Qz2feWVV7CysoK1tbWcqkoJjoc6\njDTUhFMWWU/aVlR1VbLjdbTJ3r59O9O42Ke1DjU0bBxKITrdD1LJhtAQDwBHGr9KbioFpRbxowFU\npTid2eCSgSMiweFwGK5v70TnS0cNBoMciUeSIhuFkmnUOb3M+P7oM8qTlq0+RwkiRbYqrbqzwMMR\nXAVNYZLBJkqLDpoc4VOmhVTIwt0EEgolrO3t7cxjyqlbKysrWFlZwfr6+hFJjUfKjsp36Nxr7XdO\nbEwHiUvTFQklXr/nnuiU7NRuoCOA2xOAw0IEkCsEn5rltrVokT86F1KF5ZIloRXpxlBX65TkBoPB\nEQJTYuE7SRBOhE58USclcU1C3J43NwukyC4K81GpVfOu5F+ESUkuSgvfx7YT5SkicNb93UJ2zCc9\nnox6WF9fx9raWvbJwN/Nzc0jwfwa1O9rNnqUAgUPkpxKcMBhW6dDknb7RqOBbreLdruNhYUFLCws\n5MKJipxaEUojOmaeJKUqnUp02jmYKaoiHCmazSba7XY2avCIvKtKbmrjS8HtV1GnibzDLqrrQZKM\n3Of6Toba+BHFAbJ8UgSmZej/pUguJUm6aq5SW0R0k+IkJOcqmOfR8xN56u42kgOQkcr29nZmk2OQ\n7+rqajYhnzGuUf0WecaBQ3OMalAUNFSlpZbDFYcYTzc3N5fNeGLIFWckHdeUUwrRMdHNZhMLCwtH\njNrD4TDX0VySUo9oq9XKLfWSkug8Nq9IJVTodW7jGXd/EfFQWtNO6sQbSSGRZ1jfp5+TpMvf4+8u\nknbcLhYNCj4ITApPk6pB6oknaOZQEva4PnYwbQs68L3WoXVAOypVxG63i06ng42NDWxtbaHX62X7\nPLiH1s0QWp6ss5mZGTQaDTSbzSPqp2pfXJ7J7ey+BUCqL51riY5e0GazieXlZczMzOQkHiAmGBKg\nzplrNpvZEc0kcIdFkaSSQmTQ1s9Jn6XvZaNhnrRs9NoUyRe9M6Xmero1b5pHlRg9vMQ/XYplB1A1\nw2PoxpWRpsuJSedFen2kZpzs7e3l7KnqiWeU/iSxWBcdWkc6hYvhHt1uNyM5hoeo9BQRi7ZXXXRi\nbm4O7XY7O0h2jUajMCjYByFer6vQAPk5uufWRsdG1mq1AIxUWTVyuoqojdptbupJ1RE/ZddyqWUc\n3JYTxZtNQnb+HD10xzInYn2nPyulerLRsLGMqwu/tqi8PK/aeSKHEb+ThCYtI81DtHKyS4/6fWZm\nJiO5mZmZLGDVJXuaSSYJT7jo0MFI4+boZKBUR6JjoG80OGl7UwmbZqRms4lWq4WlpaXs0GiHer1+\nhDSjEC5GS/iinEXEm0JpRHfPPfeg0WhkpKfz6ZzolFw4avDQmDgVnaepikSdfdL7UuqgVxIboRO0\n2kcUTpZOcj51LAWNNmcaSQAuxWpZpKRFJzrm150a/jx9Jq+NZojotEAPvPawpP39/YzQSHwq0dVq\ntWyJrrsBRWTHo9/vZ6sJRbYwlX5d6m61WtlagIuLi7h8+TKWl5dx+fLlnPpKKd+XJvP2MRwOcw4M\nfRffPylKG8rYEbhir47g/B0Rg8abRdv5naatRTuip9EPtznQOeFQcnAV0aWl6F6FGnjVnlIEEmFE\nakVQ1bRer+ccAepQAvIxcEWqOOvR7WcupQPISXK075LUomenENkrX6vwtqlSlJqDOLiQiCItg4MX\nHYpzc3PZYqdLS0tYXl7OLYCqdjk+250YkRnEw9DoaCRhTrokV+lEx+8s7Hq9np1ToiN0RFapxUMG\nTjPdTEeqE2nH9fmz0fUkB1cTVG0tkjpU8mVZ6vciqD3NZ6OMKweNUdR8uP1SPZ5OdJENLtqrQ9ML\n5KetUVK7dOlSNphMQlwpo/Zrley0PbiqqOEgRXO/tV/Ozs5moR/tdjuT3q5cuYKlpaXsPCMi9Lm+\nsK7bVnnoXHMl1WjueRFKIzotcG3gqvpEROLqzUkdDCeFkm+RKuYdOEVywNFYPR/ZUpJWSpLUMpwE\nHMmjTl9UDh7TxMaasr+knsHGHy0uyjJ06VbLxre21OsmaReT2lkvMrTNFkl07Ie6UGY0YFMoIbFd\nuXIF9957b/Z5+fLlnKPQNS9fpcanDfI7B9KI7NwLX4RSrbDaCKNOOk5imlQ9Oa20F51zKUkbmZKB\nG/IjR0eKeJz0o7JwiZjPJaJQnCgMx6Gqq6dFA4iddPV5kZrqC2gyHRHRqWqsqpfa7HwQ0Gt1ddu7\nAVqPjFhYWlrKhZDMzMxkDoler4dbt26FgxE3tyHJXblyBcvLy1heXsbS0hIWFxczL+v8/HwudrVW\nq2UbUkWLQujhA6BKcuMECEVpRPfcc8/hsccey1QnVzWKyMsJ8k5x7do1PP7443f8HCAv0emkfDpg\n5ubmMqJzg/r+/j6ef/55vO1tb5voPSmy9+/undVPbUCutly/fh1PPPFEktTVvsf81uv1I0TnNjCV\nJlRKUFOEq0zPPfccHn300TCsQElSZ9NE6n60VJc6sSbBNNvLneAk6WAe5+bmsLCwkAXjcmOhhYWF\nbCpYp9NBv9/PlZmuLMLwsOXlZXzzm9/EQw89hFarhVarlYV7kZi8jFVd1rhHV111EPQVhCYlOaBE\noosqiSrKuMRH6uy003ISuHrAc6riqaOAv3UEe/HFF7O0jJNYVVLRc35vSkpU547PDybRPfnkk0fe\nC+RND2yAvsySetWYV1fnXYVyryg/P/vZz+KJJ57I8qDQd7u0rJ+8NrJFHWfgvKhEp+2zXq+j3W5j\ndnY2C/1ot9tYWlrKVi5pNpvo9Xq58iKR6dSsxcVFXL9+HU8//XS4Qx7nk2vZMoIiZapRidzbhdu7\nz7VER0xTMjsPUImOleQd0yP2dd4gR83IFhm9q8j+xf+KbH+uOkaG51Q+6STgc5S4lfR0xFbJjQQT\nqd5R/lOhMvpuJzd9NyVY99YfRzK46GA+Z2dn0Ww2M0lO91lVI78vlkEypEeVzoZWq4WrV6+GziUd\nSCaB1t+0zFN3h3GiQoUKdzXKGMauYbRXbIUKFSpME88BeLzsRFSoUKFChQoVKlSoUKFChQoXCE8B\n+CqArwP44Bm/+yMAbgL4kpy7AuDTAL4G4FMAls8gHQ8A+BcA/w3gvwD8eolpmQfweQAvAPgygN8v\nMS3EDID/APDJktNyA8B/HqTl+ZLTsgzgYwC+glE9PVxSWn4Qo/Lg0cGo/ZbZXs4dZgB8A8CDAO7B\nqHM9dIbvfxTADyNPdH8I4LcOvn8QwB+cQTpeD+CHDr63AfwPRuVQRloAoHnwOQvgcwDeUWJaAOA3\nAfwVgE8c/C4rLS9h1IEVZaXlLwC8/+D7LIClEtNCXALwXYwG7rLTcq7wYwD+SX7/9sFxlngQeaL7\nKoDXHXx//cHvs8Y/AHjnOUhLE8AXALy5xLTcD+AzAJ7AoURXVlpeAnCvnSsjLUsA/jc4X3Z7eReA\n6+ckLUmUEUd3H4Bvye9vH5wrE6/DSJ3FwefrCq49DTyIkZT5+RLTcgkj6fomDlXqstLyJwA+AEDn\ncJWVliFGpPtFAL9SYlreAGAFwEcB/DuAPwPQKiktil8A8OzB97LTkkQZRHfedxIe4mzT2AbwdwB+\nA0C3xLTcxkiVvh/AT2AkTZWRlp8B8ApGtp9UnOdZlsuPYzQIvRvAr2Fk+igjLbMAfgTAhw8+t3FU\nEzrrtlsH8LMA/jb476zTUogyiO47GOnzxAMYSXVl4iZGojYAfA9GHe0scA9GJPeXGKmuZaaF6AD4\nRwA/WlJaHgHwcxipjM8C+EmMyqescvnuwecKgL8H8PaS0vLtg+MLB78/hhHhvVxCWoh3A/g3jMoG\nKL/tJlEG0X0RwPdjpK7VAbwHhwbnsvAJAO87+P4+HJLOaaIG4M8x8p79aclpuYpDD1kDwE9hJFGV\nkZbfwWjwewNGatE/A/jFktLSBLBw8L2FkT3qSyWl5WWMTD4/cPD7nRiZFz5ZQlqI9+JQbQXKKZdz\njXdj5GX8BoAPnfG7nwXwfwAGGDWcX8bIq/YZnK1b/B0YqYsv4NBN/1RJaXkLRnafFzAKpfjAwfky\n0qJ4DIeDYBlpeQNGZfICRiFAbKtllctbMZLoXgTwcYwcFGWlpQVgFYcDAUpMS4UKFSpUqFChQoUK\nFSpUqFChQoUKFSpUqFChQoUKFSpUqFChQoUKFSpUqFChQoUKFe5m/D+f1L+qCkr+EAAAAABJRU5E\nrkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x113ff8790>"
       ]
      }
     ],
     "prompt_number": 42
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We'll notice some of the clusters are empty. This is because our sampler assigns zero mass to some of the cluster; we can see this by plotting a histogram of, e.g. the last assignment vector generated by our Gibbs sampler."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.hist(usps_chain[-1], bins=40)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 52,
       "text": [
        "(array([  3.,  14.,   0.,   4.,  33.,   0.,  11.,   6.,   6.,  23.,  14.,\n",
        "          0.,  22.,  13.,   2.,   0.,   6.,   2.,   4.,  19.,   0.,   4.,\n",
        "          2.,   8.,  14.,   0.,  29.,  64.,   0.,   4.,   8.,  11.,  23.,\n",
        "          8.,  16.,  30.,  11.,  10.,   4.,  47.]),\n",
        " array([  0.   ,   0.975,   1.95 ,   2.925,   3.9  ,   4.875,   5.85 ,\n",
        "          6.825,   7.8  ,   8.775,   9.75 ,  10.725,  11.7  ,  12.675,\n",
        "         13.65 ,  14.625,  15.6  ,  16.575,  17.55 ,  18.525,  19.5  ,\n",
        "         20.475,  21.45 ,  22.425,  23.4  ,  24.375,  25.35 ,  26.325,\n",
        "         27.3  ,  28.275,  29.25 ,  30.225,  31.2  ,  32.175,  33.15 ,\n",
        "         34.125,  35.1  ,  36.075,  37.05 ,  38.025,  39.   ]),\n",
        " <a list of 40 Patch objects>)"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAHfCAYAAABwLo3rAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFzlJREFUeJzt3X2MbHddx/HPaZcGaO/QbjC3N1BsLZYCQcEHMIIyEEoo\nwVpjgpBIblD5RwQiESgkpvuXNCQGTEz8Qx5yRYI2PDTFh9hSOhGjFJGWp1IrDY2A7QXTC3cLJoId\n/zin3b17d3tnvjuPO69Xsrkzs3N2fvc3Z+e+7zln5iQAAAAAAAAAAAAAAAAAAAAAC+FpSW7f9vW9\nJG9Msp7k5iR3J7kpyfnzGiAAwKI7K8l9SS5K8q4kb+1uf1uS6+Y1KACARffSJJ/uLt+V5HB3+cLu\nOgAAu3h/kt/tLp/Ydnuz4zoAwIHWjHHfc5J8K8kzknwnbTRdsO37D6Q9TuoRl1566fCee+7Z7xgB\nAGbhniRPHfXOZ43xg69M8m9pAypJjqfdjZckR5J8+7SR3HNPhsOhrx1f11577dzHsGhf5sS8mBfz\nYk7My7y/klw6RheNFVGvTvLhbddvTHK0u3w0yQ3jPDAAwDIbNaLOTfKSJB/bdtt1Sa5I+xEHL453\n5wEAK2RtxPt9P8kTd9z2QNqwYkz9fn/eQ1g45mR35mV35mV35uV05mR35mUyxjmwvGLY7WMEAFho\nTdMkY7TROMdEAQDQEVEAAAUiCgCgQEQBABSIKACAAhEFAFAgogAACkQUAECBiAIAKBBRAAAFIgoA\noEBEAQAUiCgAgAIRBQBQIKIAAApEFABAgYgCACgQUQAABSIKAKBARAEAFIgoAIACEQUAUCCiAAAK\nRBQAQIGIAgAoEFEAAAUiCgCgQEQBABSIKACAAhEFAFAgogAACkQUAECBiAIAKBBRAJyi11tP0zRj\nf/V66/MeOsxUM+WfPxwOh1N+CAAmqWmaJJXX7iZe81lm7bo/ehvZEgUAUCCiAAAKRBQAQIGIAgAo\nEFEAAAUiCgCgQEQBABSIKACAAhEFAFAgogAACkQUAECBiAIAKBBRAAAFIgoAoEBEAQAUiCgAgAIR\nBQBQIKIAAApEFABAwagRdX6SjyT5apI7kzwvyXqSm5PcneSm7j4AACth1Ij6kyR/l+TpSX4qyV1J\nrkkbUZcluaW7DgCwEpoR7vOEJLcn+Ykdt9+V5IVJjie5MMkgyeU77jMcDof7HCIAs9Q0TZLKa3cT\nr/kss3bdH6mNkoy2JeqSJN9J8oEkn0/y50nOTXI4bUCl+/PwOAMFAFhmayPe52eS/F6Sf03ynpy+\n626YPf7bsrGx8cjlfr+ffr9fGCYAwGQNBoMMBoPy8qNssrowyb+k3SKVJC9I8va0u/delOT+JEeS\n3Bq78wCWnt15rKpp7M67P8k30h5AniQvSfKVJJ9IcrS77WiSG0YeJQDAkhu1tn46yXuTnJPkniSv\nTXJ2kuuTPCXJvUlemeS7O5azJQpgydgSxaoad0vUyHcsElEAS0ZEsaqmsTsPAIAdRBQAQIGIAgAo\nEFEAAAUiCgCgQEQBABSIKACAAhEFAFAgogAACkQUAECBiAIAKBBRAAAFIgoAoEBEAQAUiCgAgAIR\nBQBQIKIAAApEFABAgYgCACgQUQAABSIKAKBARAEAFIgoAIACEQUAUCCiAAAKRBQAQIGIAgAoEFEA\nAAUiCgCgQEQBABSIKACAAhEFAFAgogAACkQUAECBiAIAKBBRAAAFIgoAoEBEAQAUiCgAgAIRBQBQ\nIKIAAApEFABAgYgCACgQUQAABSIKAKBARAEAFIgoAIACEQUAUCCiAAAKRBQAQIGIAgAoEFEAAAUi\nCgCgQEQBABSIKACAgrUR73dvkpNJ/i/JD5M8N8l6kr9O8uPd91+Z5LsTHyEAwAIadUvUMEk/yXPS\nBlSSXJPk5iSXJbmluw4AsBLG2Z3X7Lh+VZJj3eVjSa6eyIgAAJbAOFuiPpnkc0le1912OMnx7vLx\n7joAwEoY9Zio5ye5L8mPpd2Fd9eO7w+7r9NsbGw8crnf76ff7487RgCAiRsMBhkMBuXld+6iG8W1\nSR5Mu0Wqn+T+JEeS3Jrk8h33HQ6Hu7YVAAuqaZrs8f/iMy0Zr/kss3bdH72NRtmd9/gkh7rL5yZ5\naZIvJbkxydHu9qNJbhh5lAAAS26U2rokyce7y2tJPpTknWk/4uD6JE/J3h9xYEsUwJKxJYpVNe6W\nqMruvHGIKIAlI6JYVdPYnQcAwA4iCgCgQEQBABSIKACAAhEFAFAgogAACkQUAECBiAIAKBBRAAAF\nIgoAoEBEAQAUiCgAgAIRBQBQIKIAAApEFACwkHq99TRNM/ZXr7c+k/E1U/75w+FwOOWHAGCSmqZJ\nUnntbuI1n0ma9brYPt7obWRLFABAgYgCACgQUQAABSIKAKBARAEAFIgoAIACEQUAUCCiAAAKRBQA\nQIGIAgAoEFEAAAUiCgCgQEQBABSIKACAAhEFAFAgogAACkQUAECBiAIAKBBRAAAFIgoAoEBEAQAU\niCgAgAIRBQBQIKIAAApEFABAgYgCACgQUQAABSIKAKBARAEAFIgoAIACEQUAUCCiAAAKRBQAQIGI\nAgAoEFEAAAUiCgCgQEQBABSIKACAAhEFAFAwakSdneT2JJ/orq8nuTnJ3UluSnL+5IcGALC4Ro2o\nNyW5M8mwu35N2oi6LMkt3XUAgJUxSkQ9OcnLk7w3SdPddlWSY93lY0munvzQAAAW1ygR9e4kb0ny\n0LbbDic53l0+3l0HAFgZa2f4/iuSfDvt8VD9Pe4zzNZuvtNsbGw8crnf76ff3+vHAADMzmAwyGAw\nKC/fnOH7f5TkNUl+lOSxSXpJPpbk59NG1f1JjiS5Ncnluyw/HA737CsAFlDTNHmU/xs/2pLxms8k\nzXpdbB/vjG30iDPtzntHkouSXJLkVUk+lTaqbkxytLvP0SQ3jDtQAIBlNu7nRD2cddcluSLtRxy8\nuLsOALAyRt5kVWR3HsCSsTuPRbHsu/MAANiFiAIAKBBRAAAFIgoAoEBEAQAUiCgAgAIRBQBQIKIA\nAApEFABAgYgCACgQUQAABSIKAKBARAEAFIgoAIACEQUAUCCiAAAKRBQAQIGIAgAoEFEAAAUiCgCg\nQEQBABSIKACAAhEFAFAgogAACkQUAECBiAIAKBBRAAAFIgoAoEBEAQAUiCgAgAIRBQBQIKIAAApE\nFABAgYgCACgQUQAABSIKAKBARAEAFIgoAIACEQUAUCCiCnq99TRNM/ZXr7c+76EDABPSTPnnD4fD\n4ZQfYvaapklS+Xs1OYjzARwsXuNYFLNeF9vHG72NbIkCACgQUQAABSIKAKBARAEAFIgoAIACEQUA\nUCCiAAAKRBQAQIGIAgAoEFEAAAUiCgCgQEQBABSIKACAgjNF1GOT3JbkjiR3Jnlnd/t6kpuT3J3k\npiTnT2uAAACLqBnhPo9P8oMka0n+KckfJLkqyX8neVeStyW5IMk1uyw7HA6HkxnpAmmaJknl79Xk\nIM4HcLB4jWNRzHpdbB9vpDZKMtruvB90f56T5OwkJ9JG1LHu9mNJrh59iAAAy2+UiDor7e6840lu\nTfKVJIe76+n+PDyV0QEALKi1Ee7zUJJnJ3lCkn9I8qId3x/mUba1bWxsPHK53++n3++PO0YA4FH0\neuvZ3Dwx1jKHDl2QkycfmNKIlsNgMMhgMCgvP/J+v84fJvmfJL+TpJ/k/iRH0m6hunyX+zsm6tQl\nHS8ALDyvccun9pwt/vO17MdEPTFb77x7XJIrktye5MYkR7vbjya5YaxRAgAsuTPtzjuS9sDxs7qv\nDya5JW1IXZ/kt5Pcm+SV0xsiAMDiGXd33rjszjt1yYXfdArgNW752J132pILsTsPAIBdiCgAgAIR\nBXBA9XrraZpm7C9gNI6JKnC8ALAM9vNa5TVuuTgm6rQlHRMFALCoRBQAQIGIAgAoEFEAAAUiCgCg\nQEQBABSIKACAAhEFAFAgogAACkQUAECBiAIAKBBRAAAFIgoAoEBEAQAUiCgAgAIRBQBQIKIAAApE\nFABAgYgCACgQUQAABSIKAKBARAEAFIgoAIACEQUAUCCiAAAKRBQAQIGIAgAoEFEAAAUiCgCgQEQB\nABSIKACAAhEFAFAgogAACkQUAECBiAIAKBBRAAAFIop96/XW0zTN2F+93vq8hw5wRl7j2Esz5Z8/\nHA6HU36I2WuaJknl79XEfJyy5IGcD1gU+/nd9Du9ZRle42pjXPzna9Zz3z7e6G1kSxQAQIGIAgAo\nEFEAAAUiCgCgQEQBABSIKFaCtygDMGk+4qBgGd7uOkvLMB/LMEaYNB9xMBnL8PrhIw5OW9JHHAAA\nLCoRBQBQIKIAAApEFABAgYgCACgQUQAABaNE1EVJbk3ylSRfTvLG7vb1JDcnuTvJTUnOn8YAAQAW\n0SgR9cMkv5/kmUl+Icnrkzw9yTVpI+qyJLd01wEAVsIoEXV/kju6yw8m+WqSJyW5Ksmx7vZjSa6e\n+OgAABbUuMdEXZzkOUluS3I4yfHu9uPddQCAlbA2xn3PS/LRJG9Ksrnje8Ps8bnsGxsbj1zu9/vp\n9/tjDRAAYBoGg0EGg0F5+VHPD/OYJH+T5O+TvKe77a4k/bS7+46kPfj88h3LOXfeqUsu/HmKKpZh\nPpZhjDBpzp03Gcvw+uHceactuTDnzmuSvC/JndkKqCS5McnR7vLRJDeM+qAAAMtulNp6QZJ/TPLF\nbOXg25N8Nsn1SZ6S5N4kr0zy3R3L2hJ16pILX/0VyzAfyzBGmDRboiZjGV4/bIk6bcmZbIka+Y5F\nIurUJRd+ha1YhvlYhjHCpImoyViG1w8RddqSC7M7DwCAHUQUAECBiAKmqtdbT9M0Y3/1euvzHjrA\noxrnc6IAxra5eSKVYxo2N6d9yCbA/tgSBQBQIKIAAApEFABAgYgCACgQUQAABSIKAKBARAEAFIgo\nAIACEQUAUCCiAAAKRBQAQIGIAgAoEFEAAAUiCgCgQEQBABSIKACAAhEFAFAgogBgQfR662maZuwv\n5mNt3gMAAFqbmyeSDAtLCql5sCUKAKBARAEAFIgoAIACEQUAUCCiAAAKRBQAQIGIAgAoEFEAAAUi\nCgCgQEQBABSIKACAgoWLqOrJF3u99XkPnbGteZ5ZOV7j4OCY9hkLh8PheCdSbM9GXTv54riPVbUM\nY5yl/czH+MvV5tBzNj/m/lSznI/Z/m62y3nOTllywZ+zxX++Zv360T7e6G20cFuiAACWgYgCACgQ\nUQAABSIKAKBARMGj8g5CAHa3Nu8BwGL7UcZ9Z8jm5rTf9ArAIrAlCgCgQEQBABSIKACAAhEFAFAg\nogAACkQUAECBiAIAKBBRAAAFIgoAoEBEAQAUiCgAgAIRteB6vfWxT4DrJLgAp6u+nh5c459g3b8v\np5r22jEcDsc7eWu7wo63TLdkxn2sqlmO8aDPx/jL1f5eyzDGg2oZ1uFZWpbXj4P4nC3DfMz6tWrR\nn+dZv3500TxyG42yJer9SY4n+dK229aT3Jzk7iQ3JTl/9CECACy/USLqA0letuO2a9JG1GVJbumu\nAwCsjFEi6tNJTuy47aokx7rLx5JcPclBAQAsuuqB5YfT7uJL9+fhyQwHAGA5rE3gZwzzKEd9bWxs\nPHK53++n3+9P4CEBgPlYK71r8dChC3Ly5ANTGE/dYDDIYDAoLz/qLFyc5BNJntVdvytJP8n9SY4k\nuTXJ5bss5915py7p3XmnLllYzrvzls0yrMOztCyvHwfxOVuG+ViWd+ctw3wsyrvzdnNjkqPd5aNJ\nbij+HACApTRKRH04yT8neVqSbyR5bZLrklyR9iMOXtxdBwBYGT5ss/JIS7I5fhnmw+68g28Z1uFZ\nWpbXj4P4nC3DfNidt2OJA7o7DwBgpYkoAIACEQUAUCCiAAAKRBQAQIGIAgAoEFEAAAUiCgCgYBIn\nIAZYSb3eejY3T8x7GCtpOea+dqJeloeIAihq/xGvfEo0+7Ucc/+jLP4Y2Q+78wAACkQUAECBiAIA\nKBBRAAAFIgoAoEBEAQAUiCgAgAIRBQBQIKIAAApEFABAgYgCAChw7rwDq3riy8ck+eGkB8OUVE/C\neujQBTl58oEpjAhgdYioA6ty4sukPfmlE2Yui9pJWJPNTc8ZwH7ZnQcAUCCiAAAKRBQAQIGIAgAo\nEFEAAAUiCgCgQEQBABSIKACAAhEFAFAgogAAClY+onq99TRNM9YXrKpF/32pjK9pmvR66zMdJ1uq\nzxnLaO3APc/THuFwOBzvvF7tpNXO+TbuY9Ufb3Zj3M98zPbcebN7rNk9z+3jzWqMFX5fdixhPk5d\nYg6vH4s9Rq+n83us6nLz+J0evY1WfksUAECFiAIAKBBRAAAFIgoAoGBt2g9w0UXPnPZDAADM3NQj\n6pvfvH6Me/9nkpdPaygAABMz9YhKxtkS9dipjQIAYJIcEwUAUCCiAAAKRBQAQIGIAgAomMGB5bBq\n1konzjx06IKcPPnAFMazm9oYDy7zAYxPRMHE/SiVE2Zubs7yH/HaGKd/zvJ5MR/A+OzOAwAoEFEA\nAAUiCgCgQEQBABSIKACAAhEFAFAgogAACvYbUS9LcleS/0jytv0PZzUMBoN5D2EBDeY9gIVkXdmd\nedmdednNYN4DWFCDeQ/gQNhPRJ2d5E/ThtQzkrw6ydMnMaiDzgvdbgbzHsBCsq7szrzszrzsZjDv\nASyowbwHcCDsJ6Kem+RrSe5N8sMkf5XkVycwJgCAhbef0748Kck3tl3/ZpLn7bxTr/crI//Ahx76\nfh58cB8jAgCYkf2c+OnX0+7Ke113/TfTRtQbtt3na0ku3cdjAADMyj1JnjrqnfezJepbSS7adv2i\ntFujtht5IAAAq2ItbbFdnOScJHfEgeUAACO5Msm/p91t9/Y5jwUAAACAVeRDOHd3b5IvJrk9yWfn\nO5S5en+S40m+tO229SQ3J7k7yU1Jzp/DuOZtt3nZSHus4e3d18tmP6y5uijJrUm+kuTLSd7Y3b7q\n68te87KR1V5fHpvktrSHl9yZ5J3d7au8vuw1JxtZ7XXlYWen/ft/ors+93Xl7LS79y5O8pg4Vmq7\nr6d9glbdLyV5Tk6NhXcleWt3+W1Jrpv1oBbAbvNybZI3z2c4C+HCJM/uLp+X9vCBp8f6ste8rPr6\nkiSP7/5cS/KZJC+I9WW3ObGutN6c5ENJbuyuj7WuTOPceT6E89Ht52MlDopPJzmx47arkhzrLh9L\ncvVMR7QYdpuXZLXXmfvT/kcsSR5M8tW0n1G36uvLXvOSrPb6kiQ/6P48J+1/6k/E+rLbnCTWlScn\neXmS92ZrLsZaV6YRUbt9COeT9rjvqhkm+WSSz2Xr87VoHU67Kyvdn4fnOJZF84YkX0jyvqzWboid\nLk67pe62WF+2uzjtvHymu77q68tZaQPzeLZ2ea76+rLbnCTWlXcneUuSh7bdNta6Mo2IGk7hZx4U\nz0/7Yndlkten3X3D6YaxHj3sz5JcknbXzX1J/ni+w5mb85J8NMmbkmzu+N4qry/nJflI2nl5MNaX\npP0H8dlptzL8cpIX7fj+Kq4vO+ekH+vKK5J8O+3xUHttkTvjujKNiBrlQzhX1X3dn99J8vG0uz5p\nHU97nEeSHEm7ctPOw8O/yO/Naq4zj0kbUB9MckN3m/Vla17+MlvzYn3Z8r0kf5vkZ2N9edjDc/Jz\nsa78Ytpdd19P8uEkL077GjPWujKNiPpckp/M1odw/ka2DthaZY9Pcqi7fG6Sl+bUA4hX3Y1JjnaX\nj2brH4VVd2Tb5V/L6q0zTdpdDXcmec+221d9fdlrXlZ9fXlitnZLPS7JFWm3NKzy+rLXnFy47T6r\nuK68I+1GnkuSvCrJp5K8JguyrvgQztNdknaf9B1p35K8yvPy4ST/leR/0x4/99q071r8ZFbzLcgP\n2zkvv5XkL9J+LMYX0v4yr9qxHC9Iuyvijpz6VuxVX192m5crY315VpLPp52XL6Y93iVZ7fVlrzlZ\n9XVluxdma2PPKq8rAAAAAAAAAAAAAAAAAAAAAAAAB8b/A8sEmpM2tg4gAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x113ee9210>"
       ]
      }
     ],
     "prompt_number": 52
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is a great segway into infinite dimensional models!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}