2.2-multinomial-variables.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from IPython.display import display, Math, Latex\n",
    "maths = lambda s: display(Math(s))\n",
    "latex = lambda s: display(Latex(s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A binary variable is good for representing something like a coin toss, but normally we want to have more states, such as the number rolled on a die. One way of representing such a variable is 1-of-K  representation, where we have a vector $\\mathbf{x}$ of size $K$ where state $k$ is represented by the $k^{th}$ element being 1 and all other elements being 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def one_of(K, k):\n",
    "    # If k is an array, have each output row\n",
    "    # represent an element of k\n",
    "    k = array(k).flatten()\n",
    "    x = zeros((k.shape[0],K))\n",
    "    x[range(len(k)),k] = 1\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 0., 0., 0., 0., 1., 0., 0., 0., 0.]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For example\n",
    "one_of(10,5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example we can create a die"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Die(object):\n",
    "    def __init__(self, mu=None, sides=6):\n",
    "        if mu is None:\n",
    "            mu = [1] * sides\n",
    "        # Normalise to be prob. dist.\n",
    "        self.mu = array(mu, float)/sum(mu)\n",
    "        self.sides = sides\n",
    "        \n",
    "    def roll(self, n):\n",
    "        # sample using uniform rand num to index y axis of cdf\n",
    "        rolls = argmax(cumsum(self.mu)[:,newaxis] > random.random(n),0)\n",
    "        return one_of(self.sides, rolls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And roll some numbers on a fair die"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 1. 0. 0. 0. 0.]\n",
      " [0. 0. 1. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 1.]\n",
      " [0. 0. 0. 1. 0. 0.]\n",
      " [0. 0. 0. 1. 0. 0.]]\n",
      "As (0-based) numbers: [1 2 5 3 3]\n"
     ]
    }
   ],
   "source": [
    "fairdie = Die()\n",
    "fairrolls = fairdie.roll(5)\n",
    "\n",
    "print(fairrolls)\n",
    "print(\"As (0-based) numbers:\", argmax(fairrolls, axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we have a loaded die, we can define a vector $\\boldsymbol{\\mu}$ of the probability for rolling each number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.28571429 0.14285714 0.14285714 0.14285714 0.14285714 0.14285714]\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "Probabilities, so normalised s.t. $\\sum_{k=1}^{K} \\mu_k = 1.0$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Say rolling a 1 is twice as likely as the other numbers\n",
    "nums = [2, 1, 1, 1, 1, 1]\n",
    "# normalise\n",
    "mu = array(nums, float)/sum(nums)\n",
    "print(mu)\n",
    "latex(r\"Probabilities, so normalised s.t. $\\sum_{k=1}^{K} \\mu_k = %s$\" % round(sum(mu)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([302., 150., 143., 127., 143., 135.]),\n",
       " array([0, 1, 2, 3, 4, 5, 6]),\n",
       " <a list of 6 Patch objects>)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAPHklEQVR4nO3df6xfdX3H8edLij+GOnS9kK4tu8x1i7jEQm4YC4lh4hTQrJiMBZZpY0jqH7hgZrKg/6jJTDCZsphtJNUy64ZiIxIaJU6GGOcfgLdYgVKZHXZybUevAxVmpqG+98c9jV/Lt73ffn/02/vx+Ui+Oed8zuec8z5p7uue+/mec5qqQpLUlhdMuwBJ0vgZ7pLUIMNdkhpkuEtSgwx3SWrQqmkXALB69eqanZ2ddhmStKLs2rXrB1U102/dKRHus7OzzM/PT7sMSVpRkvzXsdY5LCNJDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoOWDfckL07yQJJvJdmT5INd+7lJ7k/ynSSfTfLCrv1F3fK+bv3sZE9BknS0Qa7cfwq8vqpeC2wELktyEfBh4Kaq2gA8DVzb9b8WeLqqfge4qesnSTqJln1CtZb+N49nu8XTu08Brwf+vGvfDnwAuBnY1M0DfA74+ySpCf2vILM3fHESu52K/Te+edolSGrEQGPuSU5Lshs4BNwN/Cfww6p6ruuyAKzt5tcCTwB0638E/EaffW5JMp9kfnFxcbSzkCT9koHCvaoOV9VGYB1wIfDqft26aY6zrnefW6tqrqrmZmb6vvdGkjSkE7pbpqp+CHwVuAg4M8mRYZ11wIFufgFYD9Ct/3XgqXEUK0kazCB3y8wkObObfwnwBmAvcC/wp123zcCd3fzObplu/VcmNd4uSepvkFf+rgG2JzmNpV8GO6rqC0keBW5L8jfAN4FtXf9twD8n2cfSFfvVE6hbknQcg9wt8xBwfp/2x1kafz+6/f+Aq8ZSnSRpKD6hKkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNWjbck6xPcm+SvUn2JLm+a/9Aku8n2d19rujZ5r1J9iV5LMmbJnkCkqTnWzVAn+eA91TVg0leBuxKcne37qaq+tvezknOA64GXgP8JvBvSX63qg6Ps3BJ0rEte+VeVQer6sFu/hlgL7D2OJtsAm6rqp9W1XeBfcCF4yhWkjSYExpzTzILnA/c3zW9K8lDSW5J8oqubS3wRM9mC/T5ZZBkS5L5JPOLi4snXLgk6dgGDvckLwVuB95dVT8GbgZeBWwEDgIfOdK1z+b1vIaqrVU1V1VzMzMzJ1y4JOnYBgr3JKezFOy3VtXnAarqyao6XFU/Bz7OL4ZeFoD1PZuvAw6Mr2RJ0nIGuVsmwDZgb1V9tKd9TU+3twKPdPM7gauTvCjJucAG4IHxlSxJWs4gd8tcDLwNeDjJ7q7tfcA1STayNOSyH3gnQFXtSbIDeJSlO22u804ZSTq5lg33qvo6/cfR7zrONh8CPjRCXZKkEfiEqiQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1aNlwT7I+yb1J9ibZk+T6rv2VSe5O8p1u+oquPUk+lmRfkoeSXDDpk5Ak/bJBrtyfA95TVa8GLgKuS3IecANwT1VtAO7plgEuBzZ0ny3AzWOvWpJ0XMuGe1UdrKoHu/lngL3AWmATsL3rth24spvfBHyqltwHnJlkzdgrlyQd0wmNuSeZBc4H7gfOrqqDsPQLADir67YWeKJns4WuTZJ0kgwc7kleCtwOvLuqfny8rn3aqs/+tiSZTzK/uLg4aBmSpAEMFO5JTmcp2G+tqs93zU8eGW7ppoe69gVgfc/m64ADR++zqrZW1VxVzc3MzAxbvySpj0HulgmwDdhbVR/tWbUT2NzNbwbu7Gl/e3fXzEXAj44M30iSTo5VA/S5GHgb8HCS3V3b+4AbgR1JrgW+B1zVrbsLuALYB/wEeMdYK5YkLWvZcK+qr9N/HB3g0j79C7huxLokSSPwCVVJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNWjbck9yS5FCSR3raPpDk+0l2d58reta9N8m+JI8ledOkCpckHdsgV+6fBC7r035TVW3sPncBJDkPuBp4TbfNPyY5bVzFSpIGs2y4V9XXgKcG3N8m4Laq+mlVfRfYB1w4Qn2SpCGsGmHbdyV5OzAPvKeqngbWAvf19Fno2p4nyRZgC8A555wzQhntmL3hi9MuYWz23/jmaZcg/Uob9gvVm4FXARuBg8BHuvb06Vv9dlBVW6tqrqrmZmZmhixDktTPUOFeVU9W1eGq+jnwcX4x9LIArO/pug44MFqJkqQTNVS4J1nTs/hW4MidNDuBq5O8KMm5wAbggdFKlCSdqGXH3JN8BrgEWJ1kAXg/cEmSjSwNuewH3glQVXuS7AAeBZ4Drquqw5MpXZJ0LKnqOyR+Us3NzdX8/PxQ27b0JaROTa18OdzSz0or/yajSrKrqub6rfMJVUlqkOEuSQ0y3CWpQYa7JDVolCdUJWkq/HJ4eV65S1KDDHdJapDhLkkNMtwlqUGGuyQ1yLtlpGW0dGeGfnV45S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGrRsuCe5JcmhJI/0tL0yyd1JvtNNX9G1J8nHkuxL8lCSCyZZvCSpv0Gu3D8JXHZU2w3APVW1AbinWwa4HNjQfbYAN4+nTEnSiVg23Kvqa8BTRzVvArZ389uBK3vaP1VL7gPOTLJmXMVKkgYz7Jj72VV1EKCbntW1rwWe6Om30LU9T5ItSeaTzC8uLg5ZhiSpn3F/oZo+bdWvY1Vtraq5qpqbmZkZcxmS9Ktt2HB/8shwSzc91LUvAOt7+q0DDgxfniRpGMOG+05gcze/Gbizp/3t3V0zFwE/OjJ8I0k6eVYt1yHJZ4BLgNVJFoD3AzcCO5JcC3wPuKrrfhdwBbAP+AnwjgnULElaxrLhXlXXHGPVpX36FnDdqEVJkkbjE6qS1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KBVo2ycZD/wDHAYeK6q5pK8EvgsMAvsB/6sqp4erUxJ0okYx5X7H1XVxqqa65ZvAO6pqg3APd2yJOkkmsSwzCZgeze/HbhyAseQJB3HqOFewJeT7EqypWs7u6oOAnTTs0Y8hiTpBI005g5cXFUHkpwF3J3k24Nu2P0y2AJwzjnnjFiGJKnXSFfuVXWgmx4C7gAuBJ5Msgagmx46xrZbq2ququZmZmZGKUOSdJShwz3JGUledmQeeCPwCLAT2Nx12wzcOWqRkqQTM8qwzNnAHUmO7OfTVfWlJN8AdiS5FvgecNXoZUqSTsTQ4V5VjwOv7dP+P8CloxQlSRqNT6hKUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lq0MTCPcllSR5Lsi/JDZM6jiTp+SYS7klOA/4BuBw4D7gmyXmTOJYk6fkmdeV+IbCvqh6vqp8BtwGbJnQsSdJRVk1ov2uBJ3qWF4A/6O2QZAuwpVt8NsljQx5rNfCDIbc91Xgup6ZWzqWV84CGziUfHulcfutYKyYV7unTVr+0ULUV2DrygZL5qpobdT+nAs/l1NTKubRyHuC5DGJSwzILwPqe5XXAgQkdS5J0lEmF+zeADUnOTfJC4Gpg54SOJUk6ykSGZarquSTvAv4VOA24par2TOJYjGFo5xTiuZyaWjmXVs4DPJdlpaqW7yVJWlF8QlWSGmS4S1KDVnS4t/KKgyS3JDmU5JFp1zKKJOuT3Jtkb5I9Sa6fdk3DSvLiJA8k+VZ3Lh+cdk2jSnJakm8m+cK0axlFkv1JHk6yO8n8tOsZVpIzk3wuybe7n5k/HOv+V+qYe/eKg/8A/pilWy+/AVxTVY9OtbAhJHkd8Czwqar6/WnXM6wka4A1VfVgkpcBu4ArV+i/SYAzqurZJKcDXweur6r7plza0JL8FTAHvLyq3jLteoaVZD8wV1Ur+iGmJNuBf6+qT3R3Ff5aVf1wXPtfyVfuzbzioKq+Bjw17TpGVVUHq+rBbv4ZYC9LTyuvOLXk2W7x9O6zMq+EgCTrgDcDn5h2LYIkLwdeB2wDqKqfjTPYYWWHe79XHKzIIGlRklngfOD+6VYyvG4YYzdwCLi7qlbsuQB/B/w18PNpFzIGBXw5ya7uNSYr0W8Di8A/dUNln0hyxjgPsJLDfdlXHGg6krwUuB14d1X9eNr1DKuqDlfVRpaesL4wyYocMkvyFuBQVe2adi1jcnFVXcDSW2ev64Y1V5pVwAXAzVV1PvC/wFi/N1zJ4e4rDk5B3fj07cCtVfX5adczDt2fy18FLptyKcO6GPiTbqz6NuD1Sf5luiUNr6oOdNNDwB0sDdGuNAvAQs9fg59jKezHZiWHu684OMV0X0JuA/ZW1UenXc8okswkObObfwnwBuDb061qOFX13qpaV1WzLP2cfKWq/mLKZQ0lyRndl/V0wxhvBFbcXWZV9d/AE0l+r2u6FBjrjQeTeivkxJ3kVxxMVJLPAJcAq5MsAO+vqm3TrWooFwNvAx7uxqoB3ldVd02xpmGtAbZ3d2W9ANhRVSv6FsJGnA3csXQdwSrg01X1pemWNLS/BG7tLk4fB94xzp2v2FshJUnHtpKHZSRJx2C4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAb9PwFZ9JsTtcsrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# roll some numbers on our loaded die\n",
    "loadeddie = Die(mu)\n",
    "loadedrolls = loadeddie.roll(1000)\n",
    "\n",
    "hist(argmax(loadedrolls,1), bins=range(loadeddie.sides+1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our probability distribution is just given by $\\mu_k$ for each $k$, i.e. $p(k|\\boldsymbol{\\mu}) = \\mu_k$ (a generalisation of the Bernoulli distribution, called the Categorical distribution)\n",
    "\n",
    "But we want to use one-of-K to represent $k$ by vector $\\mathbf{x}$, so we can write it as\n",
    "\n",
    "$p(\\mathbf{x}|\\boldsymbol{\\mu}) = \\prod_{k=1}^K \\mu_k^{x_k}$ \n",
    "\n",
    "This works because $\\mu_k^{x_k}$ will be 1 for every element except that where $x_n$ = 1, so the product will give the $\\mu_k$ for the state represented by $\\mathbf{x}$. It is just a clever way of \"picking out\" the correct $\\mu_k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.         0.14285714 1.         1.         1.         1.        ]\n"
     ]
    }
   ],
   "source": [
    "print(mu**[0,1,0,0,0,0])\n",
    "\n",
    "# only kth element is non-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define our probability mass function\n",
    "Die.pmf = lambda self, x: prod(self.mu**x, axis=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some examples with our weighted die"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle p(\\mathbf{x}=[1, 0, 0, 0, 0, 0]|\\boldsymbol{\\mu}) = 0.28571428571428575$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle p(\\mathbf{x}=[0, 1, 0, 0, 0, 0]|\\boldsymbol{\\mu}) = 0.14285714285714288$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x1 = [1,0,0,0,0,0] # \"roll a one\"\n",
    "x2 = [0,1,0,0,0,0] # \"roll a two\"\n",
    "\n",
    "p1 = loadeddie.pmf(x1)\n",
    "p2 = loadeddie.pmf(x2)\n",
    "\n",
    "maths(r\"p(\\mathbf{x}=%s|\\boldsymbol{\\mu}) = %s\" % (x1,p1))\n",
    "maths(r\"p(\\mathbf{x}=%s|\\boldsymbol{\\mu}) = %s\" % (x2,p2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now consider the likelihood over N observations -- we simply need to multiply the previous result over all the samples\n",
    "\n",
    "$p(D|\\mu) = \\prod_{n=1}^{N}p(\\boldsymbol{x}_n|\\boldsymbol{\\mu})=$\n",
    "\n",
    "$= \\prod_{n=1}^{N}\\prod_{k=1}^K \\mu_k^{x_k}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Die.likelihood = lambda self, data: product(self.pmf(data))\n",
    "Die.loglikelihood = lambda self, data: sum(log(self.pmf(data)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now given some observed dice rolls, we can plot how likely it is that the number one is weighted a certain way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, '$\\\\mu_1$')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEQCAYAAABiGgneAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd5xU5fX48c+BlQWkCmikCKhYFjsbyFfFgg0rakgECxYMxhZrAogkipXoL9hFYoklCkhi3MSCUVBjQGARUIrIIqCANCk2igvn98e5G4ZhZnd2d2bulPN+veY1s3fuPM+ZZZkzz32aqCrOOedcKtQJOwDnnHO5y5OMc865lPEk45xzLmU8yTjnnEsZTzLOOedSpiDsADJNy5YttUOHDmGH4ZxzWWX69OlrVLVV9HFPMlE6dOhAaWlp2GE451xWEZElsY775TLnnHMp40nGOedcyniScc45lzKeZJxzzqWMJxnnnHMp40nGOedcyniScc45lzI+TyafqMK8efDuu7BpE3TsaLf994cGDcKOzjmXgzzJ5IP16+GWW+Bvf4NVq3Z+vnlzuPZau7Vsmf74nHM5yy+X5bp//xsOPhhGjYITT4SnnoKFC+Hrr6G0FMaMgWOOgWHDoH17GDQItmwJO2rnXI7wlkyu2rYNbrwRHnwQDjwQXnkFiot3PGe33aBLF/jlL2HuXLj7bhg+3C6njR0Le+0VSujOudzhLZlcNXSoJZhrroHp03dOMNGKiuCFFyy5zJ0LRxwBb72VnlidczkroSQjIj1FZL6IlInIoBjPF4rImOD5KSLSIeK5wcHx+SJySlVlikjHoIwFQZn1KqtDRLqKyMzgNktEzokoa7GIfBI8lz+rXj71lLVKfvUreOih6nXq/+IXdhmtdWs4/XR49dXUxemcy32qWukNqAssBPYG6gGzgKKoc64CRgaP+wBjgsdFwfmFQMegnLqVlQmMBfoEj0cCV1ZRR0OgIHi8J7Aq4ufFQMuq3mPkrUuXLprV3npLtW5d1ZNPVt2ypeblfPON6s9+prrLLqqvvZa8+JxzOQko1RifqYm0ZLoCZar6uapuAUYDvaLO6QU8GzweB5wgIhIcH62qm1V1EVAWlBezzOA1PYIyCMo8u7I6VPUHVS0PjtcHNIH3lJtWrLCWSFERvPwy7LJLzctq3BjeeAMOOQTOPdcvnTnnaiSRJNMG+DLi56XBsZjnBB/4G4AWlbw23vEWwPqIpBFZV7w6EJFuIjIH+AT4dcTrFXhLRKaLyIB4b1BEBohIqYiUrl69upJfRYYbOBA2boRx46BJk9qX16yZJZcDDoCzz4ZZs2pfpnMurySSZCTGsejWQrxzknW80jhUdYqqdgZ+CgwWkfrB80ep6hHAqcDVInJMjDJQ1VGqWqyqxa1a7bSxW3b473/huefgpptgv/2SV+5uu8H48XZ/zjmwdm3yynbO5bxEksxSoF3Ez22B5fHOEZECoCmwtpLXxju+BmgWlBFdV7w6/kdV5wHfAwcFPy8P7lcBr2CX6XLP1q02iqxtWxgyJPnl77GHTeRctgz69rX6nHMuAYkkmWlAp2DUVz2s070k6pwS4OLgcW9gQtARVAL0CUaGdQQ6AVPjlRm8ZmJQBkGZr1ZWR1BGAYCItAf2BxaLyK4i0jg4vitwMjA7sV9Llvnzn2HmTLj/fth119TU0a0bPPaYXT5LRSJzzuWkKidjqmq5iFwDjMdGhT2tqnNEZBg2mqAEeAp4XkTKsNZFn+C1c0RkLDAXKAeuVtWtALHKDKocCIwWkTuBGUHZxKsDOBoYJCI/AtuAq1R1jYjsDbxiYwkoAF5U1Tdr9mvKYN98Yx/6xx1nkypTqX9/G948fLjV17NnautzzmU9scaDq1BcXKylpVk0pWbECJvZP3Uq/PSnqa9v0yab2Ll2LcyebX01zrm8JyLTVXWnWd8+4z+bbd1qky2PPjo9CQagfn14/nlYswauuio9dTrnspYnmWz26quweDHccEN66z38cLjtNltc86WX0lu3cy6r+OWyKFl1uax7d1i6FMrKoG7d9NZdXm71f/qprXW2557prd85l1H8clmuKS2FDz6A3/wm/QkGoKAAnn3WJn+muyXlnMsanmSy1QMP2NIv/fuHF8N++9lmaGPG2IRN55yL4kkmG61YYR/sl12WnOVjamPgQEs2V11lrRrnnIvgSSYbjR1rfSJXXBF2JFBYCCNHwuefw113hR2Ncy7DeJLJRqNH2+rIBx4YdiTm+OPhoovgj3+0gQDOORfwJJNtliyByZOhT5+qz02n+++3zdFuvjnsSJxzGcSTTLYZO9buzzsv3Dii7b67bfn82ms+CMA59z8+TyZKxs+TKS6GOnVsGZlMs3kzdO5s/TSzZtkwZ+dcXvB5MrmgrAymT8+8VkyFwkK7bDZ3LowaFXY0zrkM4Ekmm4wZY/epXm25Nnr1soEAv/89rFsXdjTOuZB5kskmo0fDUUdBu3ZVnxsWEVsZeu1auOeesKNxzoXMk0y2mDfPltbPtFFlsRx6KFx4ITz8sK2t5pzLW55kssXrr9t9r17hxpGoYcNg2zZbrdk5l7c8yWSLt96CoqLMvlQWqUMHW2rmmWesFeacy0ueZLLBxo3w/vtw8slhR1I9Q4bArrvavXMuLyWUZESkp4jMF5EyERkU4/lCERkTPD9FRDpEPDc4OD5fRE6pqkwR6RiUsSAos15ldYhIVxGZGdxmicg5icadNf7zH9v2ONuSTMuW8NvfwiuvwIcfhh2Ncy4EVSYZEakLPAqcChQBfUWkKOq0/sA6Vd0XGAEMD15bBPQBOgM9gcdEpG4VZQ4HRqhqJ2BdUHbcOoDZQLGqHhbU8YSIFCQYd3YYPx7q1YNjjw07kuq74QZo1cpWA3DO5Z1EWjJdgTJV/VxVtwCjgeje517As8HjccAJIiLB8dGqullVFwFlQXkxywxe0yMog6DMsyurQ1V/UNXy4Hh9oGIJg0Tizg5vvWW7UDZsGHYk1deoEQwaBG+/bZf8nHN5JZEk0wb4MuLnpcGxmOcEH/gbgBaVvDbe8RbA+oikEVlXvDoQkW4iMgf4BPh18HwicRO8foCIlIpI6erVq+P+IkKxbJkNXT7llKrPzVRXXmnbMw8dCr6MkXN5JZEkIzGORX9SxDsnWccrjUNVp6hqZ+CnwGARqZ9g3ASvH6Wqxapa3KpVq1inhOff/7b7bOuPidSgge2g+f778M47YUfjnEujRJLMUiBy3GxbYHm8c0SkAGgKrK3ktfGOrwGaBWVE1xWvjv9R1XnA98BBCcad+d56C/bYAw4+OOxIaufyy6FtW2/NOJdnEkky04BOwaivelhHfknUOSXAxcHj3sAEteWdS4A+wciwjkAnYGq8MoPXTAzKICjz1crqCMooABCR9sD+wOIE485s27ZZS+bkk23l5WxWvz7cequNMnvjjbCjcc6lSZWfXEH/xjXAeGAeMFZV54jIMBE5KzjtKaCFiJQBNwKDgtfOAcYCc4E3gatVdWu8MoOyBgI3BmW1CMqOWwdwNDBLRGYCrwBXqeqaKurIDjNmwJo12X2pLNKll0L79nD77d6acS5P+H4yUTJqP5kRI+DGG63zv3XrsKNJjlGj4Ior4M03s3swg3NuB76fTDaaNMmWZ8mVBANwySW2NI63ZpzLC55kMpWqJZkjjww7kuSqVw8GD4bJk32kmXN5wJNMpvriC1i+PPeSDMBll0GbNt6acS4PeJLJVJMm2X0uJpnCQlsF4IMP4N13w47GOZdCnmQy1aRJtoJxts+Piefyy20VgDvuCDsS51wKeZLJVJMmQbduUFBQ9bnZqH59uPlmmDjR+meccznJk0wm+u47mDUrNy+VRbriCmjRAu66K+xInHMp4kkmE02bBlu35n6S2XVX2wrgtddg5sywo3HOpYAnmUxU0en/f/8XbhzpcPXV0KQJ3H132JE451LAk0wmmjQJOneGZs3CjiT1mjWDa66BceNg3rywo3HOJZknmUyzbZt1hOf6pbJI119v2wHce2/YkTjnksyTTKaZPx/WrcuvJNOqFQwYAC++CEuWhB2Ncy6JPMlkmorhvPnQHxPppptABO6/P+xInHNJ5Ekm03z0ETRuDJ06hR1JerVtC/36wZNPwsqVYUfjnEsSTzKZZuZMOPTQ7N+krCZ+9zvYvBkefDDsSJxzSZKHn2QZbNs2m4R52GFhRxKO/faD3r3h0Udhw4awo3HOJYEnmUyycKHN9j/88LAjCc/gwfDNN/DYY2FH4pxLAk8ymWTGDLvP5yRz+OHQsyc88ABs3Bh2NM65WkooyYhITxGZLyJlIjIoxvOFIjImeH6KiHSIeG5wcHy+iJxSVZki0jEoY0FQZr3K6hCRk0Rkuoh8Etz3iCjr3aCOmcFt9+r/itJoxgxbELOoKOxIwjV4MKxaBc88E3YkzrlaqjLJiEhd4FHgVKAI6Csi0Z+C/YF1qrovMAIYHry2COgDdAZ6Ao+JSN0qyhwOjFDVTsC6oOy4dQBrgDNV9WDgYuD5qNguUNXDgtuqKn8jYZo502b6FxaGHUm4une3Idz33Qfl5WFH45yrhURaMl2BMlX9XFW3AKOBXlHn9AKeDR6PA04QEQmOj1bVzaq6CCgLyotZZvCaHkEZBGWeXVkdqjpDVZcHx+cA9UUkOz+lZ8zI307/SCK2qdnixTBmTNjROOdqIZEk0wb4MuLnpcGxmOeoajmwAWhRyWvjHW8BrA/KiK4rXh2Rfg7MUNXNEceeCS6VDQ2SWGZascLmh+Rzf0ykM86wy4b33utbNDuXxRJJMrE+mKP/18c7J1nHq4xDRDpjl9CuiHj+guAyWvfgdlGMMhCRASJSKiKlq1evjnVK6nmn/47q1LHWzOzZ8PrrYUfjnKuhRJLMUqBdxM9tgeXxzhGRAqApsLaS18Y7vgZoFpQRXVe8OhCRtsArQD9VXVhRqKouC+6/BV7ELtPtRFVHqWqxqha3atWqkl9FClXsp3LooeHUn4n69IG99oJ77gk7EudcDSWSZKYBnYJRX/WwjvySqHNKsE53gN7ABFXV4HifYGRYR6ATMDVemcFrJgZlEJT5amV1iEgz4DVgsKr+tyIgESkQkZbB412AM4DZCbzfcMyYAXvvDU2bhh1J5thlF9ui+b//hQ8+CDsa51wNVJlkgv6Pa4DxwDxgrKrOEZFhInJWcNpTQAsRKQNuBAYFr50DjAXmAm8CV6vq1nhlBmUNBG4MymoRlB23jqCcfYGhUUOVC4HxIvIxMBNYBvy5+r+iNPFO/9j694eWLX0bAOeylKh3qu6guLhYS0tL01vpt9/a7pB33AG33preurPBnXfC0KG25M4hh4QdjXMuBhGZrqrF0cd9xn8mmDXL7r3TP7arr4ZGjeCPfww7EudcNXmSyQQVnf5+uSy25s3hiitg9GhYtCjsaJxz1eBJJhPMmWN73bduHXYkmeuGG2xYs29q5lxW8SSTCebNs4mHGTxXNHRt2timZk8/7ZuaOZdFPMlkgnnz4MADw44i8/mmZs5lHU8yYfv6a1tx2JNM1XxTM+eyjieZsM2bZ/f5vrx/ogYNsk3NRo4MOxLnXAI8yYStIsl4SyYxRxwBJ58MI0b4pmbOZQFPMmGbOxcaNrQ1ulxiBg2yzv9nn636XOdcqDzJhG3ePNh/fxue6xJz3HHQrZtNzvRNzZzLaP7JFraK4csucSK2RfOiRb6pmXMZzpNMmL77Dr74wvtjauLMM22r6nvvhW3bwo7GOReHJ5kwffqp3XuSqb46daw1M3s2/OtfYUfjnIvDk0yYfPhy7Zx3HnTsCHff7Vs0O5ehPMmEae5cKCiAffYJO5LsVFBgqwBMmQLvvht2NM65GDzJhGnePOjUyXaAdDVzySXwk5/AXXeFHYlzLgZPMmHyNctqr359uOkmeOcdmDo17Gicc1E8yYRl82YoK/P+mGS44grYbTdvzTiXgTzJhGXBAht66y2Z2mvcGK6/HkpKtu8y6pzLCAklGRHpKSLzRaRMRAbFeL5QRMYEz08RkQ4Rzw0Ojs8XkVOqKlNEOgZlLAjKrFdZHSJykohMF5FPgvseEWV1CY6XichDIhm0YYuvWZZc11xjyebuu8OOxDkXocokIyJ1gUeBU4EioK+IRF/j6Q+sU9V9gRHA8OC1RUAfoDPQE3hMROpWUeZwYISqdgLWBWXHrQNYA5ypqgcDFwPPR8T1ODAA6BTcelb5G0mXefNs5vr++4cdSW5o3twSzcsvw/z5YUfjnAsk0pLpCpSp6uequgUYDfSKOqcXULFa4TjghKDV0AsYraqbVXURUBaUF7PM4DU9gjIIyjy7sjpUdYaqLg+OzwHqB62ePYEmqjpZVRV4LqKs8C1YAO3a2eKYLjluuMEGAtxzT9iROOcCiSSZNsCXET8vDY7FPEdVy4ENQItKXhvveAtgfVBGdF3x6oj0c2CGqm4Ozl9aRdwAiMgAESkVkdLVq1fHOiX5ysp8fkyytWplgwBeeAE+/zzsaJxzJJZkYvVjRE+vjndOso5XGYeIdMYuoV2RyPk7HFQdparFqlrcqlWrWKck38KFsO++6akrn/z2tzZJ8957w47EOUdiSWYp0C7i57bA8njniEgB0BRYW8lr4x1fAzQLyoiuK14diEhb4BWgn6oujDi/bRVxh+Obb2D1am/JpELr1nD55fCXv9jio865UCWSZKYBnYJRX/WwjvySqHNKsE53gN7AhKAfpAToE/SRdMQ636fGKzN4zcSgDIIyX62sDhFpBrwGDFbV/1YEpKpfAd+KyM+Cvp5+EWWFa2GQB70lkxoDB9r98OGVn+ecS7mCqk5Q1XIRuQYYD9QFnlbVOSIyDChV1RLgKeB5ESnDWhd9gtfOEZGxwFygHLhaVbcCxCozqHIgMFpE7gRmBGUTrw7gGmBfYKiIDA2Onayqq4Argb8ADYA3glv4ysrs3lsyqdGuHVx6KTz5JNxyC7SJ2RXnIqjagMe337aFExYtstvKlVC3rl2BbNDAVkE68EA46CA45RTbbSGDJga4DCTqq9fuoLi4WEtLS1Nbyb332jL133xjcztc8i1ebJ+IV10FDz4YdjQZa84ceOIJ+NvfYHlwMbltW2tkd+wIe+5pCai8HL791kaHz5sHK1bYuR07Qq9eMGCAT/nKdyIyXVWLo49X2ZJxKVBWBrvv7gkmlTp0gIsuglGjYNAg+7R0gC008fe/w4gRMGkS1Ktne8D17Aknnmi/uqosX27b+JSUwOOPwwMPwOmnw803w7HHeuvGbefLyoTBR5alx5Ah8OOP3jcTUIW33oKuXeEXv4BVq+D++2HZMhg3zsZLJJJgwMZXDBhgiebLL+H22+0y2/HHW7Ly+bCugieZMPgcmfTYZx/o1w9Gjtx+LShPLV4Mp51m/Shr1sCzz9rGrDfdBC1b1q7sVq3g97+HJUusdfThh3DwwTb+4vvvkxK+y2KeZNJt40ZYutRbMukyZIh1KORpa2bbNnj4Yeuo/+AD+NOfrJXRr5916CdTgwa2Tulnn8GFF8If/wjFxb5mab7zJJNuixbZvbdk0mOffeDii613e9mysKNJq+XLoUcP+M1v4OijYfZsW3mnsDC19e6xBzz9tG3xs2EDdOsGjz7qO2TnK08y6eZzZNLv1lth69a8WgVgwgQ4/HCYNs0+8N94A9q3T28MPXpYK6ZHD1u79KKLbBsll188yaRbRZLxlkz6dOxo2zSPGmWXKnOYqq0PetJJto/btGk2ZSis0V6tWtnggDvugL/+FU4+GdauDScWFw5PMulWVgZNm0KL6LU9XUoNGWKfwDm8e+aWLZZQbrkFfvlLSzCZsPFqnTrWmHzxRRsUcOSRvn5pPvEkk24LF1orxicSpFeHDvCrX9kqADn4Cbd+PZx6qo0au+02+0Bv1CjsqHbUt6+tKLBqFRxzzPaFL1xu8ySTbmVl3h8TliFDbH2UYcPCjiSpvvoKuneH99+3dUH/8IfM/Q7TvTtMnAibNsFxx9m2Si63eZJJp/Jym7Dg/THhaN3alpl5/nmbJJIDvvjCWgWLFsGbb9pAukx36KE2MGHzZks0n30WdkQulTzJpNMXX1ii8ZZMeAYOtAkdt90WdiS1tnChJZhVq2wm/wknhB1R4g45xFo0P/5ogxTyfK5sTvMkk04+six8u+8O110HY8Zk9SzBhQttjbBvv7VWwZFHhh1R9R10kLW+1q61pWjWrw87IpcKnmTSqaKn01sy4br5ZmjWzPpostCXX1qrZeNGaw106RJ2RDV3xBG2WOenn8LZZ1tfjcstnmTSaeFCqF/fVwQOW/PmtjLza6/Bf/4TdjTVsmKFJZh16+wS2SGHhB1R7Z10kg1YeO8961PylQFyiyeZdFq0yIbS1vFfe+iuvdaS/aBBWfOptm6dTWZctgxefz27WzDRzj/fFmQYOzanpzLlJf+0S6clS9K/toeLrWFDG+s7aZJNSc9wmzbBWWfZ4pYlJXDUUWFHlHy/+50trDl0KLyaGRuluyTwJJNOnmQyy2WX2e6Zt9xia5tlqK1b4YILbBXl55/PrlFk1SFiK//89KeWbD75JOyIXDIklGREpKeIzBeRMhEZFOP5QhEZEzw/RUQ6RDw3ODg+X0ROqapMEekYlLEgKLNeZXWISAsRmSgi34nII1FxvRvUMTO47V69X08Sff+9beThSSZz7LIL3HmnLU/8wgthRxOTqg2Gq9jJ8pe/DDui1GrQAF55xTaNPftsW8XZZbcqk4yI1AUeBU4FioC+IhK9IlJ/YJ2q7guMAIYHry0C+gCdgZ7AYyJSt4oyhwMjVLUTsC4oO24dwCZgKHBznLdwgaoeFtxWVfV+U+aLL+zek0xm6d3bvjoPHWrDtTLMiBG2TP7NN9teLfmgTRt4+WVr+F92WdZ0mbk4EmnJdAXKVPVzVd0CjAZ6RZ3TC3g2eDwOOEFEJDg+WlU3q+oioCwoL2aZwWt6BGUQlHl2ZXWo6veq+gGWbDKXJ5nMVKcO3HefjQt+8MGwo9lBSYkll96982/PtaOOsoEAf/87PPRQ2NG42kgkybQBvoz4eWlwLOY5qloObABaVPLaeMdbAOuDMqLrildHVZ4JLpUNDZLYTkRkgIiUikjp6tWrEyiyBpYssfu99kpN+a7mjj0WzjzT1shP1b9/Nc2caSOuiott0ct8HJB400022OHmm231ZpedEvnTjfXBHN2AjXdOso4nGke0C1T1YKB7cLso1kmqOkpVi1W1uFWrVlUUWUNLlth+t61bp6Z8VzvDh8N339nGJyH76is44wybzvPqqzYQLh+J2PyZtm3hvPN8RYBslUiSWQq0i/i5LRC90tD/zhGRAqApsLaS18Y7vgZoFpQRXVe8OuJS1WXB/bfAi9hlunAsWWL/WwoKqj7Xpd+BB8Lll8Pjj4e6NPDmzXDuuTYn5p//9Hm7zZvD6NE2N+jKK71/JhslkmSmAZ2CUV/1sI78kqhzSoCK9V97AxNUVYPjfYKRYR2BTsDUeGUGr5kYlEFQ5qtV1BGTiBSISMvg8S7AGcDsBN5vavjw5cx3++1QWGgTNkKgah+kH34Izz0Hhx0WShgZp1s3W8909OiMHQToKqOqVd6A04DPgIXAkODYMOCs4HF94GWsY38qsHfEa4cEr5sPnFpZmcHxvYMyyoIyCxOoYzHWqvkOa/EUAbsC04GPgTnAg0Ddqt5rly5dNCXatVO96KLUlO2S5667VEH1nXfSXvVDD1nVt96a9qozXnm56tFHqzZurLpwYdjRuFiAUo3xmSrq7c8dFBcXa2lpaXILLS+3b8i33JIR1/xdJTZutEtnTZvCRx9ZP1oavPsunHginHYa/OMf+dnRX5UlS2wvmqIiW3IuTf80LkEiMl1Vi6OP+59yOixbBtu2+ciybNCggQ1p/vhjePrptFS5dKlNsuzUyS4HeYKJrX17eOwxmDwZ7r8/7GhcovzPOR0qhi97n0x26N0bjj7atgJI8ZTzzZutuo0bbaZ7kyYprS7r9e0LP/85/P73vuxMtvAkkw6eZLKLCDzwgC0DlOLLm9dfD1Om2FDdAw5IaVU5QcQGADZtatsC/Phj2BG5qniSSQefiJl9unSB/v1tFYB581JSxV/+AiNH2mC2n/88JVXkpFat4IknYMYM3xYgG3iSSYclS2zb3wYNwo7EVcfdd0OjRvCb3yR9gsasWTZc+fjj/YOyJs45x1ZqvvNOWx3BZS5PMunwxRfeislGrVrZp9jbb8Pf/pa0YjdssH6Y5s3hpZd8fm5NPfggtGxpi2j6ZbPM5UkmHXwiZva64gobN3vjjbZdQy2p2ofiokUwZgzssUcSYsxTu+1mo81mzPDRZpnMk0yqqVpLxpNMdioogEcesVWak3Bda8QIW1l4+HDo3j0J8eW5c8+1VuHtt8Onn4YdjYvFk0yqrV5t41M9yWSvo4+2oUz33w9z59a4mMmTYeBA24zrxhuTGF+ee+QR2HVXG6eRwRuc5i1PMqnmw5dzw3332XaNv/61TaytpjVrbMJlu3bwzDM2FNclxx572IjzSZNstJ7LLJ5kUs2TTG5o1Qr++Edbz+TZZ6s+P8K2bdCvH6xaZTs+NmuWohjz2IUXwimnwODBtoKCyxyeZFKtYkdMH12W/S691C6d/fa31jRJ0PDh8MYb9m27S5cUxpfHKiZplpfD1Vf7lgCZxJNMqi1ZYnMtmjcPOxJXW3Xq2PWYDRts28YEvP8+3Hor9OljV9pc6nTsCMOG2bbVSRxx7mrJk0yqVcyR8YvwuaFzZ+u9f+45GD++0lNXrbK1tvbZB0aN8j+BdLj+ejjiCLj2Wt9JM1N4kkm1pUttR0yXO2691RYaGzAAvv025inbtsFFF8HatdYP07hxmmPMUwUF8Oc/W4IfPDjsaBx4kkm9ZcugTZuwo3DJVL8+PPWUzZ255ZaYp9x9N7z1Fjz0kM3ldOlzxBFw3XV2ZXPy5LCjcZ5kUqm8HFau9CSTi4480q7JPPoofPDBDk9NnAh/+ANccAFcfnlI8eW5YcNsuPiAAb7kTNg8yaTSihV23cSTTG666y4bmn7ZZfDDD4B9pzj/fNhvP/I5A7MAABtWSURBVPsm7f0w4WjUyCZpzp4Nf/pT2NHkt4SSjIj0FJH5IlImIoNiPF8oImOC56eISIeI5wYHx+eLyClVlSkiHYMyFgRl1qusDhFpISITReQ7EXkkKq4uIvJJ8JqHRNL8X37ZMrv3JJObGjWy3TMXLIBBg9i61RLMhg3WD9OoUdgB5rezzrJlZ26/HT7/POxo8leVSUZE6gKPAqcCRUBfESmKOq0/sE5V9wVGAMOD1xYBfYDOQE/gMRGpW0WZw4ERqtoJWBeUHbcOYBMwFLg5RviPAwOATsGtZ1XvN6k8yeS+44+3rQAefphhly5iwgS7gnbQQWEH5sBWaq5b1+fOhCmRlkxXoExVP1fVLcBooFfUOb2AimnQ44ATglZDL2C0qm5W1UVAWVBezDKD1/QIyiAo8+zK6lDV71X1AyzZ/I+I7Ak0UdXJqqrAcxFlpcfy5XbvSSa33XMPb7W5lDueb8/Ffbdw6aVhB+QqtG1ruzW8+aa1Ll36JZJk2gBfRvy8NDgW8xxVLQc2AC0qeW284y2A9UEZ0XXFq6OyuCMXmIgVNwAiMkBESkWkdPXq1ZUUWU3LlsEuu9iSJC5nLV3bkAu+f4LOzOGxuteGHY6Lcs01ttLCddf53JkwJJJkYvVjRDc8452TrOOJxpFITDsfVB2lqsWqWtwqmQlh2TLYc0+bKe5y0o8/wnnnwabyXRh35QQavjDKvzJnmLp1bbvmVatgyJCwo8k/iXz6LQXaRfzcFlge7xwRKQCaAmsreW2842uAZkEZ0XXFq6OyuCNnQcaKO7V8jkzOGzTIVv998knY/8GroFs3GzdbsWadywhduliL5vHHYcqUsKPJL4kkmWlAp2DUVz2sI78k6pwS4OLgcW9gQtAPUgL0CUaGdcQ636fGKzN4zcSgDIIyX62ijphU9SvgWxH5WdDX0y+irPTwJJPTxo2z4bFXXWWtGXbZBf76V5sfdeGFvrlJhrnjDmjd2ufOpJ2qVnkDTgM+AxYCQ4Jjw4Czgsf1gZexjv2pwN4Rrx0SvG4+cGplZQbH9w7KKAvKLEygjsVYq+Y7rAVTFBwvBmYHdTwCSFXvtUuXLpo0jRqpXndd8spzGePTT1UbN1bt1k1106aoJ597ThVUhw0LJTYX39//bv80990XdiS5ByjVGJ+poj6ubwfFxcVaWlpa+4K++QaaNrU9SH7729qX5zLG99/bVbEVK2x/+Xbtok5QtZbM6NEwYQIce2wocbqdqdrOpG+/DXPmQIcOYUeUO0RkuqoWRx/3HulU8TkyOUnVLrfMnQsvvRQjwYBN8x85Evbd19b4X7ky7XG62ETg4Yft3ufOpIcnmVTxJJOTHn4YXnzR1sY66aRKTmzc2DptNmywZQC8fyZj7LWX9c+8/rr9E7nU8iSTKp5kcs7779teZWedFXfx5R0dfDA89phdMrvttlSH56rh2mtttebf/AbWrQs7mtzmSSZVPMnklGXL4Je/hL33tv3KEp76dMkltoDmnXfalo0uI0TuOzNop9UYXTJ5kkmVZctsy+UGDcKOxNXS5s3Qu7d1+L/yio3nqJZHHoHiYhsM8OmnKYnRVd8RR8ANN9iupf/5T9jR5C5PMqnic2RygipceSV8+CH85S9QFL00bCIaNIC//93ue/WyfhqXEW6/3UaYDRhgXyZc8nmSSRVPMjnh4YfhmWfg97+Hn/+8FgW1a2e9zJ9/bruZ+UCAjLDrrrYKwKef2m6mLvk8yaSKJ5ms9847cOONNq/iD39IQoHdu9t+zK+95nOnMkjPnnYl8+674eOPw44m93iSSYWKbZdbtw47EldDCxZYR/8BB1Szo78qV15pQ5pGjLCv0C4jjBhhXaj9+9t/X5c8nmRSYeVK33Y5i61bB2ecYRP2SkpsyktS/elPcPrpNo52/PgkF+5qomVLG59RWgoPPBB2NLnFk0wq+PDlrPXjjzaSbPFi+Mc/bMhy0tWta8sFHHQQ/OIXtjaNC90vfmHjMoYOtZasSw5PMqngSSYrqdpSIxMm2ByKo49OYWWNG8O//gXNmsGpp/om9BlAxObOFhbaZbNt28KOKDd4kkkFTzJZafhwSy6DB0O/fmmosG1bu1z2449w8sm+xlkGaN0aHnzQ5s089FDY0eQGTzKp4NsuZ52//tWSS9++Njk/bQ480BbR+uora9H4HJrQ9etnfXKDB8Nnn4UdTfbzJJMKvu1yVpkwAS69FI47zubEpP2frVs3m0Mze7Ylmm+/TXMALpKIrQLQoIGtCuRTmmrHPwVTwefIZI2PP4Zzz4VOnWzJmMLCkAI59VTbf2bqVDjzTPjhh5ACcWDfER95BCZPhvvvDzua7OZJJhVWrLC/UpfRFi60rpBGjeCNN6wPPlTnngsvvGAdAr16wcaNIQeU3/r2tVUehg71AYC14UkmFVasgJ/8JOwoXCW++sr2gykvh3//2/YYyQh9+tg1u3fesbk0330XdkR5SwSeeMK6Vi+4wBuXNZVQkhGRniIyX0TKRGSnhbFFpFBExgTPTxGRDhHPDQ6OzxeRU6oqU0Q6BmUsCMqsV4s6FovIJyIyU0SSsKdyArZsgbVrPclksLVr4ZRTbJn311+3vveM0q+fLTPw3nu25okPBghNixa2MOq8efC734UdTZZS1UpvQF1gIbA3UA+YBRRFnXMVMDJ43AcYEzwuCs4vBDoG5dStrExgLNAneDwSuLImdQTPLQZaVvUeI29dunTRWvnyS1VQfeKJ2pXjUmLdOtUuXVTr1VP997/DjqYKL7+sWlCgWlysumZN2NHktRtusP/Wr70WdiSZCyjVGJ+pibRkugJlqvq5qm4BRgO9os7pBTwbPB4HnCAiEhwfraqbVXURUBaUF7PM4DU9gjIIyjy7hnWEY8UKu/eWTMb55hvrX//4Y1t5/8QTw46oCr1722iETz6Bo46yZQhcKO6+2zY6veQSWL487GiySyJJpg3wZcTPS4NjMc9R1XJgA9CiktfGO94CWB+UEV1XdesAUOAtEZkuIgMSeK+150kmI333nXVxlJbCyy/b46xwxhnw1ls2UfPII2HWrLAjykv168PYsdYvc/75vohmdSSSZCTGMU3wnGQdr0kdAEep6hHAqcDVInJMjHMRkQEiUioipatXr451SuI8yWScDRtsFNnkyfDiizZwK6sccwx88IGteXbMMTZSwaXdAQfYwtnvvWebnbnEJJJklgLtIn5uC0Q3GP93jogUAE2BtZW8Nt7xNUCzoIzouqpbB6pacb8KeIU4l9FUdZSqFqtqcavaztKvWBpk991rV45Liq+/hhNOsBbMmDG2CGJW6twZJk2C9u3tmt+jj4YdUV666CK47DK46y5rYLqqJZJkpgGdglFf9bBO95Koc0qAi4PHvYEJQUdQCdAnGBnWEegETI1XZvCaiUEZBGW+WpM6RGRXEWkMICK7AicDsxP7tdTCihU24aJ+/ZRX5Sq3ciX06GET6V95pZY7W2aCdu3gv/+F006Da66Bq66ydc9cWj38sOX888/3brJEVJlkgv6Pa4DxwDxgrKrOEZFhInJWcNpTQAsRKQNuBAYFr52DjRabC7wJXK2qW+OVGZQ1ELgxKKtFUHa16wD2AD4QkVlYYntNVd+syS+pWnyOTEZYuND6ysvKbLHjrOmDqUrjxpYxBw60azcnnOA90WnWsKENHCkvh3PO8fkzVRFrDLgKxcXFWlpaiyk1xxxji1+9+27SYnLVM326fdnfutV2Ou7WLeyIUuSll+Dyyy3xjB5ti6+5tHn9dRuX0aePLbAqsXqH84iITFfV4ujjPuM/2bwlE6o337TP2gYN7MpSziYYsHVPpk61y7MnnGDLR/tqjmlz2mn2K3/pJV/frDKeZJLNk0woVG0fkNNPh333tT7y/fcPO6o06NwZpk2D886zRbaOOw6WLAk7qrwxeLBNZxo40C6huZ15kkmmH36wZdo9yaTVjz/Cr38N118PZ51l60u2bh12VGnUuLFdr3n+eZtHc8ghtiyNXwpPORF49llrMV9wgQ2TdzvyJJNMFcOX99gj3DjyyPLlNoJs1Ci45Rb4299sVeW8IwIXXmhJ5uCD4eKLrcPgyy+rfq2rlYYNoaTEdvc480xYsCDsiDKLJ5lk8omYaTVxIhx+uC3D/uKLNnch7/eJ69jRZgs+8IANPunc2Tau976alGrVyraLAFvT1Af8bZfv/yWTy5NMWpSXwx132NpjLVpY33ffvmFHlUHq1oXrrrM1z376U7j6aujaFaZMCTuynNapk41mXLXKxmGsWhV2RJnBk0wyVSQZv1yWMgsX2ijx3//eho5OnQpFRWFHlaH23hveftuGP61YAT/7mU1XX7Ys7MhyVrdulmiWLLEvQV9/HXZE4fMkk0wrV9q18douTeN2ogpPPgmHHQZz51o/91//mqf9L9UhYtn400/h5pvtl9apE9x6qy1L7ZLumGOsj+azz2xjvNouh5jtPMkk04oV0LIl7LJL2JHklAUL7PLDr34FxcW2VP/554cdVZZp3Bjuu8+SzdlnWwfW3nvDvff67pspcOKJtjDDvHnQvTt88UXYEYXHk0wy+RyZpNq0yT4LDzkEPvrIRpC9804GbZWcjTp2tFES06ZZP83gwdChA9xzD6xfH3Z0OeXUU20Rza++siWO5s0LO6JweJJJphUrvD8mCVRh3DjbFvnWW22C5dy51pLJ+9FjyVJcbOuifPihJZtbbrEFOG+4wSdzJlH37jbY78cf4eijYcKEsCNKP/8vm0wrV3pLppbefx+OPdaW5G/c2Fou48bl2eTKdOrWzZLNzJl2Ge2RR+wy2tlnw/jxsG1b2BFmvcMOsyWO9tjD9jV64IH8mifrSSZZVP1yWS1Mnmz/AY891lZOfvxxu0TWo0fYkeWJQw+1FQM+/9zWSJk0ySZ8dOpk48V9Tfta2WcfG0F+5pnWWLz44vzpCvMkkyzffGOdCH65LGHbttkonGOOsZ2FZ8yA//f/bJjyr38NBQVVl+GSrF0729D+yy+t76Z9exsv3rGjrYs2cqRPAKmhxo1tRYphw+CFF7a3cHKdJ5lk8YmYCfv6a/jTn6zPpVcv6wIYMQIWLYIbb7QVlF3ICgtthuuECdaKueMO+xu/8krYc09rYj74oLV8XMLq1LF1TN97z75kHXMMDBoEGzeGHVnqeJJJlop1yzzJxLRli01SO/9861+56SYb7f3ii9Zyuf56n/OSsdq3txEY8+bZ+PEhQyzhXH+9XQcqKrJrQK+9lj/XgGqpe3dbZq5/fxg+HA44wLYHz8W+Gt+0LEqNNy0bO9aWW//kEzjooOQHloV++MG+CL/6qi2DvnYtNG9u6zj+6le2jqPLYgsXwj//acnlP/+BzZvtGmeXLvYp2r27rTKw++5hR5rR3nvPVgGaNctGoN1+Oxx/fPZtghZv0zJPMlFqnGQeesj+UlatytsZ/9u2WY6dONFWM5kwwS4DNGpkS/D37Wud+/XqhR2pS7qNG22wwDvv2BDBadOs+QrWEura1VYzPfRQ64zYc8/s+xRNoa1b4Zln7FLaihU2wnzgQBvkly19k55kElTjJDNkiLV7t2zJm8kcK1fayNcpU7bfKtZq2mcf2znwzDPtunNhYbixujTbtMkSzdSp22+RI9SaN7dOuaIi211u331tJFvHjrZ2fp7atMm2ArrvPhtlufvu9uXsoovgiCMyOy/XKsmISE/gQaAu8KSq3hv1fCHwHNAF+Bo4T1UXB88NBvoDW4HfqOr4ysoUkY7AaGA34CPgIlXdksw6KlPjJNO/v+39m2OLD377rQ00WrTIlndZsMBWJvnkk+1rMonYZ0XXrjYA6fjjbZCScztYv976dGbNgjlzrI9n7lxYs2bH81q1stbPXnvZJi2tW9ttjz2231q0yOkm8datdhXyuefsiuSWLfZ/6qST7GrAUUfZryaTkk6Nk4yI1AU+A04ClgLTgL6qOjfinKuAQ1T11yLSBzhHVc8TkSLgJaAr0Bp4G9gveFnMMkVkLPB3VR0tIiOBWar6eDLrqOz91jjJnH66rR/x0UfVf22KqNpM482b7bZxI3z//fYNPL/5xm7r11t/yddf2//3FSvstnw5rFu3Y5nNmsF++1l/ysEH25IvXbpAkybhvEeXA9ats6/tCxZYa2fJErtfutRu8RbybNLEkk3z5vaH2by5Hau47bqr3Ro1stZRw4Y2dLF+/e23wkJLVoWFtuZgxa2gwLZMyIBP8XXrbOjz+PF2Gbpi9Z+WLa11c8AB1gDs0MFyccuW9mtp0iS94dcmyfwfcJuqnhL8PBhAVe+JOGd8cM5kESkAVgCtgEGR51acF7xspzKBe4HVwE9UtTyy7mTVERl3LDVNMmc1e4+FWzvAXu0Tfk28X33kcdXYt23btt+2bt1+Ky+3xFJxq44mTewP9Cc/sS+Le+65/Qtl+/Z2NaNFi4z4f+fyyXff2Re4lSu3377+evu3ovXr7bZu3Y7fnpLRFVCRbCJvdepsv69Tx/5DVNwncquQyOMoW7UO0zd1ZtrGg5ix6UA+2lTEgh/b8922XWOe30A2smudjRTKFgpkKwVspY5sQ9CI23YzVrelsEnNrm3HSzKJdCm1ASL3cF0KdIt3TpAcNgAtguMfRr22TfA4VpktgPWqWh7j/GTVsRMRGQAMANirhqsv7rM3FG7dvL0NlaB4f0/Rf3ORt8i/64ovWwUF2/8/RH4hKyzc/oWtQQP7YtewoX25q/jC17Qp7LabLx7tMlSjRvYNp1OnxF+jur3p/t131nyv+HnzZuv82LTJrkNt2WLHKr6Zbdmy4ze2yG9wFd/qKu4jv/HF+0YYeYuMr6rHMdQFuqJ05RPgE1BFFdZuacSi73dnxaZmfL2lMas3N+Hb8gb8sLWQH7bWY9PWemzVOpRrHbZqHVQrUsyOpE7yr3MnkmRifQzuFFucc+Idj9UzXtn5yaxj54Oqo4BRYC2ZWOdUZcRHx9bkZc65VBDZfoksx0d7CvZtu0XYgcSRyDCopUBkemsLRO9g/b9zgktZTYG1lbw23vE1QLOgjOi6klWHc865NEkkyUwDOolIRxGpB/QBSqLOKQEuDh73BiaodfaUAH1EpDAYNdYJmBqvzOA1E4MyCMp8NZl1JPZrcc45lwxVXi4L+j+uAcZjlwSfVtU5IjIMKFXVEuAp4HkRKcNaF32C184JRovNBcqBq1V1K0CsMoMqBwKjReROYEZQNkmuwznnXBr4ZMwoNR7C7JxzeSze6LL8mJrunHMuFJ5knHPOpYwnGeeccynjScY551zKeMd/FBFZDSyp4ctbYnN98om/5/yQb+85394v1P49t1fVnWa+epJJIhEpjTW6Ipf5e84P+fae8+39Qures18uc845lzKeZJxzzqWMJ5nkGhV2ACHw95wf8u0959v7hRS9Z++Tcc45lzLeknHOOZcynmScc86ljCeZGhCRniIyX0TKRGRQjOcLRWRM8PwUEemQ/iiTJ4H3e6OIzBWRj0XkHRFJfA/qDFXVe444r7eIqIhk/XDXRN6ziPwy+LeeIyIvpjvGZEvgb3svEZkoIjOCv+/TwogzWUTkaRFZJSKz4zwvIvJQ8Pv4WESOqHWlquq3atywbQMWAnsD9YBZQFHUOVcBI4PHfYAxYced4vd7PNAweHxlNr/fRN9zcF5j4H1s++/isONOw79zJ2z7jebBz7uHHXca3vMo4MrgcRGwOOy4a/mejwGOAGbHef404A1sw82fAVNqW6e3ZKqvK1Cmqp+r6hZgNNAr6pxewLPB43HACSISa5vobFDl+1XViar6Q/Djh9gupNkskX9jgDuAPwKb0hlciiTynn8FPKqq6wBUdVWaY0y2RN6zAk2Cx03J8t11VfV9bD+ueHoBz6n5ENupeM/a1OlJpvraAF9G/Lw0OBbzHFUtBzaQuVtwVyWR9xupP/ZNKJtV+Z5F5HCgnar+K52BpVAi/877AfuJyH9F5EMR6Zm26FIjkfd8G3ChiCwFXgeuTU9ooanu//cqVbkzpttJrBZJ9DjwRM7JFgm/FxG5ECgGjk1pRKlX6XsWkTrACOCSdAWUBon8Oxdgl8yOw1qr/xGRg1R1fYpjS5VE3nNf4C+q+v9E5P+w3XkPUtVtqQ8vFEn/7PKWTPUtBdpF/NyWnZvQ/ztHRAqwZnZlTdRMlsj7RUROBIYAZ6nq5jTFlipVvefGwEHAuyKyGLt2XZLlnf+J/l2/qqo/quoiYD6WdLJVIu+5PzAWQFUnA/WxhSRzVUL/36vDk0z1TQM6iUhHEamHdeyXRJ1TAlwcPO4NTNCgVy0LVfl+g0tHT2AJJtuv00MV71lVN6hqS1XtoKodsH6os1Q1m/ftTuTv+h/YIA9EpCV2+ezztEaZXIm85y+AEwBE5EAsyaxOa5TpVQL0C0aZ/QzYoKpf1aZAv1xWTapaLiLXAOOx0SlPq+ocERkGlKpqCfAU1qwuw1owfcKLuHYSfL/3AY2Al4PxDV+o6lmhBV1LCb7nnJLgex4PnCwic4GtwG9V9evwoq6dBN/zTcCfReQG7LLRJVn8hREReQm73Nky6Gf6A7ALgKqOxPqdTgPKgB+AS2tdZxb/vpxzzmU4v1zmnHMuZTzJOOecSxlPMs4551LGk4xzzrmU8STjnHMuZTzJOOecSxlPMs4551LGk4xzGUBEioJ9aebEeb6piGwTkVxYUcHlEU8yzmWGinXPpsd5/ghs8cKP0hOOc8nhSca5zFCRZOKtf9YluPck47KKJxnnMkNVLZmdkoyI7CsiI0VkpoiUx9tS17kw+QKZzoVMROoCh2KLTs6Ic1pFkolMQp2B04Ep2BdG/9LoMo7/UToXviKgITAvYhvr/xGRJsC+wLpgH5cK/1TVdqraG7+M5jKUJxnnwpdop/8OrZwc3p3R5RBPMs6FryLJzIzz/GnBvbdWXNbxJONc+Cr6WzZGPyEiTYHzgx/jtXScy1ieZJwLkYgUYJ3+AH1FpDDiuT2BMUCb4ND8NIfnXK15knEuXAdh+8YvBboCS0TkdRGZDCwEWgA/Buc+KSL9wgnTuZrxJONcuCr6Y94BzgG+AnoAewFPAScCj2GX0uoA00KI0bka83kyzoWroj9mmqqOB8bHOOf64OZc1vEk41y4qlpOJi4Racj2kWftgSYi0jv4eZqqLklCfM7Viqhq2DE4l5dEZBfgW+wyWBNV3VTN13cAFsV5+lJV/Utt4nMuGbwl41x4DgYKgRnVTTAAqroYm6TpXMbyloxzzrmU8dFlzjnnUsaTjHPOuZTxJOOccy5lPMk455xLGU8yzjnnUsaTjHPOuZTxJOOccy5l/j/1qNe7G2jPmAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_likelihood_1(r, *a, **kw):\n",
    "    prob_roll_1 = linspace(0,1,100)\n",
    "    ps = []\n",
    "\n",
    "    for p1 in prob_roll_1:\n",
    "        # Vary the probability of rolling a 1, and \n",
    "        # divide the remaining probability equally\n",
    "        # among the other 5 numbers\n",
    "        mus = (ones(6) * (1-p1))/5.\n",
    "        mus[0] = p1\n",
    "        ps.append(Die(mus).likelihood(r))\n",
    "    \n",
    "    plot(prob_roll_1, ps, *a, **kw)\n",
    "\n",
    "# Plot some arbitrary data (where \"0\" means rolling a one)\n",
    "plot_likelihood_1(one_of(6,[0,2,1,5,3,1,1]), 'red')\n",
    "plot_likelihood_1(one_of(6,[1,0,1,2,0,1,0,0,0,0]), 'blue')\n",
    "    \n",
    "xlabel(r\"$\\mu_1$\", size=22)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Likelihood')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAETCAYAAAA8rh0/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXSc913v8fd3RvsuWfIiW4rjLM7WxA6iTZue0oWW0JamQIEW2kvbgKGUnrIcuHAu9xYK93Au3HILHC6tS1NKoRuFLrTQwqVp06RNGidx3DSOE8dxLMXWZq3WMlrme/+YGUeRZWskzTPPos/rHB9pNKN5vo+Xj37+Pr/n9zN3R0REkicVdgEiIhIMBbyISEIp4EVEEkoBLyKSUAp4EZGEUsCLiCRU5ALezO40s0Eze7QE7/UKMzu85Nesmb0x/9yrzOyh/NfvMbMrN169iEh0WNTmwZvZy4BzwN+5+w0lfN824Diwy92nzewJ4HZ3P2pmvwK80N3fXqrjiYiELXIjeHe/GxhZ+jUzu8LMvmpmD5rZt8zsmnW89ZuAf3P36cKhgKb8583A6XUXLSISQRVhF1Ckg8Avu/uTZvYi4P8Cr1zje7wZ+LMlj38B+FczmwEmgFtKUqmISEREPuDNrAF4CfCPZlb4cnX+uZ8A3r/Ctz3r7j+y5D12AC8AvrbkNb8OvNbd7zez3yIX/r9Q+jMQEQlH5AOeXBtpzN33LX/C3f8Z+Oci3uOngc+7+zyAmXUAN7n7/fnnPwN8tUT1iohEQuR68Mu5+wTwtJn9FIDl3LTGt3kL8Kklj0eBZjO7Ov/41cDRDRcrIhIhkRvBm9mngJcD7WbWB7wP+Dngr83s94BK4NPAI0W+326gC/hm4WvuvmBmvwj8k5llyQX+O0t3FiIi4Qt0mqSZtQB/A9xAbtbKO939O4EdUEREzgt6BP/nwFfd/U1mVgXUBXw8ERHJC2wEb2ZN5Nooe7zIg7S3t/vu3bsDqUdEJIkefPDBYXfvWOm5IEfwe4Ah4GP5i6IPAu9196mlLzKzA8ABgO7ubg4dOhRgSSIiyWJmz1zsuSBn0VQANwN/7e77gSngd5a/yN0PunuPu/d0dKz4Q0hERNYhyIDvA/qWzDX/HLnAFxGRMggs4N29H+g1s735L70KeCyo44mIyPMFPYvmPcA/5GfQnADeEfDxREQkL9CAd/fDQE+QxxARkZVFfqkCERFZHwW8iEhCKeAT7Fxmgc8e6iWzsBh2KSISgsgtNialMZVZ4B0f+y4PnBxlbiHLW2+5LOySRKTMNIJPoFy4P8BDp8Zob6jms4d6wy5JREKggE+Y2flF3vm3D3DomRE++DP7ePcrruBI3ziPnZ4IuzQRKTMFfMJ85cgZ7n96hD9900382E2dvHHfTqrSKY3iRTYhBXzCPHRqlMbqCn58/04AWuureM312/j8w88yO6+LrSKbiQI+YQ73jnFjVzOp1PkNyvmZH+xifGaef39sIMTKRKTcFPAJMjO3yOP9k+zranne12+9op2dLbV89gG1aUQ2EwV8gjx6epzFrLOvq/V5X0+ljJ/q2cU9x4fpHZkOqToRKTcFfIIcPjUGcMEIHuANN3UC8K0nh8tak4iERwGfIId7x9jZUktHY/UFz+3eUk99VZpj/ZouKbJZKOAT5HDvGPu6Lxy9Q65Ns3d7I0f7J8tclYiERQGfEEOTGZ4dm2H/Cu2Zgmt2NPH4mQmC2mhdRKJFAZ8Qh3sv3n8vuHZ7IxOzC5wZny1XWSISIgV8QhzuHaUiZdyws/mir7lmRxMAj6sPL7IpKOAT4nDvGNfsaKSmMn3R1+zd3gjA0TPqw4tsBgr4BMhmnSO945dszwA01VSys6WWx3WhVWRTUMAnwFND55jMLFxwg9NKrt3RyONn1KIR2QwU8AlwpG8cgH1dF++/F1yzvYkTw1Pa5UlkE1DAJ8DJs1OkU8ZlW+pXfe01OxpZzDrHB8+VoTIRCZMCPgF6R6bZ0VxDZXr1P85rtudn0uhCq0jiKeAToHd0hq7WuqJeu3tLHdUVKU2VFNkEFPAJ0DsyTVdbbVGvrUinuHpbo2bSiGwCFUG+uZmdBCaBRWDB3XuCPN5mNDu/yOBkpugRPMA12xu569hQgFWJSBSUYwT/Cnffp3APRt9obn33rrbiA37v9kaGz2UYmswEVZaIRIBaNDHXOzIDUHSLBuDa/JIFx9SmEUm0oAPegX83swfN7MBKLzCzA2Z2yMwODQ2pbbBWvYUR/BpaNFdvyy1Z8MSAAl4kyYIO+Fvd/WbgR4F3m9nLlr/A3Q+6e4+793R0dARcTvL0jkxTXZFacZOPi2lvqKKuKn3+h4OIJFOgAe/up/MfB4HPAy8M8nibUe/IDF1tdZhZ0d9jZnS31Wl/VpGECyzgzazezBoLnwOvAR4N6nib1amRabpai++/F3S11XFKAS+SaEGO4LcB95jZI8B3ga+4+1cDPN6m1Ds6vaYZNAVdrXX0jsxodyeRBAtsHry7nwBuCur9Bcan55mcXVjTBdaC7rZaZuYXGT43t6b+vYjEh6ZJxtj5GTRrmCJZ0L0l90NBbRqR5FLAx1jhIumudYzgC6P+Ps2kEUksBXyM9a7jLtaCwg+FU2cV8CJJpYCPsVMj0zTVVNBcW7nm762tSrO1sVotGpEEU8DHWO/IzPle+np0tdXpZieRBFPAx1jv6PS6ZtAU5G52milhRSISJQr4mMpmnb7RmXX13wu62uo4PT7D3EK2hJWJSFQo4GNq6FyGuYXsuu5iLehqrcUdTo9pFC+SRAr4mCpcHN21gRF8d5vmwoskmQI+pgpz4DfUg9fNTiKJpoCPqWdHc22VXRto0WxrrKEqndJMGpGEUsDHVP/ELK11ldRUptf9HqmUsau1VssGiySUAj6mBiYybGuq2fD7aNlgkeRSwMfUwMRsSQJec+FFkksBH1O5gN/4Mr9dbbWMz8wzPj1fgqpEJEoU8DG0sJhl+FyG7SUawQO60CqSQAr4GBo+N0fWYWuJevCALrSKJJACPoYGJmYBSnaRFTQXXiSJFPAx1J8P+FK0aJpqKmmsqdByBSIJpICPocHzI/jS7KW6s6WWZ8dmS/JeIhIdCvgYGpjIkE4ZWxpKE/CdLbUawYskkAI+hvonZuloqCadspK8X2dLDafHFfAiSaOAj6GBiVm2NW+8/17Q2VLL2PQ8U5mFkr2niIRPAR9DAxOzbGssTXsGcj14gDMaxYskigI+hkq1Dk1BZz7gdaFVJFkU8DEzO7/I+Mw820vcogE4owutIokSeMCbWdrMHjazLwd9rM2gcJPT1hK2aLY1VpMybd0nkjTlGMG/FzhahuNsCgMTGYCSjuAr0im2NdWoRSOSMIEGvJntAl4H/E2Qx9lM+ku4TMFSmgsvkjxBj+A/CPw2kL3YC8zsgJkdMrNDQ0NDAZcTf+fvYm0MIOA1i0YkUQILeDN7PTDo7g9e6nXuftDde9y9p6OjI6hyEmNgYpaayhRNtRUlfd/OlhrOjM2SzXpJ31dEwhPkCP5W4A1mdhL4NPBKM/v7AI+3KfTnp0ialeYu1oKdLbXMLWYZnsqU9H1FJDyBBby7/66773L33cCbga+7+1uDOt5mUaqt+pbrbM5NlTytC60iiaF58DETWMC3FAJefXiRpChtI/ci3P0bwDfKcawkc3cGJmbZXqJlgpfqbMn90FDAiySHRvAxMjG7wOx8NpARfHNtJXVVabVoRBJEAR8j5+9iDSDgzUxz4UUSRgEfIwMl3KpvJZoLL5IsCvgY6R8v7VZ9y+1sqdEIXiRBFPAxMjiZm6MeRA8eclMlh8/NMTu/GMj7i0h5KeBjZGBilqaaCmoq04G8//llg8d1oVUkCRTwMTI4kQnkAmvBDk2VFEkUBXyMDE7OlnQd+OV2nt/ZSQEvkgQK+BgZnMwEGvCFNeY1ghdJBgV8TLh7LuADbNFUV6TpaKzmjG52EkkEBXxMTMwsMLeQDXQED5oLL5Ikl1yLxsy+B1x0gXB3v7HkFcmKBidzo+qOgAN+Z0sNj/dPBnoMESmP1RYbe33+47vzHz+R//hzwHQgFcmKhvJz4LeWeCen5Tqba/n644O4e8nXnBeR8rpkwLv7MwBmdqu737rkqd8xs3uB9wdZnDyncJPT1oDuYi3obKlldj7L6PQ8bfVVgR5LRIJVbA++3sxeWnhgZi8B6oMpSVZSaNGUowcPmkkjkgTFrgd/B3CnmTXnH48B7wymJFnJ4ESG2so0DdXBLuG/dC78DTubV3m1iERZUWmR3zj7JjNrAszdx4MtS5YbnMzQ0VgdeF9cG3+IJEdRLRozazazPwO+DvynmX1gyWheyiDou1gL2uqrqK5IKeBFEqDYHvydwCTw0/lfE8DHgipKLpS7ySn4gDczdrbUamcnkQQotqF7hbv/5JLHf2Bmh4MoSFY2NJHhZVd1lOVYnS21Wo9GJAGKHcHPLJtFcyugBCiTmblFJjMLgd/kVNCpjT9EEqHYEfy7gI/n++4GjAA/H1hV8jzP3eRUroCvZXAyQ2ZhkeqKYNaeF5HgFTuL5jDPzaLB3ScCrUqe5/wc+AAXGluqMBd+YDxD95a6shxTREpvPbNovq5ZNOU1WOYRvNaFF0kGzaKJgcGJ8tzFWqC7WUWSIbBZNGZWA9wNVOeP8zl3f9/6ytzcBiczVKSM1rryrA2zQxt/iCRCkLNoMsAr3f0mYB9wm5ndsr4yN7fByQztDdWkUuVZ3bGmMk17Q5XWhReJuWJH8L8M/N2yWTRvv9Q3uLsD5/IPK/O/Lrq2vFxcuW5yWio3F143O4nEWbGzaB5hHbNozCwNPAhcCfyVu9+/wmsOAAcAuru7iyx7cxmcmGVXa21Zj9nZXMvxoXOrv1BEIquogDezauAngd1ARWHBK3e/5Hrw7r4I7DOzFuDzZnaDuz+67DUHgYMAPT09GuGvYGgyw/7u1rIes7OllrufHNLGHyIxVmwP/ovA7cACMLXkV1HcfQz4BnDbGuvb9OYXs5ydmivbDJqCzpYapucWGZ+ZL+txRaR0iu3B73L3NYWzmXUA8+4+Zma1wA8D/2utBW52Z8/NAcHv5LTc0rnwLWWavSMipVXsCP7bZvaCNb73DuAuMzsCPAD8h7t/eY3vsek9t5NTee5iLXhuLrwutIrE1SVH8Gb2PXIzXyqAd5jZCXLTH43cRJkbL/a97n4E2F/CWjelwYny3sVaoJudROJvtRbN68tShVxUYZmCcq0kWbClvooqbfwhEmurBfyou0+YWVtZqpELFFo07Q3lDfhUyuhsrqFPAS8SW6sF/CfJjeIfJNeqWTpfzoE9AdUleYOTGdryo+ly62qro29kuuzHFZHSuGTAu/vr8x8vL085stzgRHn2Yl3JrtY6/v10fyjHFpGNW+0i682Xet7dHyptObJc/8Qs25vLO4OmYFdrLWen5pjKLFBfXeyMWhGJitX+1X7gEs858MoS1iIr6B/PcENnOEvvd7XlNvvoG51h7/bGUGoQkfVbrUXzinIVIhfK3cWaYVuZdnJariu//k3vyLQCXiSGit3Rqc7Mfs/MDuYfX2VmmkIZsKHJDO6EF/D5EXzvqC60isRRsVMzPgbMAS/JP+4D/iiQiuS8/vxOTtubw7nIuqW+itrKNH2jmiopEkfFBvwV7v4nwDyAu8/w/CmTEoCB8VzAhzWCNzN2tdbSq6mSIrFUbMDP5RcMcwAzu4LckgUSoPMj+JACHnJtml6N4EViqdiAfx/wVaDLzP4B+E/gtwOrSgAYmMhQmTba6sNbzbGrtZa+kWlyG3SJSJwUu6PTf5jZQ8At5Foz73X34UArEwYmZtnaWBPqhhtdbXVMZhYYn5nXssEiMVPsLJr3u/tZd/9KfsnfkfxIXgLUPx7eTU4Fu1qfmwsvIvFSbIum28x+F85v3/cF4MnAqhIgN4IPs/8O0NX23Fx4EYmXYgP+HcAL8iH/L8Bd7v77gVUluDv9E7OhzaApKIzgNRdeJH7WshbNnwMfBu4FvmlmN2stmuCcyywwPbcY2hz4gubaSppqKugdUYtGJG7WuhbNKHBd/utaiyZAAxPhzoFfKjdVUiN4kbjRWjQR1T+eu80gEgHfWseTg5NhlyEia7Rai+at7v73ZvYbKz3v7n8WTFkShZucCrraarnr2CDuHuqUTRFZm9VaNPX5jystJag7XwI0cH4dmvADfldrHZmFLEPnMmxtDL8eESnOai2aD+c//sHy58zs14IqSnIB31xbSU1lOuxSlkyVnFHAi8TIRjb6XLFtI6XRPz7LtqZwZ9AUdJ2/2UkXWkXiZCMBr2ZsgAYiMAe+4PxceN3sJBIrGwl49eAD1B+Bu1gLaqvStDdUa7kCkZhZbRbNJCsHuQG1gVQkLGadoclMJC6wFnS31XLy7FTYZYjIGlxyBO/uje7etMKvRndf7YdDl5ndZWZHzez7Zvbe0paeXMPnMmRD3KpvJZe3N/D0sAJeJE420qJZzQLwm+5+Lbllht9tZtcFeLzE6A95J6eV7OmoZ2Aiw7nMQtiliEiRAgt4dz9TWKvG3SeBo8DOoI6XJFG6yalgT3vuloiTGsWLxEaQI/jzzGw3sB+4f4XnDpjZITM7NDQ0VI5yIm+wsA5NyAuNLbWnowGAp4bOhVyJiBQr8IA3swbgn4Bfc/eJ5c+7+0F373H3no6OjqDLiYX+iVkqUkZ7fXQC/rItdZjBiSGN4EXiItCAN7NKcuH+D+7+z0EeK0n6xzNsbawmlYrOrQY1lWl2tdZyQi0akdgILOAttyrVR4GjWpRsbc6Mz7AtQlMkCy5vb+CEWjQisRHkCP5W4G3AK83scP7XawM8XmL0jc6cXx4gSva01/P08BTuusdNJA5WW01y3dz9HrScwZotZp3TYzO8/sYdYZdygSs66pmeW2RgIlo3YYnIysoyi0aK1z8xy0LW6WqL4Ag+P5NGbRqReFDAR0xffkGvXa3RWwni8vxc+Kd0oVUkFhTwEdObX9Arij347U011FameVpTJUViQQEfMX2j05jBjpbo9bhTKePy9npODKtFIxIHCviI6R2ZYVtjDdUV4e/ktJLLO+p1s5NITCjgI6ZvdPr8FnlRdEV7PX2j02QWFsMuRURWoYCPmL7RmfM7KEXRno4Gsg6nzmp3J5GoU8BHyPxiljPjM3RFcAZNwZ6O/EwatWlEIk8BHyH947NknUiP4AtTJXWhVST6FPAR0hvhOfAFjTWVdDRW60KrSAwo4COksKl1FO9iXWpPe73WhReJAQV8hPSOTpMyIr/Oy7U7mjjWP8liVouOiUSZAj5C+kZn2NFcS2U62n8s13U2MT23yMmzatOIRFm0k2ST6RudjnT/veD6ziYAvn/6gg26RCRCFPAR0jsS7TnwBVdtbaQybTymgBeJNAV8RGQWFhmYnI30XawFVRUprtrayPdPj4ddiohcggI+Ik6PzeIRnwO/1PWdTTx2ekK7O4lEmAI+IvpGc3Pgo3wX61LXdTZxdmqOwclM2KWIyEUo4COiMAd+V8TnwBdc39kMoDaNSIQp4COid2SaipSxvSnac+ALrt3RCKALrSIRpoCPiL7RGTpbakmn4rFPeWNNJZdtqdNUSZEIU8BHxDMj0V4HfiXXdzbx2BkFvEhUKeAjIJt1jg9MctXWxrBLWZPrdjTxzNlpJmbnwy5FRFaggI+AZ8dmmJpbZO/2eAV84ULrUbVpRCJJAR8Bx/onAbh6W9wCPrdkgdo0ItEUWMCb2Z1mNmhmjwZ1jKQ4NlAI+IaQK1mbjsZq2huqdKFVJKKCHMH/LXBbgO+fGE8MTLKzpZbGmsqwS1kTM+O6zmYefVZz4UWiKLCAd/e7gZGg3j9JjvVPxq7/XvAD3a0cG5hkfFoXWkWiJvQevJkdMLNDZnZoaGgo7HLKbn4xy4mhqdj13wtu2dOGO9z/9NmwSxGRZUIPeHc/6O497t7T0dERdjlld3J4irnFLHu3x6v/XrCvu4XqihTfOaGAF4ma0AN+s3vuAms8R/DVFWl6drdy3wl140SiRgEfsif6J0mnjCs64jmCB7jl8i0cPTPB6NRc2KWIyBJBTpP8FPAdYK+Z9ZnZHUEdK86ODUyye0sdNZXpsEtZtxdfsQVQH14kaiqCemN3f0tQ750kTwycO78yY1zduKuF2so0950Y4bYbdoRdjojkqUUTopm5RU6eje8MmoKqihQ9u1v5zlMawYtEiQI+RMcHz+EOe2Me8AC37NnCsYFJzp7TDk8iUaGAD9H5GTQxvclpqef68JpNIxIVCvgQPTEwSVVFit1b6sMuZcNesLOZuqq02jQiEaKAD9Gx/kmu2toQm12cLqUyneIHd7fx7aeGwy5FRPIU8CHJZp3vPTvOtTuawi6lZF6xt4OnhqZ4It96EpFwKeBDcrR/gpGpOV6S710nwetu7CRl8MXDz4ZdioiggA/NvcdzrYxbr2wPuZLS6Wis5tYr2/ni4dO4e9jliGx6CviQ3HP8LFdtbWBbU03YpZTUG/ftpG90hodOjYZdisimp4APQWZhke8+fTZRo/eC11y/jeqKFF94+HTYpYhsegr4EDx8aozZ+WwiA76xppIfvm4bX/neGeYXs2GXI7KpKeBDcO/xYdIp40V72sIuJRC339TJyNQc9zypKZMiYVLAh+Ce48PctKuZppjtwVqsl+/dSnNtpWbTiIRMAV9mE7PzPNI7xksT2J4pqKpI8doX7OBr3x9gRGvEi4RGAV9m9z11lqzDSxIc8AB3vHQ3mYVFPvTNp8IuRWTTUsCX2b3Hh6mtTLO/uyXsUgJ15dZG3rh/Jx//9kkGJmbDLkdkU1LAl5G7c/eTw7zw8jaqK+K7g1Oxfu1VV7OYdf7y60+GXYrIpqSAL6NvPTnM08NTvO4Fm2PXo+4tdfzMD3bx6e/20jsyHXY5IpuOAr6MPvTNp9jWVM3t+zvDLqVs3vPKq0injA/+P43iRcpNAV8mR/rG+PZTZ3nnrZdvivZMwfbmGt52y2V8/uG+8+vviEh5KODL5MPfPEFjdQU/+6LusEspu19/9dVcubWBX/3kQ2rViJSRAr4MTg5P8W+PnuHnbrmMxoTe3HQp9dUVHHxbD4tZ55c+8SAzc4thlySyKSjgy+Aj3zpBRSrFO2/dHXYpodndXs9fvGU/R/sn+O1/OkI2q+WERYKmgA/YPU8O89lDvfzkD+xka8KWBl6rl+/dym/9yF7+5ZHT3PHxBxib1l2uIkFSwAfo4VOjHPjEIfa0N/Bfb7sm7HIi4V0/dAV/ePv13HN8mNf/5T08+ux42CWJJJYCPiDH+id5+8ceoKOxmk/c8UJa6qrCLikSzIy3vXg3n/mlF7OYdX7ir7/Nf//Co5w6q4uvIqUWaMCb2W1mdszMjpvZ7wR5rKiYnJ3nI3ef4Gc/ch81lSn+/o4XbfrWzEpu7m7ly+95KT++byefeaCXl//vu3j3Jx/iK0fOaIEykRKxoPbONLM08ATwaqAPeAB4i7s/drHv6enp8UOHDgVST1Bm5xc5MTTF8aFzHD41xj8e6mUys8Ate9r4ozfewJVbG8MuMfIGJma5896n+dT9p5iYXcAMrt3exDU7GtnTXs/u9nq2NtbQWldJa30VDdUVVFekMLOwSxcJnZk96O49Kz4XYMC/GPh9d/+R/OPfBXD3P77Y96w34H/sL+9hdr40U++W/m4s/b1xIJt1sg6LWWd6boGpuUXmFp7btSidMn70hu0ceNkebtyV7MXEgrCwmOXIs+Pc++Qw9z19luOD5xiYyFz09bWVaSrTRmU6RUXaqEilMIOUGal89psZBmBQ+HFQzA8G/eiQcmqtq+Kzv/zidX3vpQK+YkNVXdpOoHfJ4z7gRctfZGYHgAMA3d3ruwnoio565kq4PZwt/ee95NO0GemUkTKjvjpNXVUFDdVpLttSz5VbG7i8vZ6ays1zl2qpVaRT3Nzdys3drbyHqwCYnlvgmbPTnD03x8j0HGPTc0xlFpmZX2RmboH5RWd+McvCorPoTtadbNZxwJ38R3/uB3cR4xkv5kUiJRTU5j9BBvxKg6AL/uW4+0HgIORG8Os50AffvH893yYxUFdVwbU7msIuQySWgrzI2gd0LXm8Czgd4PFERGSJIAP+AeAqM7vczKqANwNfCvB4IiKyRGAtGndfMLNfBb4GpIE73f37QR1PRESeL8gePO7+r8C/BnkMERFZme5kFRFJKAW8iEhCKeBFRBJKAS8iklCBLVWwHmY2BDyzzm9vBzbbpp865+TbbOcLOue1uszdO1Z6IlIBvxFmduhi6zEklc45+Tbb+YLOuZTUohERSSgFvIhIQiUp4A+GXUAIdM7Jt9nOF3TOJZOYHryIiDxfkkbwIiKyhAJeRCShYhfwq23kbWbVZvaZ/PP3m9nu8ldZOkWc72+Y2WNmdsTM/tPMLgujzlIqdrN2M3uTmbmZxX5KXTHnbGY/nf+z/r6ZfbLcNZZaEX+3u83sLjN7OP/3+7Vh1FkqZnanmQ2a2aMXed7M7C/yvx9HzOzmDR/U3WPzi9yyw08Be4Aq4BHgumWv+RXgQ/nP3wx8Juy6Az7fVwB1+c/fFefzLfac869rBO4G7gN6wq67DH/OVwEPA635x1vDrrsM53wQeFf+8+uAk2HXvcFzfhlwM/DoRZ5/LfBv5HbDuwW4f6PHjNsI/oXAcXc/4e5zwKeB25e95nbg4/nPPwe8yorZZTmaVj1fd7/L3afzD+8jt3NWnBXzZwzwh8CfALPlLC4gxZzzLwJ/5e6jAO4+WOYaS62Yc3agsF9jMzHfEc7d7wZGLvGS24G/85z7gBYz27GRY8Yt4FfayHvnxV7j7gvAOLClLNWVXjHnu9Qd5EYAcbbqOZvZfqDL3b9czsICVMyf89XA1WZ2r5ndZ2a3la26YBRzzr8PvNXM+sjtK/Ge8pQWmrX+e19VoBt+BKCYjbyL2uw7Joo+FzN7K9AD/FCgFQXvkudsZing/wBvL1dBZVDMn3MFuTbNy8n9L+1bZnaDu48FXFtQijnntwB/6+4fMLMXA5/In3M2+PJCUfLsitsIvpiNvM+/xswqyP3X7lL/LYqyojYuN6l4oU4AAAI5SURBVLMfBv4b8AZ3z5SptqCsds6NwA3AN8zsJLle5ZdifqG12L/XX3T3eXd/GjhGLvDjqphzvgP4LIC7fweoIbcoV1IV9e99LeIW8MVs5P0l4Ofzn78J+Lrnr2DE0Krnm29XfJhcuMe9LwurnLO7j7t7u7vvdvfd5K47vMHdD4VTbkkU8/f6C+QuqGNm7eRaNifKWmVpFXPOp4BXAZjZteQCfqisVZbXl4D/kp9Ncwsw7u5nNvKGsWrR+EU28jaz9wOH3P1LwEfJ/VfuOLmR+5vDq3hjijzfPwUagH/MX0s+5e5vCK3oDSrynBOlyHP+GvAaM3sMWAR+y93Phlf1xhR5zr8JfMTMfp1cq+LtMR6sYWafItdia89fV3gfUAng7h8id53htcBxYBp4x4aPGePfLxERuYS4tWhERKRICngRkYRSwIuIJJQCXkQkoRTwIiIJpYAXEUkoBbyISEIp4EVWYGbfMLO9+c+3XGwNb5EoU8CLrOxK4Mn85zcC3wuxFpF1UcCLLJPfFevZJasW3ggcMbM9ZvZRM/tciOWJFE0BL3KhfcCRJY9/ADiS35zijpBqElkzBbzIhW4it3IhZnYVuZ121KKR2FHAi1xoH5Ays0eA/wEc5bklqEViQ6tJiiyTX2p6v7tPLvv6FuB/Aq8G/sbd/ziM+kSKpYAXWcLMGoEH3f3qsGsR2SgFvIhIQqkHLyKSUAp4EZGEUsCLiCSUAl5EJKEU8CIiCaWAFxFJKAW8iEhC/X/UVRE6mCRRfAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Using some samples from our loaded die\n",
    "plot_likelihood_1(loadeddie.roll(100))\n",
    "xlabel(r\"$\\mu_1$\")\n",
    "ylabel('Likelihood')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also write the likelihood as\n",
    "\n",
    "$p(D|\\boldsymbol{\\mu}) = \\prod_{k=1}^K \\mu_k^{m_k}$\n",
    "\n",
    "Where $m_k$ is the number of data points where $x_k = 1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now to find the maximum likelihood solution for $\\boldsymbol{\\mu}$, we maximise $\\log p(D|\\boldsymbol{\\mu})$ with the constraint that $\\sum_{k=1}^{K}\\mu_k = 1$, which we can enforce using a Lagrange multiplier $\\lambda$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\frac{\\partial}{\\partial \\mu_k}[\\sum_{j=1}^{K}m_j \\log \\mu_j + \\lambda(\\sum_{j=1}^{K}\\mu_j - 1)] = 0$\n",
    "\n",
    "$\\frac{m_k}{\\mu_k} + \\lambda(1 - 1) = 0$ (The partial derivative of terms in the summations where $j \\neq k$ vanish)\n",
    "\n",
    "$\\mu_k = -\\frac{m_k}{\\lambda}$\n",
    "\n",
    "Sub into $\\sum_{k=1}^{K}\\mu_k = 1$\n",
    "\n",
    "$-\\sum_{k=1}^{K} m_k = \\lambda$\n",
    "\n",
    "so $\\lambda = -N$ (because the summation over the number of data points in each state gives the total number of data points)\n",
    "\n",
    "And $\\mu_k^{ML} = \\frac{m_k}{N}$ (which is just the fraction of the data points where $x_k=1$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now consider the distribution over the $m_k$s, which would give the number of ones, twos, threes etc. rolled, the same way the Binomial gave us the number of heads. It is called the Multinomial:\n",
    "\n",
    "$\\operatorname{Mult}(m_1,m_2,...,m_K|\\boldsymbol{\\mu},N) = \\binom{N}{m_1 m_2 ... m_K} \\prod_{k=1}^K \\mu_k^{m_k}$    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import factorial\n",
    "\n",
    "def multinomialpmf(ms, mu):\n",
    "    # ms: the number of observations of each state\n",
    "    ms = array(ms, float)\n",
    "    # We normalise our weights mu, so the input can be unnormalised\n",
    "    mu = array(mu, float)/sum(mu)\n",
    "    # The number of observations can be inferred from ms\n",
    "    N = sum(ms)\n",
    "    norm = factorial(N)/prod(factorial(ms))\n",
    "    return norm * prod(mu**ms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We cannot really plot the multinomial directly, because of the many dimensions of $\\boldsymbol{m}$.\n",
    "\n",
    "But we can show the distribution of the number of ones rolled (constraining the other numbers), and with various weighted dice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Probability of rolling $m_1$ ones')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEhCAYAAACDefxEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOydd3hUVfr4P28mlRBaQErondBCAigoSFFapIgouO6CbdW1set+Ley6Cqg/1tVdxbY2BGRVQECkKQrYEBQCRAihGANIQockQApJJuf3x72JSZhMZpKpyfk8zzy5c+4573nvncm895zznvcVpRQajUaj0ThCgLcV0Gg0Go3/oI2GRqPRaBxGGw2NRqPROIw2GhqNRqNxGG00NBqNRuMw2mhoNBqNxmG00dBoNBqNw2ijodFoNBqH0UZD41FE5LCIXOeOtiKyV0SGlK9butzdiEgXEdklIhdE5GEXyCt9HWWu35PXpdEUE+htBTS+j4gcBpoCViAbWAc8pJS66E29yqOU6u5IuXk9dyulNrhBjceAr5VSfdwguwwVXa9G4070SEPjKGOVUnWBWKAf8GT5CiKiH0KgDbDXkYr6fmn8EW00NE6hlEoHPgN6QMmUyeMishvIFpFAEekmIl+LSKY5hTKunJh+IpIsIhkiMl9EQotPiMgTIvKLOb2TLCI3OtHW5vRVuSmeRUBrYLWIXDR1X16u/qsi8rKt67d3bSKyCRgKvGbK7lyBLs7eL5vYmLr6PxHZLSJZIrKk+N6ISGypKbOPzXPP2pF7t4isF5H/mvf5oIhEi8h0EflVRM6IyERHdHTwvlWou3m+hYgsF5HTInKo9LSfeS/TzWs7ICLDndFLUwWUUvqlX3ZfwGHgOvO4FcaT9DOlziWa5WFAEJAC/A0IBoYBF4AupeonmfUbAd8Dz5bq62agBcYDzWSM6bDmDrYtrafNYxvnmpt9NDDfBwKngDgb98HutZl1vsaY+rJ3L529X45cx2Fgm3nvGgH7gPtMmUeA6WZfE4H80vfNho6vAGeB4YAFWF5OxnRgpxPfH0eu8TLdzXMBwA7gKbNteyAVGAl0AY4CLcy6bYEO3v5/qekvPdLQOMpKEckENgPfAP+v1LlXlFJHlVK5wFVAXeCfSql8pdQmYA1wa6n6r5n1zwHPlT6nlPpYKXVMKVWklFoC/Az0d6RtVVBKHQe+xTBWAKOAM0qpHTaqO3JtjuDs/XJG7jHz3qwGYkz5gea5AqXUCowfaHv0BuYopTYqpaxAMrBHKTVXKVWAYbgDAUSkvohsM0dWPSqQ58g12tIdjKnQJkqp2WbbVOAdYArGGlsIEC0iQUqpw0qpXxy+W5oqoY2GxlEmKKUaKKXaKKXuN3/wijla6rgFcFQpVVSq7AgQVUH9I2YbAERkqogkmtMYmRjTYI0daVsNFgK/N49/DyyqoJ4j1+YIzt4vRzlR6jgH44e6BZCulCqdA+Eo9umF8aNeTLSN9/tL9RMPLLMjz5FrtKU7GGtELYq/D+Z34m9AU6VUCvBnYCZwSkQWi4grvg8aO2ijoXEFpX+QjgGtRKT0d6s1kF7qfaty544BiEgbjKfIB4FIpVQDjKdaqaxtNfQFWAn0Mp+UbwA+qKCdI9fmbP+uklkRx4EoEanoHpbB/AyCgYOlimMwptSK6VX83hy9nK5Eh+pc41HgkPnAUvyKUEqNMfv/UCl1DYZxUcDzDsjUVANtNDSu5keMNYLHRCTI3EcwFlhcqs4DItJSRBphPDUuMcvDMf7xTwOIyB2YC+4OtHWGkxhz4wAopfIwnpQ/BLYppX6txrU5iztklmYrxjTOg+ai+3jKTveVpzfGVFQRgIjUw/hB3l2uzk9O6FCda9wGnDcXvMNExCIiPUSknxh7YoaJSAiQB+Sa16pxI9poaFyKUiofGAeMBs4AbwBTlVL7S1X7EPgCY0EzFXjWbJsM/Bvjh+4k0BNjsZvK2jrJHOBJc7rj/8yyhWZ/FU1NOXptTuEOmTbkTwTuAjIxpt/WAJcqaNKby0cVKUqpHABztNCjXB1HdKjSNZprKmMxRjuHzPbvAvUx1jP+aZadAK7AeJDQuBEpO9Wp0dRORKQ1xjx9M6XUeW/r405E5EfgTaXUfBfKXAC8qJRKcpVMjW+iRxqaWo/59PwIsLgmGgwRuVZEmpnTU9MwRg+fu1D+OmAE8I6I3O4quRrfRO9I1dRqRCQcYyrsCIa7bU2kC7AUwyPpF2CS6WrsEooXpTW1Az09pdFoNBqH0dNTGo1Go3EYbTQ0Go1G4zA1ek2jcePGqm3btt5WQ6PRaPyKHTt2nFFKNbF1rkYbjbZt25KQkOBtNTQajcavEJEjFZ3T01MajUajcRhtNDQajUbjMNpoaDQajcZhavSahi0KCgpIS0sjLy/P26po3EhoaCgtW7YkKCjI26poNDWKWmc00tLSiIiIoG3btpSNFq2pKSilOHv2LGlpabRr187b6mg0NYpaNz2Vl5dHZGSkNhg1GBEhMjJSjybdze6l8FIPmNnA+Lt7qbc10niAWjfSALTBqAXoz9jN7F4Kqx+GAjOBY9ZR4z1Ar1u8p5fG7dS6kYYv0LZtWw4fPsyQIUNcLnvIkCEcPnwYZzc1fvzxx3Tr1o2hQ4e6TJcxY8aQmZlpt86QIUNs7qVJTExk3bp1NtucPXuWoUOHUrduXR588EGX6Kpxko2zfzMYxRTkGuWaGk2tHGk4w8pd6byw/gDHMnNp0SCMR0d2YUKfqqRv9m3mzZvHG2+84VKjUdGPviMkJiaSkJDAmDGXB1ANDQ3lmWeeISkpiaQknb7BK2SlOVeuqTHokYYdVu5KZ8aKPaRn5qKA9MxcZqzYw8pd1Uvf3KRJEywWC40aNQKMdZY77riDnj170qdPH7766isAFixYwMSJExk1ahSdOnXiscceK5HxxRdfMGDAAGJjY7n55pu5ePEiAI0aNcJisdCkic0IAHz00Uf07NmTHj168PjjjwMwe/ZsNm/ezH333cejjz5apv7999/PqlWrALjxxhu58847AcPIPPnkkwD873//o3///sTExHDvvfditRoZN9u2bcuZM2cAeOaZZ+jatSvXX389t956Ky+++GJJHx9//DH9+/enc+fOfPfdd+Tn5/PUU0+xZMkSYmJiWLKkbEbX8PBwrrnmGkJDQ5299RpXUb+l7fI6kZ7VQ+NxHDIaIjJKRA6ISIqIPGHjfIiILDHP/ygibUudm2GWHxCRkZXJFJF5IvKTiOwWkWUiUreyPqrD5Le2XvZatPUwAP/6fD+5BWVTDucWWJm5ei8A57LzL2vrCNu3b6dVq1asWLECgNdffx2APXv28NFHHzFt2rSSRdzExESWLFnCnj17WLJkCUePHuXMmTM8++yzbNiwgZ07d9K3b1/+85//ALBixQpatWrF9u3bL+v32LFjPP7442zatInExES2b9/OypUreeqpp+jbty8ffPABL7zwQpk2gwcP5rvvvgMgPT2d5ORkADZv3sygQYPYt28fS5Ys4fvvvycxMRGLxcIHH3xQRkZCQgLLly9n165drFix4rLpqMLCQrZt28bLL7/MrFmzCA4OZvbs2UyePJnExEQmT57s0H3VeJDYP1xeFhAIOWdg80ugUy7UWCo1GiJiAV7HyO8bDdwqItHlqt0FZCilOgIvAc+bbaOBKUB3jAQ3b5iJ4e3J/ItSqrdSqhfwK/CgvT7cyfEs2943mTkFLu1n8+bN/OEPxj9h165dadOmDQcPHgRg+PDh1K9fn9DQUKKjozly5Ag//PADycnJXH311cTExLBw4UKOHKkwVEwJ27dvZ8iQITRp0oTAwEBuu+02vv32W7ttBg0axHfffUdycjLR0dE0bdqU48ePs3XrVgYOHMjGjRvZsWMH/fr1IyYmho0bN5KamnrZ9Y0fP56wsDAiIiIYO3ZsmfMTJ04EIC4ujsOHDzt62zTepCAPQupBvShAoH4rGPcq9LgJNsyEdY9CkbUyKRo/xJE1jf4YieVTAURkMTAeSC5VZzww0zxeBrwmhvvKeIwUmpeAQyKSYsqjIpnF6TbN9mGAsteHqmYWqSX3DqjwXIsGYaRn5l5WHtUgDIBG4cF22zuKvUsICQkpObZYLBQWFqKU4vrrr+ejjz5yWT8VERUVRUZGBp9//jmDBw/m3LlzLF26lLp16xIREYFSimnTpjFnzpwq91t8jcXXp/EDrnsarrwPIpqWLe81xTAkW16BC8fhpnchKMw7OmrcgiPTU1HA0VLv08wym3WUUoVAFhBpp61dmSIyHzgBdAVeraQPt/HoyC6EBVnKlIUFWXh0ZBeX9jN48OCSKZ2DBw/y66+/0qVLxX1cddVVfP/996SkpACQk5NTMjKxx5VXXsk333zDmTNnsFqtfPTRR1x77bWVthswYAAvv/wygwcPZtCgQbz44osMGjQIMEZCy5Yt49SpUwCcO3fuslHPNddcw+rVq8nLy+PixYusXbu20j4jIiK4cOFCpfU0HubCSThjfO8uMxgAAQEw4hkY/S84lgi5GZ7VT+N2HDEathzeyz86VlTH2XLjQKk7gBbAPqB4QtsRPRCRe0QkQUQSTp8+baOJ40zoE8WciT2JahCGYIww5kzs6XLvqfvvvx+r1UrPnj2ZPHkyCxYsKDPCKE+TJk1YsGABt956K7169eKqq65i//79lfbTvHlz5syZw9ChQ+nduzexsbGMHz++0naDBg2isLCQjh07Ehsby7lz50qMRnR0NM8++ywjRoygV69eXH/99Rw/Xjb9dL9+/Rg3bhy9e/dm4sSJ9O3bl/r169vtc+jQoSQnJ9tcCAdjkf2RRx5hwYIFtGzZsmStReNmPn8C3hkGeeft17vyXnhwG9RrYUxTnXdZSnKNt1FK2X0BA4D1pd7PAGaUq7MeGGAeBwJnMH7ky9QtrueITLP8WmCNvT7s6R4XF6fKk5ycfFmZxv1cuHBBKaVUdna2iouLUzt27HB7n/qzdjE/b1Dq6XpKffVP59p9OVOpf3VUKn2ne/TSuBwgQVXwu+rISGM70ElE2olIMMbC9qpydVYB08zjScAms+NVwBTT86kd0AnYVpFMMegIJWsaY4H9lfSh8QPuueceYmJiiI2N5aabbiI2NtbbKmmcoSAX1v4VIjvCNX92rm3vKRAYCvPj4ecv3aOfxmNUuhCulCoUkQcxnvQtwHtKqb0iMhvDGq0C5gGLzIXucxhGALPeUoxF80LgAaWUFaACmQHAQhGphzFS+Qn4k6mKzT40/sGHH37obRU01WHzS5BxCKZ+CoEVT53apEkXuPtL+GASfDgZxs617bKr8QukJj+s9+3bV5XfE7Bv3z66devmJY00nkR/1i7ky6fhwgmY+FbVZeSdh6VT4det8NBOqF/zIivUFERkh1Kqr61zOoyIRqOpnOtnVX/DXmg9uO1jOL77N4OhFOjgkn6FDiOi0Wgq5uB6OLzZOHbFj7slCFrGGce7/mdMV+VnV1+uxmNoo6HRaGyTmwmfPghfPuWesCDWAkj5EhbEw8VTrpevcQvaaHgBHRr9N6oSGv3LL78kLi6Onj17EhcXx6ZNm1yir6Ycm54xYknF/8c9U0h974ApH8Kp/TDv+t82DWp8Gm00KqOWZCcrDo1eHGHXFaxbt44GDRpUqa09o9G4cWNWr17Nnj17WLhwYUncLo0LSdsB2+dB/3ugRYz7+ukyGm5fC5cuwHsjYPt7teL/zZ/RRsMexdnJso4C6rfsZNX8IuvQ6NULjd6nTx9atGgBQPfu3cnLy+PSpUtOfQYaO1gLYc2foW5TGPp39/fXMg7u+hI6joAv/uby/zeNa9HeU/PjLy/rPgH6/xE2zLKdneyzx42UltlnDRfC0txReVyl4rDltkKj79+/nxEjRpTEkkpMTGTXrl2EhITQpUsXHnroIcLCwkpCo4eHh/P888/zn//8h6eeeqpEpr3Q6Dt27KBhw4aMGDGiJDT6pk2bePHFF+nbt6yXXXFo9HHjxpGenl4SImTz5s1MmTKlTGj0oKAg7r//fj744AOmTv3tvpQOjV5YWEhsbCxxcXEl54tDo69bt45Zs2axYcMGZs+eTUJCAq+99prde7l8+XL69OljN+yKxklEIOY2w8MptJ5n+ozsAEc2V5wNUKeQ9Rm00bDH+QqSLeWec2k3mzdv5qGHHgIqDo0OlIRGz8zMLAmNDpCfn8+AAZVH2y0dGh0oCY0+YcKECtsMGjSIl19+uSQ0ekZGRklo9FdeeYWFCxeWhEYHyM3N5Yorrrjs+opDowMuC42+d+9eHn/8cb744guH22gcIMACV93n+X51NkC/QBsNeyOD+i3NoXL58lbG3/BIh0YWlWFvg6UOjW6btLQ0brzxRt5//306dOjg+AVp7LPmEWg9AHrd7Pm+K/x/qyBLoMYr6DUNewx/6vJcAEFhRrkL0aHRL8deaPTMzEzi4+OZM2dOyWhL4wJ+3gAJ8yDjsHf699D/m6Z6aKNhj163wNhXzJGFmZ1s7Csun1/VodEvx15o9Ndee42UlBSeeeYZYmJiiImJKTFamipSkAvr/gqRneDqh72jQ8n/mzmyCApzy/+bpnro2FMaj3Dx4kXq1q1LTk4OgwcP5u2333Z7pFv9WTvBpmfh2xdg2mpoN9jb2sCav8D+tfCXZLDoWXRPo2NPabzOPffcQ3JyMnl5eUybNk2HRvclstJg88vQa7JvGAyAoU/CyDnaYPgg+hPReAQdGt2Hqd8SbnkfWtp8sPQO4W7N5KypBnpNQ6OpzRSamyK7joG6V9iv62l+3gBvDTZ2i2t8Bm00NJraSm4GvBoHiT46CgwOh+M/wYHPva2JphTaaGg0tZWNs40NrE17eFsT27S6EiJawN4V3tZEUwptNDSa2khaAiTMhyvvg+a9vK2NbQICjJA+KRsgL8vb2mhMtNHwAjo0+m9UJTT6tm3bSvZn9O7dm08++cQl+tYaigMSRjSHoX/ztjb26T4RrPmw3/Z3QeN5tPdUJaxNXcvcnXM5kX2CZuHNmB47nfj2NoIc+jnFodFdaTQq+tF3hMTERBISEhgzZsxl53r06EFCQgKBgYEcP36c3r17M3bsWAID9dfZLruXGlNSWWmAMkYZIRHe1so+LftCn9/rUCI+hB5p2GFt6lpmbpnJ8ezjKBTHs48zc8tM1qZWL96UDo1evdDoderUKTEQeXl5iM4xXTnlw/wD7Fzo+2HHRWD869BukLc10ZjU+kezOz6/47KykW1HMqXrFF7e8TJ51rwy5/KseczZNof49vFk5GXwyNePlDk/f9T8SvvUodGrHxr9xx9/5M477+TIkSMsWrRIjzIqY+Ns/w47nvmr4XrbtLu3Nan16P80O5zMOWmzPOuSaxfldGh050OjX3nllezdu5d9+/Yxbdo0Ro8eTWhoqENtayX+HHZcKXh/PDRsC3/Q61fexiGjISKjgLmABXhXKfXPcudDgPeBOOAsMFkpddg8NwO4C7ACDyul1tuTKSIfAH2BAmAbcK9SqkBEhgCfAofMblcopWZX7bJ/w97IoFl4M45nH7+svHl4cwAahjZ0aGRRGTo0uvOh0Yvp1q0b4eHhJCUlXTZK0pTCn8OOi0D0BPh+LmSfgfDG3taoVlPpmoaIWIDXgdFANHCriESXq3YXkKGU6gi8BDxvto0GpgDdgVHAGyJiqUTmB0BXoCcQBtxdqp/vlFIx5qvaBqMypsdOJ9RS9uk11BLK9NjpLu1Hh0a/HHuh0Q8dOlRiXI4cOcKBAwec9hardQx/CqTcv7s/hR3vMRGUFfat9rYmtR5HFsL7AylKqVSlVD6wGCgfT3s8sNA8XgYMF2N1cjywWCl1SSl1CEgx5VUoUym1TplgjDS89igU3z6emQNn0jy8OYLQPLw5MwfOdLn3lA6Nfjn2QqNv3ryZ3r17ExMTw4033sgbb7xB48b66dMuna431r+DI3BnmH+30bQHRHbUG/18AaWU3RcwCWP6qPj9H4DXytVJAlqWev8L0Bh4Dfh9qfJ5pjxHZAYBO4FB5vshGFNfPwGfAd0r0z0uLk6VJzk5+bIyjfu5cOGCUkqp7OxsFRcXp3bs2OH2PvVnXYod7yv1dD2l0tx/393GxmeVmt1YqZwMb2tS4wESVAW/q46sadjyZyw/SV1RnYrKbY1wyst8A/hWKfWd+X4n0EYpdVFExgArgU6XKStyD3APQOvWrW10o/EGOjS6l+kyBib8F1r08bYmVefKe6H/PRDWwNua1GocMRppQKtS71sCxyqokyYigUB94FwlbSuUKSJPA02Ae4vLlFLnSx2vE5E3RKSxUupMaUWUUm8Db4ORhMmB69N4AB0a3cuER0LM77ytRfXQC+A+gSNrGtuBTiLSTkSCMRa2V5WrswqYZh5PAjaZQ5xVwBQRCRGRdhgjg232ZIrI3cBI4FalVFFxByLSzFwnQUT6m7qfrcpFazS1iv1rIeE9KLJ6W5PqcywRFo6D8+WfWzWeolKjoZQqBB4E1gP7gKVKqb0iMltExpnV5gGRIpICPAI8YbbdCywFkoHPgQeUUtaKZJqy3gSaAltFJFFEit07JgFJIvIT8AowxTRMTlPFZho/Qn/Gpfj+FfjxbQiweFuT6hMcDoe+geRPva1JraXW5Qg/dOgQERERREZG6vATNRSlFGfPnuXChQu0a9fO2+p4l8yj8HIPGPYkDH608vr+wH+vgeA6cNcX3takxqJzhJeiZcuWpKWlcfr0aW+ronEjoaGhtGzpBxvX3M1ecwd194ne1cOVdJ8Am54xdrP7w+bEGkatMxpBQUH66VNTe0haBi1iIbKDtzVxHd1vNIzG3k9g4EPe1qbWUeuMhkZTa8jPhoAgYzd1TSKyA8TcBvVaeFuTWok2GhpNTSU4HP640Qj4V9OY8Ia3Nai16HwaGk1NRCnIM7c21VSHj7zzcKry8Dka16KNhkZTEzm5F17oAAdrsIfRh5NhxR+9rUWtQxsNjaYmkrTc2MznZ2FD1qauZcSyEfRa2IsRy0bYz5LZ7QY4sRvO/uIaeRqH0EZDo6lpKGUYjfbXQl3baX99EafTK0ebEZoriHzrrnTNtR29EK7R1DTSd0LmEbj2cW9r4hRzd861mV75pR0vEd8+ntSsVHad3FW2UeueXJ+0gnqDH+VgxkH2nN5TcuqlHS/ZlDd351yXpzeoTWijodHUNJKWgyUYuvrXD+OJ7BM2y0/lGAm+dp7cyayts8qetECfjGPUO/sLW09s5cWEF6vcj8YxtNHQaGoafe8w1jL8LIR4RemVm4U3A2BMuzFcE3VN2ZO5mUQWFkJkBybVa8bItiNLTt227rYSg2NLnqZq6DUNjaam0bgT9LrZ21o4zfTY6YRYymasLJ1euU5QHZqFNyv7atyVoGY9AAgPCi9z7pG4RzySrrm2oY2GRlOT2LkIfv7S21o4jbXIyuh2o5k1cJbz6ZUzDsPyu+Fkcpni0umai7mv9316PaOa6OkpjaamYC2AL5+CDsOMnOB+xPKfl7P4wGLeHfGu8z/qQXWMdZyGbaFpdJlT8e3jiW8fT0ZeBiOWjeDI+SOuU7qWokcaGk1NIfUbyD0HPW7ytiZOkVeYx1s/vUV4YDgNQxo6L6DuFdB2ECStqDBkSsPQhozvOJ41qWs4naMjXFcHbTQ0mppC0jIIqQ8dh3tbE6dYcmAJp3JP8XDsw1XPcdP9Rjj3i7HZrwKmRU8j2BLM3rN7K6yjqRxtNDSamkBBHuxbA93GQmBI5fV9hIv5F3l3z7tc3eJq+jXrV3VB3caBWIzRRgW0qteKTTdvYkirIVXvR6ONhkZTI8g4DGENoad/TU0t/3k5mZcyeSi2mnkxwiMh5ncQ3thutTpBdVBKcSb3TPX6q8XUunSvGk2NRSnjFeA/z4IFRQVsP76dgVEDPdbn7K2z2XJsC2tuXENggPYFsoW9dK/+8+3SaDS2Kcw3PKdE/MpgWIusBAUEudZgWAvtBjAEuCbqGtIvprPhyAbX9VuL8J9vmEajsU3yp/DvLnDukLc1cZgT2SeI/ySerce2ulbwJ/fCwnFQVFRhlSGthtCmXhvm751PTZ5pcRfaaGg0/k7SMggMhQZtvK2Jw7y9+21O5pykdb3WrhXc6Xo4nwbpFU9LB0gAU6Onknw2mYSTevraWbTR0Gj8mZxzkLLRcDn1k6mpX8//yic/f8LNnW8mqm6Ua4V3GW0Ea7TjRQUwrsM4GoY05OODH7u2/1qAQ98yERklIgdEJEVEnrBxPkRElpjnfxSRtqXOzTDLD4jIyMpkisgHZnmSiLwnIkFmuYjIK2b93SISW50L12hqBPvXQFGBX23oez3xdQIDArmn1z2uFx5aHzpeD8kr7U5RhQaG8ub1bzJ74GzX61DDqdRoiIgFeB0YDUQDt4pIdLlqdwEZSqmOwEvA82bbaGAK0B0YBbwhIpZKZH4AdAV6AmHA3Wb5aKCT+boH+G9VLlijqVHsWQaN2vtNhr4j54/w2aHPuK3bbTQOs+8eW2V6TIQLx+HoD3arRUdGExoYqtc1nMSRkUZ/IEUplaqUygcWA+PL1RkPLDSPlwHDxdjaOR5YrJS6pJQ6BKSY8iqUqZRap0yAbUDLUn28b576AWggIr9FItNoaiPXPg4jnjM8p/yA1hGtefO6N7mjxx3u66TzKJi2GlpdWWnV7Se2M+HTCTq0iBM4YjSigKOl3qeZZTbrKKUKgSwg0k7bSmWa01J/AD53Qg+NpnbR9mroOsbbWjiEUgoRYWDUQOqH1HdfRyF1od1gCLBUWvWKOldwKOsQH+3/yH361DAcMRq2HmHKj+cqquNseWneAL5VSn3nhB6IyD0ikiAiCadP66cHTQ1m6xtwYk/l9XyEBzc9yDu73/FMZznnYP3f4dcf7VZrU68Nw1sPZ8mBJeQU5HhGNz/HEaORBrQq9b4lcKyiOiISCNQHztlpa1emiDwNNAEecVIPlFJvK6X6KqX6NmnSxIHL02j8kMyjsH4GHFzvbU0cYuuxrXyb9i11gup4psPAUEiYDz9VPoKY1n0a5/PP80nKJx5QzP9xxGhsBzqJSDsRCcZY2F5Vrs4qYJp5PAnYZChMTo0AACAASURBVK5JrAKmmN5V7TAWsbfZkykidwMjgVuVUkXl+phqelFdBWQppS7PDanR1Ab2mi6lPSZ6Vw8HUErxys5XaB7enJs7eyijYHAd6DIK9q0ydonbIeaKGGKaxLAoeRGFRfbrahxIwqSUKhSRB4H1gAV4Tym1V0RmAwlKqVXAPGCRiKRgjDCmmG33ishSIBkoBB5QSlkBbMk0u3wTOAJsNcMkr1BKzQbWAWMwFtNzADeupGk0Pk7ScmgRa3hO+Tibjm4i6WwSswfOJtgS7LmOu99o3KdD31QaLn567HSy8rMIEP/Y6+JNdMBCjcbfOJMCr8XByP8HAx7wtjZ2UUoxafUk8q35fDL+E88GCCzIgxc6QvcJMP41z/VbA7AXsFCHeNRo/I0zB4xNbN1v9LYmlSIiPD/oec7nn/d8RNmgUIiKg72fwK7/Qf2WMPwp6HWLzerZBdnMT5rP1VFX0+cK/9j34g200dBo/I2u8fDoL2AJ8rYmDtGxYUfvdLx7KaT9CAW5xvuso7D6YePYhuGwiIWlB5ay/9x+XhuuRyYVoSfwNBp/ojDfyJnhBwZj6YGl/PXrv3rPlXXj7N8MRjEFuUa5DUIDQ7m16618k/YNqZmpHlDQP9FGQ6PxJ76eA/+92sif4cPkFuby5k9vcib3DGGBYd5RIivNuXJgctfJhFhCWJi8sMI6tR1tNDQaf0EpwxsoopnPjzQW71/M6dzTPBz7MOKtECf1WzpXDjQKbcSEjhNY/ctqnRK2ArTR0Gj8hfQdkHnE5yPaXsi/wLykeVwTdQ1xTeO8p8jwpyCo3CgnKMwot8Mfov/AoKhB5Jaf2tIAeiFco/Ef9iwzckV0u8HbmthlUfIisi5l8XCfh72rSPFi98bZxiK4WGDsKxV6TxXTpl4b5g6b6wEF/RNtNDQaf6DIariOdhphuNv6MFO6TiGqbhTdIrt5WxXDQPS6BXZ9AJ/eD407Odz0cNZhTuWcon/z/m5U0P/QRkOj8QeUgpHPQT3fD+zcKLQR4zuWz57gZbqMhpvmQaTj7r9Pb3makzknWXPjGs/vMfFh9JqGRuMPWAKh5yRoM8DbmlTIiewT3P757RzMOOhtVS6nTiPj/oVEONxkavRU0i+ms/HXjW5UzP/QRkOj8XWsBbDlVTh/WVBnn+LNn95k9+nd1A2q621VbJN9Fr6fCxmHHao+pNUQ2tRrw4KkBTq7Xym00dBofJ3Ur+GLJ+FYorc1qZDDWYdZmbKSW7rcQou6Lbytjm0KsuHLpyD5U4eqWwIsTI2eStLZJHac3OFm5fwHbTQ0Gl8nabmx+F1JpFZv8kbiGwRbgrm7593eVqViGrSGZr1g3xqHm4zrMI5m4c34JfMXNyrmX+jVHY3GlynIM37kuo+HwBBva1OGtalrmbtzLieyT6BQDG05lMZhjb2tln26jYWv/h9cOGFskqyE0MBQ1t24jiAf30zpSfRIQ6PxVXYvhZeiIf+CkaFv91Jva1TC2tS1zNwyk+PZx1Fm1uUtx7ewNnWtlzWrhK7xgIID6xxuUmww0i5UHH6kNqGNhkbji+xeakRkzTlrvM8+bbz3EcMxd+dc8qx5ZcouWS8xd6ePb4q7IhoadzbS5TrBe0nvMW7lOB1aBG00NBrfxMkIrZ7mRPYJp8p9BhH40xa47mmnmg1vPZzCokI+3PehmxTzH7TR0Gh8kSpEaPUkzcJtrwdUVO5TFK9PFBU53KRNvTYMaz2MJQeWeC/Uu4+gjYZG44tUIUKrJ5keOx2LWMqUhVpCmR473UsaOcnHt8OKPzrV5Pbut3M+/zwrU1a6Ryc/QRsNjcYX6XfX5WUORGj1FMNaDyMoIIhQSyiC0Dy8OTMHziS+fby3VXOM0Ppw8HMovORwk5grYmgV0YoXE16k18JejFg2wvcX/t2AdrnVaHyR7DOAQL3mcP54pfmtPc261HXkWfOYP3I+fZv19bY6ztP1BtixAFK/gc4jHGqyNnUtp3JOUVBkJMA6nn2cmVtmAviPsXQB2mhoNL6GtQB+Wmy4h075wNva2KR5eHMmdJzg3XwZ1aHdYAiOgP1rHDYac3fO5ZK17Mgkz5rH3J1ztdHQaDRe5OB6yDkDsVO9rUmFDIwayMCogd5Wo+oEhkCn6439GkUvQYCl0iZ+6zHmYhxa0xCRUSJyQERSROQJG+dDRGSJef5HEWlb6twMs/yAiIysTKaIPGiWKRFpXKp8iIhkiUii+fKNyV2NxtW0vgpG/ws6+GbYkPWH13Mu75y31ag+cbfDNX9xON+6X3uMuZBKjYaIWIDXgdFANHCriESXq3YXkKGU6gi8BDxvto0GpgDdgVHAGyJiqUTm98B1wBEb6nynlIoxX77hsK7RuJrwxnDlvUY4dB/j2MVjPPbtY7y/931vq1J92l8LAx6AoFCHqk+PnU6opWxdv/IYcxGOjDT6AylKqVSlVD6wGCifYWU8sNA8XgYMFyOb/HhgsVLqklLqEJBiyqtQplJql1LqcDWvS6PxT5KWQ+JHRtIlH2TxgcUIwuQuk72timvIzTAyIjpwv+PbxzNz4EyahzcvKbu9++21aj0DHDMaUUDpPfdpZpnNOkqpQiALiLTT1hGZthggIj+JyGci0t2B+hqN/6CUEUxv50Jj57KPkVuYy/KDyxnWehjN6zavvIE/sPcTY8/GqWSHqse3j+eLSV+w9dathFhCyLyU6V79fBBHjIatb295s1xRHWfL7bETaKOU6g28CtjcYSMi94hIgogknD59uhKRGo0PcfRHOJsCfX7vbU1ssjZ1Lefzz3Nbt9u8rYrr6BIPiFPh0gHqBtfl/pj7GdDCdzMpugtHJk3TgFal3rcEyqcQK66TJiKBQH3gXCVtK5NZBqXU+VLH60TkDRFprJQ6U67e28DbAH379vXNMb6mxrFyVzovrD/AscxcWjQI49GRXZjQx8l83rsWQXBdiJ7gOpkuJDUrlW6NuhF7RazXdHA5EU2hVX/D9XbI4041vbPHnW5SyrdxZKSxHegkIu1EJBhjYXtVuTqrgGnm8SRgkzLyI64CppjeVe2ATsA2B2WWQUSameskiEh/U/ezjlykRuNOVu5KZ8aKPaRn5qKA9MxcZqzYw8pd6Y4LuXQRkj6B7jdCSF3XyHQxj/V7jEVjFiE+OHVWLbrGw4ndkGHL98Y+xy8e5+ujX7teJx+m0pGGUqpQRB4E1gMW4D2l1F4RmQ0kKKVWAfOARSKSgjHCmGK23SsiS4FkoBB4QCllBcO1trxMs/xh4DGgGbBbRNYppe7GMEZ/EpFCIBeYonTiXo0P8ML6A+QWWMuU5RZYeXZtMpcKy5YHWQKYGGvEj9qScoajGUbwu4gLqQwIa8meuqMYbEfmC+sPeGW0cS7vHI1CGxFi8a1EUC6h6w1GGtjD30HDNk41fWv3W6w7tI6vb/maOkF13KSgbyE1+Xe3b9++KiEhwdtqaGo47Z5YW+mCXDERoYHsmWlsV3rgw52s3X281FlFs4hQfvj7dRXKFODQPz3rrZN+MZ0bPrmBWQNnMa7DOI/27TEyDkPDtk43235iO3euv5N/Df4Xo9uNdrla3kJEdiilbMaH8T1HcI3Gz2heP5RjWXk2y5f/qeyu6dIzO8+M78Hfx3RD8jIgIAgVXBdLgFGhRYMw0jPL5dMwyz3Nkv1LUErRv1l/j/ftMapgMADimsZxRZ0rWJe6rkYZDXvoKLcaTTXIK7ASEXZ5/uiwIAuPj+pKiwZhZV7N6//2o98oPNgoS3qb5vPiaBFmpWk9Y/PYoyO7EBZUNrSFAH++rpNbr6c8uYW5LP/ZcLOt0Tuf887Dx3cY+2ScIEACGNNuDJvTN5OZVzvcb7XR0Giqwbo9xzl48gJT+rUiqkEYAkQ1CGPOxJ6OrT1YC43NfG0GQkjdkuIJfaKYM7FniczGdYO5pW9Lbu7bqmJZbqBGutnaIiQC0rbDnmVON41vH08RRSSeTnSDYr6Hnp7SaKrBjX2i6Nw0gh5R9asmIGUDXDxhc2/GhD5RNg3PtkPnaN2oDs3qOxb+ojp8fPBjujbqWrPcbG0hYnhR7VgA+dkQHO5w0y4Nu/DVLV/RKLSR+/TzIfRIQ6NxkrwCKw99tIvkY+cRkaobDDD2ZoQ3gc4jK68L5OZbuf+Dndz6zg+csLGO4mpeH/46swfOrnlutrboGg+FeZCy0almIlJiMGqyY1Ex2mhoNE6QV2Dlj+8nsGb3MfafOF95A3vknDOyx/We8lve6koIC7bw1h9iOX3hkkcMR+OwxnSL7ObWPnyG1gMhrKGx0c9J8q353P3F3SzYu8D1evkY2mhoNA5SbDA2p5zh+Zt6ley3qDJ1GsG938KV9znVLK5NIxbe2c+thiP9Yjp3rr+TgxkHXS7bZ7EEQr+7oXFnp5sGW4LJK8xjTarzBsff0EZDo3GA0gbjXzf14hZXLUg37W6kcnWS0obj/a2HXaNLKRbvX8zOkzupF1zP5bJ9mmFPwuD/q1LTMe3GcDDjID9n/OxipXwLbTQ0GgcJDBD+dVMv13gwpe+E5X+ErKqHBYlr04hPH7yav47oUn19SpFTkMPyn5dzXZvrarabbUVYC+BMitPNRrYdiUUsrDu0zg1K+Q7aaGg0dsjNt5KVU0BokIX3bu/nOpfXne/DvtWGq2c16NDE2BB4LDOXae9t43jW5RsCnWVN6hou5F+o+W62FbHyT7BwLBQVOdUsMiySq1pcxbrUdTV6QVwbDY2mAnLzjSmpqfO3UWgtcp0HUX6OsYms+wQIdc30z6kLl9hxJIMpb/9QLcOhlOKj/R/RrVE3YprEuEQ3v6PDcLhwDI7vcrrpbV1v43fdfkdBkWMpZP0RbTQ0Ghvk5lu5+/3tfP/LGaZe1YZAiwv/VfathkvnXZo3I6ZVA96/qz9nL+ZXy3AUqkImdZ7Efb3vqx1utrboPBLE4nSODYBBLQcxrfs0gi3BblDMN9BGQ6MpR7HB2PLLWV6c1Jub4qrpJVWeXYugYTtoc7VLxca2blhtwxEUEMRt3W5jWOthLtXNr6jTCNpeDfvXVql5dkE2q39ZXWNHG9poaDTleOrTJLb8cpZ/3+wGg1FUBFGxMOABt6R0LTYcTeuFEhpoqbxBKdIvprP84HLyCt2/adDn6ToWzhyAM857Qm0/sZ2/bf4bW49tdYNi3keHRtdoypGWkcNPR7OI7+W/ebCVUogIlwqNhfwr6lUecuTF7S/yv33/4/ObPq+dXlOluXgazv4Mra6EAOeMb4G1gKEfD+XqFlfz/ODn3aSge7EXGl2PNDQajCmpd79LpahI0bJhHfcYjCIr/PKV8dfNFK9HPLZsN5Pe3MoxG2HWS5NTkMOKlBW11822PHWbGEEknTQYAEGWIEa0GcFXR78ipyDHDcp5F200NLWe3Hwrdy3cznPr9rHz1wz3dZT6FSyaAAc858d/+8C2ZGQbaxz2DEetd7O1RcZh+HwGnD9eadXyjGk3htzCXL46+pXr9fIyOsqtplayclc6L6w/wLHMXIICA8gvLOKlyb3p29aNkUp3/Q/CGkGnEe7roxx9zDWOqfO2ccOr3xFkCeDU+Uu0aBDGoyO7MKFPFEopPtz3Ye12s7VFQR788AZEdjDCizhBbNNYmoU3I/FUIvHtPZtp0d1oo6Gpdazclc6MFXtKcnDnFxYRZBEEN7qY5pwzvHH63gWBns2z3ad1Q+4e3I6XvvxtUTc9M5cZK/YAMKhrKCGBIdza9dba62ZriyZdoFEHw/XWSaMRIAEsvWEpDUMbukk576GnpzS1jhfWHygxGMUUWBUvrD/gvk73fAzWfOjjnemfpdvTLivLLbDywvoDRIZFsjh+cc3N/11VRKDbDXD4O8h1PitfscGoac5G2mhoah0Vze1XtlhcLX7+AprHQLOe7uvDDhVe84UzZF3KQkQIEP1zcBldx0JRofH5VYFXd73Knzb8ycVKeRf9LdHUOppE2J4eatEgzGa5S/jdUpjygfvkV0JF19aoxWZGLR9FdkG2hzXyE6LioEk3yMuqUvNQSyjfH/ue9ItVD0zpa2ijoalV5BcWEWy5fN4+LMjCoyNdGy22BKUM180qhEB3FY+O7EJYUDn3UclH6m3j6qirCQ9yPL1prSIgAO7fCv3/WKXmo9uNBuCzQ5+5Uiuv4pDREJFRInJARFJE5Akb50NEZIl5/kcRaVvq3Ayz/ICIjKxMpog8aJYpEWlcqlxE5BXz3G4RqeFJizXu4N9fHCAtM4+7rm5LVIMwBIhqEMaciT1t5uOuNgV58MYAI0ChF5nQJ4o5E3uWXPMVESGENNjFpaJsbu3yO6/q5vOIGIa/wPnpy5YRLYlpEsPa1KqFJPFFKvWeEhEL8DpwPZAGbBeRVUqp5FLV7gIylFIdRWQK8DwwWUSigSlAd6AFsEFEitNiVSTze2AN8HU5VUYDnczXlcB/zb8ajUN8n3KGt75N5XdXtuYfY7vzj7Hd3d/p/jVwep/hautlJvSJKjGMSimGLX6Jk1lR7DhQnzi9n69ilIK3h0CzHjD+daebx7eP57kfn+NgxkE6N3Q+K6Cv4YjLbX8gRSmVCiAii4HxQGmjMR6YaR4vA14Tw3dvPLBYKXUJOCQiKaY8KpKplNpllpXXYzzwvjJcEX4QkQYi0lwp5fzOG02tIzffyl+X/kSHJuH8Iz7acx3v+h/UbwXtrvVcnw6QfC6ZM/m/EtvgXto1qettdXwbEYjsCAc+M3bzO7lLfETbEaRfTCciqHq5U3wFR4xGFHC01Ps0Ln/CL6mjlCoUkSwg0iz/oVzb4jmAymQ6okcUoI2GplLCgi08O6EHzeqHEhbsfGiIKpH5K6R+Ddc+bsyN+xDdI7uzcvxKWka0JMTi2X0jfkm3GyBpGfz6gxEB1wkahTbir33/6ibFPI8j32Rbu33KOx5XVMfZ8urqgYjcIyIJIpJw+vTpSkRqagNZuUaI6uuim9Ijqr7nOk78yPgb41trBsX7Bjo06FBiMN7+9hee/jTJm2r5Nh2vA0uwMd1YBaxFVrYc20JqZqqLFfM8jhiNNKB0jsuWwLGK6ohIIFAfOGenrSMyq6IHSqm3lVJ9lVJ9mzRpUolITU0n5dQFrvnnJtbsruzr5QbaXwvD/g4N23i+bxusTV3LiGUj6PV+L/r9rx+rf1ldcu5cdgELtx7h86QTXtTQhwmJgPZDDKNRhc16l6yX+PNXf+aDfd5zu3YVjhiN7UAnEWknIsEYC9urytVZBUwzjycBm8y1h1XAFNO7qh3GIvY2B2WWZxUw1fSiugrI0usZGntcKrTy8EeJBAUG0N+dMaUqovVVMPhRz/drg7Wpa5m5ZSbHs41/mTxrHrO3zi7x6nnk+s70jKrPEyt2cyJL59OwycCH4LpZoJzLHQ5QJ6gOQ1oNYf2R9RRY/Ts5U6VGQylVCDwIrAf2AUuVUntFZLaIFMcdmAdEmgvdjwBPmG33AksxFs0/Bx5QSlkrkgkgIg+LSBrGSGK3iLxr9rEOSAVSgHeA+6t99ZoazYvrD5B8/DzP39TLoXwSLmX3x3Bqv2f7tMPcnXPJs5Y1BnnWPObunAtAcGAAL0+J4VJBEY8sTaSoqGaFvnAJ7QZDj4lVCpcOEN8unqxLWWw5tsXFinkWnYRJUyP59uBppr63jd9f1ZpnJ3g4dEduJvy7i5EDPP7fnu27Anot7IWysWwoCLun7S55v3jbr/ztkz0suXcA/bwxOvN1Mg5D6jcQN63SquUpKCpg2NJhDGgxgH8N/pfrdXMh9pIw6Si3mhrJL6cv0rVZBH8f40H32mKSlkNhnmE0fIRm4c1KpqbKl5dmcr9WxLZpSOemNcM91OVsfBaSPobV040d/sOfgl63ONQ0KMBIzrTl2BasRVYsVRyxeBvf8gPUaFzEHVe3Y9WD13jOvbY0uxZB0x5GgEIfYXrs9Mtca0MtoUyPnV6mTERKDMaWX86Qk1/oMR19nt1LYX+x84CCrKOw+mGj3EGmx01nzY1r/NZggDYamhrGip1pfH3gFGDM03ucE0lwbJcxyvCR3BRKKUa3G82sgbNoHt4cQWge3pyZA2dWmCDo6Lkc/jBvG7NXJ9s8XyvZONsYQZamINcod5B6wfWwBFj8Oly6np7S1Bh+PnmBGSv2cFX7SK7t3MSzCYV2LzV+PLKOAgI+tGFu06+beHvP27w67FW+mORYiO9Wjepwz+D2/PfrX7i2cxNG93RDznR/I+vynCR2yyvgu7TveO7H5/gw/kMahfrfupEeaWhqBJcKrTy8OJG6IYG8cHMvzxuM1Q+bBgNAwRd/c2rawl3kFuby/PbnKSgqcPoH6pHrO9O7ZX2eWLGH41luzDXiL1QUpdjJ6MVNw5uSfjGd9YfXu0Apz6ONhqZG8K/PD7Dv+HleuLkXV0R42L124+zLI6A6OW3hLt7Z/Q7Hs4/z9yv/TmCAcxMLQZYA5k7pQ4G1iL8sScRa291whz8FQeXykgSFGeVO0LlhZzo17MS61HUuVM5zaKOh8Xt2p2Uyb/Mhpg1ow7CuTT2vgIumLVzN4azDLNi7gLHtxxLXNK5KMto2DufZCT0Y07M5Ab6xROM9et0CY18xAlAixt9rn3DYe6o0Y9qNIfF0ImkXvPsdqQraaGj8np5R9Xlpcm9mjOnmHQVcNG3hahYmLyTEEsIjfR+plpyJsS2ZOqAtIqI3/fW6Bf6SBDMzYeRzsOFpOPy902LGtBsD+GdyJm00NH6LUoqT5/MQEW7s05LQ8pnpPEWPiZeXVWHawtXM6D+Dd0a8Q+OwxpVXdoD1e08w4Y3vyb6k3XAB6Hg91G0GXz3ndDyqFnVbcH/M/cQ29b9cctpoaPyWRT8cYfi/vyHl1EXvKZGfA3tXQvgVUK8lJdMWY1+p0rSFK8gtzCWnIIdgSzA9Gvdwmdz6YUHsSc9i1uq9LpPp1wTXgUF/hSPfGyHwneRPvf9U5WlDb6KNhsYvOXjyAs+t3Ufftg3p0MSL+a0zDoG1ACa9B4/sNaYt/pLkNYMBxuL3uJXjOJ9/3qVyr2ofyf1DOrA0IY11e3SsUMAIJ1KvJWx6tkrRb1MyUvgu7Ts3KOY+tNHQ+B15BVYe/mgXEaGBvDCpt2fda8vTtDtMT4R2g7ynQykOZx1m/t759G/Wn3rB9Vwu/8/XdaZ3qwY8sXw3xzK1Gy6BIXDtY3DmZ+MBwkle2vkSs7bOoqgKkXO9hTYaGr/j+c/3s//EBV6Y1JsmEV7aRFdUBDsXQWG+8cPhAyilmLNtDqGW0GovfldEkCWAV6bEYC1SerRRTMzv4M8/QaP2TjdtHt6ckzkn6f1+b0YsG1ESqt6X0TvCNX7Byl3pvLD+AMcycwkLtjCoU2OGdr3CewrtmA9rHzGS83Sf4D09SrHx141sObaFJ/o/4bLFb1u0iQzny0euZduhc1z9z00cy8ylRYMwHh3ZhQl9oioXUNOwBEFYQ+NB4uJJqOfY7vm1qWv5NOXTkvfHs48zc8tMgArDu/gCeqSh8XlW7kpnxoo9pGfmooCcfCsJh8+xcle6dxS6cAI2zDLyK0SP944ONvjq6Fd0btiZyV0mu72vbYfOlflM0jNzmbFij/c+E19g+Z2w6EYosjpUvbIcJ76KNhoan+eF9QfILSj7j5hbUMQL6w94R6HPHjcC193wss8EJQR49upneXfEu07v/K4Ktj8Tq/c+E1+g2zg4vQ/2fuJQ9RPZtlPrVlTuK2ijofF5Klpw9cpC7MH1kLzSSOMa2cHz/dsg/WI6J7JPICI0DG3okT596jPxFaInGCHxv/p/YK18L0v5XCaVlfsK2mhofJ6IUNtPzi0ahNksdyt1m0L3iXD1w57v2wZKKWZvnc1ta2/zaO7piu598wYejvvlSwQEwJAZcO4X2L2k0urTY6cTail7v4ItwZflOPE1tNHQ+DRrdh/jfF4hlnKzQGFBFh4d2cXzCrWIgZvn+4zHVPHi95097yTIEuSxfh8d2YUwGzvwr+3cxGM6+CRd443kW7sXV1o1vn08MwfOLMlxUi+4HgXWAjo17OQBRauOzhGu8VmUUtw+fzu5+VZu7tuSlzf87D1PneO7IWEeXDcLwhp4rl875BTkMP7T8dQLrseSG5Z4ZC2jNKU92lo0CGVY1yt4emx3Ai21/Fk08yhENDO8qpxplpfJDStvoHPDzswbMc+r+4/s5QjXRkPj01wqtJJXUET9MM89RV9GkRXevc7Il/HgdsO90geYu3Mu7+55l4WjFvpUDKNT5/NIOX2RgR3c5/brFxTkggQ4NSpdsn8Jz/74LC9c+wKj2o5yo3L2sWc0avkjgcYXOX3hEtMX7yIjO5+QQIt3DQbA9nfh2E4Y9U+fMRhKKS7kX2Bch3E+ZTAA/vFpEnctSGB3Wqa3VfEeF07C3BhImO9Us0mdJ9G1UVf+nfBvcgpy3KRc9dBGQ+NT5OQXctfC7Xyx9yRpGT7giZOVDhufgQ7DocdN3tamBBHhyaue5Jmrn/G2KpfxzIQeNAoP5s4FCRw955s/fG4noik07gTf/dsIaukglgALM/rPQCnF0QtHK2/gBRwyGiIySkQOiEiKiDxh43yIiCwxz/8oIm1LnZthlh8QkZGVyRSRdqaMn02ZwWb57SJyWkQSzdfd1blwje9hLVI8/NEuktKzePXWPvRsWd/bKhn5EooKIf7fPrMn48fjP7L3rBFpNkB877nviohQFtzRj/xCK7fP30ZWjue8unyKYU9C9inY/o5TzWKbxvLZxM/o0sgLjh4OUOk3TkQswOvAaCAauFVEostVuwvIUEp1BF4CnjfbRgNTgO7AKOANEbFUIvN54CWlVCcgw5RdzBKlVIz5erdKV6zxSZRSzFy1lw37TjFzXHeui/ZCBj5bDH8abnoHGrXztiaAsfj95PdPMmvLLHx5PbJT0wjentqXo+dyeW5dsrfV8Q6tr4KO+z4QzAAAIABJREFU18Hml+HSBaeaBlmCyLfm+2QecUceU/oDKUqpVKVUPrAYKB87YTyw0DxeBgwXY+l/PLBYKXVJKXUISDHl2ZRpthlmysCU6RuBfTRuJTOngK8OnOLewe2ZOqCtt9UxAhEqBQ1aQbex3tamhHf2vMOJ7BPMuHKGd6P7OsBV7SN5Z1pf/uatjIq+wNC/Q+45SFrhdNMVP6/g/775Pzanb3aDYlXHEaMRBZSeXEszy2zWUUoVAllApJ22FZVHApmmDFt93SQiu0VkmYi0ckB3jZ/QMDyYNQ9dw+OjunpbFYMNT8P/bnJoZ6+nOJR1iAV7FzCuwzj6XNHH2+o4xLWdm9CgTjCXCq2s/umYt9XxPFGxcM/XEDvV6aYTO02kTb02PL/teY9u3KwMR4yGrceZ8uPiiuq4qhxgNdBWKdUL2MBvI5uyiojcIyIJIpJw+vRpW1U0PsT2w+f4x8okCqxFNKgTTECADzw9H9sFP74JDduCxTcCQSulmPPjHMIsYfwl7i/eVsdpFm09wkMf7eKDH494WxXP06KPsR5WeMmpZsGWYB7v9ziHzx9m0b5FblLOeRwxGmlA6af6lkD5R4aSOiISCNQHztlpW1H5GaCBKaNMX0qps0qp4rv+DmAzT6JS6m2lVF+lVN8mTWr57lQf55fTF/nj+wl8n3KGnEuORQZ1O9ZCWD0dwpvAdU97W5sSilQRcU3j+Evfv7g17Lm7uH1gW4Z2acI/Vibx1f5T3lbH8+xZBi91h+yzTjUb1HIQQ1oO4a2f3uJUjm/cN0eMxnagk+nVFIyxsL2qXJ1VwDTzeBKwSRmrdKuAKaZ3VTugE7CtIplmm69MGZgyPwUQkdJB6scB+5y7VI0vcfrCJW6fvw2LCAvu6E/9Ol7ei1HMtrfg+E8w+nkI9QHvLRNLgIV7e9/LzZ1v9rYqVSLQEsBrv4slukU9HvhwJ0npWd5WybM06wU5Z+H7l51u+li/x+jSqAsX8y+6QTHnqdRomOsLDwLrMX6olyql9orIbBEZZ1abB0SKSArwCPCE2XYvsBRIBj4HHlBKWSuSacp6HHjElBVpygZ4WET2ishPwMPA7dW7dI23yMkv5O6F2zl94RLzbu9H68g63lbJwFoICe9Bp5FGxFIvszZ1LSOWjaDnwp4MXjzYL7K62SM8JJD3pvWjYZ1g/rwkEWuR73p/uZwmnaHnLbDtHWPjnxO0qteK90e/T/sGzmcGdAc6jIjG4+z6NYOp87bxn8kxXO8rrrXF5J03wj9EeFevtalrmbllZpkkPaGWUGYOnOnTWd0c4eeThvtpp6YRXtbEw5xLhVf7Qv8/GiNZJ8nMy2TJgSXc3fNuLAGXB4t0JTqMiMan6NO6Id89PtS3DMbpA4abbWg9rxsM8N+sbo7QqWkEnZpGoJTik11pXCr0kfUsd/P/2zvz8CiqrA+/p7MTlhC2AGGXAZRNRFZlQEQYUHAZgRmczXHFBXE+FMbRyYiOKKKDCzC4gcy4ICKijKAiiygge8CACIgQSAgQSEIga9/vj1sJnU530kk6XR247/P0k/Stqlu/qq6qc5dT58S2hcvHwZZ5cO5UhTffkLqBV7a/woc/fuh/bRXAGA1DwHj96wO8vf4gADG1wm3VUoKcDJg/EpbcY7eSYmpqVreKsCM5g4nv7+CRRYk4L5ahqoFT4M5VlYphNrTVUHo26cnL214mI9e+OSFjNAwBYVliCk8t283GA+nB9ybzyqlw5hj0vc9uJcXUjajrsTzYs7pVhO4tYpg0tAMfbz/KjC8ukjSxdZtBEyv4hdNZoU1FhMm9JpOZl8nL216uBnG+YYyGodrZfDCdiQu307NVfWaM7hYcbzInLoQXO0NCjI4N1HYQNPfoxR1QilqQk6+cTERIyZDakSGRQZ/VraKMH9iO3/Rqyaur9vPOxkN2ywkMSsGS+2DpAxXetENsB8Z0GMMHez/gh3R7DG1wvLlkuOBwTdCDQIPocF77fU8iPWR7CziJC+GTB/WEdxGH1uvyrqNtkaSUYkHSAubsmMOC4Qu4vt31iAgzt84kNTuVuOg4JvSYUOMnwd0REaaOuoyUjHMkLP2evMJCXlv7k33JtgKBiE7ktWEWXPWQjoZbAe7rfh9n889SN9xzb7S6Md5TBr+zZNsRpizeybn88xOcEaEOnr2la3A8AF7srBMquVOvBUzcFXA5+YX5PLXxKRb/uJhrW17L01c9Ta2wIHFDDhBncguYtWofb31zsMR1ExUWwjM3dwmO68afnDkOM7tBh2Hw6zftVlMK4z1lCCjTV/xQ4sYHyC1wMn1FkIxbZyRXrLwaOZVziju/uJPFPy7mrq53MWPgjIvOYADUjgjl4+1HS1035/ILg+e68Se1G0Hvu2HXh3Ds+/LX98DPmT8z5espAU/WZIyGwa8opThy2nPypKNeygOGUtrd0Rv14gMmpYgFSQvYeXwnz179LA9c/kBQ5scIFN6uD9uvm+qi3wMQURfWTq/U5qdyTvHpgU+ZmzjXz8LK5uK9Qg1+58DxM4ydu8Hr8mYxUQFU40b2SXhvnI4r1agjhEaWXB4WBYOfCJic3EIdRu3ebvfyzoh3GN52eMD2Hax4uz4a1A4i92x/UisWev4JDm3UDhkvdtbzaj7SvXF3RrYbyfyk+fycGbhAkMZoGPzC3LX7GTbza3anZDKmZzxRYSUvraiwECYNtSkTWfIWmN0X9n0B1z0N934LI1/WcxiI/nvDSwGZBFdKMf/7+dyy9BZO5ZwiLCQsaDO0BZpJQzsQ5eYoIcCJM3n8fDLbHlHVSeJC+G4uZB0FlJ5n++TBChmOiVdMJCIkguc2PVd9Ot0w3lMGv5BX4GTIpU34+w2X0rhOJH3bNSz2nrLdC6Zec2jQHoY9A0276rKuowPuKZVfmM/UDVP5aN9HDGk1pJRL7cVO0fXhet1MGNyeerXCaNUgGoBdRzK4rFnd4HDbriornyzpwQf6+8onfb42G0Y15N5u9/L85udZm7yWAfEDqkFoSYz3lKFSZOXkM33FD/Rr15BhneNQSgXXjZy2W7fihj8P1RynxxdO5Zxi4uqJbDm2hbu73s347uMv6vmLyrDrSAY3vLKOQR0aM/XGzjS3c7jTHyTEUDo1EYBAwmmfq8kvzGd+0nzGdhhL7fDafpFmvKcMfuXz71MZ8sJaFmz4mb1W8LmgMRhKwcZ/w79/Cbs/gVMH7VYEwLObni2e8L7/8vuNwagEHePq8NjwTqzff5IhL6zhzXU/1exIuV4dLxSsme5z1siwkDDu6HIHtcNrByTagulpGHwmLTOHvy/9ns92pdIxrg7TbulK9xYxdss6z5k0WDJez120vw5GvQq1G9sqyamcOMRBek46R88cpXPDzrbquRA4nH6Wvy3ZxZq9x+ndJpb37uoTPI2WiuDpJdPQKIjrrBs/t6+oUObIpJNJPLXhKV4c+CJNoqsWdLOsnoaZ0zD4zLf7T7JyTxqThnbgrgFtCQsJotayUvDub+DYLj0kdeUd+s1b2+Qo3k56m3VH1jHr2lnERsYSGxlrm54LiRaxtZj3pytZuuMoWTkFiAhKKfIKnUSE2j8U6TNF8xYrn9TvCNWL1x58XUdD7hltMLJPaOPS++5yh1nrhtdl76m9zNgyg+cGVN/EuOlpGMpkX1oWe4+dYXiXpiilSMnIsdd11p28s/pmCo2AI1u162zjTgGXsezAshIhP5pFN2NL2haGtBrCP6/6J5HuLr4Gv7IsMYXnVuzhnzd1of8lNS8drle+fQU+fwzie8GNs6HhJWWu/ur2V5mzYw4NIhuQnpNe6fAzZk7DUGFyCwr515d7GT5zHU99mkRegRMRCS6DkZIIcwfClwn6e/MethmMhG8TSMlOQaFIyU5hS9oWBrcYzPO/fN4YjADQsHY4Aox7fSP/98EOTmXn2S3JP/S9D25+DU7shTlXwYbZZUbHbRbdDICTOSeLr8WEbxP8mvXRDE8ZSgQXbBYTxS1XNOd/O1PZl3aGkd2a8fj1lxIeanP7InGhSze+ObToA0kfQ60G0H6IrdI8JUwCSEpPMhPeAaJ32wYsf2gAL638kblrD7BqTxpP39SFHCsMSVC4flcGET1c1fpq/WLq8sk6gdOgv3pcffaO2aXKipJ3+SvYpTEaFznuwQWPnD7HSyv3ERMVxlt/upJBHeydSAZKTxhmJEPGImjaHW5bDNENAi5JKcWe9D2sOryKlOwUj+tcSAmTagKRYSE8MqwjN3RrxuQPE/l2/wk+2Jxc4tqesngnQM0yHAB1m8Jv34cd78Il1+qyc6cgMqbE3F0gkncZo3GR8+zyPaWCxAFEhYcEh8EAzy9BAZw9GXCDcSz7GG/seoPVh1eTkp2CIIQ5wsh35pda90JKmFST6NS0LovH92fAc195DYBY44wGaOPQ/bf6f2chvDNGz+GNfAViWgD6mvPUiPHntWj6zhchu45kMO2zPdw06xtSMkoPqwCkeikPCE4nHEuCTW/Ah3d6DmMOAYlKm5WXxfKflvPt0W8BCHGEsGTfEjrEduDJfk+yavQqpvafSmRIyXmLCzFhUk0ixCEcPe35Gj5y+hz/+nJvgBX5GXFAt7FweBPM6gtb3walmNBjQrVfi6ancYFzKjuPTQfT2fhTOn/q35r4+rVITM7gjXUH6BofQ+2IUM7kln6JKKAT3gW52lW2KHPeR3fBzg/0/7WbaN/1Ag89jSpEpXX3dnL1MEnNTmXV4VWsOrSKTcc2UeAsYEirIfRr1o+GUQ1ZN3Yd4SHng+gVbXehJ0yqaTSLifIYcTki1EFalg4YqZRi7NwNtG0UTa82sfRu0yC4nD28IQI9b4d218DH9+ssgLs/YcSNsyF+GDMPfESqA+KcMKHVML9ei8bltgbiPnHtPrl39PQ55qzZz8YD6fxgvbEdEepg9m09uKZjE87mFSAIUeEhHhMm+SXxTYmJaxf/c4CcDB3Z89B6OLQBjmyBwlz4y16o0wT2fwVZqdCyD9Rvow2I+0tQYVGVDjJY5O3kOnkdERLBP/r9gxFtRzDuf+NIPJ5I67qtGdRiEINaDqJrw66EBEE4EoPv+HJtZ+cW8OC72/juYDpZObrxFF8/ioeH/IKbe8QXv2Fd9PJgefeeLTidOmXx1reh9z3w2aQq3ytludz6ZDREZBgwEwgBXldKTXNbHgG8DVwBnATGKKUOWsumAH8GCoEHlVIryqpTRNoA7wGxwFbgd0qpvLL24Y1gMRr+vNCWbDvCuo9m8RDv0UxOcFQ1ZIZzDKktb+DWni24uUc8aVk5DJq+mh6t6tO7TSy92zaga3w9ry8+bVr6b1psnU5jdZw0acThHpO4cuTdlT/gxIUs+3ISM+vWIjU0hLiCQiZknGFE/79Bn7th23/h4/HgCNWT2S37QMu+utUU7jkB0bLVj5dsPbW9iREDp1ZK3uAPBpN2Nq1UeVytOL649Qt2HN9BnfA6tK3XtlL1G4IHX++9QqdiT2om3/2UzsYD6dzaM57BnZqw60gGd8zfTK82sUSEOVi6/Si5BeddXqvawPKrEXIW6myAfshKWSWjISIhwF5gCJAMbAJ+o5RKcllnPNBVKXWPiIwFblJKjRGRS4F3gV5AM+BL4BfWZh7rFJGFwGKl1HsiMgfYoZSa7W0fZWmvjNF46YOH+ThjBcdDhUYFilH1hvLgrS9UqA5Xlmw7wsIv/sqx2M3FdTZJ78noIf8sdXG4Bv3beyyLE2dyyTxXQGZOPpnn8mlUJ4Jtn86la/ibzImtXfxAvif9DGtO30b7a2/n/mt0vuGCQiehvryx7ekBn3mWEddO994ycRZC5hHIzYKcTMjN1H8bd4K4zix75VISooUcx/n9RzqdJGTmMmLCPp3q8vhuPRwVHl2uRE89g8iQSJ7o+wRXN7+akzknSc9Jp31Me2IiY0g6mcSivYtIz0k//zmXzlvD3qJDbAe6zO/icT+CkPiHxPLPmeGiYXdKJrNW72fjgZPFQ1ru1IkI5YHBl1A3Moy6UWHUjQyjZ+v6RIaFkJNfSIhDPEZPqJZevp+CIFbVaPQFEpRSQ63vUwCUUs+4rLPCWme9iIQCqUAjYLLrukXrFR2ee53ANOA4EKeUKnDdt7d9qDIOoKJG46UPHmbBmRUlHnbhTsU1zjZc0XIY4GDskHsBWLn5fySn7ceJQimFUhAWGs7vh+nlS77+gJT0ZL776TMS6x0lz3HeLS7S6aTzqdZER1xJTr6TnIICcvILia0Ty9z7HwXgT/96jpNZJVvD8Q3i6HNmLi83CnPT6GTM6RxuGvY87TvdBMDOzf8mLy9Lh9dQTlCF1K/fjradtZ3d/vXTFOTnsD7pPeZFh5Pn9oD/3Zlz9I7pgMrLhvxzqPxs4jqMos11z1CYfZINr1yK4vzlqYAWl99Om+ueYcjrnUgNKz1d1qiggK/+vJvMvEw+2f8JBc6C8x9VwIDmA+jSqAup2am8sfMNClQBhc5Clh9czjlPcxpuvDr4VQbED2Bt8loe/+ZxYiNjaRDZQIfwiIrltk63EV8n3mtPo2l0Uz7/9efl7sdw8aGUou2U/3l8HHti89+upWHtCGZ8/gMvf7WPWuEhllEJpW5kGP+5ozeDZ6zxOOdSLyqMHX+/DoD1+09yOP0sDofgED3BHxEawrDO2htq26FTnDiTR4hDD6H1+PBq6uWVdq89G9WUWo/u8fl4qxp7qjng2t9JBnp7W8d62GcADazyDW7bFplQT3U2AE4rpQo8rO9tHyd8OAaf+DhjBTluyYPyHMJyx0GWH51DuFMxFm0U3t/6NOsjMkusW6/Qye+t5UuSprMl8hzUB51K5jw5DgeJ9X8mz3GoRHm7AgG00ThX9x1+ql/SXbCeCuM/4SElDIbW6GBBbC0OfDeNOZbReDjxZVJDSu53iCOGFyyjMX7fu2Q5BOqUfls5x+HgtTq1eM2ZrK+QUCAqhHGFJ5gM5IdFck9caXfcO2NieRA45mUY7HiILs/IyWDad9NKLY+NjKVLoy5k5WXx2cHPCJVQQh2hZRqMyb0mF8d16hjbEYAB8QNYM2aN120evuJhjz0X4+1k8EZRNARPD/nmMZGsmPhLMs/lW6MCBcREhQHQ/5KGhIU4SizLys0nItThNY1txrnz7tvvbTrEx9uPlljesHZ4sdF4ddU+vtx9vgE00nEz08Jep5acfyP+rArnufwxxa31quKL0fAU9c3d4Hpbx1u5p3GTstb3VQcichdwF0DLli09bOKd46FeAtwpRUKLOwhxeVjf0edRbjj1EyJ6WEOA8LDzSXUe7P8Y6VmHmfjDHI+B8/IF5na+31U4tSLqFX99rO/jnMlJL1oIQJ2oBvx2w+NeNU4Y8FTx12evfoa8wjxr3w5wOIiNblq8/KVBL1GIkztWP+Q1sN+bQ9/Ux2YtbxylDUV4WBQLfrXAkl109NC4ll4eF16PlPzMUvU1CdfH16x2M74e8zWhjlBCHCGEOkIJldDi/bSv3551Y9cVb3fdous8+p43jW7KuE7jPJ+PMjDeTobKMGloB4/DSZOGdqR2RCi1I0JpRknPqz5tG9Cnred3ibwZoab1zjfkEm64jElDO6CUnndxuuWteeL6y5gw+Bc4laJQKW6ZBeTDI6ELaSYnOaoa8FzBaD7J7RVQo5EMtHD5Hg8c9bJOsjV0VA9IL2dbT+UngBgRCbV6G67re9tHCZRSc4G5oIenfDi+YhoVKNLCSj9AGxcobhn8UImyXl1HlllXj86jAIhLmkOqh7PcpBD6XuF9srlLp1s8lsdtnUFKfkap8qbhMXRqez6cRo/2N5Spr2fra4q381bflXFXetzWIQ66N+7ute4Jff5KwrrHyVHnW0yREsbEPjr0QYgjhJhI30OqT+gxwe89gxFtRxgjYagQnjILVmXi2psRenRYx+Lv9aPDqR/tPUd6ywYlHUeaxUSx9PRVLM27qkS5PxNW+fJy3yagvYi0EZFwYCyw1G2dpcAfrP9/DXxlzTUsBcaKSITlFdUe+M5bndY2q6w6sOr8uJx9+I1R9YYS6RYMLNLpZFS9oZWu86F2NxHpligm0ql4qN1NlapvQp8pREpYyfokjAl9pnjZIrD1gX4gJ1w1labRTRGEptFNSbhqaqUf0iPajiChX0LJ+volmIe+IeDceHlzvpl8DT9NG8E3k6+pkrvtjZc355mbu9A8JgpBP9ir6uruKc+67g35Lw+9ry63w4F/od1j31RKPS0iTwKblVJLRSQSWABcjm79j1VKHbC2fQy4HSgAHlJKfeatTqu8LeddbrcBtymlcsvahzeCwXsK/OsuCmW/mBYM9RkMBvvwhxtvld/TqKkEy3saBoPBUJMw+TQMBoPB4BeM0TAYDAaDzxijYTAYDAafMUbDYDAYDD5zQU+Ei8hx4OdKbt4QP75tXk0YjVUn2PVB8GsMdn0Q/BqDTV8rpVQjTwsuaKNRFURkszfvgWDBaKw6wa4Pgl9jsOuD4NcY7PpcMcNTBoPBYPAZYzQMBoPB4DPGaHhnrt0CfMBorDrBrg+CX2Ow64Pg1xjs+ooxcxoGg8Fg8BnT0zAYDAaDzxijYTAYDAafMUbDAyIyTER+EJF9IjLZbj3uiEgLEVklIrtF5HsRCcqUcyISIiLbRORTu7V4QkRiRGSRiOyxzmVfuzW5IiITrd93l4i8a0V6tlvTmyKSJiK7XMpiReQLEfnR+ls/CDVOt37nRBH5SER8T+gSAH0uy/5PRJSINLRDmy8Yo+GGiIQArwK/Ai4FfiMil9qrqhQFwF+UUp2APsB9QagRYAKw224RZTATWK6U6gh0I4i0ikhz4EGgp1KqMzqFwFh7VQEwDxjmVjYZWKmUag+stL7byTxKa/wC6KyU6grsBSqfMKbqzKO0PkSkBTAEOOS+LJgwRqM0vYB9SqkDSqk8dG6PUTZrKoFSKkUptdX6Pwv9sKt85pZqQETigRHA63Zr8YSI1AUGAG8AKKXylFKn7VVVilAgyspUWYvSGTMDjlJqLaUzZo4C5lv/zwduDKgoNzxpVEp9bmUDBdiAzgpqC17OIcCLwCN4SGMdTBijUZrmwGGX78kE2QPZFRFpjU5MtdFeJaX4F/oGcJa3ok20BY4Db1lDaK+LSLTdoopQSh0Bnke3OlOADKXU5/aq8koTpVQK6AYN0NhmPeVxO/CZ3SJcEZGRwBGl1A67tZSHMRqlKZ0kPEgtv4jUBj5EZ0TMtFtPESJyPZCmlNpit5YyCAV6ALOVUpcD2dg/rFKMNS8wCmgDNAOiReQ2e1XVfKxMogXAf+3WUoSI1AIeA56wW4svGKNRmmSghcv3eIJgWMAdEQlDG4z/KqUW263Hjf7ASBE5iB7eu0ZE/mOvpFIkA8lKqaIe2iK0EQkWrgV+UkodV0rlA4uBfjZr8sYxEWkKYP1Ns1mPR0TkD8D1wDgVXC+otUM3DnZY90w8sFVE4mxV5QVjNEqzCWgvIm1EJBw9+bjUZk0lEBFBj8XvVkpVLYF5NaCUmqKUildKtUafv6+UUkHVSlZKpQKHRaSDVTQYSLJRkjuHgD4iUsv6vQcTRBP1biwF/mD9/wfgYxu1eEREhgGPAiOVUmft1uOKUmqnUqqxUqq1dc8kAz2sazToMEbDDWuy7H5gBfomXaiU+t5eVaXoD/wO3YLfbn2G2y2qBvIA8F8RSQS6A/+0WU8xVg9oEbAV2Im+V20PNSEi7wLrgQ4ikiwifwamAUNE5Ee098+0INT4ClAH+MK6X+YEmb4agwkjYjAYDAafMT0Ng8FgMPiMMRoGg8Fg8BljNAwGg8HgM8ZoGAwGg8FnjNEwGAwGg88Yo2EwGAwGnzFGw+AV0eyw3qQtKhstIn/0sO48EdkcUIEBJtiPUUSeEJEjIuIUkXl267ED99/IX7+ZiHS2QpYPdCl7VUTeqGrdNQ1jNAxlMRqoD7zjVvZHW9QYvCIiPYF/oF9i6w9MtVfRRcF0YJyIXGK3kEBijIahLB4EFlixjwzViJWwKrwKVXS0/r6qlFqvlNrvD11244fzUm0opQ4C64B7bZYSUIzRqGGIiENEskVkgoi8YGUAOyUik6zlvxORJBE5IyKLRSSqkvu5BB0gb5FL2TzgFuCXVlddiUiC23ZDrOxo2SKyTkQuc1t+lYisEZGzInJSRF4TkTrlaJknIpt9qHu1iCxyKxto6ezsVtcI6zydFZFlorPPXSI6I2K2tU5XL3puFJ0FLsfScanb8nKP0UXHjSLyPZAD9C7jHIwWkZ0ikisih0XkadF5Nop+lwXWqhnuwygVqauC59uX47xMRJaLSLpVz24Rua8MbV7PS3m6fcVH3eOtfWSLyCdAUy/VfYjubVw8z1KllPnUoA9wCTpU+0HgGXQ01P9YZa8Ay9DJj+5H57IYX8n93A2cARwuZe2Ar9DxkPpYn3hr2Tx0dNPtwBhgJDpD2vecD1fTH8gF3geGo+NnHQEWlaOl3Lqt9Va71wUMtM5NZ7e6tgA3A7cBp9DGcbN13L+y9pXkVv88dA6OA8A4a/ud6PwrkRU5RquuE9Zx3IaO2RTv5fivs45hPjrj2yPWPua4/C5TrXUGWb9L3crU5e/fEtiPviaHo4Mujgcml/Nblzovvuh22X5zGd/L1Y0OSa+A2cBQdEyyw1bZQLf9dbPKu9n9bAjUx3YB5lPBH0w/qBTwoEtZK6tspdtDLhGYWcn9zAU2eShfBKz2UD4PnaegvUvZjZaujtb3r4FVbttdg8tD3YuWcuu2ylZ7eGgNpLTRKADauazznLXO713Khltlndx0KKCf27kvAO6pyDG61NXdh99ig4c6HwEKOW+0/2jVV9sPdfnltwQaWt+7VOC683hefNHtsn1ZRsMX3d8Bn7mt8xqejUaoda7urMx9VhM/F0+X6sKhC3Aa3Qoqoijj3DRlXcku5ekAIjJbtGeNrxEq49A2QoLeAAAEF0lEQVQtvopwUCn1o8v3olDj8aITzfQFFopIaNEHPSacD1xR2borqLGoLtcx/33W3688lLlnbUxTSn1b9EUp9TO619KrEsd4RCm1vSyhonPW9wA+cFv0Pnp4uW9Z21ehLn/8lunoFvocERkjIr5m9CtxXvx1DnzRbe3rckqHd/eYs0bpqNin0ffLRYExGjWPLsA6VXJyuiu6tbO2qMC6QVoDu6yid6lYkqFIdDe+Irjn2M5zqas+EALMQt+gRZ9cIIySia8qWndF8VbXaQ9l7vV7SjCUhh7zrugxHvNBa0NrW/d1i77H+lBHZeqq8m+plHKih5VSgTeBVBH5WkQuL0enuz5/nQNfdDdC9x7cf+eyEkvlUrnrsEZS4Ukkg+10RbewXOkG7FFKuT7ku6AbBYlQnMweEU/ZbD2Sjn9bT6fR3fsE4H8elvsjO2IO4O5pU5GHqi94ai03Ro/3V/QYfen1nUA/2Nz328T6m+5DHf6uy+fjVErtAW4RnWnyauBZYJmIxFtGxRPu5yWQuo+jG2Du+yqrlxRTAQ01HmM0ahCiPaHaAe7J57t6KctGT0RWhh/w3O3PoxKtKqVUtohsADoopZ6spKbySAYGuJUN8fM+GotIv6IhKhFpie7BvVUdx6iUKhSRLcCtlBySHI12dFgf6Loqc5xWz/grEXkB/d6Pzw/aQOsWke3oyXDXRE03e1m3EVALPXF/UWCMRs3iMnTvwZOBeMlD2fdltObK4xvgCRFppJQ67lK+BxglIjeiH9JHlVK+9hIeAVaKiBM9oZ4FtER7ez2mlKrqjfcR8GcReRHtsTMI7f3iT04AC0TkceAc8CR66GKetbw6jvHvwAoReQudc70L2lvqNaVUsk11lXuclsvy8+ie8QH08NCjwA6lVEVb5gHTjfaWWiwis9HX1C/RHlue6InuvXzrZfkFh5nTqFl0Ac7i0nsQkfroyeBEt3W7eiirCKvRLUH3m2UW8Dl6jHoTcJevFSql1qF7Ao3Q7xV8gr6JD+Pb+H559S8D/gr8Gn2ztwIeqmq9bvwMTEIPcbwHZAJDlVI5lga/H6NS6nN0rvWeVn0PATPQbtW21OXjcaZa/z8GfIa+dnajXXiDVrdS6iN0KuAbgCXoiXFvKVmHAWuUUicrdkQ1F5Pu9SJDRJRSyqeJDRGZCVyilBpRzbIMhhqH5Wn1M/q9k//YrSdQmJ7GRYKIvC4iydb/ySLyug+bTQcGisgvqledwVAjuRU9RPme3UICielpGMpERMYCKUqpNXZrMRiCCRH5DfqdkrXlrnwBYYyGwWAwGHzGDE8ZDAaDwWeM0TAYDAaDzxijYTAYDAafMUbDYDAYDD5jjIbBYDAYfMYYDYPBYDD4jDEaBoPBYPAZYzQMBoPB4DP/D0EKGOhzeSB0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 15 # max number of ones\n",
    "m1s = arange(n+1)\n",
    "N = 35\n",
    "\n",
    "for w1 in [1,2,3]:\n",
    "    probs = []\n",
    "    for m1 in m1s:\n",
    "        # constrain so always sums to N\n",
    "        ms = [m1,5,5,5,5,n-m1]\n",
    "        assert sum(ms) == N\n",
    "        probs.append(multinomialpmf(ms,[w1,1,1,1,1,1]))\n",
    "    plot(m1s, probs, 'o--', label='\"one\" of weight %s' % w1)\n",
    "    \n",
    "xlabel(\"$m_1$ (the number of ones rolled)\", size=15)\n",
    "legend(loc='upper left')\n",
    "title(\"Probability of rolling $m_1$ ones\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can plot the distribution over a magical 3-sided die in 2D with a colourmap (only showing two variables on the axes, since the third can be inferred as they sum to N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'pmf of \"three sided\" die with weights [2,1,1]')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUgAAAEWCAYAAAAEkA60AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2de5xdRZXvvz8SQkAkCFEeCRKcBMfgoAxcYC5cB0Ex4CM4gwP4RBm5zoCC6FUYFSPKvXCvAio4TBQQQQzIgEaNPBQigyKQAAIJRGII0hKBQIg8TLCT3/2jqpNzdp/H7s7p7nNOr+/nsz+9q3ZV7bXP2b1OVa2qtWSbIAiCoD+bjbQAQRAE7UooyCAIgjqEggyCIKhDKMggCII6hIIMgiCoQyjIIAiCOrSVglTiEkmrJN0xiPrvlPSopOck7bUJcsyX9M+DrT8USHplfq4xda7PknT5INsedN1abTSTtRVIulDS58rIM1xI+qmkD5QsOyTvmKSDJK3Pn/+MVrc/CHl+J+nF4f4uWkVbKUjgQODNwGTb+w6i/peBE21vbfvu4kVJzn/nSzoonw/7P9JgsP37/FzrhvI+kqZIWp7Pl0uaMtA2hkNW2x+x/UXYoBR6hupeA5DpMNuXbmo7+TuwpLGDbOKx/Plfl9t7q6RbJT0j6Y+SvinppQ3u/0VJ90nqlTSriaxvlHSzpNV9700ltv8K+N+DfI4Rp90U5K7ActvPb0L9RS2Upymb8BIHwXAxAfgSsDPwGmAy8P8alF8KfAr4SYm2nwcuBv7XJsrYlmyygsy9jNMkLc5D40skjc/XDpLUI+lTkp6QtELSEZIOl/RbSU9L+rdc9jjgW8Df5eHBF2rcazNJn5X0SG7vO5ImSNpC0nPAGOA3kn5XUvYZwL8BR+V7/qbi8q6SfinpWUk3SJqY6/T9uh8n6ffATTl/f0m/yr/Sv+nroeZrEyRdlJ//D5K+1GCovK+kBZL+JOlxSecU7js2p3eT9Iss343AxEI7jeRpWLfkZ1e3jRqylnp+SeMl/bnis/5s7sVsk9NfknRePv92Tr8E+Cmwc/4On5O0c25yXH5HnpW0SNI+dZ7lC5K+ns83l/S8pP+b01tKWiPpZSU+1w3DZkljJH1F0kpJD0s6Uf17hTXfMeCW/PeZ/Dx/J2lq/rxX5zavLPtd2b7C9nW2X7C9CvgmcECD8pfa/inwbIm277B9GbCsrDwdhe1NOoDlwP3ALsB2wC+BL+VrBwG9wOnA5sCHgSeBK4CXAnsAa4BX5fLHArc2uNeHSL9urwK2Bq4BLqu4bmDqAOWfBVxeyJsP/A7YHdgyp8/K16bk+3wHeEm+Pgl4Cjic9KPz5px+ea7zA+A/cvlXAHcA/7OOPLcB78vnWwP7F+47tqLcOcAWwBtIL/Pl+VozeerWHcDn1uj+RVkH8vy3AP+Yz2/I38NhFdfemc+/XXjPemp8r2vyZzAG+D/Ar+vc82Dgvnz+3/M9b6+49puSn+t84J/z+UeAxaTe2suAnxU+k/k0f8fGVsj4PeAz+b7jgQPrPEu/z6JGmfOAOSW+48uBWSXfhzeRRn+l/sc65WjVEPt824/afho4Ezim4tpfgDNt/wWYQ+ppfNX2s7YXkYbEe5a8z3uAc2wvs/0ccBpwtIZmmHuJ7d/a/jNwFfD6wvVZtp/P198LzLM9z/Z62zcCC4DDJe0AHAacnMs/AZwLHF3nvn8BpkqaaPs5278uFpD0SuC/AZ+zvdb2LcCPKoo0kqdZ3aYMpI1BPP8vgL/P3+mewNdyeny+538NQNRb82ewDrgMeF2dcrcB0yRtT1L2FwGTJG0N/H2WCRp8rjXa/CfSe97j1Gs7q0aZZu9YJX8hTSHtbHuN7VsbPXg9JL0Z+ACp0xI0oVUK8tGK80dIcx19POWNk/V/zn8fr7j+Z1JPqQw75/Yr7zUW2KG8qKX5Y8X5C/SXsfKZdwXelYddz0h6hmRw2ilf2xxYUXHtP0g9qVocR+pVPCjpTklvq1FmZ2CVq+dqKz+XRvI0q1uGgbQx0Of/BakX9LfAfcCNJCW1P7DU9soByFn8DsfX+jHNCmpBvs8bsgy/Ig1DKxVko8+1yM5UvyOP1ijT7B2r5FOAgDvydMGHGpStiaT9SaO3I23/dqD1RyOt6nntUnH+SuCxFrVb5DHSS1p5r16qFe5AGaw7o8p6j5KG+h8uFpK0E7AWmGi7t2mj9kPAMZI2A/4BuDr3bCpZAbxM0ksqlNQrK2RqJM+uTeqWodn9K3mUATw/STG9Gngn8Avbi3OP9a1sVFRFWuGS6hek4fRewJ05/RZgXzbOCdb9XGuwgjS87mOXegVr0O95bP+RNEWFpAOBn0m6xfbSMg0qLXubC3zI9s8HIMuoplU9yBMkTZa0HcnoUXoCeYB8D/h4NhBsTVo+cGXJf7x6PA5MyQppsFwOvF3SW/Lk/HglA9Vk2ytIc2lfkbSNkqHpryT9fa2GJL1X0sttrweeydlVy2VsP0Lq8XxB0rj8D/P2kvI0q9uUgbQx0Oe3/QKwEDiBjQrxV8D/pL6CfBzYXtKEgTxHgV8A7wcW236RPJ8IPGz7yVym7udao72rgJMkTZK0LfDpAcjyJLCeNNcOgKR3VdxnFUmJllpGJem1wHXAR233mwpRWuo2vyK9eZ7S2AwYm59zTL7WZ4CbktOb5bKbp6TGSxo3gGdta1qlIK8g/RMsy8eXWtRukYtJc0m3AA+TJuE/uoltfj//fUrSXYNpwPajwEzSj8OTpJ7G/2Lj5/t+YBxp0n4VcDW1h2UAM4BFSlb5rwJH215To9y7gf2Ap4HPk4xGZeWpW3cADKSNgTw/JGW1OcmY05d+KRt7clXYfpD047ksD313rlWuCb8iGUv67rGY9H5tuGeJz7WSb5L+J+4F7gbmkUY7TZVa/pE4E/hlfp79SfOvt+f3Yi5wku2HSz7bJ4CXAxdpo6W/cjncLiTjaqXsfybZEj6Tz99XUfYR4A85/YZ8fR5pFPHn/NxdgbKVafANpMWh/2z7Zy2RKAi6EEmHARfa3rVp4U27zxuA60nTGkfZvr5EnXuAQ2w/VaLsZ4Enbf9HSXmWkKz/V9ke8LzpSBMKMgiGAElbAm8k9aZ2AP6TtMzo5BEVLBgQ7baTJgi6BQFfIE0p3A08QCyt6Tg2uQcZBEHQrQzpPuI8/H6WNDHda3ufbOm+krRbYDnwT3khbRAEQVsxpD3IrCD3qVzcq7TH9WnbZ0k6FXiZ7YZLIKStDNtuSO+992CMlEEQ1GPhwoUrbb98U9qYKvmFkmVXwPW2R9wdWzNGwhPNTNJOCYBLSevNmqwR2xY4fkNqwYJZQyBWEIxeJA10N1U/XiAtVi3DrEE4SBkJhtpIY+AGSQsl9Wm4HfLi4b5FxDW3nEk6XsmrzYL00QdB0M6I1OMqc3QKQy3rAbYfk/QK4EZJD5ataHs2MBtA2jksSUHQ5mxGWmnfTQypgrT9WP77hKRrSftaH5e0k+0VeZ/yEwNtt5aTY7t/XhAEw4dI25+6iSEbYkt6ibJbdyWnpoeS/EbOJblbIv/94VDJEATB8BFD7IGxA3CtpL77XGH7Okl3AlcpeRD/PfCuIZQhCIJhoht7kEOmIG0vo4aD0rzf85Chum8QBCNDXw+ym+i25wmCYISIHmQbUzTchNEmCIaXsGIHQRDUIXqQQRAEDeg2hdJtzxMEwQgRPcgOIuYkg2B4CSt2EARBHbrRSBMexYMgaAl9Q+wyR6n2pBmSlkhaml0jFq9vIenKfP32ikiL20u6OQcnO79QZ5yk2ZJ+K+lBSf/YSIboQQZB0BJaOcTOYWYvAN4M9AB3Sppre3FFseOAVbanSjoaOBs4ihSN8nPAa/NRyWeAJ2zvnkM9b9dIjlGjIGNOMgiGlhYbafYFluYdeUiaQ/IlW6kgZwKz8vnVwPmSZPt54FZJU2u0+yHgrwFy7PmVNcpsIIbYQRC0hAE6q5jY5+81H8cXmptEijveR0/Oq1nGdi+wGti+rnxSX1iCL0q6S9L3Je3Q6JlGTQ8yCIKhZYA9yJW292nSXJGiX9gyZSoZC0wGfmn7FEmnAF8G3teoQhAEwSYjWmrF7gF2qUhPBh6rU6ZH0lhgAvB0gzafIoUnuDanv0+ax6zLqFWQ4XQ3CFqLgM3LapTepiXuBKZJ2g34A3A08O5CmT7fsrcBRwI3uUEUQtuW9CNSTKybSF7FFtcrD6NYQQZB0FokGNsiBWm7V9KJwPXAGOBi24sknQEssD0XuAi4TNJSUs/x6I2yaDmwDTBO0hHAodkC/ulc5zzgSeCDjeQIBRkEQUuQYPMxrWvP9jxgXiHv9IrzNdRxuG17Sp38R4A3lJUhFGQQBC1hQD3IDqHLHicIgpFCgs23GGkpWksoyApiMXkQbAJd6K2iyx4nCIIRIxRkEARBA7pMo3TZ4wRBMGKItCCniwgF2YCYkwyCARBD7CAIgjoICCt2EARBDaIHGQRBUIdQkKObmJMMgiaEkSYIgqAG0YMMgiCoQyjIIAiCOoQVO6gknO4GQQXRgwyCIKhDFyrIiGoYBEFr6NtqWOYo05w0Q9ISSUslnVrj+haSrszXb5c0JedvL+lmSc9JOr9O23Ml3d9MhlCQQRC0hgHGfW3YlDQGuAA4DJgOHCNpeqHYccAq21OBc4Gzc/4a4HPAJ+u0/Q/Ac2UeacgVpKQxku6W9OOc3i1r+4ey9h831DIEQTAM9BlpyhzN2RdYanuZ7ReBOcDMQpmZwKX5/GrgEEmy/bztW0mKslpEaWvgFOBLZYQYjh7kScADFemzgXNtTwNW0STsYqchzao6gmDUMLAe5ERJCyqO4wutTQIerUj35LyaZWz3AquB7ZtI+UXgK6Twr00ZUgUpaTLwVuBbOS3gYJK2h6T9jxhKGYIgGCYGpiBX2t6n4phdo7UixZCuZcpsLCy9Hphq+9p6ZYoMtc3pPOBTwEtzenvgmaztofavAgD5FyX/qkwYUiGDIGgRrdMoPcAuFenJwGN1yvRIGktSFE83aPPvgL1zSNixwCskzbd9UL0KQ9aDlPQ24AnbCyuzaxStqfFtz+77dYGthkTGIAhaSGut2HcC07LNYhwp5vXcQpm5wAfy+ZHATbbr9iBt/7vtnXNI2AOB3zZSjjC0PcgDgHdIOhwYTwrifR6wraSxuRdZ61ehqwgHF8GooYXrIG33SjoRuJ6kUi+2vUjSGcAC23OBi4DLJC0l9RyP3iBK6iVuA4yTdARwqO3FA5VjyBSk7dOA0wAkHQR80vZ7JH2fpO3nkLT/D4dKhiAIhpEWbzW0PQ+YV8g7veJ8DfCuOnWnNGl7OfDaZjKMxDrITwOnZK2/PelXIAiCTqeF6yDbhWER1fZ8YH4+X0Za4xQEQTfRhVsNu+xx2p9wcBF0LaEggyAIGhAexYMgCGoQPcggCII6hMPcIAiCOkQPMhgKYjF50BWEggyCIKhD31bDLiIUZBAErSF6kEEQBHUQyetCFxEKsg2JOcmgI4khdhAEQR1iiB0EQdCALtMoXfY4QRCMGDHEDkaCmJMMOoIuHGJHXOwgCFpDa8O+ImmGpCWSlko6tcb1LXLo6KU5lPSUnL+9pJslPSfp/IryW0n6iaQHJS2SdFYzGUJBBkHQGlroMFfSGOAC4DBgOnCMpOmFYscBq2xPBc4lhZSGFA/7c8AnazT9Zdt/DewFHCDpsEZyhIIMgqA1tNaj+L7AUtvLbL9ICtEys1BmJil0NKRQ0odIku3nbd9KUpQbsP2C7Zvz+YvAXaS4WHXpshmD0UE43Q3akoHNQU6UtKAiPbsQG3sS8GhFugfYr9DGhjI5yNdqUhiXlU1FlbYF3g58tVG5UJBBELSO8lbslSmkc13KhIguHUa6qlKKof094Gs5BExdQkEGQdAaWmvF7gF2qUjXChHdV6YnK70JpPCvzZgNPGT7vGYFYw4yCILW0For9p3ANEm7SRpHink9t1BmLil0NKRQ0jfZbtiDlPQlkiI9uYwQTfW9pL8CemyvzfGt9wS+Y/uZMjcIgmCU0MIeZJ5TPBG4njRwv9j2IklnAAtszyWFjL4sh5B+mqREkyjScmAbYJykI4BDgT8BnwEeBO6SBHC+7W/VfaQmChdJ9wD7AFOysHOBV9s+fBDPPSiknQ3HD9ftuoIw2gQDQdLCJnOCTdnndfKC60reb2c2+X7DQRl9vz5r83cC59n+uqS7h1qwIAg6jC7cSVPmcf4i6RjSWP/tOW/zoRMpCIJOxaNwL/YHgY8AZ9p+WNJuwOVDK1YQBJ2GN4MXR5vDXNuLJX0S2F3Sa4EltpvuYQxGlnBwEQw3FvSOKbswZv2QytIqylixDyJt51lOmmXYRdIHbN8ytKIFQdBJWGLd2LKTkC8OqSytoszTfAU41PYSAEm7k1ah7z2UggVB0HmsG9Ndk5BlFOTmfcoRwPZvJYWRJgiCKoxY12Uec8soyAWSLgIuy+n3AAuHTqRgKIg5yWCoMaJ3FCrIfwFOAD5GmoO8BfjGUAoVBEHnYcSLZb3hdghlFOQBwIW2zxlqYYIg6FxG6xD7WOBCSU8B/5WPW22vGkrBgiDoPEadgrT9fgBJO5M8ZlwA7NysrqTxpOH4Frns1bY/nxeazwG2I3n0fV/27hsMI+F0N2g1o3IOUtJ7gf8B/A3JU+/5pF5kM9YCB9t+Llu9b5X0U+AU4FzbcyRdSIor8e+DfYAgCNqDNMTurs3YZZ7mPOB3wIXAzbaXl2k4+2V7Lic3z4eBg4F35/xLgVmEggyCjicZacaNtBgtpem+INsTgQ8B44EzJd0h6bIm1YAUmSy7S3sCuJGkaJ+x3ZuL9JDiStSqe7ykBSluxQtlbhcEwQhioJcxpY5OocwQexvglcCuJJ+QEygR9wHA9jrg9TlAzrXAa2oVq1N3Nsk1evYHGQRBe9N9Q+wyO8tvJbk5uxc4yvar+ww3Zcnex+cD+wPb5vgRUDvORDBCSLOqjiAYCH3LfMocZZA0Q9ISSUslnVrj+haSrszXb5c0JedvL+lmSc9JOr9QZ29J9+U6X1N2K16PMgryS7b/1fYVtnvyTd5V4uFennuOSNoSeBPwAHAzyRoOycfkD0vIEARBB9AqBSlpDGnFzGHAdOAYSdMLxY4DVtmeCpwLnJ3z1wCfAz5Zo+l/J4UnmJaPGY3kKKMg+2lu4LQS9XYCbpZ0LykAz422fwx8Gjglx5HYnhRXIgiCDqfFPch9gaW2l+VlgHOAmYUyM0mGXoCrgUMkyfbztm8lKcoNSNoJ2Mb2bdmI/B3giEZC1J0wkHQYcDgwSdLXKi5tA/TWrrUR2/cCe9XIX0Z6+CAIuggj1pbfajgxGWA3MDvbHfqYBDxake4B9iu0saFMDguzmtTpWlnnnpNyO5Vt1jQS99FoRvUxYAHwDqqdUzwLfLxRo0F3EA4ugoEwwK2GK5sE7ao1N1g01pYpsynl6ytI278BfiPpCtt/adRIEARBi/di9wC7VKRrGXT7yvRkw+8EUvjXRm1ObtJmFWXWQYZyDIKgFC1cB3knME3SbpLGkWJezy2UmUsy9EIy/N7kBnGsba8AnpW0f7Zev58mRuLuWrQUBMGI0cqthnlO8UTgemAMcLHtRZLOABbYnksy8F6WDb5Pk5QoAJKWk+wl4yQdQYqKsJjkvvHbwJbAT/NRl0ZGmstsv0/SSba/OvhHDbqFcHARNKLV7s5szwPmFfJOrzhfA9Rccmh7Sp38BcBry8rQSN3vLWlX4EOSvkNhgtN2o7F+EASjjGTF7q692I0U5IXAdcCrSFbsSgXpnB8EQQCMMm8+tr8GfE3Sv9v+l2GUKQiCDmU0Osz9F0mvI/mEBLglLwIPgiDYwKgMuSDpY6S9i9fkrO9Kmm3760MqWdARxGLyoI9RqSCBfwb2s/08gKSzgduAUJBBEGxggFsNO4IyClLAuor0Ompv2QmCYBQzWnuQlwC3S7o2p48gPPAEQVCDUacgbZ8jaT5wIKnn+EHbdw+1YEFnEnOSo5dRGdUQwPZdpBCtQRAENRlV6yCDIAgGyqgbYgdBEJShG8O+NlSQOS7E9bbfNEzyBF1GzEmOHkbdHKTtdZJekDTB9urhEioIgs5jtM5BrgHuk3Qj8Hxfpu2PDZlUQRB0JKNxDvIn+QiCIKjLqFwobvvSHNf6lbaXDINMQVM2H6J2hz66Rjjd7V5aPQcpaQbwVZJH8W/ZPqtwfQtS6Na9gaeAo2wvz9dOI8XNXgd8zPb1Of/jpO3TBu4jreuuCg9bSdOYNJLeDtxD8g2JpNdLKsaGCIJglJOs2FuUOpqRDcQXAIcB04FjJE0vFDsOWGV7KnAucHauO50UfmEPYAbwDUljJE0CPgbsY/u1JMV7NA1oqiCBWaQ41s8A2L4H2K1EvSAIRhF9Q+wyRwn2BZbaXmb7RWAOMLNQZiZwaT6/GjgkB+OaCcyxvdb2w8DS3B6kUfOWOQriVmxqVEOgt4YFu2Es2SAIRicDUJATJS2oOI4vNDUJeLQi3ZPzapax3QusBravV9f2H4AvA78HVgCrbd/Q6HnKGGnul/RuYIykaaQu6q9K1AuCYBQxwDnIlbb3aXC9lsewYsesXpma+ZJeRupd7kYaEX9f0nttX15PiDIK8qPAZ4C1wPdIYRi/WKJe0DKKRpkyX9tgDDnFdnsL6aEx4sRi8u6gxesge4BdKtKT6T8c7ivTk4fME0jhX+vVfRPwsO0nASRdA/x3YPAK0vYLwGeyo1zbfrZZnSAIRh8t3mp4JzBN0m7AH0jGlHcXyswFPkBy4H0kcJNtZyPyFZLOAXYGpgF3AOuB/SVtBfwZOARY0EiIMiEX/htwMfDSnF4NfMj2wpIPGgTBKKCVy3xs90o6kTRiHQNcbHuRpDOABbbnkvzSXiZpKanneHSuu0jSVcBi0jDoBNvrSH5tryZ5JusF7gZmN5KjTH/4IuBfbf8XgKQDSU509xzoQwdB0N20cquh7XnAvELe6RXna4B31al7JnBmjfzPA58vK0OZp3m2TznmG9wqKYbZI0qt+cUtm6TLfNXFOcc/N0nXytv0ecqYk+xMRtVOGkl/m0/vkPQfJAONgaOA+UMvWhAEncSoUpDAVwrpym5prIMMgqAfo8bdme03bkrDknYh7ZPckWQ9mm37q5K2A64EpgDLgX+yvWpT7hUEwcizns1KbSPsJMpYsbcF3k9SaBvKl3B31gt8wvZdkl4KLMwu044Ffm77LEmnAqcCnx6c+N1Ks3WPxflFyIsMKti+yfVabRTnIP9USD9do04xr1hn09dSxpxk5zCahth9zAN+TfJ8sb5sw7ZXkLbzYPtZSQ+QtgDNBA7KxS4lzWeGggyCDme0zUH2Md72KZtyE0lTgL2A24EdsvLE9gpJr6hT53gg78+csCm3D4JgGDCjaA6ygsskfRj4MWm7IQC2a423+iFpa+A/gZNt/yk522iO7dnkRZzSzmEUCoK2Z3SGXHgR+H+k/dh9isrAq5pVlLQ5STl+1/Y1OftxSTvl3uNOwBMDFzsIgnZjtA6xTwGm2l45kIazX7aLgAdsn1NxqW//5Fn57w8H0u7opIyziqJRZtdCekohvdUg5Hi8Rt5DhfTyQro40Ki1x2BghpvwSt6eGLF2NIV9zSwCXhhE2wcA7yMF/Lon5/0bSTFeJek4kl+2mluFgiDoLEZrVMN1wD2SbqZ6DrLhMh/bt1LbLxskLxpBEHQZo3GI/YN8BEEQ1GVUzkHavrRZmWC4KbNQfHJ1csfCnONrazQxuUZeJT079M9bUMh7ZlqhwF2F9PIaDRcXl9dyitGYWEw+8hixbv0oU5CSHqbG3mvbTa3YQRCMHrxerF0zyrYaApVxI8aTjCrbDY04QRB0KrZY19tdPcimUQ1tP1Vx/MH2ecDBwyBbEASdhGFd75hSRxkkzZC0RNLS7LeheH0LSVfm67fnHXt9107L+UskvaUif1tJV0t6UNIDkv6ukQxlhth/W5HcjNSjLE54BSNOcV5ym+pkcc7xbTWaeFN1cvzk6jWMa56p8bXfWlij+YPCnOTVhxUq3F7jxg8U0sX1ljEn2QnYovcvrelBShoDXAC8mRSE605Jc20vrih2HLDK9lRJRwNnA0dJmk4Kv7AHKSbNzyTtnsMufBW4zvaRksbRZEFwmSF2pV/IXrKLsjIPGQTBaEKsX9eydZD7AkttLwOQNIfk6KZSQc4EZuXzq4Hz8waVmcAc22uBh3PMmn0lLQLeQPIohu0XSTsF61LGir1JfiGDIBglGCg/BzlRUmVEwdnZ/0Ifk4BHK9I9wH6FNjaUyUG+VpO2lE0ieSCrrDuJNBR5ErhE0uuAhcBJtp+vJ2SZIfYWwD/S3x/kGc3qBkEwilgvWFO6B7nS9j4NrtfaZFJcTVOvTL38scDfAh+1fbukr5L80X6unhBlnuaHwGqStl3bpGzQLhS/2R0L6QP7Vzl0j+pt8QcVQg9NnNB/O/6SXV9dlb72Pe+sSi+7YI/qCucVOwHA0uKiiOI85R8K6U2fk4SYlxwSiv6RB08PsEtFejLwWJ0yPZLGkvwiPt2gbg/QY7vvBbuapCDrUkZBTrY9o0S5IAhGM8khZKu4E5gmaTfSL+TRwLsLZfoc39wGHAncZNuS5gJXSDqHZKSZBtxhe52kRyW92vYS0pbnxTSgjIL8laS/sX3fQJ4uCIJRRgsVZJ5TPBG4HhgDXGx7kaQzgAW255K8hV2WjTBPk5QoudxVJOXXC5yQLdgAHwW+my3Yy4APNpKjjII8EDg276hZSxrf2/aeA3vkIAi6GtOKsOgbm7PnkUK+VOadXnG+hjrewGyfCZxZI/8eqje/NKSMgiwuZAuCIOiP6TorRZllPo8MhyBBJcWf4eLXVMtIUchrNtSZ2P+nfm8WVqU/yCVV6R1vWN2kUfj4oedWpc86oToe2/n7fKp/pU8WHFzcWlz0/stCenmNO8di8hGntXOQbUF3ebcMgmDkCAUZBEFQh1CQQRAEdQgFGYwMxbeu1ltYDIZVcEL7x4LzipXFQGD0i/O140OFOccL+1f5y03V6Um7VTu4+Prp1XOQe76z/2qx4y+/uDrj2ILn3vkHFGRVUDwAAA+4SURBVO/aX5AhWEwec5KDIBRkEARBDdYDa0ZaiNYSCjIIgtYQQ+wgCII6hIIM2oNa82tPFdI91cn7p1enb+3fwsK9965K/3HahKr0jvRfB3lNIeuBe6rT//QP1ekPn395/xufUJ08/vzLqjOKc5ILaji86LdWMpzuDjuhIIMgCBoQCjIIgqAG0YMMgiCow3oGM5PR1oSC7Aia7c2G/usgC1vo/zilOv3j/rGKbnjTzKr0JXtUe4I67SPn9avzD4V1kGcW5iSvKpT//On048P7VM9L3rvf31Slzz+5sH+7uHcb4I9PFzKKn9kT/esM0PVMON1tgoF1TUt1FKEggyBoHTHEDoIgqEEXzkFuNtICBEHQJfQpyDJHCSTNkLRE0lJJ/WLHSNpC0pX5+u2SplRcOy3nL5H0lkK9MZLulvTjZjJEDzIIgtbQwq2GksYAFwBvJi3qvVPSXNuVMWSOA1bZnirpaOBs4ChJ00nhF/YgxaT5maTdK8IunAQ8ABQcFPQnFGRHUusnuOCcot9C6buqkz+rEdbwy9XJM86rjoY58dD+UQ0/fFW1gWVWtW8KFhcWjtfkF9XJd+z3o6r0DwqREnvm1zDSfGvfQkbRaFPLvFo0bA08XkAsJi/QuiH2vsBS28sAJM0BZlIdZGsmMCufXw2cL0k5f47ttcDDOWbNvsBtkiYDbyWFYzilmRAxxA6CoDUMbIg9UdKCiuP4QmuTgEcr0j05r2YZ272k8NTbN6l7HvApUn+3KdGDDIKgNQwsaNdK242CZ6nOHcqUqZkv6W3AE7YXSjqojJBD1oOUdLGkJyTdX5G3naQbJT2U/75sqO4fBMEw07cOsszRnB5gl4r0ZOCxemUkjQUmkOZW6tU9AHiHpOXAHOBgSTWcA2xkKHuQ3wbOB75TkXcq8HPbZ2Wr1KnAp2vUDRpS5me6OAe3vJAuBsYCvl3trGIN21Wlj59VcGwLLDn01VXpYtCu6dcW5Ci+4tA0COeUguw9+9SYg7yu0GnoKXj/7TcnC/3nJTc9ZumonpNs7TKfO4FpknYjeUM+Gnh3ocxc4APAbcCRwE22LWkucIWkc0hGmmnAHbZvA04DyD3IT9p+byMhhkxB2r6l0uyemQkclM8vBeYTCjIIugPTsq2GtnslnQhcD4wBLra9SNIZwALbc4GLgMuyEeZpkhIll7uKZNDpBU6osGAPiOGeg9zB9goA2yskvaJewTxpmyduJ9QrFgRBu9DirYa25wHzCnmnV5yvAd5Vp+6ZJEt1vbbnkzpoDWlbI43t2cBsAGnn4uRsEATtRhfupBluBfm4pJ1y73EnansQCAZFs/mzMh914e3+dsExbU//QF9fOfazVelrC2sWD3/nT6rSe9I/aFeRe6l2VrGcKU3rsG0h3bN9IWM7+lNcO1r87279nCR08bxkFyrI4V4H2TepSv77w2G+fxAEQ0XfMp8yR4cwZD1ISd8jGWQmSuoBPg+cBVwl6Tjg99SZPwiCoEMJd2flsH1MnUuHDNU9gyAYQSLsaxAEQR0GtpOmIwgF2bUMxmhTnGEvGDF+tmf/Kgt2qEou+8EeVenzD6pO11wUvmNB1ucKxqCiw4sFNdp4pkZeFf0NTP3ziv8OQ/Pf3rWLycOjeBAEQQO6zIodCjIIgtbQhct8QkEGQdAawkgTdC5l5tOKDh2KG2uLDjCAZwou+q4uOJK4unqOkqk1brtjYS6w+FY+V0j31Gjjj8XNVk8V0u1rPeiaOcnoQQZBEDQgFGQQBEENYplPEARBHWKZT9A91PqpL+YVx0u1nP01c8xbcBqxtIbTiH55NZz5DliOZmlo1+5Ox85JxhxkEARBHdbTMoe57UJENQyCoHW0LiYNkmZIWiJpaQ7RUry+haQr8/XbKyMYSDot5y+R9Jact4ukmyU9IGmRpJOayRAKMgiC1uGSRxMkjQEuAA4DpgPHSJpeKHYcsMr2VOBc4Oxcdzop/MIewAzgG7m9XuATtl8D7A+cUKPNKmKIHTSgOEdXa86uOKYqOqEtzv3V2hNdnHMsvpbFOrXkaObsttbYr5jXnhNoo8rp7kb2BZbaXgYgaQ4pptXiijIzgVn5/GrgfEnK+XNsrwUezjFr9s1Bu/pCvjwr6QFSvOzKNquIHmQQBCPBREkLKo7jC9cnAY9WpHtyXs0ytnuB1SSrYNO6eTi+F3B7IyGjBxkEwUiw0najgL+qkVccnNcr07CupK2B/wROtl0c8lQRCjIIghbRUjN2D7BLRXoy/aOq95XpkTSWFP706UZ1JW1OUo7ftX1NMyFiiB0EQYtoaVCaO4FpknaTNI5kdJlbKFMZ4+pI4CbbzvlHZyv3bsA04I48P3kR8IDtc8oIET3IYBNpZsgpY/gYjtewjBztuXC8Fu25mLx1K8Vt90o6EbgeGANcbHuRpDOABbbnkpTdZdkI8zRJiZLLXUUyvvQCJ9heJ+lA4H3AfZL63DD/W46/XRMlhdvepLjYxTncoDOoZbUuEgpyU9lUBSlpYZM5wRJtvM5wQ8nSO27y/YaD6EEGQdAius9bRSjIIAhahOm2vYahIIMhpkyPort6HSNBe8xJdp+3ilCQQRC0iBhiB0EQ1CF6kEEQBHWIHmQQBB3AyMxJRg8yCIKgDt3nMTcUZBAELSKG2EEQBA2IIXYQBB3G8DjdjR5kEARBHUJBBkEQ1CGs2EEQBHXoPiv2iDjMbRbOMQiCTqSlDnPbgmHvQVaEc3wzyTX6nZLm2q4bWSwIgtZTbbjZae9Nb7H7htgj0YPcEM7R9otAXzjHIAg6mtb2IJuNNHNIhSvz9dtzpMK+a6fl/CWS3lK2zSIjoSDLhHNE0vF9ISHhhWETLgiCwdLXgyxzNKZipHkYMB04RtL0QrHjgFW2pwLnAmfnutNJ4Rf2AGYA35A0pmSbVYyEgiwTzhHbs23vk9yybzUMYgVBsGn0GWnKHE0pM9KcCVyaz68GDsmBuWYCc2yvtf0wsDS3N+DR60hYscuEcyywYiV84RFgIrByyCRrLZ0ia6fICZ0ja6fICRtl3XXTm1pxPcyaWLLw+DQ63MBs27Mr0rVGmvsV2thQJgf5Wg1sn/N/XajbN0pt1mYVI6EgN4RzBP5A6gq/u1EF2y8HkLSgEwL9QOfI2ilyQufI2ilyQmtltT2jFe1kyow065Wpl19rxNwwauGwK8h64RyHW44gCNqaMiPNvjI9ksYCE0jhXxvVHdDodUTWQdqeZ3t3239l+8yRkCEIgrZmw0hT0jjSSHNuocxc4AP5/EjgJqc41nOBo7OVezdgGnBHyTar6LSdNLObF2kbOkXWTpETOkfWTpET2lTWeiNNSWcAC2zPBS4CLpO0lNRzPDrXXSTpKmAxyWR+gu11AAMdvSop3CAIgqDIiAyxgyAIOoFQkEEQBHXoCAXZzs4tJF0s6QlJ91fkbSfpRkkP5b8vG0kZ+5C0i6SbJT0gaZGkk3J+W8krabykOyT9Jsv5hZy/W95S9lDeYjZuJOWsJO/UuFvSj3O6LWWVtFzSfZLu6VuH2G7ffzvR9gpyMNuDhplvk7YzVXIq8HPb04Cf53Q70At8wvZrgP2BE/Jn2W7yrgUOtv064PXADEn7k7aSnZvlXEXaatYunAQ8UJFuZ1nfaPv1Fesf2+37bxvaXkHS5s4tbN9CsqBVUrkF6lLgiGEVqg62V9i+K58/S/qHnkSbyevEczm5eT4MHEzaUgZtIGcfkiYDbwW+ldOiTWWtQ1t9/+1EJyjIUs4t2owdbK+ApJSAV4ywPP3Ink/2Am6nDeXNQ9Z7gCeAG4HfAc/Y7vN00E7vwXnAp0ibkSFtd2tXWQ3cIGmhpONzXtt9/+1CJ6yDLOXcIiiPpK2B/wROtv2n1OFpL/K6tddL2ha4FnhNrWLDK1V/JL0NeML2QkkH9WXXKDrismYOsP2YpFcAN0p6cKQFamc6oQc5COcWI87jknYCyH+fGGF5NiBpc5Jy/K7ta3J228pr+xlgPmnOdNu8pQza5z04AHiHpOWk6Z+DST3KdpQV24/lv0+Qfnj2pY2//5GmExTkgLcHtQGVW6A+APxwBGXZQJ4buwh4wPY5FZfaSl5JL889RyRtCbyJNF96M2lLGbSBnAC2T7M92fYU0rt5k+330IaySnqJpJf2nQOHAvfTZt9/W2G77Q/gcOC3pHmoz4y0PAXZvgesILlJ7iFZK7cnWQMfyn+3G2k5s6wHkoZ69wL35OPwdpMX2BO4O8t5P3B6zn8VaU/tUuD7wBYj/ZkW5D4I+HG7yppl+k0+FvX9L7Xb999OR2w1DIIgqEMnDLGDIAhGhFCQQRAEdQgFGQRBUIdQkEEQBHUIBRkEQVCHUJBBPyTNlzTkQackfSx7FvruUN8rCAZDJ2w1DDoISWO9cQ9yM/4VOMwpdnEQtB3Rg+xQJE3Jva9vZp+JN+RdJ1U9QEkT8zY4JB0r6QeSfiTpYUknSjol+zH8taTtKm7xXkm/knS/pH1z/Zdk/5d35jozK9r9vqQfATfUkPWU3M79kk7OeReSFi7PlfTxQvnxki7JfgvvlvTGivtcI+m67Lvw/1bUOVTSbZLuyrJsnfPPkrRY0r2Svtyijz8YLYz0SvU4BncAU0j+HV+f01cB783n84F98vlEYHk+P5a0s+OlwMuB1cBH8rVzSc4r+up/M5+/Abg/n//vintsS9rd9JLcbg81dmAAewP35XJbk3Zw7JWvLQcm1qjzCeCSfP7XwO+B8fk+y0jhPccDj5D26U8EbgFekut8Gjgd2A5YwsbYS9uO9PcWR2cdMcTubB62fU8+X0hSms242ckX5LOSVgM/yvn3kbb49fE9SP4uJW2T90YfSnLM8MlcZjzwynx+o+2iX0xI2xuvtf08gKRrgP9B2kpYjwOBr+f7PyjpEWD3fO3ntlfnthYDu5KU9XTgl9kz0TjgNuBPwBrgW5J+Avy44ScTBAVCQXY2ayvO1wFb5vNeNk6fjG9QZ31Fej3V70NxD6pJbrz+0faSyguS9gOeryPjYHypNapTfOaxufyNto/p11CaHjiE5EjiRJK3nSAoRcxBdifLSUNb2OhRZqAcBSDpQGB17rVdD3w0ewVC0l4l2rkFOELSVtmDzDuB/ypR5z35HruTeqlLGpT/NXCApKm5zlaSds/zkBNszwNOJoVvCILSRA+yO/kycJWk9wE3DbKNVZJ+BWwDfCjnfZHk6/DerCSXA29r1IjtuyR9m+TZBuBbthsNrwG+AVwo6T5Sb/hY22vrOfa1/aSkY4HvSdoiZ38WeBb4oaTxpF7mx2s2EAR1CG8+QRAEdYghdhAEQR1CQQZBENQhFGQQBEEdQkEGQRDUIRRkEARBHUJBBkEQ1CEUZBAEQR3+P5iZ+bG9azMvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 50\n",
    "ms = arange(N+1)\n",
    "\n",
    "# default to NaN (so white background for region where sum>N)\n",
    "o = ones((N+1,N+1)) * NaN\n",
    "\n",
    "for m1 in ms:\n",
    "    for m2 in arange((N-m1)+1):\n",
    "        m3 = N - m1 - m2\n",
    "        assert (m1 + m2 + m3) == N\n",
    "        assert (m3 >= 0)\n",
    "        o[m2,m1] = multinomialpmf([m1,m2,m3],[2,1,1])\n",
    "                \n",
    "jet()\n",
    "imshow(o, origin='lower')\n",
    "xlabel(\"number of ones\")\n",
    "ylabel(\"number of twos\")\n",
    "colorbar()\n",
    "title('pmf of \"three sided\" die with weights [2,1,1]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}