GenSamples.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mudigonda/Tools/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.\n",
      "  warnings.warn(self.msg_depr % (key, alt_key))\n"
     ]
    }
   ],
   "source": [
    "from matplotlib.patches import Ellipse\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import pymc3 as pm\n",
    "import scipy as sp\n",
    "import seaborn as sns\n",
    "from theano import tensor as T\n",
    "import theano"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.random.seed(3264602) # from random.org\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 100\n",
    "\n",
    "mu_actual = sp.stats.uniform.rvs(-5, 10, size=2)\n",
    "\n",
    "cov_actual_sqrt = sp.stats.uniform.rvs(0, 2, size=(2, 2))\n",
    "cov_actual = np.dot(cov_actual_sqrt.T, cov_actual_sqrt)\n",
    "\n",
    "x = sp.stats.multivariate_normal.rvs(mu_actual, cov_actual, size=N)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "var, U = np.linalg.eig(cov_actual)\n",
    "angle = 180. / np.pi * np.arccos(np.abs(U[0, 0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mudigonda/Tools/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py:892: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.\n",
      "  warnings.warn(self.msg_depr % (key, alt_key))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAFzCAYAAADrOKo/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8XHd95//X6DIjaXQZ2VJ8ix3Ft6PYsePYjh07jpPY\nCQEatiy7LeXXdmkTGuAHpC3Q7hZaQtsQCi3l2tJCE+jSdlm6sGzJAiUkJI4Tx9c4vupYdizfLcuW\n5ippNJezf4ykyLauo5k558y8n49HHnHiuXzOHFvv+d49lmUhIiIizlBmdwEiIiLyJgWziIiIgyiY\nRUREHETBLCIi4iAKZhEREQdRMIuIiDhI1sFsGMY9hmF0Gobx9lwWJCIiUsqyCmbDMBYCvw9sz205\nIiIipS3bFvN54D8C4RzWIiIiUvIqsnmSaZr9AIZh5LYaERGREjdhMBuG8QjwPsACPIP/ftw0zWfz\nXJuIiEjJmTCYTdN8Cnhqum9kWZbl8Xim+zIiIiJukVXoZdWVnc0bezweuroiOXg7+zQ317n+GqA4\nrqMYrgF0HU5SDNcAxXEdxXANkLmObGQ7K/vthmH8AngQeNIwjJ9m9e4iIiJylWwnf/0Y+HGOaxER\nESl52vlLRETEQRTMIiIiDqJgFhERcRAFs4iIiIMomEVERBzEMcGcTqfp7r6S03/S6fS472lZFp//\n/Gf44Acf5rHHPsDp06cAePLJP+W97/01HnvsAzz22AfYseNlkskkjz76KO9//29z6NDB4df4oz/6\nGF1dl0Z9/RdffD53H1AOvPHGCT7ykfcD8Ed/9HEAPvKR93Py5BtXPe611/byx3/8X/NaS3v7MZ5+\n+ht5fQ8RETfKxQYjOREM9rBt2y+ora3NyetFo1E2b76PGTNmjvmYl156kVgsxte//jTnzp3ly1/+\nAp///BcB+MAHPsyGDZuGH7tjx8usWbOGTZu28jd/82VuvXUFO3a8zOLFS2luvuG6175w4TzPPvvv\n3HPPlpxcT64M7b722c/+1aQely9LlixlyZKleX0PkVJiWRaxWIxw+BInT54jkRggmUwBFpZlXfU4\ny7IoKyujoqKSiooKKioqqK9vYPbsOfh8PvsuQgAHBTNAbW0t9fX1BXu/s2dPs2zZcgDmzbuRzs6L\nV/0BHikSidDU1MSMGTOJRCKk02n+9V//B08+OXrAffGLn+fo0SN8+9v/QDqd5vz5c1y8eIHf/u3f\n4X//7//FE098DoCHHrqfZ575OSdPvsGXvvSXeDweamr8fPKTj+P3v/klJZlM8pnPfJqLFy/g8/n4\n4z/+U3bv3smrr77MlStX+PSnn2Tbtl/w7LM/pby8nLvvvod3v/vX6eq6xJ/8yX/D6/WyaNGS4dcb\nel+AZ575IceOmcTjcf78z//iqut48cXn+e53/5mKigpaW5fxoQ/97lW//+STf0pFRQXhcJg/+7PP\n8vnPf4YLF86TTCZ55JH3s3r1Wnbv3slXv/rXzJzZxPz5CwgEGrn99jV8//vf44knPsdzzz3L9773\nL1RUVGAYrTz22Md4+ulvEI1GOH36FBcunOexxz7G+vUbpniHRYqLZVn09HRz6dIl+vp6GRgYYGAg\nPvhPgrKyMpqaGkgmM1+uy8ogsznj6F+0LStFIpFiYKCfnp4rHDlykIqKCqqrq/H761i8eAl1dYX7\nmSwZjunKtsPChYvZufNV0uk0p093cP78OYLBIADf//73+N3f/SCf/vQnCYWCzJo1i9OnT3PmzGnm\nzJnD//2//8b99z/Id77zLf7iL/6c9nbzqtd+z3t+k1WrVvNbv/U+IBOsX/vaNygrK7umNZr59Ze+\n9Ff84R9+ki996W+54471fP/737vq9X7602eYObOJr3/9Kd7xjv/I9u3bALh0qZOvfe0bJBIDvPDC\nc3z960/xta99gxdeeI5Llzr513/9Lvff/xa+8pW/o6mp6br3BZgxYyZf/erf8+CDb+df//W7w/+/\nr6+Pf/zHp/nqV/+er3717+nsvMihQweu+xwbGgI88cTnePbZn9LU1MyXv/x1nnzyL/nyl78AwN/9\n3df41Kee4K//+mu0t5vD1+/xeOjr6+Ob3/xbvvKVr/M3f/NNzp8/x759ewDo6urir/7qKzz22Mf4\n4Q+/P9nbKlI00uk0586dYd++3bz00gv8+7//X/bt20VPz2USiTgej4XP56Wuro6ZM2fQ2Bigqqpq\nyj1eHo+HqqoqZs6cSUNDA16vl4GBfnbs2M7zz/+MvXt30dfXl6erlGs5qsVcaHfeuZFDhw7w4Q8/\nyqJFS2hpaQEs3vrWX6K+voHFi5fwT//0bb71rW/ye7/3Bzz33E/40pf+ig996DH+4R/+jt/4jd/i\nwoXz/P7v/yGf+tR/43Of++KY7zXUMh/L0aOH+dznnsCyLJLJBK2tVz/eNE3Wrl0HwNatDwDwk588\nM/y4I0cOc/bsGR577ANYlkVfXz8XLlygo+MkW7ZkHn/77WvYuXPHde+9evXa4Rp37doB3A3AyZMn\nuHTpIh/96IeHu8kuXrzArbeuvOr5t9ySqeHgwQMcPLifAwf2Y1nWYFdaks7OCyxevGTwM7/rqrH/\nM2dOMX/+Any+KgBWrVo9/CVn5cpVADQ330Bvb2zcz0+kGCSTSc6ePcPly1309kaJxXqprq6ipqYG\nn8+LzzejYLV4PB4CgQAAicQA27Y9T319gBUrVlJbm90e0DI5JR3MAO973weGf/3ud7+TxsYZNDa+\n+Yd/06Z7+MIXMt27n/3sZ+nqivD009/gPe/5TS5evDg8JtPb2zvu+1RUVALXj90mk0kAqqqq+MpX\n/m7M55eXl2FZ109mq6ysGPx3JRs3buLjH/+jq37/n//525SVZd4znR69m35kTSN/XVnpxTCW8YUv\nfGXMuq6t4b/8l4fZuvUt177DqK8/9N8j60okklRVZUK6vLx8+P+PNcQg4mapVIqOjpNcuXKZ3t4o\nfX19+P1+qqqqhv9xgrKyMhobG7Esi5df3sbs2XNZuXJV3ueilKqS7so+frydz372zwB49dVXMIxb\nAPjjP/5Dzp8/B2RmKC9cuGj4OZcvd3Hu3Bluv30NM2bM4OLFC8Tj/Xi9V0+YKCsrI5VKXfeefr+f\nK1cuD79/X18m0BcvXsqrr74CwHPP/Wy4O3dIa+sy9u7N/L9XXtnOd77zrat+3zBa2bdvD/F4P5Zl\n8eUvf4GBgQHmz7+Jo0ePDF7LyNd8M+hef/01AA4fPshNN908/P8XLFjAqVMnh7v3n3rq77l8+fLo\nHyaZFve2bS8A0NPTzd///d8AMHPmTE6fPkUqlWLXrleves78+Tdx7tyZ4W6y/fv30dp6y5jvIeI0\nyWSSPXt2sWfPruEv2hM93jTbePnll3j22Z9w/vwZIE1NTQ0zZ850TBiPxuPx0NjYSDQa5vnnf0Yw\n2GN3SUXJUS3maDRa0NdatGgxlmXxO7/zXnw+H48//gQA73rXr/L4439EVVU11dU1fOITjw8/59vf\nfoqHH84sOVq1ajXf+96/8NhjH+S9733kqte+6aabOXasja9+9YtXzTRfvHgpVVVVfPCDj7BixUpm\nz54DwO/+7sf4/Oc/wz//8z8O1vKZq17v/vsfZM+eXXz4w49SWVnJJz/5p+ze/WbIzZo1m1/91ffw\noQ89Ojz5y+v18iu/8h4+9an/xrZtL7Bo0eIRr/jmOG9PTw8f//hjRKNRnnjic5w5cxoAn6+Kj3zk\no3z844/h9XpZutS4Zpz6alu2PMC+fXv54AcfJp22ePjhR4FMr8QnPvEHzJ07j5aWhZSVvfl9MPNZ\nPMZHP/physvLWblyFStW3Mbu3TvHv3kiDpBMJnn66W/Q3d0NwL59e3j44UepqLj6R2sikeD48Xau\nXOkiGo1QW1uL1+tlxozCdU3nUmVlJZWVlezatYNVq9Zyww3Xr0yR7HkK2EVojXe+Zjqdzvm3r0Cg\n8aoQmK5iOiO0kNexe/erzJ/fwuzZs/nLv3yS229fw/33Pzit19S9cJZiuI5srmHPnl38/Oc/G/45\nk06nuf/+t7B27ToGBgZobzfp7r5Mb2+M2to6Kisr81H6Vfx+H7FYPO/vMyQYDLFixW3MmTM3Z69Z\nDH+eAJqb67Lq63dMi7msrGzcNcfiXpYFn/jEx6ip8TNjxkzuvXer3SWJ5EUymeTEiXZ6e2P09sZo\naGjA5/MV9drgQKCBAwdeo66uPmf7UJQ6xwSzFK916+5k3bo77S5DJOdWrVrN3r176Oh4g1AohM/n\n495778Xr9VJVVbxhfK1AIMDOna+wZcsDmhCWAwpmEZEsxGIx2tqOcNNNCygrA5/Py5IlxlWrCUqF\nx+OhutrH/v37uP32NXaX43oKZhGRSbIsi46Ok5w9e5pYLEYg0MCMGTNcO4krlyoqKunsvEA8Hi/q\nrvtCUDCLiEwgGo3S1naEnp4reL1eqqur8fm8dpflOIFAgAMH9nPHHevtLsXVFMwiIqOwLIuTJ09w\n9uwZ+vt7aWgI0NDQYHdZjlZWVsaVK12k0+mcrogpNQpmEZER4vE4Bw++Tnf3Zaqrq6mq8pXURK7p\n8vv9nDx54qpDc2RqFMwiIkAoFGLHjh2Ew0EaGhqG94mWqfH5fFy61KlgngYFs4iUtM7Oixw7ZuLx\nJKisrKaxsdHuklwvEgnbXYKrKZhFpORYlsWpUx10dLxBOp2ktrYOv7+2oDtmFbNUKk1/f7+j9/12\nMgWziJSMdDrNsWNtnDt3hoqKCmpqqu0uqSj5/TVcvHiBlpabJ36wXEfBLCJFL5FIcPjwQS5duojf\n76eubvLnCadSqeEzwkt1A5GpqqqqoqenW8GcJQWziBStWCzG4cMH6e6+nNWErlQqxTPP/IhwOASA\naR7joYfeoXCehNGOvZXJUTCLSNGJRqMcPLifUChEY2Mg65252ttNwuHQ8JrccDhEe7tJa+uyXJZb\nlNJpBXO2FMwiUjR6e3s5cOA1QqHQ4HaZmmFtF7WYs6dgFhHXi8fj7N+/j2Cwm0AgQGPj2F3WUxkz\nXrLEwDSPDXdl19c3sGSJkdviRa6hYBYR10okErz++j6uXMmMIU+0BnmqY8bl5eU89NA7NPkrCx6P\ntuTMloJZRFwnmUxy8ODrdHZemFQgD8lmzLi8vFxjylmoqHgzXpLJJPv37wMyZ1iP/D25nj4dEXGN\nVCrFoUMHuHjxPPX19Tpu0cGGwjeZTPL009+gu7sbgH379vDww48qnMeR1SdjGEY58BSwCCgHPm6a\n5iu5LExEZEg6nebIkUOcP3+W2trarLfN1JhxYViWNdzlv3//Prq7u4d7Kbq7u9m/fx9r166zs0RH\ny/Yry28CUdM07zYMYxnwLUAHcIpITlmWxdGjRzh79hR+f820D5bQmHFhJJMJqqtr7C7DtbIN5u8A\n/zL46y5A/UkikjOWZXHixHFOnjxOdXV1Tk960phx/sXjcerq6oHMmPK+fXuGu7JnzJjBqlWr7SzP\n8bIKZtM0U8DQIrXf482QFhGZlosXL3DkyEEqKsppaGiwuxzJQjw+QFNTE5AZa3744Uc1+WsKJvx0\nDMN4BHgfYAGewX8/bprms4ZhfAi4HXhHXqsUkaIXDofYv38ficTAlPayFueprq65aoigoqJCY8pT\n4LEsK6snDgb2fwJ+2TTNxCSekt0biUhRGxgYYOfOnVy+fJlAIIDH47G7JJmmiooKNm3aZHcZTpDV\nH+ZsZ2UvBN4PbJ5kKAPQ1RXJ5u0co7m5zvXXAMVxHcVwDVDa12FZFocPH+TcuTPU19fj9dbQ2zuQ\npwon5vf7iuI8ZruvI5EYoLl5zrT+XBfT34tsZNvR/wiZCV8/NgxjqHv7LaZpJrN8PREpIadOneLY\nsSPU1NRkvfRJnCkSibJ2bYvdZbhatpO/Pgl8Mse1iEiR6+np5sCB/VhWShO7ilRlpRefz2d3Ga6m\nqXEikncDAwO89toewuEwDQ31dpcjeVRT47e7BNdTMItI3liWxZEjhzh79jQNDQ0K5SKXSqWoq1NP\nyHQpmEUkL86ePUNb22Gqqqo0jlwiIpEIt956m91luJ6CWURyqre3l717d5NIxKmvVwu5lJSXlw/v\n+CXZUzCLSE4MLX86f/4MDQ0BfL5Ku0uSAkqn0zQ2zrS7jKKgYBaRabt48SIvvPAyXm8lgYC6rUtR\nKBRi06bb7S6jKCiYRSRriUSCPXt2YVlx6upq7S5HbFRdXYPfrxnZuaBgFpGsHD9+jBMn2qmvr6eh\noaEods2S7CQSA8yaNdvuMoqGgllEpiQUCrJv3x7KyjyabS1AZrevO++82+4yioaCWUQmJZ1O89pr\ne7lypSun5yOL+zU0BK46TUqmR8EsIhM6ffoUpnkEv9+vUJarxGIxli69xe4yioqCWUTGlFmTvItE\nIqG9rWVUqVSauXPn2V1GUVEwi8h1Rh7JGAgE8Pm8dpckDmRZFoFAo87QzjEFs4hc5cqVK+zfvwef\nz6fJXeNIpVK0t5sALFlilOQYaygUZOPGzXaXUXQUzCICZFo/+/fv4/LlS+q2nkAqleKZZ35EOBwC\nwDSP8dBD7yi5cPb766itrbO7jKJTZncBImK/y5cv8/Of/zt9fVFXhHIqlaKt7QhtbUdIpVIFf//2\ndpNwOERZWRllZWWEw6Hh1nOpCIfDGMYyu8soSmoxi5SwoSVQ3d1XXHMko1qrzuD1+mhqarK7jKKk\nFrNIierqusTzz/+M/v6YI0N5rFaxE1qrS5YY1Nc3kE6nSafT1Nc3sGSJUdAa7BSLxVi8eIndZRQt\ntZhFSkymlbyHYLDbsccyOr1VXF5ezkMPvaNkJ395PB7mzr3R7jKKllrMIiWks7NzsJXc5+hzc8dr\nFTultVpeXk5r6zJaW5flNZTtHk+/Vm9vLy0ti+wuo6ipxSxSAlKpFK+9todQqMexreTJKqXWqhN7\nDtJpi5aWm217/1KgFrNIkbt48TzPPfczBgbijm4ljzRRq7hQrVW7OWE8faTe3l4WLtTYcr6pxSxS\npFKpFPv27SEUChIIOH8J1Eil1Cp2k1QqTUtLi91lFD0Fs0gROn/+HIcPH6C2tpb6enduADHUKi5l\nS5YYmOax4a5sO2d/x2IxFi1aast7lxoFs0gRSaVS7Nmzi2g07IqNQmR8Tuo5sCxLreUCUTCLFInL\nly/z2mu78fv91NW5s5Us13NCz0EkEmHZslttraGUKJhFXM6yLA4efJ3OzgtqJU+CDp+YGsuy8Hp9\nzJmjox0LRcEs4mLRaJRdu16hsrJSoTwJhVx+VCxfAILBIHfddY/dZZQUBbOIS5040c6JE+00NDTo\nPNxJGrn8CBhefpTrrmInrj/ORjKZ5IYbZlNbW2t3KSVF65hFXCaRSLB9+4ucOdNBIBBQKDuQ09Yf\nZysSibJy5Sq7yyg5CmYRFzl//hzPP/8s5eVl1NT47S7HdZyynacb9Pb2smTJUte18ouBurJFXCCd\nTrNv325CoR4aGwN2l+NahVp+5KT1x9myLFi4cLHdZZQkBbOIwwWDPezevRO/v9o1W2o6WSGWHzlp\n/XE2QqEwt9++1u4ySlZWwWwYRjPwj0AVUAl81DTN3bksTETg6NHDnDlzikBArWS3ccL642yk02lq\na2tpamqyu5SSle0Y828A/900zS3AJ4EncleSiPT19fHCC8/R1XVRoZwnTjtO0SlCoRCrV99hdxkl\nLasWs2maXxzxnwuAM7kpR0ROneqgre2wZlzn0WjLmd797v9sc1X2GxgYYN68+fh8PrtLKWlZjzEb\nhjEL+BFQC2zJWUUiJSqdTrNr16v09kZpbGy0u5yiNtp65ra2NlpaSnuyU19fP3fdtcLuMkrehMFs\nGMYjwPsAC/AM/vtx0zSfBdYZhvFWMuPND+azUJFiFgoF2blzB7W1fm3mILYIhyPceutK9dI4gMey\nrCk/yTCMzcAB0zSDg//dZZpm8wRPm/obiZSAY8eOcfjwYbWSCyiVSvHDH/6QUCjTld3Q0MA73/lO\nV82cnozMOHobAK2trWNeXyqVwufzsWnTpkKWVwqy+paTbVf2u4DbgS8bhrECOD2ZJ3V1RbJ8O2do\nbq5z/TVAcVxHMVxDOp2mrW0/nZ1X8Pv9xGJxu0vKmt/vc139DzzwtuuWM7ntGkYzdC+uHUd//fVD\nY24LGgwG2br1Qcf8nSqGv9+QuY5sZBvMfw78o2EY7wK8wAezfB2RkjTUdT1r1gz8fu3gZQe3Lmea\nrMnuC57pwr6t6HoL3CzbWdlXgIdyXItISThxop3jx9tpbAxQUVFBPO6spTrFciqSTCyVSlFXV8/c\nuTrS0Um085dIgaTTaXbufIX+/j7HbqtZLKciyeS2BY1EImzdqnm7TqNgFimAkbOundx1XahjESX/\nJtoWNBJRF7ZTKZhF8mxk17VIIY01jp5KpaitVRe2U+nYR5E8SaVS7NixndOnO1wTyvk+FlHbYDpD\nOBxh7dp1dpchY1CLWSQPgsEedu3aQW1traO7rq+Vz1ORNH7tDJFIhBUr1IXtZApmkRw7fvwYJ04c\nd00r+Vr5Wkak8Wv7qQvbHRTMIjmSSqXYtWsH8Xi/a0NZils4HOH++zUL2+k0xiySA+FwiOee+xlg\nUVNTY3c5thprHHnhwsX09/fS1XWJVCqV8/FrGZ+6sN1DLWaRaTp9+hRHjx7SXteMPY4M8JOf/Biv\nt4poNEZ/fz+/9mtvV0gUSDKZoKEhoC5sl1Awi2TJsiwOHNhPV9clhfKgscaRh35dUVHBDTfMIp1O\n88YbxzW+PEnT2Y3Nsiz6+uJs2LA5X+VJjimYRbKQTCZ55ZWX8Higvj67jepFRnNtCAPTms0eCoXY\nsOFuHefoIgpmkSkKBnvYuXMH9fV16oq9xnjbQE60PaSMPhSwePGirGezR6NRli5tpa5OXx7dRMEs\nMgUdHSc5duyorbOunXzIxHjroPO1PrqYjDYUcObMmaxeK5lM4vfX0dKyMJclSgEomEUmwbIsXntt\nL8HgFQIBe0N5upt05DvYx1oHXezHLObL/Pk3EY8PTLm3obe3ly1b3pLv8iQPFMwiE0gkEmzfvo2K\nijJqa+3rEkylUjz77L/zxhvHaWpqpqysbMqbdGj3LWcbbSigtbWV1tbWKX2Z6ukJcuedG4db3uIu\nCmaRcVy5coU9e3bS0FBv6w+5oUA9efIEFy5coKenO6sxWu2+5WzjDQVM9h7FYjEWLVpMQ4M2uXEr\nBbPIGE6caHfM1ppDgdrU1Ex3dzf9/f1cvtzFwoWLNYmqyEynyz+zcYufxYuX5rgqKST1c4hcw7Is\n9uzZyalTJwkEGib1nEKdmuTxeFi61GDu3LnccsuyKXdD5/v0qInodKn8ikaj3H333XaXIdOkFrPI\nCP39/bzyyjZ8Ph+1tbWTek4hxm2vHXu8+eZFPPDAg1N+j3yeHjURjW/nVygUYs2a9fo8i4CCWWRQ\nV9cl9u3bTSAQmNJmDIUYt81loNo1O1rj2/nT19fPggUtzJgxw+5SJAcUzCJAe/sxOjqOO3prTS03\nktGkUinKysoxjFvsLkVyRGPMUtIsy2Lv3t2cPXsq61msdo/busV0PyeNT48uHI6wYcNddpchOaQW\ns5SsVCrF9u0vUlbmwe/3Z/06do3bOnkHsNFM9DmNdz0anx5dMBhi/foNVFToR3kx0d2UkhSLxXj5\n5Repq8vNfteF7mZ2a1CN9TlNdD0an75eJBKltXUZgYBzh18kO+rKlpLT1XWJ7dtfoKGhwfFBNpaR\nQTVyBzC3Krbrybd4PE5z8w0sWHCT3aVIHiiYpaS88cYJ9u/fS2Njo47BcxGN478plUphWR5Wrlxl\ndymSJ+rKlpJx4MB+Ll/upKFhcpuGONloeyovXLiYtrYjw7/vpN6AicbDxzsuEuxdf+0klmURiUTZ\nulWHUxQzBbMUvXQ6zY4d20mlkrYeQpFL1wbVwoWL+clPfuzIMefJjIdPJni1XCyziciGDZsccV8l\nf9SVLUWtr6+P559/Fo8Hqqqq7C4np4aCqrV1GW+8cdyxY7STHT8eeT0KnutFImFuvfU26urq7S5F\n8kwtZilamZOhXp3yTl4iTtPX18+cOTcyd+48u0uRAlCLWYrS6dOn2Lt3V8lM8nLy5Cgn1+YGyWQS\nr9fLsmW32l2KFIhazFJ0Dh06QGfnhUmfDFUMnDw5ysm1OZ1lWcRivZrsVWIUzFI0LMvi1VdfZmBg\ngLq64pjkNRVOnhw13sYiCuyxBYNB7r77vuGNVaQ0TCuYDcOYBRwF3mma5rbclCQydQMDA7z00gtU\nV/uori6uSV654rQQdOvuZYUSDodZtWo1NTU1dpciBTbdFvPngRO5KEQkW6FQkFdffZmGhga1LMbg\nxBDUNptj6+3tZcGCFmbNmmN3KWKDrH+KGYZxHxAGDuauHJGp6ezsZOfOV2hsbFQoj2O6W17qZKfC\nicfj1NbWsXRpq92liE2yajEbhlEJfAr4ZeDLOa1IZJJOnerg2LEjBALZHdcok5Ov1vZEu32Vosx2\nmxZr1qyzuxSx0YTBbBjGI8D7AAvwDP77p8A3TdMMG4bB4P8XKZjDhw9z/LiZ9RnKpWY6IZivLmfN\n1r6aZVlEozG2bHmgJJb4ydg8lmVN+UmGYWwn0w3uARYBl4BfMU3z6DhPm/obiYxi165dXLp0idra\nWrtLcZVMd3QbAK2trZMOwcOHD7Nz587hYE6n06xfv57ly5fnrdZSFAwGeetb31p0O9SVuKy+YWUV\nzCMZhvEt4FuTmJVtdXVFpvVedmtursPt1wDuvQ7LstixYzvJZIKZMxuIxeK21JHL2c1+v8+265is\na7uy6+sbruvKdst1jHff7LyGYDDI+vUbc9ID5Na/3yMVwzUANDfXZRXMuVjHrJaw5F0qlWLbtl/g\n9Vba2qJw4uzmfCuGLmcn37dwOMRtt92uYRkZNu1gNk3z4VwUIjKWvr4+tm9/kdpav+0/SEt1iY+T\nNy+ZDKfet0gkyuLFhpZFyVW085c4WjgcYseOlwkEGjQhRopKLNbL7NlzaGlZaHcp4jBa+CmO1dnZ\nyauvvkyTH1JEAAAgAElEQVRjo3NOh9KBDO7ktPs2MDCA31/L8uUrbKtBnEstZnGk06dPYZqHHbdG\nuRjGW0uRk+5bKpUilUpzxx3rbXl/cT4FsziOaR7l7NnTjp0MY8d4q9P2uXajXNy36d4HrVWWyVAw\ni6Ps37+Pnp4rJXk61FicPKO4lOTiPgSDIe65Z4vunYxLY8ziCENrlMPhIH6/3+5ybHXtvtTT3ed6\nrNeVqZnufQgGg9x550ZtICITUotZbOeUNcpOMFqrbPHiJXl53WJsdTu1y19rlWUqFMxiq4GBAV58\n8XlHrFF2gtHW20JmFvF0DntobzcJBnvo6ekBMttqOmEdby7l+8tHtvuNRyIRFi3SWmWZPAWz2Ka3\nt5ft21+gvr5eRzaOIxczilOpFMePtzMwMADAlSuXWbt2bc5rtVO+NxHJ5j5Eo1FuvPEmbr5Za5Vl\n8hTMYotIJMwrr2zXxiHXGKtVNv0ZxaN98dGXoamayn3o7e3lhhtmYxg6V1mmRsEsBdfd3c3u3Tto\nbGy0uxTHydd62/JyD0uWLKWnpxuAxsYZlJcX1xciJ53v3N/fRyAwg1tvXWnL+4u7KZiloDo7O3n9\n9b0K5XHkY530UGgN9U7YvfNVPjhlE5F4PE5VlZ9Vq1YX/L2lOCiYpWDOnz/HoUMHHLebl1Pkc0bx\nUGi1tbVx5swp5s+fP63Xu7ZWp7D7sI1EIkF5eYV29ZJpUTBLQXR0dNDe3kYg0GB3KVMynbCcynML\ntZypvd2ko+MkbW1HaW9v5x3v+OWsJpJdW+u73/2fc1qnGyWTSZLJFJs3b9a8CZkWBbPk3fHjxzh1\n6g0aGurtLmVKphOWU31uIY4lbGtrY+/ePSNmZl9hyRKD5cuXT+l1Rqu1ra2NlpbFOavVbVKpFP39\n/dx77/0KZZk2TcuUvDLNo5w61UFdnbtCGaa301OuduuC3O3YdebMKeLxOB6PB4/HQzwe58yZU1m/\n3kRKZaexdDpNNBpj8+YtWvYnOaEWs+TNoUMHuHTpInV1tXaX4nhjzSjOZRf3/Pnz8fl8wy1mn8+X\n1VjzaLW2trbS358cfkyp7DRmWRbhcJh7772figr9OJXc0J8kyYsDB/bT3X2Z2lr3hvJ0lt9M9blj\nzShuazuSsy7u1tZlrF69lo6ONwBoaVmYdVf54sVLhieRtbYuGwzcN4O5EF3zI9mxFadlWQSDQTZv\n3oLX6837+0npUDBLzu3du4toNEpNTY3dpUzLdJbfZPPcfM8oLi8v5z/8h1+eVE1jBd21LeF4fCAn\nNU93kp0drfOeniB33bWZ6urqvL6PlB4Fs+SMZVns3LmDgYF+qqszh1E49VCByZpOWOYiaHO9acZk\nahov6MZqCa9Zc3vWdU83WAvdOofMSVHr12/Q8aSSFwpmyQnLsnjlle2k08nhE6LcPs7ohC8Vdmya\nkYugm0rddgTrdASDIVavvoPGxhl2lyJFSsEs02ZZFi+99CLl5R58Pt/w/3fbD9yRnPSlwu5NM0aa\nSku4UHUXcivOUCjEihW30dTUnJfXFwEFs0yTZVls2/YLKivLqaiotLucnLn2S0Uw2MOzz/47LS0t\nruySn4rxgi4fLfjpBmuhehXC4RDLlt3KnDlzc/7aIiMpmCVrE4Wykw4VmI50Os3x4+309HTT1dXl\nui75qZoo6HLdEs5FsOa7dR4Oh7jllluZN296W5mKTIaCWbKS6b4ev6XslEMFsjHyS8Xly10ANDU1\n4/F4XNUln61Cd587qbv+WgplKTQFs0zZUChXVEzcfe3kH7jjGfmloqOjg0uXLmmrxRKkUBY7aP84\nmZKphLLbDX2peOCBB2loCJBOp0mn067tkpepUSiLXdRilkkrpVAeyc1d8qUgH8vaFMpiJwWzTEqp\nhvIQt3bJF7t8LGtTKIvd1JUtEyr1UJbpyecpU7k8xQsUyuIMajHLuDKh/ELJhLITdvsqJk7aqGUi\noVAIw1iuUBbbqcUsY3ozlMtyFspOPqN3KER2797N7t27eeaZHzmuRrfJdYv2WkuWGNTXN0x7Yl44\nHOa2225j/vwFOatNJFtqMcuo8hXKTm49uXkL0VKVi4l54XCmpdzS0kJXVyQfZYpMSVbBbBjGe4E/\nB44P/q9nTdP8bM6qElvlI5RBwVeKCrH723Qm5oVCYVpbl6ulLI4ynRbzd03T/MOcVSKOkK9QdgO3\nbSHqhvFwJy81y4TyMoWyOI66smVYvkM52+ArVAA5OUSu5fRhgZGcuNRMoSxONp1gvtcwjB8DlcAf\nmKa5P0c1iQ0sy2L79m15bSlnE3xjBVC+ODFERqNhgeyFwyF1X4ujTRjMhmE8ArwPsADP4L//B/C4\naZo/MQzjTuC/Aysneq3m5rrpVesAxXANcP11vPjii9TWevF6vXl/7zVrbp/0Yw8fPkx/f4yqqkxd\n/f0xzpw5yfLly/H7fRM82x2yuY6qKi9eb8VwMKfTaaqqvLZ+Jm64H6FQiE2bNjB//uhLoor177cb\nFcM1ZGvCYDZN8yngqXF+/1XDMJoMw/CYpmmN91pun/HY3Fzn+muA669j795d9PbG8Pl8JBJxGyu7\nXn//AAMDyasCqL9/AIBYrLC15qNL3e/3ZXUd8+ffzOuvH7pqWGD+/JsL/pkMyfY6CikYDHHbbbdT\nVRUY9e9xsf79dqNiuAbI/stFtrOy/wA4Y5rmdw3DuBXomiiUxZkOHNhPJBKhpqba7lJG5ZQJWUNd\n6qFQkO7uK7z00ou8972PFKSHYTRuGg93gmAwyOrV62hqarK7FJEJZTvG/C/AdwzD+ABQDjySu5Kk\nUI4ePUx3dxc1NX67SxmTUwKovd0kFApy/Hg78Xgcy7L49ref5pFHfifn9Uy2Ze6W8XC79fT0sG7d\nBhobZ9hdisikZBXMpmmeA7bkuBYpoPb2Y1y4cJbaWueP4zglgLq7rxCPx4fPZY7FojmfcOWm2dZO\nZ1kWwWCQjRvvpq6u3u5yRCZNW3K6TDKZZM+eXezZs4tkMpnVa5w4cYLTp99wRSg7xZIlBn5/DZZl\nYVkWPp+XQKCRjo6OnG4vmu8tLEvFUChv2nSvQllcR+uYXSSZTPL009+gu7sbgH379vDww49SUTH5\n23ju3FlOnTqmH1ZTVF5eznvf+wjf/vbTxGJRAoFGzp07R3l5OV1dXa5v2bphs5LJSqfThMNh7rln\nK1VVVXaXIzJlajG7yP79++ju7h5uTXV3d7N//75JP7+zs5PDhw/S0NCQxyqLl9fr5ZFHfoe3ve3t\nzJkzhwULFlBeXp7Tlu1kD2XI5WEgxXR4RyqVIhqNcd99DyiUxbXUYi4R3d3dvP76XgKBgN2luNrI\n8e6urq68vP5Ek91yPQ5dLJuVJBIJ4vEB7rvvfle3+EXUYnaRVatWM2PGjOHW1IwZM1i1avWEzwuH\nQ+zevUOhnEO5Om5wNEPh39q6bNSA0Tj09QYGBkil0tx771aFsrieWswuUlFRwcMPPzrcfb1q1eoJ\nx5djsRg7dmynsbGxECXm3dBYaFWVl/nzb7bth7BTlnHlglPWimcrHo9TUVHJXXfdNTxjXsTNFMwu\nU1FRwdq16yb12P7+frZvf7FoWsoju3C93gpef/2QrROu7FrGlesgdfOXjN7ePurq6lizZnJ/J0Tc\nQMFcpJLJJNu2PU8g0JD3VkShZvSO1YXr9LHQXH8++QhSp6wVn4pYLMbMmc2sXLnK7lJEckrBXITS\n6TQvvvg8dXV1BQllbYgxtnx9Pm4M0lyKRiPMm7cAw7jF7lJEck6Tv4pM5vjGF6muripIOBZyIlI+\nJ1zliyZq5V4oFOammxYqlKVoqcVcZPbu3YXHw5Q2HXGLkV24dk/+EnsEgyGWL1/BvHk32l2KSN6o\nxVxEjhw5RDQaxucr3Lm4hW7FDnXhLl++3BWh7MZWvhNZlkVPTw9r1qxTKEvRK75mVYnq6DjJhQvn\nqKsr7P7XTpjRW6jJZ9m8jxM+H7dLp9OEQmHuuuseamtr7S5HJO8UzEWgs7OT9vY227batHMiUqEm\nn03nfUp9otZ0JJNJ+vr6uO+++207+1qk0NSV7XLhcIjXX99bsvtfF2pylSZxFV48HieZTHHffQ8o\nlKWkqMXsYgMDA+zY8TKNjcWxgYjIkN7ePmpq/Nxxx3rt5iUlRy1ml0qn02zb9gsCgdJsKQ8p1OQq\nTeIqnGg0QlNTM+vW3alQlpKkFrMLZdYqb6Omprrkf3AVanLVRO9TTOcZ2ykUCrNw4SIWLVpidyki\ntlEwu9DevbvxeKyiXKucjUJNrhrrfbT7WW5ojbJIhrqyXebo0cNEo6G8rlVOpVK0tR2hre0IqVQq\nb+9TLDQxbHq0RlnkampyuUhHRwfnz5/N61plt7b+1JXsTlqjLHI9tZhdorOzk+PHj+Z9AxE3tv6G\nvkzs3r2b3bt388wzPypoS18Tw7KTTCaJxWLcd9/9CmWREdRidoFIJMz+/XtobGy0uxRHGvllAij4\ncZDa3Wvq4vE4lgX33ffA8H0TkQwFs8Mlk0l27NhesFBessTANI8Nd2Wr9Tc52t1r8np7e6mpqdUa\nZZExKJgdzLIsXnnlpYLuf+3G1p++TLhHJBJh3rz5+hIjMg4Fs4O9/vpreDwUPBjd1vpz45eJUhQM\nBgeXQ823uxQRR1MwO1RHxxt0d3dRW1vY06Lcym1fJkqJZVkEgyHWr99AIKB5EiITUTA7UHd3N8eO\nma7dblNLl2RIIpEgHo/rdCiRKVAwO8zAwAC7d7/q2oMp7FgHrS8CztTX14fX6+Pee+/XzGuRKVAw\nO4hlWbz88jYaGurtLiVrhV665NYNUYpdJBJl1qzZ3HrrSrtLEXEdBbOD7N27i8rKCrUupqC93SQU\nCtLT0w1kvtwUcg2zXC8YDLF06S20tLTYXYqIKykBHOL48WNEImHXj8MVehesVMqivf0Y58+f5/z5\n87S3HyOVsvL2fjK2kXteK5RFspd1i9kwjI8Dvw4MAP+/aZp7c1ZVienqukRHxwnq69052WukXC5d\nmtzYcXqS/0/yKZlMEolE2bx5C9XV1XaXI+JqWQWzYRjLgF8FVgO3Ab8MKJiz0NfXx759u/O2s5cd\nE6NysXRpsmPH5eXlLF68hJ6eHgAaGxtzdo2aVDY58Xicqio/W7Zoe02RXMi2xfwQ8D3TNC1g/+A/\nMkXpdJpXXnmJQCA/M7DHCjc3mOwksqFdv4Yel6uuc00qm5xYLMqMGc1s3bqZrq6I3eWIFIVsg7kF\nSBmG8ZPB1/iYaZoHclZVidi161Wqqnx52y94rHBbs+b2vLyfHfK165fdB2O4QTgcZuHCxSxcuNju\nUkSKyoTBbBjGI8D7gKEZNR7gBuCnpmm+zTCMu4B/ANblrcoiZJpH6e/vpaamxu5SHGkq+19r16/C\nyuzkFWTVqrXMmjXL7nJEio7HsqY+g9UwjMeBNtM0/+fgf3eapjnR31BNlR3U2dnJjh078taFPSSV\nSvHDH/6QUCgTbg0NDbzzne90TXdsKpWira0NgNbW1oLWbfdnZ+e1jyeVShGNRrn//vvx+/12lyPi\ndFl1h2YbzOuBD5im+duGYbQC/2Sa5toJnma5fQyqublu2uNoyWSS5577WcF29hptApPf7yMWixfk\n/fOlENdQiMlfo13HtePb9fUNjhjf7uvro6LCy513bryullz83bBbMVwDFMd1FMM1ADQ312UVzFmN\nMZumudMwjLcZhvEKmZbwh7J5nVL06qsvU19f2GMc1c2bHbs+OyeOb0ciEebMmceyZbfaVoNIqch6\nHbNpmp8GPp2zSkqAaR4lmUzi9eZvnaeW+EguDZ0MtXLlKubMmWt3OSIlQVtyFkh3dzenT5+koSF/\nXdha4lMcpjLxLZ+SySSxWC+bNt2j8WSRAlIwF0AqlWLPnp15n+zlxC5Qmbp8LQGbit7eXny+arZu\nfYs2DREpMAVzAezatYO6ulq7y5ACycVwgp1zA0KhEAsWtGAYt9jy/iKlTsGcZydOtNPf31eQrkCn\ndIGWMjcPJ2h9sogzKJjzKBwOceLEcQKBwhxO4YQu0Kkoxolqbh1OSCYT9PXFdQiFiAMomPMknU6z\nc+eOgoXyELcsj3Jzy7LYxGIx/P46tmzZnLftYUVk8jSrI0/27NlFTY1aHmMZ2bIsKysbblm6XaHP\no56uUCjIjTfexLp1dyqURRxCLeY86Og4STQaprZWE75KjVuGE4bGk9esWUdTU7Pd5YjICArmHItG\no5jmkbydr1wsinmimtOHExKJBPH4APfcs5Wqqiq7yxGRayiYc8iyLHbteiXv65WLgVtalsUmGo1S\nX9/Axo0aTxZxKgVzDh05cojKykr9wJukQrYsi3EG+FRYlkUoFGLp0lZaWhbaXY6IjEPBnCPhcIjz\n58/kdcvNUjXdUC31GeCJRILe3j42bNhEXV293eWIyAQUzDlgWRZ79uxyXSi7oRU5VqhOhVvXFudC\nNBqltraerVvv1taaIi6hYM6BTBe2uz5Kt7QixwrVNWtut7kyZ1PXtYh76Sv0NIXDIc6dO4vP57O7\nlCkp1nXEo3Hb2uLpSiQSRCJRNmzYpFAWcSF3NfMcZqgLu9C7e5WSXCyrKqUZ4ENd1xs2qOtaxK0U\nzNPgxi7sIW5ZR5yrUHX62uLpGtowxDBuUStZxOXcmSoOMNSF7dbWsptakcUeqtOVTCbp7e1j48a7\nNetapAgomLNgWRZ797q/CzvfgeeGWd+Fkq/PIhqNUFcXYMuWTeq6FikSCuYsHD16mIoKfXTjccus\n70LIx2fxZtf1Mlpabs5VqSLiAPqKPUWRSJizZ8+4bhZ2oZXSrO+J5PqzSCYThMMRNm68W6EsUoTU\n7JuCzCzsnTQ0aBxP7BGNRqivD7Bhg/a6FilWajFPweHDmS5s/UCcWKmtHR5PLj4Ly7Lo7u5m4UKD\ntWvX68+gSBFTi3mSEokE7e3t+Hw1dpfiCm6a9Z1v0/0s4vE4iUSSu+++D7/fn68yRcQhFMyTtH//\nPurr6+nvT9pdimtomdObsv0sQqEQs2fP5dZbV6qVLFIiFMyTEAz20NNzhXnzZgEKZsm/ZDJJLBZj\n9ep1zJw50+5yRKSAFMyT8PrrrxEIuOvkKHGvaDRCbW0DW7dqbbJIKVIwT+DUqVNYVsruMmQMxbSJ\nSTqdJhQKsXz5Cm68cYHd5YiITRTM40in05jmEdfv8FWsimkTk97ePsrKyrjvvgfwer12lyMiNlI/\n2TgOHTpAba1mYTtVsWxi0tMTZM6cuWzefJ9CWUTUYh5LX18fFy6co7Gx0e5SpEglEgmi0QQbNtxF\nfb16ZUQkQy3mMbz22l5N+HI4N29iEomEqamp5Zd+6ZcUyiJyFbWYR3HpUie9vTHq6+tsraOYJjbl\ngxs3MUmlUkQiEW67bTWzZs3W2mQRuU5WwWwYxieABwALKAdmmabZmsvC7GJZFocOve6IUC6WiU35\n5KZNTGKxGFVV1WzZ8hadTiYiY8rqp4Npmk8CTwIYhvFfgOZcFmWnEyfaHfFDc+TEJmB4YpNbQkje\n9OYRjbfQ0rLQ7nJExOGmlUCGYZQDHwTuzUk1NrMsi46Ok7a3lqV49Pb2Ul5ewebNW6iurra7HBFx\ngelO/noX8FPTNOO5KMZu7e0mPp8zlqu4eWKTZL7k9fT0sGDBzdx9970KZRGZtAlbzIZhPAK8j8x4\nsmfw34+bpvks8AjwaF4rLBDLsjh9uoP6emectezGiU2SEYvFqKz0cs89W6mqqrK7HBFxGY9lWVk9\n0TCMGmCnaZorJvmU7N6oQI4cOcLp06f1g1SyNrSl5ooVK1i8eLHd5YiI/bJadjGdMebbgLapPKGr\nKzKNt8sfy7LYv/8wDQ31xGJj98r7/b5xf98tiuE6nHYNsVgMr9fHunWb8Xq9k/6z3txc59i/F1NR\nDNdRDNcAxXEdxXANkLmObEwnmOcAl6bxfMdobzepqvLZXYYjXLt2GlB3+jjS6TThcBjDWMZNN7XY\nXY6IFIGsg9k0zR8AP8hhLbY5c+Y0dXW1dpdhu2vXTh892obHA5FI5pur1lJfLRqNUl1dw333PUBl\nZaXd5YhIkbB/wa7NTp8+RXl5ce1MOtUdw4Ye39HRQSgUHH58R8cbeDwemptvALSWesjQ7l233HIr\n8+freEYRya2SD+Y33jhOTU3xnCA11R3DRj7+8uUuuru7MYxWbRU5hlgsSlWVWskikj/F1VScoosX\nL5BKJe0uI6emehTiyMc3NWU2cLt8uYt0Ok1Ly0JaWm7WWmoyX2CCwSBLl97Chg2bFMoikjcl3WI2\nzaOcO3cG0MQmgLKyMhYvXsKsWbNoaWnR5K9BkUiE2to67XEtIgVRsj9lwuEwP/zhD0gmE0D+JjYV\n+oSoJUsMTPPYcFf2RK3cax8fCDTywAMPXlVnqY4pJ5MJenv7Wb78VubMmWd3OSJSIko2mH/0ox+S\nSAwMB1A+JjbZcULUVHcM0w5j17Msi1AozOzZc9iwYbPG20WkoEo2mCORUN5/4Np1QtRUj0J009GJ\n+dbbG6OsrIKNG++mtlZL6ESk8EoymC9dusTChQs5c+bMpLt8pbgNLYFauvQWWlputrscESlhJRnM\nJ0+eIBBozHsX7lTHe8UeoVCYxsZGtm59sOS78UXEfiUXzOl0mmCwh8bGQN67cDV+62zxeJxEIsnq\n1Xcwc+ZMu8sREQFKMJjfeOMENTWFOxtX47fOM3QKVEvLQpYu1WYqIuIsJRfMFy+ex+fTgRWlKhKJ\nUFPj5777HsDr9dpdjojIdUoqmOPxONFohBkzZthdihRYIpGgv7+f5ctXMHv2XLvLEREZU0kFc1vb\nERoaGuwuQwrIsiyCwRBz5sxl40atSRYR5yupYA6Fgjp3uYREo1EqK73cdddmrUkWEdcomWBOp9P0\n9sYUzCUgHo8Tjw9wyy3LmTfvRrvLERGZkpIJ5s7Oi5r0VeRSqRThcIQFC26itXWZuq1FxJVKJpgv\nXDhfVOcuy9VCoSCBwEy2bNE5ySLibiUTzL29Mbxe/cAuNtFojIqKCtat20hDQ8DuckREpq0kgtmy\nLGKxGF6vfnAXi3g8Pri3dSvz599kdzkiIjlTEsEcCgUpK5v8eGOhz1CWyRs6bOK225Zz++0LNI4s\nIkWnJIL59OlT1NXVTeqxdpyhLBMbWo/c1NTMli0bmDOnka6uiN1liYjkXJndBRRCLBYdPhN5IiPP\nUC4rKxs+Q1nsEw6HSKXSbNp0D2vW3EFFRUl8nxSRElUSP+F6e3upq9MGE24Ti8UoKytj1ao7aGpq\nsrscEZGCKIlgTqWSk36szlC2XyzWi2VZLF3ayo03zre7HBGRgir6YLYsi2Ry8sGsM5Tt09/fRyKR\nYtGipbS0tNhdjoiILYo+mOPxOB7P1IbSdYZyYfX39zMwkGDhwsXcfPNCzbQWkZJW9MHc19dLRYVa\nvE4Uj8fp74/T0nIzixcvVSCLiFACwRyJRPB6vXaXISMkEglisRgLFrRgGLcokEVERij6YM5sxalg\nzlYuN1tJJBJEo1FuvHEBGzbcPeklbCIipaTog7miooJ0Oq0JXFnI1WYryWSScDjC3LnzFMgiIhMo\n+mD2en0kk0mdOJSFkZutAMObrUx2YlwmkMPMmjWHdes2amMQEZFJKPqflD5f1ZSWS8n0DQwMEIvF\nmD17rgJZRGSKsvqJaRjGHOBpwEdmW8/fN03ztVwWlit+fw3JZMruMlxpqput9Pb2kkqlufHG+dx1\nV6smdYmIZCHbpsxHgR+YpvlNwzA2AE8Cb8tdWbnj99diWWm7y3ClyW62Eg6Hqaz0snDhEhYs0BGM\nIiLTkW0wdwEzB389Y/C/Hcnj8VBVVW13Ga411mYrlmURCoWoqfGzcuVqbrjhBhuqExEpPtkG85eA\nXYZhvBeoAzblrqTc8/l8dpdQNFKpFKFQiECgkY0b76a2dnLHaYqIyORMGMyGYTwCvA+wAM/gv38K\n/E/TND9rGMbbgS8A/ymfhU5HIDCDnp7LmoQ0DUNrkJuabuC++9bpy46ISJ54LMua8pMMw/gx8EnT\nNF8zDMMLtJumOdHg4tTfKEcSiQT/9m//xsyZMyd+sFwlFouRTqeZP38+K1as0HpwEZHJy2oGbLZN\nyHbgTuA1YB1wbDJP6uqKZPl2uVBJLBaf1iv4/b5pv4YTTHQdQ93VdXX1LFjQwo03zsfj8dDd3VvA\nKsfX3Fxn85+n3NB1OEcxXAMUx3UUwzVA5jqykW0wfxZ4yjCMXyXTEn4sy9cpmAULbuLUqZNUV2si\n2FiGWsdNTTewerW6q0VE7JBVMJumeRH4pRzXklctLQt5440TdpfhOEOt4/r6BgxjGXPnzrO7JBGR\nklZSs6EWLVrKyZPHqKnx212K7aLRKJYFzc3NrFmzXgd9iIg4REkF80033cS5c2dIpVIlOYlpaO9q\nv38Oy5evYNasOXaXJCIi1yi5Y37Wr99AJBK1u4yCSSaTdHd3098fp74+wNatD3LvvfcqlEVEHKqk\nWsyQ2cnqjjvWs2vXDhobG+0uJy8SiQSRSBi/v44ZM5pYu/ZOdVWLiLhEyQUzQGPjDDZs2MSOHdsJ\nBAJFcdjCwMAA0WiU2to6mppuYN26jTrqUkTEhUoymCFzUtJdd93Dq6++TFWVz5Utyng8TiwWo66u\ngVmz5rBhw+KSHDsXESkmJRvMALW1tWzd+hYOHnydzs4LNDQ02F3SuBKJBOFwBK/Xi9/vZ968+bS0\nLKSsrOSmCoiIFK2SDmbInD61cuUqgsGbOHz4ELFYhIaGBkd0b8fjcaLRKFVV1fj9fpqaZrFmzXqq\nqqrsLk1ERPKk5IN5SCDQyF133U00GuHQoYMEgz34fD5qamoKVkNfXy+9vX3U1Pjx+/3ceOMCbrxx\ngcaKRURKiIL5GrW1ddx550Ysy+L8+XOcO3eGSCTCwEAcjyeAZZVl3Zq2LIv+/n76+/uwLPB6vXi9\nvs0/nusAAAckSURBVMF/vMyZcyNz587TOLGISAlTMI/B4/Ewb96NzJt3I5BZDzwwEKa9/RQDAwMk\nkwmSySTJZJJUKnXV8zK57cHj8VBWVobX66Wy0jscvk1NzVRXVzuiu1xERJxFwTxJFRUVzJlzEzU1\nM+wuRUREipim84qIiDiIgllERMRBFMwiIiIOomAWERFxEAWziIiIgyiYRUREHETBLCIi4iAKZhER\nEQdRMIuIiDiIgllERMRBFMwiIiIOomAWERFxEAWziIiIgyiYRUREHETBLCIi4iAKZhEREQdRMIuI\niDiIgllERMRBFMwiIiIOomAWERFxEAWziIiIg1Rk8yTDMGqAfwRmAVHgt0zTvJTLwkREREpRti3m\nR4HjpmluBj4D/HnuShIRESld2QbzEmAXgGmaLwObclaRiIhICcs2mA8CbwcwDOMeYEHOKhIRESlh\nE44xG4bxCPA+wAI8g/9+HFhhGMY24EVA48siIiI54LEsa1ovYBiGH9hhmubK3JQkIiJSurLqyjYM\n422GYfzZ4H/+JvCT3JUkIiJSurJaLgX8AviQYRg7gCvAe3JXkoiISOmadle2iIiI5I52/hIREXEQ\nBbOIiIiDKJhFREQcJNvJX+MyDKOZzF7aVUAl8FHTNHdf85gE8BJvro3eapqmowa8J3kdvw78LpAC\nvmma5tMFL3QchmGUA08Bi4By4OOmab5yzWPccC8mcx2OvhdDBjfl+R7w26Zp/niU33fD/ZjoGhx/\nLwzDqAC+DdwEJMlcS8c1j3HsvTAM46+BO4E08Humae4Z8Xv3k9kuOQn8xDTNJ+ypcmITXMdJ4PTg\n71nAr5umecGWQidgGMatwA+BvzZN82+v+b0p3Y+8BDPwG8B/N03zu4ZhbAaeAB685jE9pmluydP7\n58q41zF4mMefAGvJfOC7DcP4gWmaQVuqHd1vAlHTNO82DGMZ8C1g/TWPccO9GPc6XHIvMAxjIfD7\nwPZxHubo+zHRNbjlXgD/H5nP+jcMw3gA+Avg1655jCPvxeDPo8WmaW40DKMVeBrYOOIhXwYeAC4A\nLxqG8b9M02yzodRxTeI6LOCtpmn22VLgJA3+mf8K8PMxHjKl+5GXrmzTNL9omuZ3B/9zAXBmlId5\n8vHeuTSJ61gP7DJNM2qaZj+ZH1R3FbLGSfgO8NHBX3cBM0Z5jOPvBRNfhxvuBcB54D8C4XEe4/T7\nMdE1uOVebAX+9+Cvf87oNTr1Xmwl0zpj8Ad8wDCMWgDDMG4GrpimeX6wdf/jwcc70ZjXMciDc+/B\nSP3A28gE71WyuR/5ajFjGMYs4EdALTDaN84qwzD+iUw30g9M0/xivmqZjgmuYzaZkBjSBcwpUGmT\nYppmikx3IsDvAf8yysMcfy8mcR2OvxcAg0GFYRjjPczR92MS1+CKe8GIOk3TtAzDSBuGUWGaZnLE\nY5x6L2YDe0b89+XB/3ec6z//S8DCwpU2JeNdx5C/Gwy3l0zT/EQhi5ss0zTTQHyMvxNTvh/TDuax\n9tI2TfNZYJ1hGG8lM057bVf2x4B/Gvz1NsMwXjRNc99068nWNK5jJFu/2Y13DYZhfAi4HXjHKE91\nzb2Y4DpGsv1b9gR/psbjmPsxjWsYyWn3AjI1rbvmYaP1IDrmXkxgvM/Y9s9/Cq6t9U+AnwLdwP8x\nDONdpmn+oPBl5dSE92PawWya5lNkJuUMMwxjs2EYAdM0g6Zp/tQwjO+M8rxvjHj8c8AKwLY/8Fle\nx3mubgnMA3bkudQxjXYNMPxD6ZeAXx5seV77PMffC5jwOhx1L2Ds65jE8xxzP7K8BlfcC8MwnibT\nmjk4OBGMa1rLjroX1zhPpvYhc3mzG3W0z/98geqaqvGuA9M0h74UYRjGj8l8/m4L5infj3wtl3oX\n8F4AwzBWkJlVN8wwjKWGYfzz4K8ryIztHM5TLdMx7nUAO4G1hmHUD46LbCQzg9MxBifqvB94l2ma\niVF+3xX3YqLrwAX3YhTXfXN2y/0YYbRv/265F88CvzL46/9AZqvhYQ6/Fz8D/jOAYRirgXOmacYA\nTNM8BdQZhrFgsO6HBh/vRGNex+Cfn58ahlE5+Nh7gEP2lDklV/2dyOZ+5GVLTsMwZpLp9q0DvMBj\npmnuNgzjvwIvmKa50zCMvyAzZpsC/o9pmn+R80KmaZLX8S7gD8lM5//KiMlijmAYxmeAd5P5UjHU\nFfkWMl10broXk7kOR98LAMMw3g78AWCQGXe6YJrmW930d2OS1+CGe1EG/AOwhMzknd8yTfOcW+6F\nYRhPkgmrFPAhYDUQNE3z/xiGsQn4PJm/J//LQWPj15ngOj4C/BbQC7xmmuZjthU6jsEvFV8gMxch\nAZwD/g04mc390F7ZIiIiDqKdv0RERBxEwSwiIuIgCmYREREHUTCLiIg4iIJZRETEQRTMIiIiDqJg\nFhERcRAFs4iIiIP8PxOytdOdkZDiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc7376ea290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "blue = sns.color_palette()[0]\n",
    "\n",
    "e = Ellipse(mu_actual, 2 * np.sqrt(5.991 * var[0]), 2 * np.sqrt(5.991 * var[1]), angle=-angle)\n",
    "e.set_alpha(0.5)\n",
    "e.set_facecolor('gray')\n",
    "e.set_zorder(10);\n",
    "ax.add_artist(e);\n",
    "\n",
    "ax.scatter(x[:, 0], x[:, 1], c='k', alpha=0.5, zorder=11);\n",
    "\n",
    "rect = plt.Rectangle((0, 0), 1, 1, fc='gray', alpha=0.5)\n",
    "ax.legend([rect], ['95% true credible region'], loc=2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "ufunc 'isfinite' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-8-52ea76b2b3db>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      2\u001b[0m     \u001b[0mmu\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mDeterministic\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'mu'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mtheano\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshared\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmu_actual\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m     \u001b[0mtau\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mDeterministic\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'tau'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mtheano\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshared\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0minv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcov_actual\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 4\u001b[1;33m     \u001b[0mx_\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mMvNormal\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'x'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mmu\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mmu\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mtau\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mtau\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      5\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/mudigonda/Tools/anaconda/lib/python2.7/site-packages/pymc3/distributions/distribution.pyc\u001b[0m in \u001b[0;36m__new__\u001b[1;34m(cls, name, *args, **kwargs)\u001b[0m\n\u001b[0;32m     23\u001b[0m             \u001b[0mdata\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpop\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'observed'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     24\u001b[0m             \u001b[0mdist\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcls\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 25\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mmodel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mVar\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdist\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     26\u001b[0m         \u001b[1;32melif\u001b[0m \u001b[0mname\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     27\u001b[0m             \u001b[1;32mreturn\u001b[0m \u001b[0mobject\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__new__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcls\u001b[0m\u001b[1;33m)\u001b[0m  \u001b[1;31m# for pickle\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/mudigonda/Tools/anaconda/lib/python2.7/site-packages/pymc3/model.pyc\u001b[0m in \u001b[0;36mVar\u001b[1;34m(self, name, dist, data)\u001b[0m\n\u001b[0;32m    254\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mdata\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    255\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdist\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"transform\"\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 256\u001b[1;33m                 \u001b[0mvar\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mFreeRV\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdistribution\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdist\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmodel\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    257\u001b[0m                 \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfree_RVs\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvar\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    258\u001b[0m             \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/mudigonda/Tools/anaconda/lib/python2.7/site-packages/pymc3/model.pyc\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, type, owner, index, name, distribution, model)\u001b[0m\n\u001b[0;32m    477\u001b[0m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdistribution\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdistribution\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    478\u001b[0m             self.tag.test_value = np.ones(\n\u001b[1;32m--> 479\u001b[1;33m                 distribution.shape, distribution.dtype) * distribution.default()\n\u001b[0m\u001b[0;32m    480\u001b[0m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlogp_elemwiset\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdistribution\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlogp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    481\u001b[0m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmodel\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmodel\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/mudigonda/Tools/anaconda/lib/python2.7/site-packages/pymc3/distributions/distribution.pyc\u001b[0m in \u001b[0;36mdefault\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m     49\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     50\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mdefault\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 51\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_test_val\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtestval\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdefaults\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     52\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     53\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/mudigonda/Tools/anaconda/lib/python2.7/site-packages/pymc3/distributions/distribution.pyc\u001b[0m in \u001b[0;36mget_test_val\u001b[1;34m(self, val, defaults)\u001b[0m\n\u001b[0;32m     55\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mval\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     56\u001b[0m             \u001b[1;32mfor\u001b[0m \u001b[0mv\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mdefaults\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 57\u001b[1;33m                 \u001b[1;32mif\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mall\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0misfinite\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetattr_value\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     58\u001b[0m                     \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetattr_value\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     59\u001b[0m         \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mTypeError\u001b[0m: ufunc 'isfinite' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model:\n",
    "    mu = pm.Deterministic('mu',theano.shared(mu_actual))\n",
    "    tau = pm.Deterministic('tau',theano.shared(np.linalg.inv(cov_actual)))\n",
    "    x_ = pm.MvNormal('x',mu=mu,tau=tau)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_samples = 4000\n",
    "n_burn = 2000\n",
    "n_thin = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with model:\n",
    "    step = pm.Metropolis()\n",
    "    trace_ = pm.sample(n_samples, step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}