4_Model-Building-with-PyMC.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:a46499385e3028806dd39835f649859ba5d1c43178162374a89aad7840169f37"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Building Models in PyMC\n",
      "\n",
      "Bayesian inference begins with specification of a probability model\n",
      "relating unknown variables to data. PyMC provides three basic building\n",
      "blocks for Bayesian probability models: `Stochastic`, `Deterministic`\n",
      "and `Potential`.\n",
      "\n",
      "A `Stochastic` object represents a variable whose value is not\n",
      "completely determined by its parents, and a `Deterministic` object\n",
      "represents a variable that is entirely determined by its parents. In\n",
      "object-oriented programming parlance, `Stochastic` and `Deterministic`\n",
      "are subclasses of the `Variable` class, which only serves as a template\n",
      "for other classes and is never actually implemented in models.\n",
      "\n",
      "The third basic class, `Potential`, represents 'factor potentials', which are *not* variables but simply\n",
      "log-likelihood terms and/or constraints that are multiplied into joint\n",
      "distributions to modify them. `Potential` and `Variable` are subclasses\n",
      "of `Node`.\n",
      "\n",
      "## The Stochastic class\n",
      "\n",
      "A stochastic variable has the following primary attributes:\n",
      "\n",
      "`value`\n",
      ":   The variable's current value.\n",
      "\n",
      "`logp`\n",
      ":   The log-probability of the variable's current value given the values\n",
      "    of its parents.\n",
      "\n",
      "A stochastic variable can optionally be endowed with a method called\n",
      "`random`, which draws a value for the variable given the values of its\n",
      "parents. \n",
      "\n",
      "### Creation of stochastic variables\n",
      "\n",
      "There are three main ways to create stochastic variables, called the\n",
      "**automatic**, **decorator**, and **direct** interfaces.\n",
      "\n",
      "**Automatic**\n",
      "\n",
      "Stochastic variables with standard distributions provided by PyMC can be created in a\n",
      "single line using special subclasses of `Stochastic`. For example, the uniformly-distributed discrete variable $switchpoint$ in the coal mining disasters model is created using the automatic interface as follows:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import pymc as pm\n",
      "import numpy as np\n",
      "from pymc.examples import disaster_model\n",
      "\n",
      "switchpoint = pm.DiscreteUniform('switchpoint', lower=0, upper=110)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Similarly, the rate parameters can automatically be given exponential priors:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "early_mean = pm.Exponential('early_mean', beta=1., value=4)\n",
      "late_mean = pm.Exponential('late_mean', beta=1., value=3)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Decorator**\n",
      "\n",
      "Uniformly-distributed discrete stochastic variable $switchpoint$ in the disasters model could alternatively be created from a function that computes its log-probability as follows:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "@pm.stochastic\n",
      "def switchpoint(value=1900, t_l=1851, t_h=1962):\n",
      "    \"\"\"The switchpoint for the rate of disaster occurrence.\"\"\"\n",
      "    if value > t_h or value < t_l:\n",
      "        # Invalid values\n",
      "        return -np.inf\n",
      "    else:\n",
      "        # Uniform log-likelihood\n",
      "        return -np.log(t_h - t_l + 1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that this is a simple Python function preceded by a Python\n",
      "expression called a **decorator**, here called\n",
      "`@stochastic`. Generally, decorators enhance functions with\n",
      "additional properties or functionality. The `Stochastic` object\n",
      "produced by the `@stochastic` decorator will evaluate its\n",
      "log-probability using the function $switchpoint$. The `value`\n",
      "argument, which is required, provides an initial value for the\n",
      "variable. The remaining arguments will be assigned as parents of\n",
      "$switchpoint$ (*i.e.* they will populate the `parents` dictionary).\n",
      "\n",
      "To emphasize, the Python function decorated by `@stochastic` should\n",
      "compute the *log*-density or *log*-probability of the variable. That\n",
      "is why the return value in the example above is $-\\log(t_h-t_l+1)$\n",
      "rather than $1/(t_h-t_l+1)$.\n",
      "\n",
      "**Direct**\n",
      "\n",
      "Its also possible to instantiate `Stochastic` directly:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def switchpoint_logp(value, t_l, t_h):\n",
      "    if value > t_h or value < t_l:\n",
      "        return -np.inf\n",
      "    else:\n",
      "        return -np.log(t_h - t_l + 1)\n",
      "\n",
      "def switchpoint_rand(t_l, t_h):\n",
      "    return np.round( (t_l - t_h) * np.random.random() ) + t_l\n",
      "\n",
      "switchpoint = pm.Stochastic( logp = switchpoint_logp,\n",
      "                doc = 'The switchpoint for the rate of disaster occurrence.',\n",
      "                name = 'switchpoint',\n",
      "                parents = {'t_l': 1851, 't_h': 1962},\n",
      "                random = switchpoint_rand,\n",
      "                trace = True,\n",
      "                value = 1900,\n",
      "                dtype=int,\n",
      "                rseed = 1.,\n",
      "                observed = False,\n",
      "                cache_depth = 2,\n",
      "                plot=True,\n",
      "                verbose = 0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Notice that the log-probability and random variate functions are\n",
      "specified externally and passed to `Stochastic` as arguments. This\n",
      "is a rather awkward way to instantiate a stochastic variable;\n",
      "consequently, such implementations should be rare.\n",
      "\n",
      "## Data Stochastics\n",
      "\n",
      "Data are represented by `Stochastic` objects whose `observed` attribute\n",
      "is set to `True`. If a stochastic variable's `observed` flag is `True`,\n",
      "its value cannot be changed, and it won't be sampled by the fitting\n",
      "method.\n",
      "\n",
      "In each interface, an optional keyword argument `observed` can be set to\n",
      "`True`. In the decorator interface, the\n",
      "`@observed` decorator is used instead of `@stochastic`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats.distributions import poisson\n",
      "\n",
      "@pm.observed\n",
      "def likelihood(value=[1, 2, 1, 5], parameter=3):\n",
      "    return poisson.logpmf(value, parameter).sum()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the other interfaces, the `observed=True` argument is added to the\n",
      "instantiation of the `Stochastic`, or its subclass:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "disasters = pm.Poisson('disasters', mu=2, \n",
      "                       value=disaster_model.disasters_array, \n",
      "                       observed=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## The Deterministic class\n",
      "\n",
      "The `Deterministic` class represents variables whose values are\n",
      "completely determined by the values of their parents. For example, in\n",
      "our disasters model, $rate$ is a `deterministic` variable."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "@pm.deterministic\n",
      "def rate(s=switchpoint, e=early_mean, l=late_mean):\n",
      "    ''' Concatenate Poisson means '''\n",
      "    out = np.empty(len(disaster_model.disasters_array))\n",
      "    out[:s] = e\n",
      "    out[s:] = l\n",
      "    return out"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "so `rate`'s value can be computed exactly from the values of its parents\n",
      "`early_mean`, `late_mean` and `switchpoint`.\n",
      "\n",
      "A `Deterministic` variable's most important attribute is `value`, which\n",
      "gives the current value of the variable given the values of its parents.\n",
      "Like `Stochastic`'s `logp` attribute, this attribute is computed\n",
      "on-demand and cached for efficiency.\n",
      "\n",
      "A Deterministic variable has the following additional attributes:\n",
      "\n",
      "`parents`\n",
      ":   A dictionary containing the variable's parents. The keys of the\n",
      "    dictionary correspond to the names assigned to the variable's\n",
      "    parents by the variable, and the values correspond to the actual\n",
      "    parents.\n",
      "\n",
      "`children`\n",
      ":   A set containing the variable's children, which must be nodes.\n",
      "\n",
      "Deterministic variables have no methods.\n",
      "\n",
      "### Creation of deterministic variables\n",
      "\n",
      "Deterministic variables are less complicated than stochastic variables,\n",
      "and have similar **automatic**, **decorator**, and **direct**\n",
      "interfaces:\n",
      "\n",
      "**Automatic**\n",
      "\n",
      "A handful of common functions have been wrapped in Deterministic\n",
      "objects. These are brief enough to list:\n",
      "\n",
      "`LinearCombination`:\n",
      ":   Has two parents $x$ and $y$, both of which must be iterable (*i.e.* vector-valued). This function returns:\n",
      "\n",
      "`Index`:\n",
      ":   Has two parents $x$ and `index`. $x$ must be iterables, `index`\n",
      "must be valued as an integer. `Index` is useful for implementing dynamic models, in which the parent-child connections change.\n",
      "\n",
      "`Lambda`:\n",
      ":   Converts an anonymous function (in Python, called **lambda functions**) to a `Deterministic` instance on a single line.\n",
      "\n",
      "`CompletedDirichlet`:\n",
      ":   PyMC represents Dirichlet variables of length $k$ by the first\n",
      "$k-1$ elements; since they must sum to 1, the $k^{th}$ element\n",
      "is determined by the others. `CompletedDirichlet` appends the\n",
      "$k^{th}$ element to the value of its parent $D$.\n",
      "\n",
      "`Logit`, `InvLogit`, `StukelLogit`, `StukelInvLogit`:\n",
      ":   Common link functions for generalized linear models, and their inverses.\n",
      "\n",
      "It\u2019s a good idea to use these classes when feasible in order to give\n",
      "hints to step methods.\n",
      "\n",
      "Certain elementary operations on variables create deterministic variables. For example:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pymc import Lambda\n",
      "x = pm.MvNormal('x', np.ones(3), np.eye(3))\n",
      "y = pm.MvNormal('y', np.ones(3), np.eye(3))\n",
      "sum_xy = Lambda('sum_xy', lambda x=x, y=y: x+y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(x[0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "x[0]\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(x[0]+y[2])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(x[0]_add_y[2])\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "All the objects thus created have `trace=False` and `plot=False` by default.\n",
      "\n",
      "**Decorator**\n",
      "\n",
      "We have seen in the disasters example how the decorator interface is used to create a deterministic variable. Notice that rather than returning the log-probability, as is the\n",
      "case for `Stochastic` objects, the function returns the value of the deterministic object, given its parents. Also notice that, unlike for `Stochastic` objects, there is no `value` argument\n",
      "passed, since the value is calculated deterministically by the\n",
      "function itself. \n",
      "\n",
      "\n",
      "**Direct**\n",
      "\n",
      "`Deterministic` objects can also be instantiated directly:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def rate_eval(switchpoint=switchpoint, early_mean=early_mean, late_mean=late_mean):\n",
      "    value = np.zeros(111)\n",
      "    value[:switchpoint] = early_mean\n",
      "    value[switchpoint:] = late_mean\n",
      "    return value\n",
      "\n",
      "rate = pm.Deterministic(eval = rate_eval,\n",
      "                  name = 'rate',\n",
      "                  parents = {'switchpoint': switchpoint, \n",
      "                          'early_mean': early_mean, \n",
      "                          'late_mean': late_mean},\n",
      "                  doc = 'The rate of disaster occurrence.',\n",
      "                  trace = True,\n",
      "                  verbose = 0,\n",
      "                  dtype=float,\n",
      "                  plot=False,\n",
      "                  cache_depth = 2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Containers\n",
      "\n",
      "In some situations it would be inconvenient to assign a unique label to\n",
      "each parent of some variable. Consider $y$ in the following model:\n",
      "\n",
      "$$\\begin{align*}\n",
      "x_0 &\\sim N (0,\\tau_x)\\\\\n",
      "x_{i+1}|x_i &\\sim \\text{N}(x_i, \\tau_x)\\\\\n",
      "&i=0,\\ldots, N-2\\\\\n",
      "y|x &\\sim N \\left(\\sum_{i=0}^{N-1}x_i^2,\\tau_y\\right)\n",
      "\\end{align*}$$\n",
      "\n",
      "Here, $y$ depends on every element of the Markov chain $x$, but we\n",
      "wouldn't want to manually enter $N$ parent labels `` `x_0' ``,\n",
      "`` `x_1' ``, etc.\n",
      "\n",
      "This situation can be handled naturally in PyMC:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N = 10\n",
      "x_0 = pm.Normal('x_0', mu=0, tau=1)\n",
      "\n",
      "x = np.empty(N, dtype=object)\n",
      "x[0] = x_0\n",
      "\n",
      "for i in range(1, N):\n",
      "\n",
      "    x[i] = pm.Normal('x_%i' % i, mu=x[i-1], tau=1)\n",
      "\n",
      "@pm.observed\n",
      "def y(value=1, mu=x, tau=100):\n",
      "    return pm.normal_like(value, np.sum(mu**2), tau)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Or, more succinctly using a list comprehension,"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = [pm.Normal('x_%i' % i, mu=x[i-1], tau=1) for i in range(1,N)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "PyMC automatically wraps array $x$ in an appropriate `Container` class.\n",
      "The expression `'x_%i' % i` labels each `Normal` object in the container\n",
      "with the appropriate index $i$. For example, if `i=1`, the name of the\n",
      "corresponding element becomes `` `x_1' ``.\n",
      "\n",
      "Containers, like variables, have an attribute called `value`. This\n",
      "attribute returns a copy of the (possibly nested) iterable that was\n",
      "passed into the container function, but with each variable inside\n",
      "replaced with its corresponding value.\n",
      "\n",
      "The Potential class\n",
      "-------------------\n",
      "\n",
      "For some applications, we want to be able to modify the joint density by\n",
      "incorporating terms that don't correspond to probabilities of variables\n",
      "conditional on parents, for example:\n",
      "\n",
      "$$\\begin{eqnarray*}\n",
      "p(x_0, x_2, \\ldots x_{N-1}) \\propto \\prod_{i=0}^{N-2} \\psi_i(x_i, x_{i+1}).\n",
      "\\end{eqnarray*}$$\n",
      "\n",
      "In other cases we may want to add probability terms to existing models.\n",
      "For example, suppose we want to constrain the difference between the early and late means in the disaster model to be less than 1, so that the joint density becomes:\n",
      "\n",
      "$$p(y,\\tau,\\lambda_1,\\lambda_2) \\propto p(y|\\tau,\\lambda_1,\\lambda_2) p(\\tau) p(\\lambda_1) p(\\lambda_2) I(|\\lambda_2-\\lambda_1| \\lt 1)$$\n",
      "\n",
      "Arbitrary factors are implemented by objects of class `Potential`. Bayesian\n",
      "hierarchical notation doesn't accomodate these potentials. \n",
      "\n",
      "Potentials have one important attribute, `logp`, the log of their\n",
      "current probability or probability density value given the values of\n",
      "their parents. The only other additional attribute of interest is\n",
      "`parents`, a dictionary containing the potential's parents. Potentials\n",
      "have no methods. They have no `trace` attribute, because they are not\n",
      "variables. They cannot serve as parents of variables (for the same\n",
      "reason), so they have no `children` attribute.\n",
      "\n",
      "### Creation of Potentials\n",
      "\n",
      "There are two ways to create potentials:\n",
      "\n",
      "**Decorator**\n",
      "\n",
      "A potential can be created via a decorator in a way very similar to\n",
      "`Deterministic`'s decorator interface:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "@pm.potential\n",
      "def rate_constraint(l1=early_mean, l2=late_mean):\n",
      "    if np.abs(l2 - l1) > 1:\n",
      "        return -np.inf\n",
      "    return 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The function supplied should return the potential's current\n",
      "*log*-probability or *log*-density as a Numpy `float`. The\n",
      "`potential` decorator can take `verbose` and `cache_depth` arguments\n",
      "like the `stochastic` decorator.\n",
      "\n",
      "**Direct**\n",
      "\n",
      "The same potential could be created directly as follows:\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def rate_constraint_logp(l1=early_mean, l2=late_mean):\n",
      "    if np.abs(l2 - l1) > 1:\n",
      "        return -np.inf\n",
      "    return 0\n",
      "\n",
      "rate_constraint = pm.Potential(logp = rate_constraint_logp,\n",
      "                    name = 'rate_constraint',\n",
      "                    parents = {'l1': early_mean, 'l2': late_mean},\n",
      "                    doc = 'Constraint on rate differences',\n",
      "                    verbose = 0,\n",
      "                    cache_depth = 2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Example: Bioassay model\n",
      "\n",
      "Recall from a previous lecture the bioassay example, where the number of deaths in a toxicity experiment was modeled as a binomial response, with the probability of death being a linear function of dose:\n",
      "\n",
      "$$\\begin{aligned}\n",
      "y_i &\\sim \\text{Bin}(n_i, p_i) \\\\\n",
      "\\text{logit}(p_i) &= a + b x_i\n",
      "\\end{aligned}$$\n",
      "\n",
      "Implement this model in PyMC (we will show you how to fit the model later!)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Log dose in each group\n",
      "log_dose = [-.86, -.3, -.05, .73]\n",
      "\n",
      "# Sample size in each group\n",
      "n = 5\n",
      "\n",
      "# Outcomes\n",
      "deaths = [0, 1, 3, 5]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pymc import invlogit\n",
      "\n",
      "def bioassay():\n",
      "    \n",
      "    a = Normal('a', 0, 0.01, value=0)\n",
      "    b = Normal('b', 0, 0.01, value=0)\n",
      "    \n",
      "    p = Lambda('p', lambda a=a, b=b: invlogit(a + b*log_dose))\n",
      "    \n",
      "    y = Binomial('y', n, p=p, value=deaths, observed=True)\n",
      "    \n",
      "    return(locals())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Fitting Models\n",
      "\n",
      "PyMC provides three objects that fit models:\n",
      "\n",
      "- `MCMC`, which coordinates Markov chain Monte Carlo algorithms. The actual work of updating stochastic variables conditional on the rest of the model is done by `StepMethod` objects.\n",
      "\n",
      "- `MAP`, which computes maximum *a posteriori* estimates.\n",
      "\n",
      "- `NormApprox`, the joint distribution of all stochastic variables in a model is approximated as normal using local information at the maximum *a posteriori* estimate.\n",
      "\n",
      "All three objects are subclasses of `Model`, which is PyMC's base class\n",
      "for fitting methods. `MCMC` and `NormApprox`, both of which can produce\n",
      "samples from the posterior, are subclasses of `Sampler`, which is PyMC's\n",
      "base class for Monte Carlo fitting methods. `Sampler` provides a generic\n",
      "sampling loop method and database support for storing large sets of\n",
      "joint samples. These base classes implement some basic methods that are\n",
      "inherited by the three implemented fitting methods, so they are\n",
      "documented at the end of this section.\n",
      "\n",
      "### Maximum a posteriori estimates\n",
      "\n",
      "The `MAP` class sets all stochastic variables to their maximum *a\n",
      "posteriori* values using functions in SciPy's `optimize` package; hence,\n",
      "SciPy must be installed to use it. `MAP` can only handle variables whose\n",
      "dtype is `float`, so it will not work, for example, on the disaster model example. \n",
      "\n",
      "We can fit the bioassay example using `MAP`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pymc.examples import gelman_bioassay\n",
      "M = pm.MAP(gelman_bioassay)\n",
      "M.fit()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This call will cause $M$ to fit the model using Nelder-Mead\n",
      "optimization, which does not require derivatives. The variables in\n",
      "`DisasterModel` have now been set to their maximum *a posteriori*\n",
      "values:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M.alpha.value"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 19,
       "text": [
        "array(0.6523226889885603)"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M.beta.value"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 20,
       "text": [
        "array(6.493563631375887)"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can also calculate model selection statistics, AIC and BIC:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M.AIC"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "8.041324910757734"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M.BIC"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "6.8139136329975152"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`MAP` has two useful methods:\n",
      "\n",
      "`fit(method ='fmin', iterlim=1000, tol=.0001)`\n",
      ":   The optimization method may be `fmin`, `fmin_l_bfgs_b`, `fmin_ncg`,\n",
      "    `fmin_cg`, or `fmin_powell`. See the documentation of SciPy's\n",
      "    `optimize` package for the details of these methods. The `tol` and\n",
      "    `iterlim` parameters are passed to the optimization function under\n",
      "    the appropriate names.\n",
      "\n",
      "`revert_to_max()`\n",
      ":   If the values of the constituent stochastic variables change after\n",
      "    fitting, this function will reset them to their maximum *a\n",
      "    posteriori* values.\n",
      "\n",
      "\n",
      "The useful attributes of `MAP` are:\n",
      "\n",
      "`logp`\n",
      ":   The joint log-probability of the model.\n",
      "\n",
      "`logp_at_max`\n",
      ":   The maximum joint log-probability of the model.\n",
      "\n",
      "`AIC`\n",
      ":   Akaike's information criterion for this model.\n",
      "\n",
      "`BIC`\n",
      ":   The Bayesian information criterion for this model.\n",
      "\n",
      "One use of the `MAP` class is finding reasonable initial states for MCMC\n",
      "chains. Note that multiple `Model` subclasses can handle the same\n",
      "collection of nodes."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Normal approximations\n",
      "\n",
      "The `NormApprox` class extends the `MAP` class by approximating the\n",
      "posterior covariance of the model using the Fisher information matrix,\n",
      "or the Hessian of the joint log probability at the maximum."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N = pm.NormApprox(gelman_bioassay)\n",
      "N.fit()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The approximate joint posterior mean and covariance of the variables are\n",
      "available via the attributes `mu` and `C`, which the the approximate posterior mean and variance/covariance, respectively:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N.mu[N.alpha]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 24,
       "text": [
        "array([ 0.65232269])"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N.C[N.alpha, N.beta]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 25,
       "text": [
        "matrix([[  0.77931584,   2.00872419],\n",
        "        [  2.00872419,  13.03545444]])"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As with `MAP`, the variables have been set to their maximum *a\n",
      "posteriori* values (which are also in the `mu` attribute) and the AIC\n",
      "and BIC of the model are available.\n",
      "\n",
      "We can also generate samples from the posterior:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "\n",
      "N.sample(100)\n",
      "plt.hist(N.trace('alpha')[:])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 26,
       "text": [
        "(array([  1.,   0.,   0.,   8.,  12.,  18.,  18.,  23.,  14.,   6.]),\n",
        " array([-2.61665925, -2.10355726, -1.59045526, -1.07735326, -0.56425127,\n",
        "        -0.05114927,  0.46195272,  0.97505472,  1.48815672,  2.00125871,\n",
        "         2.51436071]),\n",
        " <a list of 10 Patch objects>)"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEACAYAAACTXJylAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADZBJREFUeJzt3X+sZPVdxvHnYXcbQRq3LWYXcZvtH7XVhAhGq0lrmEQx\ni38U+oc0GBOs0TRagRg1XUrSvYkxljZUQjSNiQtZqcU0VNtuFcpKmBYSBYu7sPwq1ZRYKrukllY2\nxMqPxz/mLFxv752ZnTlzzn7mvl/JJGfOnHO+n5M797nf+z2/nEQAgDrO6LsAAMCpIbgBoBiCGwCK\nIbgBoBiCGwCKIbgBoJixwW17l+17bD9q+xHbVzfzV2w/bftw89rTTbkAAI87j9v2Tkk7kxyxfbak\nByVdJulySc8n+Xg3ZQIATto67sMkxyQda6ZP2H5c0nnNx15wbQCAdUw9xm17t6QLJf1zM+sq2w/Z\n3m97+wJqAwCsY6rgboZJbpd0TZITkj4h6S2SLpD0jKQbFlYhAOD/GTvGLUm2t0n6gqQ7kty4zue7\nJR1Mcv6a+dwEBQBmkGTsUPSks0osab+kx1aHtu1zVy32HklHN2h8aV/79u3rvQb2j/3bjPu3zPuW\nTNffHXtwUtI7Jf2apIdtH27mfUjSFbYvkBRJX5f0/qlaAwDMbdJZJfdp/V75HYspBwAwCVdOzmgw\nGPRdwkKxf7Ut8/4t875Na+LByZk3bGdR2waAZWVbmefgJADg9ENwA0AxBDcAFENwA0AxBDcAFENw\nA0AxBDcAFENwA0AxBDcAFENwA0AxBDcAFDPptq4AChndQr973JeoWwQ3sHS6DlGeG941hkoAoBiC\nGwCKIbgBoBiCGwCKIbgBoBiCGwCKIbgBoBiCGwCKIbgBoBiCGwCKIbgBoBiCGwCKIbgBoBiCGwCK\nIbgBoBiCGwCKIbgBoBiCGwCKIbgBoBiCGwCKIbgBoJixwW17l+17bD9q+xHbVzfz32j7kO0nbd9l\ne3s35QIAnGTjD+2dknYmOWL7bEkPSrpM0vskfSvJR21/UNIbkuxds27GbRtA+2xL6vr3zuJ3vT22\nlcTjlhnb405yLMmRZvqEpMclnSfp3ZIONIsd0CjMAQAdmHqM2/ZuSRdKul/SjiTHm4+OS9rRemUA\ngHVtnWahZpjkM5KuSfL86N+xkSSxve7/SSsrK69ODwYDDQaDeWoFgKUzHA41HA5PaZ2xY9ySZHub\npC9IuiPJjc28JyQNkhyzfa6ke5K8fc16jHEDHWOMu765x7g9+hbsl/TYydBufF7Slc30lZI+O0+h\nAIDpTTqr5F2SvizpYb32Z/xaSQ9I+rSkN0t6StLlSb6zZl163EDH6HHXN02Pe+JQyRyNE9xAxwju\n+uYeKgEAnH4IbgAohuAGgGIIbgAohuAGgGIIbgAohuAGgGIIbgAohuAGgGIIbgAohuAGgGIIbgAo\nhuAGgGKmegIOsAxWP7kJqIzgxibT/S1Pu22TP06bAUMlAFAMwQ0AxRDcAFAMwQ0AxRDcAFAMwQ0A\nxRDcAFAMwQ0AxRDcAFAMwQ0AxRDcAFAMwQ0AxRDcAFAMwQ0AxRDcAFAMwQ0AxRDcAFAMwQ0AxRDc\nAFAMwQ0AxUwMbts32z5u++iqeSu2n7Z9uHntWWyZAICTpulx3yJpbTBH0seTXNi87my/NADAeiYG\nd5J7JT23zkduvxwAwCTzjHFfZfsh2/ttb2+tIgDAWLMG9yckvUXSBZKekXRDaxUBAMbaOstKSZ49\nOW37LyUdXG+5lZWVV6cHg4EGg8EszQHA0hoOhxoOh6e0jpNMXsjeLelgkvOb9+cmeaaZ/j1JP5Pk\nV9esk2m2DXTFtkbH1TttteM2+9lHftfbY1tJxh5DnNjjtn2bpIsknWP7G5L2SRrYvkCjb8jXJb2/\nhXoBAFOYqsc904bpceM0Q497cW3yu96eaXrcXDkJAMUQ3ABQDMENAMUQ3ABQDMENAMXMdAEOAKw2\nOmOnW5v5TBaCG0AL+jjNcvNiqAQAiiG4AaAYghsAiiG4AaAYghsAiiG4AaAYghsAiiG4AaAYghsA\niiG4AaAYghsAiiG4AaAYghsAiiG4AaAYghsAiiG4AaAYHqSAXvTxxBRgWRDc6BFPTQFmwVAJABRD\ncANAMQQ3ABRDcANAMQQ3ABRDcANAMQQ3ABRDcANAMQQ3ABRDcANAMQQ3ABQzMbht32z7uO2jq+a9\n0fYh20/avsv29sWWCQA4aZoe9y2S9qyZt1fSoSQ/Junu5j0AoAMTgzvJvZKeWzP73ZIONNMHJF3W\ncl0AgA3MOsa9I8nxZvq4pB0t1QMAmGDug5NJou5vrAwAm9asD1I4bntnkmO2z5X07HoLraysvDo9\nGAw0GAxmbA4AltNwONRwODyldTzqME9YyN4t6WCS85v3H5X0X0mut71X0vYke9esk2m2jc1p9Oiy\nPp6As+xtboZ9HLW5rPliW0nGPq5pYnDbvk3SRZLO0Wg8+8OSPifp05LeLOkpSZcn+c6a9QhubIjg\nXpb2+mtzWfOlleCeo3GCGxsiuJelvf7aXNZ8mSa4uXISAIohuAGgGIIbAIohuAGgGIIbAIohuAGg\nGIIbAIohuAGgGIIbAIohuAGgGIIbAIohuAGgGIIbAIqZ9UEKWDKju/UBqIDgxipd334UwCwYKgGA\nYghuACiG4AaAYghuACiG4AaAYghuACiG4AaAYghuACiG4AaAYghuACiG4AaAYghuACiG4AaAYghu\nACiG4AaAYghuACiGBykAKKnrpzYlXT5oZDyCG0BRm/eJTQyVAEAxBDcAFENwA0Axc41x235K0n9L\nelnSi0ne0UZRAICNzXtwMpIGSb7dRjEAgMnaGCo5vQ63AsCSmze4I+kfbX/F9m+1URAAYLx5h0re\nmeQZ2z8s6ZDtJ5Lce/LDlZWVVxccDAYaDAZzNgcAy2U4HGo4HJ7SOm7raiDb+ySdSHJD8z6n05VG\nGG90FVrXFzR0/f3YDG1uhn3so013duWkbSUZOwQ981CJ7bNsv76Z/kFJvyTp6KzbAwBMZ56hkh2S\n/q65X8BWSX+d5K5WqgIAbKi1oZLv2zBDJaUwVLIsbW6GfeyjzSUZKgEA9IPgBoBiCG4AKIbgBoBi\nCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4A\nKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbgBoBiCG4AKIbg\nBoBitvZdwKm46aY/15e+9ECnbZ511lbdeuv+TtsEgHGcZDEbttP2ti+55L268843SfrZVre7sRe1\nZcsH9NJL3+uovf7YlrSY78IGLXbc3mZpczPsYx9tWovKyu9ryVYSj1umVI975CJJ7+2ore9J+kBH\nbQHAdBjjBoBiCG4AKGbm4La9x/YTtr9m+4NtFgUA2NhMwW17i6Q/k7RH0k9IusL2j7dZ2OluOBz2\nXcKCDfsuYMGGfRewYMO+C1igYd8F9G7WHvc7JP1bkqeSvCjpbyRd2l5Zpz+Cu7ph3wUs2LDvAhZo\n2HcBvZs1uM+T9I1V759u5gEAFmzW0wG7PmlTkrRli3Tmmddr27ZPdtTiy3rhhY6aAoApzXQBju2f\nk7SSZE/z/lpJryS5ftUyvYQ7AFQ36QKcWYN7q6SvSvoFSf8p6QFJVyR5fJYiAQDTm2moJMlLtn9X\n0hclbZG0n9AGgG4s7F4lAIDFWOiVk7b/yPZDto/Yvtv2rkW21zXbH7P9eLOPf2v7h/quqU22f8X2\no7Zftv1TfdfThmW+cMz2zbaP2z7ady2LYHuX7Xua7+Qjtq/uu6Y22f4B2/c3efmY7T/ZcNlF9rht\nvz7J8830VZJ+MslvLqzBjtm+WNLdSV6x/RFJSrK357JaY/vtkl6R9BeSfj/Jv/Zc0lyaC8e+KukX\nJX1T0r9oiY7N2P55SSck/VWS8/uup222d0rameSI7bMlPSjpsmX5+UmS7bOSvNAcR7xP0h8kuW/t\ncgvtcZ8M7cbZkr61yPa6luRQkleat/dL+tE+62lbkieSPNl3HS1a6gvHktwr6bm+61iUJMeSHGmm\nT0h6XNKP9FtVu5KcPAH5dRodP/z2esst/CZTtv/Y9n9IulLSRxbdXo9+Q9I/9F0ExuLCsSVhe7ek\nCzXqMC0N22fYPiLpuKR7kjy23nJz34/b9iFJO9f56ENJDia5TtJ1tvdK+lNJ75u3zS5N2r9mmesk\n/W+ST3VaXAum2b8lwpH4JdAMk9wu6Zqm5700mv/gL2iOl33R9iDJcO1ycwd3kounXPRTKtgjnbR/\ntn9d0i9rdE57Oafw81sG35S0+gD5Lo163SjC9jZJn5H0ySSf7bueRUnyXdt/L+mntc7NWRZ9Vslb\nV729VNLhRbbXNdt7JP2hpEuT/E/f9SzY2Cu5iviKpLfa3m37dRo9SunzPdeEKXn0fL39kh5LcmPf\n9bTN9jm2tzfTZ0q6WBtk5qLPKrld0tskvSzp3yX9dpJnF9Zgx2x/TaODCCcPIPxTkt/psaRW2X6P\npJsknSPpu5IOJ7mk36rmY/sSSTfqtQvHNjzlqhrbt2n0bL83SXpW0oeT3NJvVe2x/S5JX5b0sF4b\n9ro2yZ39VdUe2+dLOqBRh/oMSbcm+di6y3IBDgDUwqPLAKAYghsAiiG4AaAYghsAiiG4AaAYghsA\niiG4AaAYghsAivk/gA1y6umW0h8AAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x110c15210>"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In addition to the methods and attributes of `MAP`, `NormApprox`\n",
      "provides the following methods:\n",
      "\n",
      "`sample(iter)`\n",
      ":   Samples from the approximate posterior distribution are drawn and stored.\n",
      "\n",
      "`isample(iter)`\n",
      ":   An 'interactive' version of `sample()`: sampling can be paused, returning control to the user.\n",
      "\n",
      "`draw`\n",
      ":   Sets all variables to random values drawn from the approximate posterior.\n",
      "    \n",
      "\n",
      "### MCMC\n",
      "\n",
      "The `MCMC` class implements PyMC's core business: producing Markov chain Monte Carlo samples for\n",
      "a model's variables. Its primary job is to create and coordinate a collection of 'step\n",
      "methods', each of which is responsible for updating one or more\n",
      "variables. \n",
      "\n",
      "`MCMC` provides the following useful methods:\n",
      "\n",
      "`sample(iter, burn, thin, tune_interval, tune_throughout, save_interval, ...)`:   Runs the MCMC algorithm and produces the traces. The `iter` argument\n",
      "controls the total number of MCMC iterations. No tallying will be\n",
      "done during the first `burn` iterations; these samples will be\n",
      "forgotten. After this burn-in period, tallying will be done each\n",
      "`thin` iterations. Tuning will be done each `tune_interval`\n",
      "iterations. If `tune_throughout=False`, no more tuning will be done\n",
      "after the burnin period. The model state will be saved every\n",
      "`save_interval` iterations, if given.\n",
      "\n",
      "`isample(iter, burn, thin, tune_interval, tune_throughout, save_interval, ...)`: An interactive version of `sample`. The sampling loop may be paused\n",
      "at any time, returning control to the user.\n",
      "\n",
      "`use_step_method(method, *args, **kwargs)`: Creates an instance of step method class `method` to handle some\n",
      "stochastic variables. The extra arguments are passed to the `init`\n",
      "method of `method`. Assigning a step method to a variable manually\n",
      "will prevent the `MCMC` instance from automatically assigning one.\n",
      "However, you may handle a variable with multiple step methods.\n",
      "\n",
      "`stats()`: Generate summary statistics for all nodes in the model.\n",
      "\n",
      "`summary()`: A pretty-printed summary of posterior quantities\n",
      "\n",
      "The sampler's MCMC algorithms can be accessed via the `step_method_dict`\n",
      "attribute. `M.step_method_dict[x]` returns a list of the step methods\n",
      "`M` will use to handle the stochastic variable `x`.\n",
      "\n",
      "After sampling, the information tallied by `M` can be queried via\n",
      "`M.db.trace_names`. In addition to the values of variables, tuning\n",
      "information for adaptive step methods is generally tallied. These\n",
      "\u2018traces\u2019 can be plotted to verify that tuning has in fact terminated. After sampling ends you can retrieve the trace as\n",
      "`M.trace[\u2019var_name\u2019]`.\n",
      "\n",
      "We can instantiate a MCMC sampler for the bioassay example as follows:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M = pm.MCMC(gelman_bioassay)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 27
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Step methods\n",
      "\n",
      "Step method objects handle individual stochastic variables, or sometimes groups \n",
      "of them. They are responsible for making the variables they handle take single \n",
      "MCMC steps conditional on the rest of the model. Each subclass of \n",
      "``StepMethod`` implements a method called ``step()``, which is called by \n",
      "``MCMC``. Step methods with adaptive tuning parameters can optionally implement \n",
      "a method called ``tune()``, which causes them to assess performance (based on \n",
      "the acceptance rates of proposed values for the variable) so far and adjust.\n",
      "\n",
      "The major subclasses of ``StepMethod`` are ``Metropolis`` and\n",
      "``AdaptiveMetropolis``. PyMC provides several flavors of the \n",
      "basic Metropolis steps.\n",
      "\n",
      "### Metropolis\n",
      "\n",
      "``Metropolis`` and subclasses implement Metropolis-Hastings steps. To tell an \n",
      "``MCMC`` object :math:`M` to handle a variable :math:`x` with a Metropolis step \n",
      "method, you might do the following:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M.use_step_method(pm.Metropolis, M.alpha, proposal_sd=1., proposal_distribution='Normal')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 28
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`Metropolis` itself handles float-valued variables, and subclasses\n",
      "`DiscreteMetropolis` and `BinaryMetropolis` handle integer- and\n",
      "boolean-valued variables, respectively.\n",
      "\n",
      "`Metropolis`' `__init__` method takes the following arguments:\n",
      "\n",
      "`stochastic`\n",
      ":   The variable to handle.\n",
      "\n",
      "`proposal_sd`\n",
      ":   A float or array of floats. This sets the proposal standard deviation if the proposal distribution is normal.\n",
      "\n",
      "`scale`\n",
      ":   A float, defaulting to 1. If the `scale` argument is provided but not `proposal_sd`, `proposal_sd` is computed as follows:\n",
      "\n",
      "    if all(self.stochastic.value != 0.):\n",
      "        self.proposal_sd = ones(shape(self.stochastic.value)) * \\\n",
      "                            abs(self.stochastic.value) * scale\n",
      "    else:\n",
      "        self.proposal_sd = ones(shape(self.stochastic.value)) * scale\n",
      "\n",
      "`proposal_distribution`\n",
      ":   A string indicating which distribution should be used for proposals.\n",
      "Current options are `'Normal'` and `'Prior'`. If\n",
      "`proposal_distribution=None`, the proposal distribution is chosen\n",
      "automatically. It is set to `'Prior'` if the variable has no\n",
      "children and has a random method, and to `'Normal'` otherwise.\n",
      "\n",
      "Alhough the `proposal_sd` attribute is fixed at creation, Metropolis\n",
      "step methods adjust their initial proposal standard deviations using an\n",
      "attribute called `adaptive_scale_factor`. During tuning, the\n",
      "acceptance ratio of the step method is examined, and this scale factor\n",
      "is updated accordingly. If the proposal distribution is normal,\n",
      "proposals will have standard deviation\n",
      "`self.proposal_sd * self.adaptive_scale_factor`.\n",
      "\n",
      "By default, tuning will continue throughout the sampling loop, even\n",
      "after the burnin period is over. This can be changed via the\n",
      "`tune_throughout` argument to `MCMC.sample`. If an adaptive step\n",
      "method's `tally` flag is set (the default for `Metropolis`), a trace of\n",
      "its tuning parameters will be kept. If you allow tuning to continue\n",
      "throughout the sampling loop, it is important to verify that the\n",
      "'Diminishing Tuning' condition of [Roberts and Rosenthal (2007)](http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.jap/1183667414) is satisfied: the\n",
      "amount of tuning should decrease to zero, or tuning should become very\n",
      "infrequent.\n",
      "\n",
      "If a Metropolis step method handles an array-valued variable, it\n",
      "proposes all elements independently but simultaneously. That is, it\n",
      "decides whether to accept or reject all elements together but it does\n",
      "not attempt to take the posterior correlation between elements into\n",
      "account. The `AdaptiveMetropolis` class (see below), on the other hand,\n",
      "does make correlated proposals."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### AdaptiveMetropolis\n",
      "\n",
      "The `AdaptativeMetropolis` (AM) step method works like a regular\n",
      "Metropolis step method, with the exception that its variables are\n",
      "block-updated using a multivariate jump distribution whose covariance is\n",
      "tuned during sampling. Although the chain is non-Markovian, it has\n",
      "correct ergodic properties ([Haario et al., 2001](http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.bj/1080222083)).\n",
      "\n",
      "`AdaptiveMetropolis` works on vector-valued, continuous stochastics:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from pymc.examples import disaster_model_linear\n",
      "M = pm.MCMC(disaster_model_linear)\n",
      "M.use_step_method(pm.AdaptiveMetropolis, M.params_of_mean)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 29
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`AdaptativeMetropolis`'s init method takes the following arguments:\n",
      "\n",
      "`stochastics`\n",
      ":   The stochastic variables to handle. These will be updated jointly.\n",
      "\n",
      "`cov` (optional)\n",
      ":   An initial covariance matrix. Defaults to the identity matrix,\n",
      "adjusted according to the `scales` argument.\n",
      "\n",
      "`delay` (optional)\n",
      ":   The number of iterations to delay before computing the empirical\n",
      "covariance matrix.\n",
      "\n",
      "`scales` (optional):\n",
      ":   The initial covariance matrix will be diagonal, and its diagonal\n",
      "elements will be set to `scales` times the stochastics' values,\n",
      "squared.\n",
      "\n",
      "`interval` (optional):\n",
      ":   The number of iterations between updates of the covariance matrix.\n",
      "Defaults to 1000.\n",
      "\n",
      "`greedy` (optional):\n",
      ":   If `True`, only accepted jumps will be counted toward the delay\n",
      "before the covariance is first computed. Defaults to `True`.\n",
      "\n",
      "`shrink_if_necessary` (optional):\n",
      ":   Whether the proposal covariance should be shrunk if the acceptance\n",
      "rate becomes extremely small.\n",
      "\n",
      "In this algorithm, jumps are proposed from a multivariate normal\n",
      "distribution with covariance matrix $\\Sigma$. The algorithm first\n",
      "iterates until `delay` samples have been drawn (if `greedy` is true,\n",
      "until `delay` jumps have been accepted). At this point, $\\Sigma$ is\n",
      "given the value of the empirical covariance of the trace so far and\n",
      "sampling resumes. The covariance is then updated each `interval`\n",
      "iterations throughout the entire sampling run [^1]. It is this constant\n",
      "adaptation of the proposal distribution that makes the chain\n",
      "non-Markovian.\n",
      "\n",
      "### DiscreteMetropolis\n",
      "\n",
      "This class is just like `Metropolis`, but specialized to handle\n",
      "`Stochastic` instances with dtype `int`. The jump proposal distribution\n",
      "can either be `'Normal'`, `'Prior'` or `'Poisson'` (the default). In the\n",
      "normal case, the proposed value is drawn from a normal distribution\n",
      "centered at the current value and then rounded to the nearest integer.\n",
      "\n",
      "### BinaryMetropolis\n",
      "\n",
      "This class is specialized to handle `Stochastic` instances with dtype\n",
      "`bool`.\n",
      "\n",
      "For array-valued variables, `BinaryMetropolis` can be set to propose\n",
      "from the prior by passing in `dist=\"Prior\"`. Otherwise, the argument\n",
      "`p_jump` of the init method specifies how probable a change is. Like\n",
      "`Metropolis`' attribute `proposal_sd`, `p_jump` is tuned throughout the\n",
      "sampling loop via `adaptive_scale_factor`.\n",
      "\n",
      "### Automatic assignment of step methods\n",
      "\n",
      "Every step method subclass (including user-defined ones) that does not\n",
      "require any `__init__` arguments other than the stochastic variable to\n",
      "be handled adds itself to a list called `StepMethodRegistry` in the PyMC\n",
      "namespace. If a stochastic variable in an `MCMC` object has not been\n",
      "explicitly assigned a step method, each class in `StepMethodRegistry` is\n",
      "allowed to examine the variable.\n",
      "\n",
      "To do so, each step method implements a class method called\n",
      "`competence(stochastic)`, whose only argument is a single stochastic\n",
      "variable. These methods return values from 0 to 3; 0 meaning the step\n",
      "method cannot safely handle the variable and 3 meaning it will most\n",
      "likely perform well for variables like this. The `MCMC` object assigns\n",
      "the step method that returns the highest competence value to each of its\n",
      "stochastic variables."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Running MCMC Samplers\n",
      "\n",
      "We can carry out Markov chain Monte Carlo sampling by calling the `sample` method (or in the terminal, `isample`) with the appropriate arguments."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M = pm.MCMC(gelman_bioassay)\n",
      "M.sample(10000, burn=5000)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------57%--                ] 5798 of 10000 complete in 0.5 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------99%----------------- ] 9929 of 10000 complete in 1.0 sec"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [-----------------100%-----------------] 10000 of 10000 complete in 1.0 sec"
       ]
      }
     ],
     "prompt_number": 30
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pm.Matplot.plot(M.LD50)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Plotting LD50\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAFwCAYAAABpb3VdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xm8VVX9//HXh0lBEEUMZXIINDAzwwE19KppiANmZZmW\nkppmmn3j68wvL5WWVk4phqakmXOmlBNON4dSJNFUQAXlK4JijohKMnx+f+yzOfucu8+4z93nnnvf\nz8eDx91nnz2sve7l7M9Z67PXMndHRERERNpel3oXQERERKSzUOAlIiIikhIFXiIiIiIpUeAlIiIi\nkhIFXiIiIiIpUeAlIiIikpLEgZeZXW1mS83s2SLbXGJmL5nZM2a2fdJziogUY2ZDzOwhM3vezJ4z\nsx9m1vczs/vM7EUzm2FmG0T2OSPzOTXPzPaNrB9lZs9m3ru4HtcjIh1HLVq8pgFjC71pZuOAYe4+\nHPgecHkNzikiUsxK4H/cfRtgNPADMxsBnA7c5+5bAQ9kXmNmI4FvACMJPs+mmJlljnU5cHTmM2y4\nmRX8vBMRKSVx4OXujwDvFtnkIOCazLZPABuY2YCk5xURKcTd33D3pzPLy4G5wCAin0eZnwdnlscD\nN7j7SndfCMwHdjazTYE+7j4zs921kX1ERCqWRo7XIGBR5PVrwOAUzisigpltDmwPPAEMcPelmbeW\nAuGXwIEEn02h1wg+u/LXL86sFxGpSlrJ9Zb3WvMUiUibM7PewJ+Bk939g+h7HsyXps8iEUlVtxTO\nsRgYEnk9OLMuh5npA1Ckk3H3/C9lNWNm3QmCrj+6++2Z1UvNbBN3fyPTjfhmZn3c59RrmfWD89br\n80tEgOo+w9Jo8ZoOfAfAzEYD70Wa+nO4e4f4d/bZZ9e9DLoOXUt7/9eWMonxVwFz3P2ivM+jIzPL\nRwK3R9Z/08x6mNkWwHBgpru/ASwzs50zx/x2ZJ8cnflvoprz16reGvHaO8r5O/O1u1f/GZa4xcvM\nbgD2APqb2SLgbKB75j/UVHe/y8zGmdl84ENgQtJzioiUsBtwBPBvM5udWXcG8EvgZjM7GlgIHArg\n7nPM7GZgDrAKOMGzn6wnAH8AegJ3ufs9aV1EW5s8eTIAZ599dp1LUl+qB0lT4sDL3Q8rY5sTk55H\nRKRc7v4ohVv0v1Rgn3OBc2PW/wvYtnalaz8UaARUD5ImjVzfBpqamupdhJroKNcBuhbp2Or9N1HP\n83fma6/3+TvztSdhSfopa8nMvL2URUTanpnhbZhcnyZ9flUuHJ9W9SaNqtrPMLV4iYh0UpMnT16b\n39SZqR4kTWrxEpG6UItX56YWL2l0avESERERaecUeImIiIikRIGXiEgnpdymgOpB0qQcr07uRz+C\n738ftt663iWRzkY5Xp2bcryk0SnHS6py8cVw4431LoWIiEjnoMBLREREJCWJAy8zG2tm88zsJTM7\nLeb9/mZ2j5k9bWbPmdlRSc8pIiLJKbcpoHqQNCXK8TKzrsALBHOfLQaeBA5z97mRbZqBddz9DDPr\nn9l+gLuvyjuWciTqwAyam0FTlUnalOPVuSnHSxpdvXK8dgLmu/tCd18J3AiMz9vmdWD9zPL6wNv5\nQVc5WlpaOOWUU3LWHXXUUey0007svffe7LPPPjzyyCMALFy4kI033pg999yTPffck7feeguARx99\nlN12240xY8bw3HPPtTrHlVdeWWmxRERERMrWLeH+g4BFkdevATvnbXMl8KCZLQH6AIdWc6Lw21H+\nuj/84Q+MHDmSV199lbFjx/LQQw8BweSZt9xyS872kyZN4q677mLZsmUcf/zx3HnnnTnvX3HFFRx7\n7LE569w99twiIiIilUoaeJXTRnwm8LS7N5nZp4H7zGw7d/8g4bmDAmSaqYcOHcrXv/51ZsyYwe67\n785jjz3G7rvvzpgxYzjnnHP4+OOP6dq1K3379qVv37688847Ocf5y1/+wgsvvMBee+3Fsccey733\n3kvv3r158cUXue666zj88MNZuXIlPXr04M9//jN9+vRh2rRpXHHFFay77rpMmjSJ3XbbjWOOOYbX\nX3+d3r17c91119GnT59aXKaIdABHHXUUBx54IF/96lfXrlu4cCEjRoxgxIgRrFixgj59+nDCCSdw\n5JFHAkFr//jx49lyyy0B+OpXv8qkSZMAuOeee/jRj37E6tWrOeaYYzjttFZptrHnBFiyZAn77LMP\nhx56KGfH5Bq8//77XH/99Xz/+9+v2fWH/vd/T1+7/N3v/qDq4+y1124cccS3EpcnzO+KqweRWksa\neC0GhkReDyFo9YraFTgHwN0XmNkrwNbArPyDNTc3r11uamqiqamposIMHDiQ119/nU033ZQFCxbQ\ns2dPjj32WG677TZ22WUX1l9//bXbduvWjVWrVtGtW1AFX/nKVzj33HN58MEHAZgxYwajRo3i0ksv\nBWD69On07NmTiy66iJtuuonx48dz5ZVX8sgjj9CtWzfcncsuu4y9996bCRMmcNNNN3HFFVcwceLE\niq5BpKNqaWmhpaWl3sWoKzOLbUEfNmwYTz31FACvvPIKhxxyCO7OUUcdBcAee+zB9OnTc/ZZvXo1\nJ554Ivfffz+DBg1ixx135KCDDmLEiBGtzhln4MCBPP/88wXL+u677zJlypQ2CbyuvPIPa5enTRtZ\n5VEe5ZNPWmoSeCngkjQlDbxmAcPNbHNgCfAN4LC8beYRJN8/ZmYDCIKul+MOFg28yhX9UFm8eDFb\nb701PXr0WLvukEMO4fHHH2e//fZj2bJla9dHg65CdtxxRwCWL1/Occcdx+LFi3nnnXf42te+xiuv\nvMKoUaPWHsPMmDNnDrNmzeLaa69l5cqV7L777hVfj0hHlf9lqrM+RVYqmXyLLbbgggsuYOLEiWsD\nr7h9Zs6cybBhw9h8880B+OY3v8kdd9zRKvACePjhh7ngggt44403OP/88/nqV7/KwoULOfDAA3n2\n2Wd5/vnn+e53v8snn3yCu3PrrbcyadIkFixYwPbbb8++++7Leeedl/ja41Xb4tWdmO/vIu1eosDL\n3VeZ2YnAvUBX4Cp3n2tmx2XenwqcC0wzs2cIkvlPdfd3Ch608jIAsGjRIm677TYeeughli9fTu/e\nvYHgA2ebbbahZ8+erFq1ivfff59ly5bRr1+/VsfK/2YYvp4xYwZbbrklf/rTn7jgggv44IMP+PSn\nP81TTz21NoBzd0aMGMGuu+7KEUccAQTBnYhIpbbffnvmzZu39vU//vEPtttuOwYNGsSvf/1rRo4c\nyeLFixkyJNvhMHjwYJ544olWx3J33njjDR577DHmzp3LQQcd1Krb8Xe/+x0nn3wy3/rWt1i1ahWr\nVq3ivPPO4/nnn2f27Nltd6EinVDSFi/c/W7g7rx1UyPLbwEHJj0PwC233MLTTz8NsDa4mTBhAn36\n9KFr165MnTqVjTfemHvuuYdJkybRq1cvttxyS8455xwAfv7znzNu3Di6dOnClClTWh1/zz335OCD\nD2bChAlANvAaPXo05557LrNnz2bAgAFsttlmbLTRRhxzzDHsttturLfeepx11ll873vf43vf+x7T\npk0DYOLEiYwbN64Wl96m9DS3SPsSbeEaNWoUixYtolevXtx9990cfPDBvPjii2Ufy8w4+OCDARgx\nYgRLly5d+96bb77J5MmT2XXXXTnnnHN47bXXOOSQQxg2bFinGuZBOV6SpsSBV1r22GMPFi5cmLMu\nTD7NN3bsWMaOHdtq/ZgxY3jssccKniPalD5+fHZUjIEDBzJrVusm7QkTJqwN0kLXXHNNweOLiJTz\nlPTs2bMZOTLIfYo+oLPffvtxwgkn8M477zB48GAWLco+VL5o0SIGDx4ce7xo+kU0oPrUpz61NtgY\nPXo0f/vb3xg3bhxTp05liy22qOzCGpgCLkmTpgwSEUlRqZakhQsXcsopp3DSSScBsHTp0rX7zJw5\nE3enX79+7LDDDrz00kssXLiQTz75hJtuuomDDjqoqjK9/PLLbLHFFpx00kmMHz+eZ599lvXXX58P\nPqjJw+ciEtEwLV7SdjRMmUh6jjvuOH70ox8BwTA4119/PQsWLOALX/jC2uEkTj75ZL7zne8AcOut\nt3L55ZfTrVs3evXqxY2ZWe27devGpZdeype//GVWr17N0UcfHZtYD7mtbHHLN998M9dddx3du3dn\n00035ayzzmKDDTZgt912Y9ttt2XcuHFtmFwv0rkkmjKoljTlRn1oyiCpF00ZVH/1ym3q23cTli0L\nc82qrbcrOPzwWVx33RWJy6McL6lGtZ9havESEemkFGgEVA+SJuV4iZ5qFBERSYlavEREpCxvvvlm\nzpOU1Vq9emUNSiPSmBR4iYh0UpXmNt166638+Mdns+66QxOdd/XqzYGajaOdmHK8JE0KvERPNYp0\nUtUEGmZf5/33Ww9AXbn288GjgEvSpBwvERERkZQkDrzMbKyZzTOzl8zstALbNJnZbDN7zsxakp5T\nakvJ9SIiIulI1NVoZl2BS4EvAYuBJ81survPjWyzAXAZ8GV3f83M+ic5p4iI1IZymwKqB0lT0hyv\nnYD57r4QwMxuBMYDcyPbfAv4s7u/BmsnzRYRkTpToBFQPUiaknY1DgKizxa/llkXNRzoZ2YPmdks\nM/t2wnOKiIiINKSkLV7lZAd1B74A7A30Av5pZo+7+0sJzy01oqcaRURE0pE08FoMDIm8HkLQ6hW1\nCHjL3T8GPjazh4HtgFaBV3Nz89rlpqYmmpqaEhZPRNqLlpYWWlpa6l0MiVBuU0D1IGlKGnjNAoab\n2ebAEuAbwGF529wBXJpJxF8H2Bm4IO5g0cBL0qOnGiUN+V+mwpud1I8CjYDqQdKUKPBy91VmdiJw\nL9AVuMrd55rZcZn3p7r7PDO7B/g3sAa40t3nJC24iIiISKNJPHK9u98N3J23bmre618Dv056LhER\nEZFGppHrRerklltgzZp6l0I6s8mTJ6vLF9WDpEtzNYqeaqyTQw+FhQths83qXRLprJTbFFA9SJrU\n4iVKrhcREUmJAi8RERGRlCjwEhHppJTbFFA9SJqU4yUi0kkptymgepA0qcVLpI70YIOISOeiwEt0\n868jPdggItK5KPAS3fxFOinlNgVUD5Im5XiJ1JFaG6WelNsUUD1ImhK3eJnZWDObZ2YvmdlpRbbb\n0cxWmdkhSc8pIiIi0ogSBV5m1hW4FBgLjAQOM7MRBbY7D7gH0Hd8ERER6ZSStnjtBMx394XuvhK4\nERgfs91JwK3AfxKeT0REakS5TQHVg6QpaY7XIGBR5PVrwM7RDcxsEEEwthewI6BUbhGRdkC5TQHV\ng6QpaYtXOUHURcDp7u4E3YzqapSSVq6EKVPqXQoREZHaStritRgYEnk9hKDVK2oUcKMFj2/1B/Yz\ns5XuPj3/YM3NzWuXm5qaaGpqSlg8aVTPPQc/+AGccEK9SyK10tLSQktLS72LISJSV0kDr1nAcDPb\nHFgCfAM4LLqBu28ZLpvZNOCvcUEX5AZeItKx5H+ZUk5N/YW/g87e1aZ6kDQlCrzcfZWZnQjcC3QF\nrnL3uWZ2XOb9qTUoo4iItAEFGgHVg6Qp8QCq7n43cHfeutiAy90nJD2fiDQ+DRwrIp2VpgwSERER\nSYkCLxGRTkrjVwVUD5ImzdVYxP77wymngB6uTF9nmbhbXW5ST8ptCqgeJE1q8SrirrvgL3+pdymk\nI+voAeZHH8Hrr9e7FCIi7YcCLxFpM8cdBwMH1rsUIiLthwKvEkq1SPzf/8HLL5d/vGefhTlzkpVJ\nOo6O3tX4xhv1LoEUo9ymgOpB0qQcrxJKBV6f/zwsWwarV5d3vM99Dnr0gP/+N3nZRNq7jh5YNjrl\nNgVUD5KmTtfiNWsWLF9eu+N99BGsWVN6u2efzQZxXbvW7vwiIiLSODpd4LXjjvDzn5e/fakWr3K/\n0X/uc/Dvf1e2j0ij09+6iEiuThd4AaxcWbtjVXJj+eST4GeXTlnrtfPWW3DOOfUuhUjjU25TQPUg\naUqc42VmY4GLCOZq/L27n5f3/uHAqYABHwDfd/d/Jz1ve1FJ4BW2nqkVIJm//AUmTYKzzkp+rOXL\nYZ11oHv35Meq1pNPwoABMHRo/crQVvS33r4ptymgepA0JWp7MbOuwKXAWGAkcJiZjcjb7GVgd3f/\nHPAz4Iok5+wI1OJV2Mcfl96mljfzPn3ghz+s3fEq9Y9/wE47wde+Vr8ytCUFXiIiuZKGADsB8919\nobuvBG4Exkc3cPd/uvv7mZdPAIMTnjOxSgatrFWOV9J9OoteveDRR4tvU+v6e/HF1uvOOqv2A3+u\nWNF63VNPBT9XrartuWppzZryn9oVEZHikgZeg4BFkdevZdYVcjRwV8JzpqotRhZXi1dxS5bUuwRw\n7rkwfXptj9mzJ/z1r7nrwidia5l3WGtHHAHDh1e3r75ktG/KbQqoHiRNSXO8yg5LzGxP4LvAbgnP\nGWvNmmCYhnICpVreDBq1xWvePFh//XqXImjdevNNOOSQ8vepdf0VOl5bBN133gljx2ZzysoZiqSW\nVq+ufDiTf/4TFi5sk+JInSm3KaB6kDQlDbwWA0Mir4cQtHrlMLPPAVcCY9393UIHa25uXrvc1NRE\nUwWzU5fTVfPww8XfHzoUnngCNt00u66WN9/wWHEtXitXwi23wLe+VbvzFTNiBGy3XTrnKuaII4LR\n/6P17J59/de/woEH5u6TVuDaFoHX1KnBoLvHH597jjSu6Y474OCD05kfcuVK+PDD7HUtXw6zZrXQ\n0tLS9icXEWnHkgZes4DhZrY5sAT4BnBYdAMzGwrcBhzh7vOLHSwaeFUqbDlwL3wT22OP4sdYtCjI\n91l3Xdhww/LOG57r6KPhqqvK2ycu8PrHP+Dww9MLvKA+3VsXXBAEyaeeGryOCwKi6w46qH4TSZ91\nFnz/+7U51jvvZJejuV5pBl6VTG0VVU39n3oqXHQR7L9/8Pq//839MhX06qhrR0Q6n0TZRu6+CjgR\nuBeYA9zk7nPN7DgzOy6z2U+ADYHLzWy2mc1MVOICwsCrnK6bYjeS+fOhX7/yzxveMK++uvx92kuO\nV1iONAObU0+F007Lvi507h12KHyMtg5S5s4Nfr77bnCuRx5JfsyLL45fn2bg1a3Kr1nV/H0sWJD7\nuj10r0trym0KqB4kTYnH8XL3u4G789ZNjSwfAxxTzrH++99gTKVqhE9dhble1froo8q2r+aGEhd4\ntWXws2pV/E231lMXPfccfOYzxW/w+fVVTY5TW+d4jRyZ+3rhQhgzJtk5Cv1+08zxqvb3XYu/zXp8\n2TCzq4H9gTfdfdvMumaCz6P/ZDY7M/MZhpmdQZCHuhr4obvPyKwfBfwBWBe4y91PTvEy2pRymwKq\nB0lTO2l7CVx4YfX7RgOvJPKDhrYIiNL+9t+9e5BblC+8EdeqPNtuW7q7Nf8GXE39tvfWk0cfDeYE\njSp0nWm2eFV7jriyX3ZZ8THX8vepUyvvNIIxBqMcuMDdt8/8C4OukQSpEiMz+0wxW1tjlwNHu/tw\ngtSK/GOKiJStXQVe5QyeWUglXY3FtOUI5sWS69va8ce3Hq8qDLwmTw5aGwsxg2eeKe88H35Y/P1y\nAq+2GDstiVWrKgsQx4yBvffOXZf/8ED+cnsOvOKceGK2C9YdXngh9/3wutJ+ajO3DP4IEPcwT1xN\njAducPeV7r4QmA/sbGabAn3cPUyRuBY4uC3KKyKdQ7sKvJK0LlXS4lXsBhTX4vXuu4UH0wznX6xE\n3PnTyLPaeuvc19GupzfeKL7vc8/Vpgz5194ILV7f/S6ccAI8+2z5+xS7ruhTm+HfbbXX9PHHxYPm\nqGoD/lLdpPfdF3Qxx1m2LPf17rvD3XfHb5uik8zsGTO7ysw2yKwbSO4T2eGYhPnrF1N8rMKGotym\ngOpB0tSuAq/QyJFw6aWV7VNJi9f8Is9WxuUn7bcfDBwYv/2AAaXPl6+9dJVVkvNTboBUarvw5j9l\nCvz859W1eIUqTXq/7jo46aRg+dprYc6c8vf93e/gc5+r7HzliPt7XbKk/ID+s5+FceOC1sxSw6XU\nsqsxuj58OvbBB1u/F33aGILfWalytrHLgS2AzwOvA7+pa2nq7Oyzz1Z+E6oHSVfi5PpaCm8Mc+fC\njBlBd0a5Kmnxmj49mGj55Zdh4sT4MoQ3CndYurTwsaptRVi+PBi2opInzQolyVcrGniVuil/8EFt\nzhme56yz4L334FOfar1NuV2Nu+9eeFuz4PccHQPsoovgX/+C3/4WjjwyGNMq3Pbxx4PcrP/93+Ln\nXrEi+L1VKr+c0b+vsAyhQYPgzDPhnHOy655/HrbZpvVxX34Z3n47mOvx2Wfj62PlyuDvptDv+Be/\ngMMOg803L1722bNh++2z68P/az16BD8vugj22it+3/bC3d8Ml83s90A4l0D+mISDCVq6FpM7zdng\nzLpWkoxDKCLtX0tLbcYibFctXoXyYIr51KeClpNKk+t/9rPiN9lyz19J4BW9yfbpAyef3Pq9Qp58\nMjf/7LjjgtHnk6gk8KpVXlqayfX5eUfFnqg85xw45ZTSx/z978s7d6FAK1/YTZhftsWLs++tWBG0\nbBUad23VquKtlzvuGPwfKVRvZ54J11zTev2bb0Lv3tmy77JL7vv5gVe0fOE+jz+e+7reMjlboa8A\nYQfydOCbZtbDzLYAhgMz3f0NYJmZ7ZxJtv82cHvcsZubm9f+U9Al0vE0NTXl/D+vVrsKvKrxn//A\nY4+17tIopdCwFdGE4XJU84h+eAN86aXy91mc9x37iivgz3+u/NyQHeX/oYfK36fcwCust9694fbI\n7WnFChg/PnucsA7STL7Ov4boxM/lBnO1CEDffz+7HA5fkv8U5MqVQcvaww+3Drhm5o2Et2pV8XI9\n80zwd11pV+PChbkPSxQKJsP/S9FW0fxtn3iisnPXgpndAPwD2NrMFpnZd4HzzOzfZvYMsAfwPwDu\nPge4mWA8wruBE9zXXsUJwO+Bl4D57n5PypfSZpTbFFA9SJraVVdjVCXfkN0rb/Eq1HUWDruQ3xVU\nSJIbcXTfSs4TtoQUOrd7kCN34onxyew331x+GcNy/e//wrHHlr/fhx8G+VRhd97rrwddf+HgtGG5\n3n67/GNCEICWG0CUGtIgHE0+ery33ip+zGJdvV/5SuEhUf7zn+yyWev8KICWFggbSsL1Z5wRzGoA\n2X123hl69Qq6aiG3xeull4IWqM02yz1/sRkdoseOCuur0Jea8HXY9dq1a5AmMGJE6/9fY8emOytD\nUD4/LGZ1waGO3f1c4NyY9f8Ctq1h0doN5TUFVA+SpoZv8YLgxlDpcBLPP59d/vDD1q1H0ZtMGNRN\nmVJ9GaPCG2AlQVv0pjl9eut1y5Zlb9bPPQc//CE89VTr46xZE99l9cAD8ecNrz3/6bRCoknhr0We\nBQsDlvwWrzjFgtDBg4OJpit1442tn0qMO8/GGxc/TrHf2e23Z1sRwycXw662aBdl9Nqjc4yGgRRk\n6zH6dxr10UfB07aQO/H1VlvB6NGtt1+zpvIWr/Baw7+BQoFXeO41a4IHY1asCObfzHf99fHHERHp\nTNpt4FXqw3nlymzrxDPPZFuB8gOv0aNzu3biXH99kJxc6PyLFgU/H3us9b6FbmaPPdb6ybSvf730\nvoWuO+6GH92/b99sjlI4F2BcELpqVfx5J0yIP2+5gWzYFVposvLw5lyLrrpiQWC0BS2sy48+CpLH\nw26z8JrCenj++dp1NYbH+fDDoHttl11a18nq1blP+eWXF7KtcdEWtjvuyN3ml7/Mni/a5R0XWOcH\nXu5BfuTTTxe+lnIDr/zX7rDnnoWPKyLSmSW+DZrZWDObZ2YvmdlpBba5JPP+M2a2fdw2+UoFXj/5\nSbZ14s034corg+VooLBiRXDze++93PXl3PzDuRej5agkn+uLXwyGL4gKg5Pw/HfdBfdkskXCoKBL\nl/iWquhNM1z++OPcoDKcBLlYsFTu5N+h/GMtWBA/1EG0tSZO/ij50evp2TN321K/+0KD3M6eDf37\nt14fzeWC1q2b+flzxZQbeEGQtA6tJ6eOBmLRIDJ63WFAFv2be+qp3G1efTVbnuh2cX+n0a7GRx+F\nX/0q6P6cNi37/tFH504KHm5fKPA65JDc19HtSo0LJ+2DcpsCqgdJU6LAy8y6ApcSTLExEjjMzEbk\nbTMOGJaZbuN7BOPoxKqkCyK/K2P58uBnNFAIb3zu8Le/Zdfvs0/p4996a+syxeUjFWspyb/hx+0T\ndvFFhw6IG2csTGCeMSO7bvJk2GCD1tsWq8dKZweI1udhh8GwYXD++YW3z0+VCAOL/GE6oiodVDX/\nad7hw4MnRAuNnp5/vPBvp5our1LBd/Rawjyn6O82f5uoaHnClq7o+T78MHebaG5XXODV1AT//Gew\n/PTT2fPus092ovLweMuWBV82fve7IMB7883WLV6FAvrwGGGA+dRTsN568dtK+6LxqwKqB0lT0hav\nnQie8lno7iuBGwmm3og6CLgGwN2fADYws5LDjsZ1Y4RdaNA6Jynu5vDqq9l10WTfuBalUueHIH+m\nEoVu7NGgMby5RUcfX7Kk9T5h4HXAAYXP9+9/55630A2+VLfam28GwQzkBo833hj8fDduEpaI6HWH\n5c0fqb1YeebOzX2d38IWdsNBENTMnx88/Zc/wXmhaWvCLtkw/yquDIWUasmJtoiF5c7/e+vVq3QX\nblj2aFfjb3+bW7crV8Z34YbLf/97du7M99+Pr49ogn5o1Kig9Ss8drkTx7+ZGSFrzJjCQ1+AcrxE\npHNLGngNAhZFXofTbJTaZjAl5H84T52a2yUVfUoM4gOv8GaS39UYvfkXGtIh3P7vf8+u69OnVKnL\nEw2y4rqu1qwJxm2KTq0SBgbFhgcKyxwGmaWmeilk3rwgmHngAVh//dbvh8ct1HoWDZDDbrxiU+Tk\nBz2//GXu60JDfwBMmpRdzg/qzj03fn0Sp55a/P0w2Ifstec/zTd8ePyYYXEPdCxalLtNNEBauTJb\n14W6GqNLseGaAAAgAElEQVSTluc/HQnwxz+2Pm74upzBepcty/19h155pfS+IiKdUdLAq9zvrvnt\nCRV/5w3nClyxIviXH1TEBV7ht+5Ro3LXR2/kd95ZvLVjwYLsctx2xfYtJ8+lUOD1/PNw//2tr3Oj\njUq3zuR3ueWLthjF+de/gp8LF8a/v2ZN8Dvo1St7w46WM9q6uGxZ7k09Lggq1joSVWwGgbhjf/BB\n0AVd67HCosF4vrPOyi6HgUtcV+oll7Tet5yWoNtuyy5/5zvZfaLBT6Hu0PD3EFffcQ9FFAq8oi1g\nG28M3/xmeceT9ke5TQHVg6Qp6The+dNsDCF3Qtm4bQpOufH3vzcTDgZ7771NvPtuExtumHvjGj06\neIoxX3jTffVV2HLLYDnaRRW9KUdvTKtXxx8vblaASp/IO/vs4CGAYuKGV8h/ECA6Qn05wzD88IfF\nzzlqVPH3f/zj4Oef/hT//r/+la3bDz8MnrY78sjs+9HA6623gvG8wnqIa/Eqd5Ln008v/J5Z9mYf\nzef77W8LP7EZdccd5ZUBgrG69tgjKHf37sHvKO7hgkKBVKGnMqOTZxdy+OHZ5eiDEtFg6r334kfY\nL9byl18m98J/7+GDLBD8HcS1bsUH0y1Ay9qx8qT+lNcUUD1ImpK2eM0ChpvZ5mbWA/gGwdQbUdOB\n7wCY2WjgPXePbbvYffdmhg5tBpqBprWDbUbFBUkQPKkFwWTBELT63H9/9v3oTSd6Q1m1qvwJuaPB\nwttvw667VvZEXJywLEcdlV2X/xRZ/g0z6STb5QaQhUa2f/TRbHC4fHm27kP5LUzz5xfvaixXfvdy\nvvDY0fkZV6+ufYvXHXcEOVPrrgs//WmwLm7IklLJ6HHry8k/jDvOJptkl999N36w22KtUPmD6v73\nv/H/1/7zn/Ja5vKfVA00Ac28/34zwf9xEZHOJ1Hg5e6rgBOBewmm2rjJ3eea2XFmdlxmm7uAl81s\nPjCVYPqNWK++GiT1ViMMgMKb3ZNP5r5fKPCqJP8nGizMnx88MVZqGIVSLrssGNU92noRXkPYauCe\nHUMsLmgKA6TX8toaC3XN1SK5Oay3uMAmv6uzS5fWgVdcgn6xct12W+mBU+N+l2vWBK1TtXbxxcHP\nhx8Ofsb9Xgo95FAsEAyfzq1UOfmHlXxJaGnJzjoQdeml5c0RWmjIDxGRzi7xlEHufjfB3GbRdVPz\nXp9YzrGSth5BcFN75JHsKNnR9VFduwY36kpyUX7yE/j0p4Nk6bib/EknBbk7lbTmvP12EFT07t26\nrNHAK3z8P+7Ye+2VXY7mDo0fH1/OWgRe0Wlkognc0Hq4iWjgFdZ3/pOKL76YDWLi3HRT6TLFXWtb\n9SCExx0wIHg6csiQ1tuEw4/kd5EWawmrZG7lsHU33LeUSubmLOSOOwq3OkfV8oEGaTthXlNn72pT\nPUia2tVcjQ8+2HpdpUHCnXcGwUd+q0v0RrDhhkEe2IMPtp5wuJQf/KBw4HXppcGQEyedlF23alXp\nLrLu3XOvMwxArrkm+Ln//rnbFwvsTj45u1zoBn/LLfHrC41qHycMkuOmGoqOxg5BkFvq93jaacEY\nYYWUM79kqWEu2sJOOwWj099TwbTJhQKgSv/WoxNPpzVEQzlBFyjwahQKNAKqB0lTu50yKFTpB/hd\nd+U+iRh3nDVrsl1D1cz7B9mR5vNNnJj7+uKLYeDA4sfq0SP3xhl2N4VPqkWHKFi+HC64oLKy5iuU\n07bVVuUNIQCwfWb+gXK6ht9+G4YOLb7N7bfDr39d3rnjfPIJ/OIX1e9fqTD/MGyVHDu2/H3DwXnz\nVRo8RYPRWbMq27etxc1uICIiHTDwKrRPtEvxmWeqnzPwvfeC0eP32y/+/ZUr4fOfz74uZ0iJbt1y\nb7rRgT3z3X5762T2YioJ0tpq7KVq85YqUWgy6bbSq1fws5YtTUmCp3LyrtIU9+VHREQaIPCqpkUq\nLqjKz+VKMllzoUFXQ9HumHICxw8+qP6hglLyW+DqIRyDrSMJA67f/a52x7zwwtodS6QcGr8qoHqQ\nNLWrHK84X/1q6W2OOSZ33KK4PKVaBl6VDE9QTpfLCQWf85T2KsxxKzTIrEgjUG5TQPUgaWr3LV7l\n+MpXcl/HPan4//5f7uu33qr+fKUGKI267LLqzyMiIiIdS4cIvHr0KH/b3XYLflb6NKOIiIhIUh0i\n8Co2ibKIiMRTblNA9SBpavc5XuWopMUrSW6XFHfAAbnzJErnNnp08Sd0pf6U2xRQPUiaOkQYUkmL\nlwKvtrPeevUugbQXV12VHQBYRESyEoUhZtbPzO4zsxfNbIaZbRCzzRAze8jMnjez58ysgtT08lQy\nllLSCaalsM98pt4lkPbiG99IbzR9EZFGkrT953TgPnffCngg8zrfSuB/3H0bYDTwAzMbkfC8OSoZ\nJTsaeJ15Zi1LIQceWO8SVK4tJtAWtX42CuU2BVQPkqakOV4HAeGt6xqghbzgy93fAN7ILC83s7nA\nQGBuwnOz997BRNA77lj+PtHAa6utkp1/vfUKTx1U7PxxLQF77AF//3uy8vztb0GeVb00YmviX/8K\n669f3rY33gjf/GbblqcjaMS/g85KuU0B1YOkKWmL1wB3X5pZXgoMKLaxmW0ObA88UWy70KBBxd/f\nYYdgQupK8rbCm8Lvfgff/nb5+8UJp42pRKGybrZZ9eX48peDn/mTaactOqckwOWXJzteGkFO797l\nb7v55m1WjA7la18LfioAExFprWTIksnhejbm30HR7dzdgYJZHWbWG7gVONndy5q975hjir9/113l\nHCW/HMHP3XdPnmj/xS9Wvk+hcyYJUk49tbogMF/cBNkbb1z9/pXsG/Wb3wQ/v/Wt6vYvJn8mhEqC\ng3Jylm64IX79RReVf5620FZTUoV22im7fP31wU/leImItFYy9HD3fdx925h/04GlZrYJgJltCrwZ\ndwwz6w78GbjO3W8vdK4hQ5qB8F8L665bvGyV3tjdszfacm+4w4cXfq/crspoF2LXrq3fHzo0WeC0\n886Vd3nGiZsKqWfP4vtEJwTPv9EWa8U74ojC74UBV1vcuMePr37fUuV5/fXCrXT1fpo26TyQT5Ro\no54+PbscF8AHWsj+/25OViCpCeU2BVQPkqakOV7TgSOB8zI/WwVVZmbAVcAcdy/6vX/rrZtZtCj7\nutSciDvskF3u0qW8ORTDG2AtukEGDy5vuwGRDti4wCuptuzSKfXgQt++hd/bYQc4/XT45S9bv/e5\nzxXeL/wdFQt0CuXKFdKtWzCVVL9+rd/r2xfef7/4/nvtVXq8uE02Kfxe0t+RezDh+QUXlLf9F74A\nTz2V7JxRn/1sLY7SlPkX0o2u3pTbFFA9SJqSfg//JbCPmb0I7JV5jZkNNLM7M9vsBhwB7GlmszP/\nxsYdLP/mVOpmFb3xrl4Ny5aVLnClLV7F5OeITZkSv10YSPTvn9tCFAqvIy65/pBDgp9NTa3fS9J6\nE4oGgnGBTKnPo2g9xu0fV24onb8H8YFM+CDFO++U3j8qbIUZNw6WLAmWw2t/773S+z/wQBDMAPzs\nZ7nvldOaVYsWr0q+kOd/CSkUpN5yS+vcvHzu2RbZ6NOKm25a+vgiIpIr0e3A3d9x9y+5+1buvq+7\nv5dZv8Td988sP+ruXdz98+6+febfPbGFyStNeFPfdttC58993adP6TKHwUq5gVexG0r+e9//fu7r\nT386+Ble1zPPwP33Fz7e7ru3XhcGap/6VO76Cy+ESy4JlpO0or32WnY57lqPP774/vmBV/5NPO6Y\nH3+cnTMTcsv/xz9m99l559b7hgFAsZa2OGHrpFk2YMgPoEoJrzUacAB0715631oEXoW78ALR/yfu\nuXVf6O842v2e79vfhldeqayMcZ58svW6IUOSH1dEpBG1q3Hcw5tT2CIRPkU2Zkxwg8h/Aq2cFpqJ\nE3Nff+97wfABG2SGes2/af7hD+WXN7rvXnu1fj8sX3hdvXoRm7c2Nrb9LxAGJX/5S+76/fcPcsM+\n+KD4yP2zZ8evDwORaKvSN75R+Dih6DUffnjpBPi439G668LAgcETqZBtSYIgCBgwAOaWGGyk0hbL\nbbapXatMNa2l4d/Ad75TetuDD45fXyrAjg4lkn+txf5GCl3PmjWtn+TccMPscpKnPCdNqn5fqR3l\nNgVUD5KmdhV4hS0Jm24a3DhGjQpehzeGXXbJ3T5u+ITm5uDn1VcHP/feO/f9Ll2CfJ4wMT8/ybtY\n68X22+e+jna7rF7devv8wKuQK67ILuc/dRfmsUVvpE88AcOGBculhkOoZDqlYcNat6zlu+++7PLx\nx8NGG2Vf77tv6+1HFBgqt3v3oMXOPT5XLn8U/FNOCX4WCp5KXWe9nyoM/wa6dMkO2vruu/HbRgOs\n7baLXx8Kf/9Dh+Z+yYjW04knxgf84Xb5gdeBBwZPDOfn5r3yCvz2t8HywoXBNr/4RfxxIfi9xrXi\nQvb/ttTX2WefrfwmVA+SrnYVeP3oR/DGG9nX4Q2h0P+HPfcsfKzwacRSrQT5N/L8IKnUk5Onnw5b\nbx2fQF7omMVcdVXu6zAY/dKXsut22ql4q8tuuwXXddJJsOWW5Z87mhc0ahT85Cett8kf6T1af3FP\nQJbTKlJO/Zx/fuvzRUUHQc3vCnRv2zG4KmkBM8tuv0GrCbaybr+99bHz68k9aPF89FG4++7cIDju\nQZMXXoiftDq//NOnw377tQ6IN988+2Vjs82C8he7hvXWSz4osIhIR9OuAq9tt819AjAUBj/VzMk4\ncGDws9A8gvm5QvlTnQwdml2OS/b+xS9g3jwYPbr1e6NGBbkslTxJmV+eMHCMlqNcl1xSuCUoLthZ\nvTpbxw8/3DqZO/9G3BZdd//8Z/Ftyzln3N9QNaKtTfkuuQSmTQuWTzwxCHzK0dyce73nnRe/Xdh1\nfckl2dbb0IUXwqxZ2de77QYjRwbLYQ5ctJ7C8221VfyTmeHfwgcfwNtvFy//mDEwdWrxbfJFr/fe\neyvbV0Sko2lXgVe+/EDlggvKH7Q0bH1Zf/3gJlQoZ+jcc+HFF7Ov994brrkmftubbiodAEXHO7rh\nBpg/P9sqUM5wF/nCOrjkkvIeHoDyRsGPCwJXrcresOMSuaNDfdRStCylJtouVIc77xyU7403gjy+\n0P/9X3llyA/WRo2Kbx2CoJv0pJPgqKOC1717F8/Tixo6NPd6Tz219TbR98eMgQkTct8fNqx0V12h\neor7vYfr1lsvfriNqHXXza3fSvPd4rqjpX6U2xRQPUiako7j1aY22SQ3R2jbbeGRR2DlSlheYuz7\nsKWrVDdWr165g6R26ZIds+iii+A//8m+16dP0I316quFk8qjI3iHwUvYypCkhWiddYIA7oMPim/3\n1lutW+0mT27dXRt3w9xvP7j22mC51JhVhSRJPIfqn9AcMCDbIhc9Rrkthb/8ZTAAajhx+jrrxOdF\nPfZY7tygRx3VOi+vlELBz223wVe+ErxO2poY1+IVNWkS/Pzn1Q0qHFXtsCzlPAkqbU95TQHVg6Sp\nXbd4rbcezJnTen337rlPV+Vzz3YLVnNjCG/cJ59c+Abxpz9VftxKbqbR4RZC5bSYbbRR64AhLlcr\nv14eeyzII6v2hh8mvxcqUzHRssT9vqLDSnTpUnpS6/AYpbq1og9LHHUUnHFG8e0Bdt01929i2rTW\ng4u+9FLrgU6LBULf/jYceWTwNONZZwV/d9UKzxP9W4mr35/9LJggvNATlOWq5v/X008XHiJGRKSj\na9ctXrVQ6Y3BLOiq22KL2pel0pHW87tIp00LWrRqIb9ewgFNyynjO+8EXVJ9+mRbBMPk99Dzz2eX\n99wTbr218PFK1fWhh2aX3UuPMh8q1q1V6jrDp0arMWwY/M//BOefMKH1OFb5dR+2MkLQCgWFr3GX\nXXKH3yikVIsXZIefSGvw07D7v1junIhIR9fhA69KB9qEoEvv5ZeD5fybUpKbVH6L1W67wX//G79t\n3M1yn32qP3e+gQPjr7Gc69twwyCwGjkSFiyI3yZM9o5ad1248srW63/60yDX6TOfye2qjVPJtFBJ\nROfOrLY7bZtt4ltMzzijugnWAf7xj/K2u/nm6o5frWJ/N+FTph99lE5ZpHxhXlNn72pTPUiaOnTg\nVUmQ9NZbwZQ+pW6yF14YPMVYStxUOflBw49+BF/7WtlFTGybbYKWrRkzglaZRx5pvU25dRYXWBUS\nHvOss+Inx+7WLQh2o0OJFFJOblAtpoOq1fyXEye2Hn9uzz2LD4WSRFjX0W7UUteS5FrL2XfgQE0p\n1F4p0AioHiRNVbcNmFk/M7vPzF40sxlmVnBEHzPrmpmj8a/Vnq+t5Sekh/Lzk3bdFb773eLH6tUr\n98mvUH7g1ZaTW8c5++xs3lOhc8fdIGv1JFq1o5VHc7byxzlr7w45BH7968oDjz594Mc/rk0Zor/r\nUuPSJTm2iIiUlqRT5nTgPnffCngg87qQk4E5QLv93rvuusGglflP851wQuXH+vBDOOyw1uvD4SDe\nfz94am6bbSo/dnvzxS/CN7+Zfb3ppvCb37TNufbdt3DuVXQYhKTBwG9+k/tAQj2Ciy5daleP0fKr\n9UlEpL6SdDUeBITjmF8DtBATfJnZYGAccA5Qo+/wbSOcQDsqyQTUUdGb3frrw4oVxbcPb5ZtdZOs\npMWrmAEDgvHKQl27tm6paevA5ZVXcge3TXq+Aw+MHyy3UZU7/ls11OLV2JTbFFA9SJqSBF4D3H1p\nZnkpUGi88AuBU4ASgwC0X/vtV/7I5LXSlje0m27KHaIhqi0CvaTHjBvMFYLR3R98sPV0QEnrrtTw\nF+3ZRhsFT52G5s0rPV2SgqfOS4FGQPUgaSoaeJnZfUDcd/+zoi/c3c2s1e3VzA4A3nT32WbWlKSg\n9fSnP+XezBpddHiGctSza+qxxwqP0n7AAUHglS9JINHo3XCzZuVO2L711qX36d8fWlqqO5+CNhGR\nyhQNvNy94AAGZrbUzDZx9zfMbFPgzZjNdgUOMrNxwLrA+mZ2rbt/J+6Yzc3Na5ebmppoins0sA42\n3LD4gK0dSXsLPHbdtbp9oqPLJ1WL4CKteq22izR/8vNyVVI3LS0ttFQb4YmIdBBJuhqnA0cC52V+\n3p6/gbufCZwJYGZ7AP9bKOiC3MBL6qO9BV7VGDYMZs6sdyk6lzAAK9Yln/9lSnPj1Z9ymwKqB0lT\nksDrl8DNZnY0sBA4FMDMBgJXuvv+Mft0gNt6OurVhdMec7wKOeCA1qPCt1d77NF6PK+OIPw7HTAg\nCLrKnSxc2gcFGgHVg6Sp6sDL3d8BvhSzfgnQKuhy978Df6/2fJKORmrxGj4crr++3qUoz2c+U/6o\n843ITEGXiEg52vUk2Z1ZWw8nUUhbnK/RE7AbvfxtSXUjIlKZDj1lUKM6/vjajRbfHjRSK5pURoFX\nY1NuU0D1IGlS4NUOXX55drktB7+MExckKXAS6ZgUaARUD5ImdTW2c0OHwrJl6Z1PQVZratUpTHUj\nIlIZtXg1gLRbvfINHw7331/9/t/6Vu6gntJxKPASEamMWrwkR1yL18UXJ2t1++pXgwnIpeP52tfg\nzjvrXQqp1uTJkzWeGqoHSZdavCTHmjWt13XvHvwTybfOOjBuXL1LIdVSblNA9SBpUuAlOdZbD/r2\nrXcp2hd1p4mISK0o8JIcc+aodUtERKStKPDqRMp5YnHQoLYvRyPZeGPo3bvepRBpGxq/KqB6kDRV\nHXiZWT/gJmAzMnM1uvt7MdttAPwe2IZgrsbvuvvj1Z5XqvfRR/UuQeN5/nnookdQpINSoBFQPUia\nktxSTgfuc/etgAcyr+NcDNzl7iOAzwFzE5yzIbS0tNS7CK0MGwaf/3xl+7TH66hWtdey8caw0Ua1\nLUtSHen3IiLS2SQJvA4CrsksXwMcnL+BmfUFxrj71QDuvsrd309wzobQHm+ML72kwKuj6EjXIiLS\n2SQJvAa4+9LM8lJgQMw2WwD/MbNpZvaUmV1pZr0SnFNERGpE41cFVA+SpqI5XmZ2H7BJzFtnRV+4\nu5tZXOp2N+ALwInu/qSZXUTQJfmTKssrIp2UmV0L3ODud9e7LB2FcpsCqgdJk3mVk/OZ2Tygyd3f\nMLNNgYfc/TN522wC/NPdt8i8/iJwursfEHM8zRIo0sm4e9mjpJnZOsA3gP2BfwC/d/cP26pslTAz\nr/aztJFMmTKFiROfY8WKKTU4Wvirr7beruDww2dx3XVX1KAsIpUzs4o+w0JJhpOYDhwJnJf52WpS\nmExQtsjMtnL3F4EvAc/HHayawotIp7IRsCXwPkF6w9UEgZiISMNIEnj9ErjZzI4mM5wEgJkNBK50\n9/0z250E/MnMegALgAkJzikinddEYIq7LwAws0V1Lk/D0/hVAdWDpKnqwMvd3yFowcpfv4SgKyB8\n/QywY7XnERHJaIkEXfu7u6bnTkiBRkD1IGmq+9CQZjbWzOaZ2Utmdlq9yxPHzK42s6Vm9mxkXT8z\nu8/MXjSzGZmBYsP3zshczzwz2zeyfpSZPZt57+I6XMcQM3vIzJ43s+fM7IcNfC3rmtkTZva0mc0x\ns1806rVEytHVzGab2V8zrxvyWsxsoZn9O3MtM2t4LXtETjMmvSsSEamdugZeZtYVuBQYC4wEDjOz\nEfUsUwHTCMoYFTuArJmNJMg7GZnZZ4rZ2mmWLweOdvfhwHAzyz9mW1sJ/I+7bwOMBn6Qqe+GuxZ3\nXwHs6e6fJxiYd8/w4Y1Gu5aIk4E5ZLONG/VanODBm+3dfafMusTXAmxnZnub2V7ED18jItLu1bvF\naydgvrsvdPeVwI3A+DqXqRV3fwR4N291oQFkxxM88r7S3RcC84GdLXjys4+7z8xsdy0xg862JXd/\nw92fziwvJ5hFYBANeC0A7h5OgtQD6ErwO2rIazGzwcA4gum1wsCjIa8lI/9hmVpcyyJgK+AzwI/a\nsOydhsavCqgeJE31niR7EMGHaeg1YOc6laVShQaQHQhE56J8jeA6V2aWQ4sz6+vCzDYHtgeeoEGv\nxcy6AE8BnwYud/fnzawhrwW4EDgFWD+yrlGvxYH7zWw1MNXdr6Q217Il8CKwDkHr4E/b7Ao6CeU2\nBVQPkqZ6B14dYuCbIgPItktm1hv4M3Cyu3+Q7dlprGtx9zXA5y2YmupeM9sz7/2GuBYzOwB4091n\nm1lT3DaNci0Zu7n762a2MXCfBWP+rZXgWrYE/kYQlImINKR6B16LgSGR10PI/Zbbni01s00iA8i+\nmVmff02DCa5pcWY5un5xKiWNMLPuBEHXH909HHutIa8l5O7vm9mdwCga81p2BQ4ys3HAusD6ZvZH\nGvNacPfXMz//Y2Z/IUgpqMW1vObuz7X5BYiItKF653jNIkgA3tyCcb6+QTAwayMIB5CF3AFkpwPf\nNLMeZrYFQVLwTHd/A1hmZjtnkoe/Tcygs20pc96rgDnuflHkrUa8lv7hk3Fm1hPYB5hNA16Lu5/p\n7kMyMzx8E3jQ3b/diNdiZr3MrE9meT1gX+BZanMtXc3sr2Z2i5ndkuJldVjKbQqoHiRV7l7Xf8B+\nwAsESbVn1Ls8Bcp4A7AE+IQgJ20C0A+4nyDnZAawQWT7MzPXMw/4cmT9KIKb0HzgkjpcxxeBNcDT\nBEHKbIInyRrxWrYlyO96Gvg3cEpmfcNdS9517QFMb9RrAbbI/E6eBp4L/0/X4lqA3sCOmfcGl1GW\nqwnyyZ6NrOsH3FegHGcAL2XKsW9MOV4CLi5wLu8MLrvsMl933e87eA3+kflX7f5T/fDDj613lUgn\nlvl/X/HnZNVzNYqIpMnMrgQ+cfcfmNkUdz+hxPZjgOXAte6+bWbd+cBb7n6+BeMGbujup2eGtbie\nYLDnQQRB4nB398xYZCe6+0wzu4sgoL0n71zeGT5LNVejSJZVOVdjvbsaRUTKtZygBQvg41IbewcZ\nBkZEOhYFXiLSKN4CdjWz3xB0mVej2LAW0Qd7wmEt8tfXdRiYWlNuU0D1IGmq91ONIiJlcfdzzOwz\nQBd3n1OD49V0iI7m5ua1y01NTTQ1NdXq0G1G41cFVA9SjpaWFlpaWhIfR4GXiDQEM7shs9gzk1tR\nTZdfmw3REQ28RKTjyf9CVW0rqboaRaQhuPth7n4Y8BXg4SoP03BDdIhIx6IWLxFpCGa2DcEjcN2B\nbcrY/gaC4Tn6m9ki4CfAL4GbzexoYCFwKIC7zzGzmwkmKV8FnBB5TPEE4A9AT+Cu/CcaG1n4jb2z\nd7WpHiRNCrxEpFF8LfPzvwTjehWVaR2L86UC258LnBuz/l8E48Z1OAo0AqoHSZMCLxFpFLMiy4PN\nbLC731m30oiIVEGBl4g0imOAxwi6G7+Icq1EpAEp8BKRRjHP3X8NYGYbu/s1pXaQ4pTbFFA9SJoU\neIlIwzCzqwhavJaW2lZKU6ARUD1ImhR4iUijOItgHK33CBLsRUQajsbxEpFGcRFwtrsvA35b78KI\niFRDgZeINIo1wP9llt+rZ0E6Cs1RGFA9SJrU1SgijeK/wEgzOwnYsN6F6QiU2xRQPUiaFHiJSLuX\nma7nVqA/YMCU+pZIRKQ6CrxEpN1zdzezPd39/HqXRUQkCQVeItLumdl4YLyZfRl4B8Ddv17fUjU+\njV8VUD1ImtpN4GVmXnorEelI3N3K3HSsu+9mZpe7+/fbtFCdiAKNgOpB0tSunmp09w7x7+yzz657\nGXQdupb2/q9CQ81s/8zPcWY2rg0+gkRE2ly7afESESniFoLE+puBjetcFhGRqinwEpF2z93/UO8y\ndETKbQqoHiRNCrzaQFNTU72LUBMd5TpA1yISR4FGQPUgaWpXOV4dRUe5MXaU6wBdi4iItA8KvERE\nRERSkjjwMrOrzWypmT1bZJtLzOwlM3vGzLZPek4REUlOcxQGVA+SplrkeE0DfgtcG/dm5rHvYe4+\n3BHeCd0AABwWSURBVMx2Bi4HRtfgvCIikoBymwKqB0lT4hYvd38EeLfIJgcB12S2fQLYwMwGJD2v\niIiISKNJI8drELAo8vo1YHAK5xURERFpV9JKrs+fFkTTA4mI1JlymwKqB0lTGuN4LQaGRF4Pzqxr\npbm5ee1yU1OTHpsX6UBaWlpoaWmpdzEkQrlNAdWDpCmNwGs6cCJwo5mNBt5z96VxG0YDLxHpWPK/\nTKmFQUQ6o8SBl5ndAOwB9DezRcDZQHcAd5/q7ndlJrWdD3wITEh6ThEREZFGlDjwcvfDytjmxKTn\nERGR2tIchQHVg6SpYeZqbGlp4c477+RXv/rV2nVHHXUUc+bMoU+fPnTp0oWf/OQnjBkzhoULF7Lj\njjvy2c9+FoBbbrmF/v378+ijj3LaaafRpUsXLr/88rXvi4h0Rgo0AqoHSVPDBF5m+Q9GBuv+8Ic/\nMHLkSF599VXGjh3LQw89BAT5JLfcckvO9pMmTeKuu+5i2bJlHH/88dx55501L6e7ry1rdFlERESk\nYQKvQtyDkSmGDh3K17/+dWbMmMHuu+/OY489xu67786YMWM455xz+Pjjj+natSt9+/alb9++vPPO\nO62ONXHiRP71r3/x8ccfc8UVV7Dddtsxc+ZMJk6cSLdu3TjggAOYOHEiEydOZObMmfTo0YOrr76a\nzTbbjJEjRzJ69Gj69u3Le++9x3rrrceLL77I9ddfT//+/dOuFhEREWmHGj7wiho4cCCvv/46m266\nKQsWLKBnz54ce+yx3Hbbbeyyyy6sv/76a7ft1q0bq1atolu3bBX8/Oc/p2fPnsyePZtf/epXXHfd\ndfz4xz/mpptuYtCgQbg7s2bNYsmSJTzyyCM8+uij/PSnP+Wqq65i8eLFXHjhhfTt25cJEyYwatQo\nLr300npUg4hIWZTbFFA9SJoaPvCKduUtXryYrbfemh49eqxdd8ghh/D444+z3377sWzZsrXr84Mu\ngPPPP58HHngAgO7duwPwySefMGjQoLXnWrBgATvuuCMAO+ywA2eeeSYAw4YNo2/fvmuPFW4jItJe\nKdAIqB4kTWmNXN9mwq7GRYsWcdttt7HvvvuyfPnyte8//PDDDB8+nJ49e7Jq1Sref/99Fi1aRL9+\n/XKO8/bbb3P//ffz8MMPc+GFF7JmzRoA1llnHZYsWbL2XMOGDePJJ58E4Mknn2SrrbYCoEuX3KpU\nbpeIiIjka6gWr1tuuYWnn34agCOOOAKACRMm0KdPH7p27crUqVPZeOONueeee5g0aRK9evViyy23\n5JxzzgGCrsRx48bRpUsXpkyZknPsfv360a9fP/bcc09Gjx69NnC64IILOPTQQ+nevfvaHK9NN92U\nMWPG0L17d6ZNmxZbVgVeIiIiks/CFqN6MzNvL2URkbZnZrh7h/iG0qifX5XmNk2ZMoWJE59jxYop\npTcuKfzVV1tvV3D44bO47rorEpdEOV5SjWo/wxqqxUtERGpHgUZA9SBpavgcLxEREZFGocBLRERE\nJCUKvEREOqnJkyevzW/qzFQPkibleImIdFLKbQqoHiRNavESERERSYkCLxEREZGUKPASEemklNsU\nUD1ImpTjJSLSSSm3KaB6kDQlbvEys7FmNs/MXjKz02Le729m95jZ02b2nJkdlfScIiIiIo0oUeBl\nZl2BS4GxwEjgMDMbkbfZicBsd/880AT8xszU0iYiIiKdTtIWr52A+e6+0N1XAjcC4/O2eR1YP7O8\nPvC2u69KeF4REUlIuU0B1YOkKWnL0yBgUeT1a8DOedtcCTxoZkuAPsChCc8pIiI1oNymgOpB0pS0\nxaucaeXPBJ5294HA54HLzKxPwvOKiIiINJykLV6LgSGR10MIWr2idgXOAXD3BWb2CrA1MCv/YM3N\nzWuXm5qaaGpqSlg8EWkvWlpaaGlpqXcxpAO577572WOPgxIdo0cPuOGGq+nfv3+NSiVSnLmX02hV\nYOcgSf4FYG9gCTATOMzd50a2uQB4390nm9kA4F/A59z9nbxjeZKyiEhjMTPc3epdjlpo1M+vMK+p\n3K62KVOmMHHic6xYMaUGZw9/9dXW20Lg34lLsc46EzjjjB8C6nKUylT7GZaoxcvdV5nZicC9QFfg\nKnefa2bHZd6fCpwLTDOzZwi6Nk/ND7pERCR9jR1obJ75l0zXrj05+uijGTx4cOJjiZQj8bAO7n43\ncHfeuqmR5beAA5OeR0RERKTRacogERERkZQo8BIR6aQ0flXgqquuUj1IajSCvIhIJ9XYOV61oxwv\nSZNavERERERSosBLREREJCUKvEREOinleAWU4yVpUo6XiEgnpRyvgHK8JE1q8RIRERFJiQIvERER\nkZQo8BIR6aSU4xVQjpekSTleIiKdlHK8AsrxkjSpxUtEREQkJQq8RERERFKiwEtEpJNSjldAOV6S\nJuV4iYh0UsrxCijHS9KkFi8RERGRlCQOvMxsrJnNM7OXzOy0Ats0mdlsM3vOzFqSnlNERESkESXq\najSzrsClwJeAxcCTZjbd3edGttkAuAz4sru/Zmb9k5xTRERqI8xr6uxdjldddRWgepB0JM3x2gmY\n7+4LAczsRmA8MDeyzbeAP7v7awDu/lbCc4qISA0o0Agox0vSlLSrcRCwKPL6tcy6qOFAPzN7yMxm\nmdm3E55TREREpCElbfHyMrbpDnwB2BvoBfzTzB5395cSnltERESkoSQNvBYDQyKvhxC0ekUtAt5y\n94+Bj83sYWA7oFXg1dzcvHa5qamJpqamhMUTkfaipaWFlpaWehcDADNbCCwDVgMr3X0nM+sH3ARs\nBiwEDnX39zLbnwF8N7P9D919Rj3KXWvK8Qoox0vSZO7lNFoV2NmsG/ACQWvWEmAmcFhecv1nCBLw\nvwysAzwBfMPd5+Qdy5OURUQai5nh7lanc78CjHL3dyLrzif4knh+5gntDd39dDMbCVwP7EiQSnE/\nsJW7r4ns2yk+v6ZMmcLEic+xYsWUGhwt/NXXt9569RrMCy88rhwvqVi1n2GJcrzcfRVwInAvMAe4\nyd3nmtlxZnZcZpt5wD3AvwmCrivzg66oaKuXiEgbyv/APAi4JrN8DXBwZnk8cIO7r8w8SDSf4MEi\nEZGKJR7Hy93vdvet3X2Yu/8is26qu0+NbPNrd9/G3bd190uKHU/TNohIChy4P/PAz7GZdQPcfWlm\neSkwILM8kNwUiriHiEREyqIpg0SkM9rN3V83s42B+8xsXvRNd3czK9YH1iH6FZXjFVCOl6RJgZeI\ndDru/nrm53/M7C8EXYdLzWwTd3/DzDYF3sxsnv8Q0eDMuhyN+HCQAo2AxvGSctTqAaFEyfW1FCan\nZpLV6l0cEWlj9UquN7NeQFd3/8DM1gNmAJMJZuB4293PM7PTgQ3ykut3IptcPyyaTa/k+moouV4a\nW7WfYWrxEpHOZgDwFzOD4DPwT+4+w8xmATeb2dFkhpMAcPc5ZnYzwQNEq4ATOkWUJSJtQoGXiHQq\n7v4K8PmY9e8QtHrF7XMucG4bFy11yvEKKMdL0qTAS0Skk1KgEVCOl6Qp8XASIiIiIlIeBV4iIiIi\nKVHgJSLSSU2ePFmDVhPkeKkeJC3K8RIR6aSU4xVQjpekSS1eIiIiIilR4CUiIiKSEgVeIiKdlHK8\nAsrxkjQpx0tEpJNSjldAOV6SJrV4iYiIiKREgZeIiIhIShR4iYh0UsrxCijHS9Jk7p7sAGZjgYuA\nrsDv3f28AtvtCPwTONTdb4t5390dMyNpmUSk/cv8X7d6l6MWws+v9uxXv/oVH3/8caJjzJw5kwce\nGMqKFVNqUKLwV1/feuvVazAvvPC4crykYtV+hiVKrjezrsClwJeAxcCTZjbd3efGbHcecA/Z/20i\nIpKSyZPP48MPjwJ6JTjKF4Ada1MgkU4q6VONOwHz3X0hgJndCIwH5uZtdxJwK/ofKyJSR6cD/etd\nCJFOLWmO1yBgUeT1a5l1a5nZIIJg7PLMqvbdHi8i0kk0N0+muVm5TcrxkjQlbfEqJ4i6CDjd3d3M\njCJdjc3NzWt/NjU10dTUlLB4ItJetLS00NLSUu9iSERzs8bxAo3jJelKlFxvZqOBZncfm3l9BrAm\nmmBvZi+TDbb6Ax8Bx7r79LxjKblepBNRcn26evfuz4cfzqP9dDUquV4aW12S64FZwHAz2xxYAv+/\nvfuPkaM+7zj+/mAbiDFgESMgYGTSUAsqkUIoNpA0jiDJgVpoQxGYENpAqVuVhgq1BlOpXtRUpJWa\nogQ14YdB+UFDgUBrWhpiQk8tAWMo2CEBE5yUxkBwCKRAJUiP8PSPmQ3j89z+mJmb2bn7vKTR7czO\n7vP97s5+97nZZ7/LWcDK7A4R8c7uZUk3AndOTrrMzMzMZoNSiVdEvCHpIuBukukk1kXEE5JWpddf\nU0EbzcxsGnTru2b7R47r1q0D/BNKVo/S83hVxR81ms0u/qixXv6oMZ8/arSiio5hIz1zfbfY3szM\nzGwmGOnEy1/vNTMzs5lkpBMvMzObPp7HK+F5vKxOI13j5Xovs5nLNV71co1XPtd4WVEzssbLzMzM\nbCZx4mVmZmZWEydeZmazlGu8Eq7xsjq5xsvMGuEar3q5xiufa7ysKNd4mZmZmY04J15mZmZmNXHi\nZWY2S7nGK+EaL6uTa7zMrBGu8aqXa7zyucbLinKNl5mZmdmIc+JlZmZmVhMnXmZms5RrvBKu8bI6\ntaLGq9Pp0Ol0mm2gmVXKNV71co1XPtd4WVGN1XhJGpO0VdJTki7Nuf6jkrZI+pakb0o6atgY/k/E\nzMzMZoJSiZekOcDVwBhwJLBS0hGTdvs+8KsRcRTwF8C1ZWKamZmZtVXZM17HAdsi4umImABuBk7P\n7hARD0TEy+nqg4DP55qZjQDXeCVc42V1mlvy9gcD2zPrzwDLeux/AXBXyZhmZlaBTmdt000YCRdc\ncIFrvKw2ZROvgasiJX0AOB84sWRMMzMzs1Yqm3g9CyzOrC8mOeu1k7Sg/jpgLCJ+MtWddb+52Ol0\nWLFiRcmmmdkoGR8fZ3x8vOlmmJk1qtR0EpLmAk8CJwHPAZuAlRHxRGafQ4F7gXMjYmOP+5pyOgn/\ndJDZzOPpJOqVN51Et76rmY8cR2c6idWrLwRg7Vp/9GqDKzqGlZ7HS9IpwFXAHGBdRFwpaRVARFwj\n6XrgN4EfpDeZiIjjcu7HiZfZLOLEq16exyuf5/Gyohqbxysi/jUilkbEuyLiynTbNRFxTXr5dyPi\n7RFxdLrsknQNwxOpmpmZWVu17ieD/JVfMzMza6vWJV5mZlYNz+OV8DxeVqdW/FZj3jYzazfXeNXL\nNV75XONlRTVW49Uk13uZmZlZm7Q68fKpYTMzM2uTVideZmZWnGu8Eq7xsjq1usbL9V5m7eUar3q5\nxiufa7ysqFlZ49XlWi8zMzNrgxmRePkUsZmZmbXBjEi8srI/tD15m5mZvcU1XgnXeFmdZkSN16Db\nOp2OkzCzEeEar3q5xiufa7ysqMZ+JLsqdSRe2W1Owsya5cSrXk688s2ffwgbNtzCQQcdVOp+FixY\nwP77719Rq6wNnHgV3GZmzXDiVS8nXvn22ON45s17vtR9TEy8yhlnnMZNN91QUausDYqOYXOnozFm\nZjb6uvVdnc7ahlvSnJ/+9AHWrCn7ONzAxMR91TXKZjQnXmZms9RsTriy/DhYnWbctxqH5TovMzMz\nq0vpxEvSmKStkp6SdOkU+3wmvX6LpKPLxqySv0LcPlUly4NOOeLk3MzMqlIq8ZI0B7gaGAOOBFZK\nOmLSPqcC74qIw4HfAz5XJqa1x3QlLGWT5W67svfTvZyXjDk5t5nK83gl/DhYrSKi8AIcD3wts34Z\ncNmkfT4PnJVZ3wockHNfkX4tKLp6bet3fZFta9eu3env5MujaJTb2n1s+7Vr0HZ39+v3nPWLUcUx\nYuWlj2upMWhUluwxMqr22uvtAS8ExIgspEvT7ahiWRdnnvnxpp9iq1nRMazsYPNbwHWZ9XOBz07a\n507ghMz6PcB7cu4r25GdLteVePXaNmpvuHlJyCgkCP0SnLxEadA3rUGfs163HeZ+yhwPdT/+o3Z8\nDsKJV72ceE3n4sRrNmoq8TpjwMTrxMz6PcAxOfeV7chOl0ch8ZqOpKbM/Qz62OSdESsTN++2UyVR\nRZ7TXm0t8/xM5/GQp0jCXub5acMb/2ROvOrlxGs6Fydes1HRMazUBKqSlgOdiBhL19cAb0bEX2X2\n+TwwHhE3p+tbgfdHxI5J9xWwNrNlRbqY2cwwni5dVxCeQLU2eROoNjuP12hMoApVPA43cOaZ93HL\nLZ5AdTYpPAl0kWytu5DMA/Y9YAmwO7AZOGLSPqcCd6WXlwMbp7iv6tPRGuWdWcr2qXu53xmjUaoz\nG7WP02aSYc9uZY+lYW+Tdxz2u77I8dxrvx7bSo1Bo7K0YfzyGa/pXHzGazYqOoZVMeCcAjwJbAPW\npNtWAasy+1ydXr+FnI8ZoyUD17CK1DGZ5SmS5PZL4vslZr3a0Ks9TrxGkxOv6VyceM1GjSVeVS1t\nGLjK8NkhGzXDJl6DGvQMrhOvejnxms7FiddsVHQMG7kfyTazenQ6nZ2W7LY6+Eey6+Uar6m5xsuK\nKDqGOfEys0Y48RrMHXfcwZVX/h1l7/6RR8Z5880fkk28mjU6iVd5Trxmo6JjmH8k28xshG3fvp3N\nm/diYuIPS97TamCfKppkZiU48TIzG3G77XYo8MGmm2E9fOMbd7Ns2YdL3cfuu8Ott36BAw88sKJW\n2Shy4mVm1oekMeAqYA5wfWTmKmyzZmu8Rkf5x+FDvPTSwWzaVK4de+55Hq+//nq5O7GR58RrGoyP\nj7NixYqmm1HaTOkHuC9WnKQ5JFPinAw8CzwkaX1EPNFsy7LGKTLhdHUJV7H4oxK73OPQjX9IqTYA\nzJnztuGjNzgeND0WNR2/qN2absBMND4+3nQTKjFT+gHui5VyHLAtIp6OiAngZuD0QW64dOlRLFq0\npNRy+eUdfvazfpHGy/WwtCbjNxm7+fhNjgdNj0VNxy/KZ7zMzHo7GNieWX8GWDbIDbdv385rr40D\n+5Zswt4lb29t8clPfop99x38eLn//vt49dXXdtp27LHHsHLlWVU3zSrixMvMrLfC8x1IsPfeq5F2\nr7I9u3j99SfZc8//HPp2l1xyLACf/vTDtcd/5ZXk7z77/HrtsScr8zhUEb/rlVf+m3Xrrhn6dhs3\n3r/LtnPOObuKJvV1xRVX1BKnivgXXvj7XHvt56axNYMZqXm8mm6DmdWrDfN4SVoOdCJiLF1fA7yZ\nLbD3+GU2O7V6AlUzs1EkaS7J79GeBDwHbAJWjlZxvZm1hT9qNDPrISLekHQRcDfJdBLrnHSZWVE+\n42VmZmZWk8ank5A0JmmrpKckXdp0e4YhabGkf5P0HUnflvSJdPt+kjZI+q6kr0ta2HRbByVpjqRH\nJd2ZrreuL5IWSrpN0hOSHpe0rI39gKSeKD2+HpP095L2aEtfJN0gaYekxzLbpmx72ten0vHgQ820\nejCDPgc5x+LyOuOn++70mq4j9lRjY8m4fd8rJH0mvX6LpKPLxhwmvqSPpnG/Jembko6qK3Zmv1+R\n9Iakj1QVe9D4klakx9m3JY3XGV/SIklfk7Q5jf87FcbeZRzL2We44y4iGltITttvA5YA84DNwBFN\ntmnI9h8I/HJ6eQFJHcgRwF8Dq9PtlwKfarqtQ/TpEuAmYH263rq+AF8Azk8vzyX5Ln8b+7EE+D6w\nR7r+D8Bvt6UvwPuAo4HHMtty2w4cmb7+56X93gbs1nQfevRtoOcg71isM356/U6v6TpiTzU2lojZ\n970COBW4K728DNhY4fM9SPzju88vMFZV/EFiZ/a7F/hn4Iya+74Q+A5wSLq+qOb4HeDKbmzgRWBu\nRfF3GcfKHndNn/EqPDHhKIiI5yNic3r5f4EnSOb8OY1kwCX9+xvNtHA4kg4hOYiuB7rf1GhVXyTt\nC7wvIm6ApD4nIl6mZf1IvQJMAPOVFHjPJynubkVfIuI/gJ9M2jxV208HvhIRExHxNMlAe1wd7Syo\n73PQ41isJX7ahrzX9LTHnmJsfEeJmIO8V/y8XRHxILBQ0gElYg4VPyIeyDy/D1LFVPYDxk79EXAb\n8EJFcYeJfw7w1Yh4BiAiflxz/B/y1i/A7wO8GBFvVBF8inEsa+jjrunEK29iwoMbakspkpaQZMUP\nAgdExI70qh1AVS/+6fa3wJ8Cb2a2ta0vhwEvSLpR0iOSrpO0F+3rBxHxEvA3wA9IEq7/iYgNtLAv\nGVO1/R0kr/+uUR8LBnkO8o7F+TXGh/zXdF2xgV3GxqIGea/I26eq5GfY96oLgLvqii3pYJJkpDtJ\nVZXF24P0/XBgv/Tj5Yclfazm+NcBvyTpOWALcHGF8fsZ+rhr+luNM6KyX9IC4KvAxRHxqvTWP5YR\nEWrBHD+Sfg34UUQ8KmlF3j4t6ctc4Bjgooh4SNJVwGXZHVrSDyT9AvDHJKfYXwZulXRudp+29CXP\nAG1vtF+SNpB8ZDbZn2VXevRjqmPxz+uIP8hrerpiZ+5nAclZmIvTM19FDXosTD6rV9UxNPD9SPoA\ncD5wYo2xrwIuS58PUd3ZzUHjzyM51k8iOTP/gKSNEfFUTfEvBzZHxIp03Nwg6d0R8WoF8Qcx1HHX\ndOL1LLA4s76Ynf/rHXmS5pEkXV+KiH9MN++QdGBEPC/pIOBHzbVwYCcAp0k6FdgT2EfSl2hfX54B\nnomIh9L124A1wPMt6wfAscD9EfEigKTbSepI2tiXrqmOp8ljwSHptsZExAenui4ttu33HOQdi5fl\n7Ddd8fNe01+MiPNqiJ0dG7+cGRuLGuS9YjqPoYHeq9KC+uuAsYjo9fFU1bHfA9yc/tO/CDhF0kRE\nrK8p/nbgxxHxGvCapH8H3g1UkXgNEv8E4C8BIuJ7kv4LWAqU+0mGYu3re9w1/VHjw8DhkpYo+U2N\ns4AqDpRapP9ZrAMej4irMletJymCJv1bdtCZdhFxeUQsjojDgLOBeyPiY7SsLxHxPLBd0i+mm04m\nKfq8kxb1I7UVWC7pbemxdjLwOO3sS9dUx9N64GxJu0s6jOSji00NtG9QfV8XPY7FuuLnvab7Jl1V\nxO4xNhY1yHvFeuC8NP5yko/md1CNvvElHQrcDpwbEdsqijtQ7Ih4Z0Qclj7XtwF/UFHSNVB84J+A\n9yr5Bu18kiLzx2uMv5Xk9UVaX7WU5ItJdRj+uOtXfT/dC3AKyTdetgFrmm7PkG1/L0ntxGbg0XQZ\nA/YD7gG+C3wdWNh0W4fs1/t561uNresLyX9aD5F81n87ybcaW9ePtC+rSd6sHyMp4JzXlr4AXyGp\nTfs/kv+IP96r7SQfF2wjGUQ/3HT7+/Qttx8ktWr/0utYrDN+Zv+fv6briD3V2Fgy7i7vFcAqYFVm\nn6vT67cAx1T8nPeMT/IFhhcz/d1UV+xJ+94IfKTOvqfrf5IZqz5R82O/iOQf0i1p/HMqjD15HDu/\n7HHnCVTNzMzMatL0R41mZmZms4YTLzMzM7OaOPEyMzMzq4kTLzMzM7OaOPEyMzMzq4kTLzMzM7Oa\nOPEyMzMzq4kTLzMzM7Oa/D8+OH5vNe70LwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x110c09c90>"
       ]
      }
     ],
     "prompt_number": 31
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.core.display import HTML\n",
      "def css_styling():\n",
      "    styles = open(\"styles/custom.css\", \"r\").read()\n",
      "    return HTML(styles)\n",
      "css_styling()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<style>\n",
        "    @font-face {\n",
        "        font-family: \"Computer Modern\";\n",
        "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
        "    }\n",
        "    div.cell{\n",
        "        width: 90%;\n",
        "/*        margin-left:auto;*/\n",
        "/*        margin-right:auto;*/\n",
        "    }\n",
        "    ul {\n",
        "        line-height: 145%;\n",
        "        font-size: 90%;\n",
        "    }\n",
        "    li {\n",
        "        margin-bottom: 1em;\n",
        "    }\n",
        "    h1 {\n",
        "        font-family: Helvetica, serif;\n",
        "    }\n",
        "    h4{\n",
        "        margin-top: 12px;\n",
        "        margin-bottom: 3px;\n",
        "       }\n",
        "    div.text_cell_render{\n",
        "        font-family: Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
        "        line-height: 145%;\n",
        "        font-size: 130%;\n",
        "        width: 90%;\n",
        "        margin-left:auto;\n",
        "        margin-right:auto;\n",
        "    }\n",
        "    .CodeMirror{\n",
        "            font-family: \"Source Code Pro\", source-code-pro,Consolas, monospace;\n",
        "    }\n",
        "/*    .prompt{\n",
        "        display: None;\n",
        "    }*/\n",
        "    .text_cell_render h5 {\n",
        "        font-weight: 300;\n",
        "        font-size: 16pt;\n",
        "        color: #4057A1;\n",
        "        font-style: italic;\n",
        "        margin-bottom: 0.5em;\n",
        "        margin-top: 0.5em;\n",
        "        display: block;\n",
        "    }\n",
        "\n",
        "    .warning{\n",
        "        color: rgb( 240, 20, 20 )\n",
        "        }\n",
        "</style>\n",
        "<script>\n",
        "    MathJax.Hub.Config({\n",
        "                        TeX: {\n",
        "                           extensions: [\"AMSmath.js\"]\n",
        "                           },\n",
        "                tex2jax: {\n",
        "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
        "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
        "                },\n",
        "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
        "                \"HTML-CSS\": {\n",
        "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
        "                }\n",
        "        });\n",
        "</script>\n"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 32,
       "text": [
        "<IPython.core.display.HTML object>"
       ]
      }
     ],
     "prompt_number": 32
    }
   ],
   "metadata": {}
  }
 ]
}