proportion_of_var_CI.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence interval for the proportion of variance : assume a beta distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.stats as sst\n",
    "from __future__ import print_function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# example of proportion of variance\n",
    "propE = .42\n",
    "propStd = .24\n",
    "propVar = propStd*propStd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Match the two coeff of a beta:\n",
    "http://en.wikipedia.org/wiki/Beta_distribution\n",
    "\n",
    "$$\n",
    "V(X) = \\frac{\\alpha\\beta}{(\\alpha + \\beta)^2(\\alpha + \\beta +1)}\n",
    "\\\\\n",
    "E(X) = \\frac{\\alpha}{(\\alpha + \\beta)}\n",
    "$$\n",
    "\n",
    "which gives us:\n",
    "\n",
    "$$\n",
    "\\beta = \\frac{\\phi - a}{\\psi a}\n",
    "\\\\\n",
    "\\alpha = \\phi \\beta\n",
    "\\\\\n",
    "\\textrm{With} \\;\\; \\phi = \\frac{E(x)}{1-E(x)} \\textrm{ , }\\;\\; \\psi = \\phi + 1\n",
    "\\;\\;  \\textrm{, and} \\;\\;\\; a = \\psi^2 V(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem is that we may end up with negative parameters for the beta distrib. That's not good ...\n",
    "So, in that case - I fudge the parameters to keep the mean of the distribution correct  ... but the variance can be way off of course - so see the check "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def beta_param(E,V, fudge=True):\n",
    "    \"\"\"\n",
    "    E: expected value (the mean)\n",
    "    V: variance\n",
    "    \"\"\"\n",
    "    \n",
    "    phi = E/(1-E)\n",
    "    psi = phi + 1\n",
    "    a = psi**2 * V\n",
    "    if fudge and (phi <= a):\n",
    "        beta = .001    \n",
    "        print(\"\\n Possible heavy fudging \\n\")\n",
    "    else:\n",
    "        beta = (phi - a)/(psi * a)\n",
    "    alpha = phi * beta\n",
    "    return alpha, beta\n",
    "\n",
    "def mean_var(alpha, beta):\n",
    "    \"\"\"\n",
    "    mean and var of beta\n",
    "    \"\"\"\n",
    "    E = alpha/(alpha + beta)\n",
    "    alpbet = alpha + beta\n",
    "    V = alpha*beta/(alpbet**2 * (alpbet +1))\n",
    "    return E,V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(propE, propV) = (.42, .24**2)\n",
    "\n",
    "# compute parameters of the beta distrib\n",
    "al, bet = beta_param(propE, propV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parameters alpha and betat of the distrib: 1.35625 1.87291666667\n"
     ]
    }
   ],
   "source": [
    "print(\"parameters alpha and betat of the distrib:\",al,bet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " approx mean and variance of the distrib: 0.4200,0.0576\n",
      " to be compared to: 0.4200,0.0576\n"
     ]
    }
   ],
   "source": [
    "print(\" approx mean and variance of the distrib: {:.4f},{:.4f}\".format(*(mean_var(al, bet))))\n",
    "print(\" to be compared to: {:.4f},{:.4f}\".format(propE, propV))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "now, what's a 95\\% CI for this beta distrib ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "betadist = sst.beta(al,bet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% Confidence Interval = [0.04 - 0.89]\n",
      "80% Confidence Interval = [0.11 - 0.76]\n"
     ]
    }
   ],
   "source": [
    "print(\"95% Confidence Interval = [{:3.2f} - {:3.2f}]\".format(betadist.ppf(.025),betadist.ppf(.975)))\n",
    "print(\"80% Confidence Interval = [{:3.2f} - {:3.2f}]\".format(betadist.ppf(.1),betadist.ppf(.9)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x3901d90>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD9CAYAAABHnDf0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlclOX6P/DPCGPuK8bRmemQgAIimKFIpo1bECnmUkI7\nEtJiafU7x05lQn0zqe/5nkw6SaaWa5ZLaOHUcRl3wS0wRcSC40Anjgikacn2/P64I4lYBpiZe5bP\n+/V6XoLz8DwXj3hxz3VvKkVRQEREzqWd7ACIiMjymNyJiJwQkzsRkRNicicickJM7kRETojJnYjI\nCTWb3GfOnLnC09OzePDgwScbO8doNOpvueWWE4GBgd/o9XqjRSMkIqIWUzU3zn3fvn2junTp8tPD\nDz+86uTJk4Prv15eXt5j5MiRB7788stwrVZbWFJS4uHh4VFitYiJiKhZzbbcR40ata9nz55ljb2+\nbt26+6dNm7ZJq9UWAgATOxGRfG2uuefl5fmWlpb2GjNmzO6QkJCjq1evfsgSgRERUeu5t/UClZWV\n6uPHjw/duXPnuKtXr3YKCws7NGLEiMO+vr55dc9TqVRc54CIqBUURVG19Gva3HLX6XSmO++886uO\nHTv+3Lt374ujR4/em5WVFdxIgDwUBQsWLJAeg70cfBZ8FnwWTR+t1ebkPnny5LT9+/ffXl1d7Xb1\n6tVOGRkZoQEBAafbel0iImq9ZssyMTEx6/fs2XNHSUmJh06nMyUlJS2orKxUA0BCQkKqn5/fmYiI\nCENQUFB2u3btauLj45cxuRMRydXsUEiL3UilUmx1L3tnNBqh1+tlh2EX+Cyu47O4js/iOpVKBaUV\nNXcmdyIiO9ba5M7lB4iInBCTOxGRE2JyJyJyQkzuREROiMmdiMgJMbkTETkhJnciIifE5E5E5ISY\n3ImInBCTOxGRE2JyJyJyQkzuREROiMmdiMgJMbkTETkhJnciIifE5E5E5ISY3ImInBCTOxGRE2o2\nuc+cOXOFp6dn8eDBg082dd6RI0eGubu7V23evHmq5cIjIqLWaDa5x8bGrjQYDBFNnVNdXe02b968\n5IiICENr9vojIiLLcm/uhFGjRu0rKCjwauqcJUuWPD19+vSNR44cGWaxyMhsNTXAf/4DXLwIXL4M\nXLok/qw9fvoJaNcOUKvF0b799Y87dAA8PABPT+DGG8XHbm6yvyMiaqtmk3tzioqKNGlpaZN37do1\n9siRI8NUKpViicDo92pqgNxcICcH+O47ceTniz/Pnwd69BCJuWtXoFs38Wft0aWL+PrKyutHRYX4\n8+efgZISoLgY+O9/gbIyoFcvkez79QN8fYGBA4EBA8SfOh2TP5EjaHNynzt37tuLFi16QaVSKYqi\nqJoqyyQmJv72sV6vh16vb+vtnVZpKZCRARw+LI7MTJF0Bw8G+vcH/P2Bu+8Gbr4Z8PICOnWyzH2r\nqq4n+8JCIC9P/EL57DPg7Fnx7sDbGxg0CBg2TBxDh4pfIETUdkajEUajsc3XUSlK8w3tgoICr0mT\nJm07efLk4Pqv9e/f/7vahF5SUuLRqVOnq8uWLYuPiora+rsbieTf5oCdVWUlsGcPsGULsHMn8P33\nQEgIEBYGjBgBhIaKsolsV66IhJ+dDRw5Io6TJ8UvmeHDRbIfOVL8ElKx94WozVQqFVrTl9nm5F5X\nbGzsykmTJm2bOnXq5gYCZHKv5+pV4MsvRUL/4gvRIp4yBYiMBAIDHaf8UVEBfPONeHdx5Aiwd6+o\n+48dC4wbJ46bb5YdJZFjam1yb7YsExMTs37Pnj13lJSUeOh0OlNSUtKCyspKNQAkJCSktiZYV1ZT\nA6SnAytWiBZ6SIhI6AsXAlqt7Ohap317UZoZOhR4/HHxd//+t/j+du4EXn5ZlI3GjxelpAkTLFdG\nIqKGmdVyt8iNXLzlXl4OrFwJvPsu0LMn8MQTwOTJQO/esiOzPkUBTp0CduwAtm4Fjh0Trfl77gEm\nThR9CUTUMKuWZSzBVZN7Tg6wZAnw8cdARATwzDOifu7K9eiLF4HPPxedtDt3ijr9lCnAvfeKUTpE\ndB2Tu53JyADmzxcdjwkJ4ujXT3ZU9ufqVeCrr4BNm4Bt24DbbgMeeki8q2HphojJ3W7k5wMvvig6\nFZOSRKK64QbZUTmGK1dE5/Lq1aJz9p57xPPT68UkLCJX1Nrkzv8yFlJeDvzlL6KD1N9fjAl/7DEm\n9pbo3Bl48EExguj0aTFi6LnnxDj+114Ts3CJyDxM7m1UUQG8846YvVleLoYEvvKKSFTUen37As8/\nD3z9NZCWJiZUBQSIuvyuXaKTlogax7JMG2RmAg8/LMZwv/mmmLhD1vPjj8CaNcB774mZtI8/Djzy\niBh9ROSsWHO3oaoqYNEiMQomJUW0Jsl2FAXYv18keYNB1OXnzuVEKXJOrLnbSH6+6OAzGsV4bSZ2\n21OpgFGjgHXrxNIHHTqIvo7oaODoUdnREdkHJnczKQqwapVYP2XqVDF8z1FnlDoTjQZITha/dEND\nxb/NmDFiOYeaGtnREcnDsowZSkvFjNJTp4C1a4HgYNkRUWMqK4FPPwXeeguorhad21OnciglOS6W\nZawkN1e85f/Tn8SiWEzs9k2tBu6/Hzh+XPSLvPmm+Df79FO25Mm1sOXehAMHgGnTgDfeAGJjZUdD\nraEootM1MVFMknrlFWD6dLbkyXFwtIyFbdkCzJolZktGNLmDLDmC2iSflCS2HUxKEuUaV17jhxwD\nk7sFvfsu8PrrYq2TW2+VHQ1ZUm2Sf/FFMXv4zTeB0aNlR0XUOCZ3C6ipEf/pt2wBtm8X29mRc6qp\nESt1vvSS2DLwjTc4CY3sEztU26iiQsw23btX1NqZ2J1bu3ai4/XMGbF5yPjxol/l/HnZkRFZBpM7\nxIzTe+8VtdgdOwAPD9kRka3ccAMwZ45Y6E2rBW65Rbx7++kn2ZERtY3LJ3dFEWPYr10Tw+W4hrhr\n6t5drDx58qRYpMzPT3Smc/gkOSqXr7kvWCD2NN29G+jSRXY0ZC8OHRK7Zrm7i1U/hw2THRG5KqvV\n3GfOnLnC09OzePDgwScben3t2rUPBAcHZwUFBWWPHDnyQHZ2dlBLg5Bl6VKxPskXXzCx0++FhYnd\ntGbNErtCxcYCP/wgOyoi8zWb3GNjY1caDIZGR3r379//u717947Ozs4Omj9//muzZs1637IhWseW\nLcCrr4phcTfeKDsaskft2omkfuYM0KeP2Dzk7bdFHw2RvTOrLFNQUOA1adKkbSdPnmxysFhZWVnP\nwYMHnywsLPzDklr2VJbZv19MYNm+nePYyXxnzoj+mUuXgNRUsSwFkbW1tizjbskgli9fHhcZGZne\n2OuJiYm/fazX66HX6y15e7OcOiWWFFi7lomdWsbPT+wCtXo1MHGiGGH1P/8jOmOJLMVoNMJoNLb9\nQoqiNHvk5+d7BQYGnmzqnF27do3x9/c/XVpa2rOh18Wt5DKZFEWnU5Q1a2RHQo6upERRHntMUTQa\nRfnkE0WpqZEdETmrX3OnWbm67mGRoZDZ2dlB8fHxy7Zu3RrVs2fPMktc09Kqq8WklYQE4IEHZEdD\njq53b2DZMjHLNSkJuPtuToAi+9Lm5H7+/Pmbpk6dunnNmjUP+vj4nLNEUNaQnCyGtf3tb7IjIWdy\n++1ieeGRI0WZLzWVm3eTfWi2QzUmJmb9nj177igpKfHw9PQsTkpKWlBZWakGgISEhNTHHnvsgy1b\ntky56aabzgOAWq2uzMzMHP6HG0nsUM3MFDXSY8cAnU5KCOQCTp0CZs4Uw2qXLeMSFmQZXDisET/9\nJKaUL1zI/U7J+qqqxHDJ5GSxdvxTT3HteGobJvdGxMeLrdc+/NDmtyYXdvasaMWrVMDy5cCAAbIj\nIkfFVSEbsHmzGLr2zjuyIyFXM2CAWGH03ntFPf7dd1mLJ9ty2pZ7UREwdCjw2WdiKjmRLGfPAg8+\nCPTqBaxYAfTrJzsiciRsuddRUwM88gjw5JNM7CTfgAFij4ARI0T/z8aNsiMiV+CULff/+z/xH2jv\nXjH8kcheZGQADz0kEv2SJZzdSs1jy/1XOTliy7Q1a5jYyf6EhgInTojhksHBwJ49siMiZ+V0LfdJ\nkwC9Hnj+eavfiqhN0tOBuDgxa/rll9kYoYZxKCQAo/H6Eq033GDVWxFZxH/+I8o0FRViMTtOsqP6\nXL4sU1MD/OUvYrISEzs5ir59ga++Au66SywhvHWr7IjIWThNcv/kEzGOeMYM2ZEQtUy7dmLNoy1b\nxNZ+zzwD/PKL7KjI0TlFcr92TexY/9ZbnOpNjuu220Rn6/ffiyG8Z8/KjogcmVOkwvfeAwICgDFj\nZEdC1DY9ewKffir2bh05kmPiqfUcvkO1vFxMEtm1S+xxSeQsjh4VyxdMmSIWIlOrZUdEMrhsh+ob\nbwBRUUzs5HxCQsQy1bm54l1pUZHsiMiROHRyP38e+OADsRMOkTPq1QvYtk2Mphk2TLxDJTKHQ5dl\nHn4Y+POfgddes+hliezSzp1iAbKnnwZeeIGDB1yFy01iOnFCtGbOngW6dbPYZYnsWlERcN99QJ8+\nwKpV/Nl3BS5Xc583D5g/nz/c5Fo0GmD3buBPfxKLj+XlyY6I7JVDJvejR0WLfdYs2ZEQ2V779sDS\npcCcOWKD7u3bZUdE9qjZ5D5z5swVnp6exYMHDz7Z2DnPPPPMO76+vnnBwcFZJ06cuMWyIf5RaqpI\n7BwaRq4sIQHYtEksPrZoEXd6ot9rNrnHxsauNBgMEY29np6eHnnu3DmfvLw83/fff3/WE0888Z5l\nQ/y9H38UEztmzrTmXYgcw+23A5mZYkvJ6GjgyhXZEZG9aDa5jxo1al/Pnj3LGnt969atUY888shH\nABAaGppRXl7eo7i42NOSQda1di0wfryoORIRoNWKjWk6dhRLGBQUyI6I7EGbV5AuKirS6HQ6U+3n\nWq22sLCwUOvp6Vlc/9zExMTfPtbr9dDr9S26l6KIkszf/976eImcUYcOwMqVwNtviwS/aRO3mHRU\nRqMRRqOxzdexyPYA9YfpqFSqBqt/dZN7a2RkiLedY8e26TJETkmlAp59VizHERUFLF4M3H+/7Kio\npeo3fJNaOUuzzcldo9EUmUym37YYKCws1Go0GqtMlK7tSOXkDaLG3X23mMk6aZLYuCYxkf9nXFGb\n/8mjoqK2rlq16mEAOHz48IgePXqUN1SSaavycuCzz8ROS0TUtMGDxTvdf/0LiIkBfv5ZdkRka822\n3GNiYtbv2bPnjpKSEg+dTmdKSkpaUFlZqQaAhISE1MjIyPT09PRIHx+fc507d76ycuVKq6Tf1auB\niAgxM4+ImufpKSY8xcWJfYU/+0zs/ESuwSGWH1AUserju++KH1IiMp+iAK+/DixbBnz+uWjVk+Nw\n6uUHDhwAqquBO+6QHQmR41GpgJdfFhOdxo3jypKuwiGSe21HqqrFv7uIqFZMjNjlKSZGlDnJudl9\nWebiRcDHBzh3Dujd2wqBEbmYU6fEiJr4eLH3MBtN9s1pyzKrVgETJzKxE1nKoEHAwYNiolNCAlBV\nJTsisga7brkrCuDvL3Zbuv12KwVG5KIuXxZrw6tUwCefAF26yI6IGuKULfc9ewA3N7ELPBFZVteu\nwNatQL9+YrDCf/8rOyKyJLtO7qmp4m0ja4JE1qFWiyGSEyeKd8f5+bIjIkux27LMtWtiwlJBgdgk\nmIisKyVFDJdMTweCgmRHQ7VaW5axyMJh1nDoEBAQwMROZCuzZ4sG1YQJYs+EUaNkR0RtYbdlmR07\nxIQLIrKdGTOANWuAadOAbdtkR0NtYdfJffx42VEQuZ4JE4AvvhATB1eulB0NtZZd1tx//FHsLnPh\ngtiEgIhsLzcXCA8HnnoK+MtfZEfjupyq5m40il1kmNiJ5Bk4UKzrNH68GBOflMSRa47ELssyO3ey\nJENkDzQaMd9k2zbguefExEJyDHaZ3NmZSmQ/brxRrAt/+LCow1dXy46IzGF3yb2oCCguBoYMkR0J\nEdXq0UPs6vTdd8CDDwKVlbIjoubYXXLftUtsgO3mJjsSIqqrSxcxiuann4Dp04FffpEdETXF7pI7\nSzJE9qtDB2DzZvHnxIki0ZN9sqvkrijsTCWyd2o1sG4dcNNNwF13iZE0ZH+aTe4GgyHCz8/vjK+v\nb15ycvK8+q+XlJR4REREGIYMGfJ1YGDgNx9++OGjrQ0mN1eUY7y9W3sFIrIFNzexFHdAgNi4/tIl\n2RFRfU1OYqqurnYbOHBg7o4dO8ZrNJqiYcOGHVm/fn2Mv79/Tu05iYmJideuXbvhjTfe+FtJSYnH\nwIEDc4uLiz3d3d1/twWAOZOYUlKAEyeA5cvb+m0RkS3U1Ig1aY4fB778EujeXXZEzscq67lnZmYO\n9/HxOefl5VWgVqsro6OjP05LS5tc95y+ffv+59KlS90A4NKlS9169+59sX5iNxdLMkSOpV074N13\ngWHDxLIF5eWyI6JaTc5QLSoq0uh0OlPt51qttjAjIyO07jnx8fHLxo4du6tfv37fX758uesnn3xy\nX2PXS0xM/O1jvV4PvV7/2+dVVWJm6tKlLf4eiEgilQp45x3g2WdF4+yrr7iaa1sYjUYYjcY2X6fJ\n5K5SqZqdj7Zw4cIXhwwZ8rXRaNR/++233hMmTPhXVlZWcNeuXf/QzVI3udd3/Dig0wGenuaETUT2\nRKUC/vEPsQbNuHFi1Bv3PW6d+g3fpKSkVl2nybKMRqMpMplMutrPTSaTTqvVFtY95+DBg7fde++9\nnwKAt7f3tzfffHN+bm7uwJYGwiGQRI5NpQLeegu4804xV+XCBdkRubYmk3tISMjRvLw834KCAq+K\nior2GzZsmBEVFbW17jl+fn5nduzYMR4AiouLPXNzcwf279//u5YGwiV+iRyfSiV2c4qMFDX40lLZ\nEbmuJssy7u7uVSkpKbPDw8O/rK6udouLi1vu7++fk5qamgAACQkJqS+++OLC2NjYlcHBwVk1NTXt\n3nzzzb/26tWrRf+kV68CR44Ao0e35VshInugUgELF4olCu68UzTcevSQHZXrsYv13P/1L7Gc6P79\nNgmFiGxAUYC5c4GMDNHJ2q2b7Igck1WGQtoKh0ASOR+VCnj7bWDoUFGm4VIFtmUXyZ2dqUTOSaUS\nkxP9/MRaNFeuyI7IdUgvy5SWAl5eQEkJ0L69TUIhIhurqQFiY8WS3tu2AR07yo7IcThsWWb3bmDU\nKCZ2ImfWrh2wYoWYxzJlCnDtmuyInJ/05M6SDJFrcHMDPvpIrAsfEyNmpZP1SE/uOTncdYnIVbi7\nA2vXiuHPcXGiXEPWIT25l5ZymjKRK7nhBrHhx3ffAXPmcNNta5Ge3MvKgJ49ZUdBRLbUqRPw+efA\nwYPAyy/LjsY5NTlD1RbKyriCHJEr6t4dMBjEzPTu3YG//lV2RM5FanKvqBBH584yoyAiWfr0EYMq\nRo0SM1gff1x2RM5DanIvKxNrTqhaPIKTiJyFRiOWILnjDqBrV+CBB2RH5BykJ3eWZIjI21ts0zdu\nnOiDi4yUHZHjk9qhWlrKzlQiEgYNArZsAR59FDh8WHY0jk9qcudIGSKqKywM+PBD4J57xBwYaj3p\nyZ1lGSKqKzJS7OgUEQGYTM2fTw2TWnNnWYaIGvLQQ8B//ysS/L59bAS2hvSWO5M7ETXk+edFK37i\nRLFcAbWM9A5V/kYmosYkJwO+vsB994lt+8h8bLkTkd1q1w744AOx/kx8PNehaYlmk7vBYIjw8/M7\n4+vrm5ecnDyvoXOMRqP+lltuOREYGPiNXq83mntzdqgSUXPUauCTT8TomcRE2dE4EEVRGj2qqqrc\nvL29z+Xn53tVVFSog4ODvz59+rR/3XPKysp6BAQEnDKZTFpFUXDhwgWPhq4lbvV7t92mKPv2/eGv\niYj+oLhYUfr3V5Tly2VHYlu/5s4mc3VDR5Mt98zMzOE+Pj7nvLy8CtRqdWV0dPTHaWlpk+ues27d\nuvunTZu2SavVFgKAh4dHibm/WFiWISJz3XgjkJ4OvPiimM1KTWtyKGRRUZFGp9P9NtJUq9UWZmRk\nhNY9Jy8vz7eyslI9ZsyY3ZcvX+46Z86cxQ899NDqhq6XWOc9lV6vR1mZnmUZIjLbwIHAxo3A1KnA\nV18550Y/RqMRRqOxzddpMrmrVKpmuy8qKyvVx48fH7pz585xV69e7RQWFnZoxIgRh319ffPqn1s3\nuSsKx7kTUcvdfjvw7rvApEliPXidTnZElqXX66HX63/7PCkpqVXXaTK5azSaIpPJ9NujM5lMutry\nSy2dTmfy8PAo6dix488dO3b8efTo0XuzsrKCG0rudf38s+gJ79ChVXETkQu7914xezUyEti/X6wH\nT7/XZM09JCTkaF5enm9BQYFXRUVF+w0bNsyIioraWvecyZMnp+3fv//26upqt6tXr3bKyMgIDQgI\nON3cjTlShoja4tlngTFjgGnTxL4Q9HtNJnd3d/eqlJSU2eHh4V8GBAScnjFjxgZ/f/+c1NTUhNTU\n1AQA8PPzOxMREWEICgrKDg0NzYiPj19mTnJnSYaI2kKlAv7xD7EG/KxZHANfn0qx0RNRqVRK3Xvt\n3Qu89JJYN4KIqLWuXBEbfUyfDrzwguxoLE+lUkFRlBZvaSRt4TAuPUBEltC5M5CWBowYAQwYIEbS\nkMTkzjHuRGQpGo1I8OHhwJ//DNx6q+yI5JO2tgyTOxFZ0tChwPvvA5MnA0VFsqORj2UZInIaU6YA\nZ8+KMfD79omSjatiy52InMpf/woEBQEPPgjU1MiORh6pyZ0tdyKyNJUKSE0FLl4UI/JclbTkznHu\nRGQtN9wAbN4slgr+6CPZ0cjB0TJE5JQ8PIBt28QYeD8/IDS0+a9xJlJb7izLEJE1BQQAy5eLJQq+\n/152NLbFDlUicmpRUcDjj4sEf+2a7GhsR8ryAzU1QPv2YmVItdomtyciF6YoYiXJbt1ES17V4sn8\n8rR2+QEpLffLl4FOnZjYicg2VCrgww+Bo0eBlBTZ0diGlA5VlmSIyNa6dBFLFISFAYMGAWPHyo7I\nuqS03JnciUiGm28G1q0D7r8fyM+XHY11SUnuHClDRLKMHSs22b7nHrFcsLNiy52IXM7TT4uFxuLj\nnXeTD2nJnS13IpJFpQL++U/gzBnn7WCV0qHKpQeISLaOHYGNG0UH69ChwMiRsiOyLJZliMhl9e8P\nrFgBzJgBFBfLjsaymk3uBoMhws/P74yvr29ecnLyvMbOO3LkyDB3d/eqzZs3N7vJFTtUiche3H03\nMHOmSPBVVbKjsZwmk3t1dbXb7NmzUwwGQ8Tp06cD1q9fH5OTk+Pf0Hnz5s1LjoiIMJgzk4otdyKy\nJwsWAB06AH/7m+xILKfJ5J6ZmTncx8fnnJeXV4Fara6Mjo7+OC0tbXL985YsWfL09OnTN/bp0+eC\nOTdlcicie+LmBqxdK2rwmzbJjsYymuxQLSoq0uh0OlPt51qttjAjIyO0/jlpaWmTd+3aNfbIkSPD\nVCpVowOLEhMTAQCnTgH5+XoA+jaETkRkOb17A59+CkRGihmsfn5y4jAajTAajW2+TpPJvalEXWvu\n3LlvL1q06IVfFwZTNVWWqU3uH30EjB/f0lCJiKwrJAR4/XWxgmRGhliywNb0ej30ev1vnyclJbXq\nOk0md41GU2QymXS1n5tMJp1Wqy2se86xY8dujY6O/hgASkpKPLZv336XWq2ujIqK2trYdVmWISJ7\n9dhjwIEDwOzZYrExR9Xkkr9VVVXuAwcOzN25c+e4fv36fT98+PDM9evXx/j7++c0dH5sbOzKSZMm\nbZs6dermP9zo1yV/q6vFFlgVFUA7aavJExE17soVYNgwYN484JFH5MbS2iV/m2y5u7u7V6WkpMwO\nDw//srq62i0uLm65v79/TmpqagIAJCQkpLb0huXlYk1lJnYisledO4v9V8eMEUk+IEB2RC1n8806\nzp0DwsOBb7+1yW2JiFrtgw+AxYtF/b1TJzkxOMxmHZzARESOIi4OCAoC5syRHUnL2Ty5szOViByF\nSgUsXQrs2SPWgXckbLkTETWha1dRf58zBzh7VnY05mPLnYioGUOGAK++Ctx3H/DLL7KjMQ+TOxGR\nGR5/HBgwAHjuOdmRmIdlGSIiM6hUwLJlgMEgNtq2d2y5ExGZqXt3YPVqICEB+OEH2dE0jcmdiKgF\nRo4EZs0CHn0UqKmRHU3jWJYhImqh+fNFQ9We91+1+R6qbLkTkaNTq8X672FhwNixQGCg7Ij+iGUZ\nIqJW8PEBFi0CHngAuHZNdjR/ZPO1ZTp3FhvRylgnmYjIkhQFmD4d8PIC/v5369yjtWvL2DS5//KL\ngq5dxW85VYtDJSKyPxcvAsHBYu13a2xC5BALh9WWZJjYichZ9O4NrFwpRs9cvCg7muukJHciImcy\nYYJYmiAhQZRq7AGTOxGRBSxcCJw+LRYZswc2Te4c405EzqpDB1F3nzNHDBqRjS13IiILGT5c1N6f\nfFJ+ecas5G4wGCL8/PzO+Pr65iUnJ8+r//ratWsfCA4OzgoKCsoeOXLkgezs7KCGrsPkTkTOLjER\nyMmRX55pNrlXV1e7zZ49O8VgMEScPn06YP369TE5OTn+dc/p37//d3v37h2dnZ0dNH/+/NdmzZr1\nfkPXYlmGiJydvZRnmk3umZmZw318fM55eXkVqNXqyujo6I/T0tIm1z0nLCzsUPfu3X8EgNDQ0IzC\nwkJtQ9diy52IXIE9lGeaTe5FRUUanU5nqv1cq9UWFhUVaRo7f/ny5XGRkZHpDb1WWsrkTkSuQXZ5\nptmFw1Qqldm/d3bv3j1mxYoVMw8cODCyodczMhJRXQ3k5wN6vR56vb4FoRIROY7a8kxUFKDXA56e\n5n2d0WiE0Whs8/2bTe4ajabIZDLpaj83mUw6rVZbWP+87OzsoPj4+GUGgyGiZ8+eZQ1dq0+fRDz5\nJHD77W0LmojIEdQtz2zcaN7s/PoN36SkpFbdu9myTEhIyNG8vDzfgoICr4qKivYbNmyYERUVtbXu\nOefPn78gFWnuAAAIKklEQVRp6tSpm9esWfOgj4/PucauxbIMEbkaWeWZZlvu7u7uVSkpKbPDw8O/\nrK6udouLi1vu7++fk5qamgAACQkJqa+++uorZWVlPZ944on3AECtVldmZmYOr3+tsjKOliEi19Kh\nA7BiBTBlChAeDvToYZv72nRVyPbtFfz4o/hmiYhcSUKC2OSjpbs3OcSSvx07Krh61Sa3IyKyK6Wl\nQEAA8PnnQEiI+V/nEEv+st5ORK6qVy8gORl4/HGgutr692NyJyKykYcfBjp3BpYutf69bJrc2ZlK\nRK5MpQL++U8xguaHH6x7L7bciYhsaNAgIC4OeP55696HyZ2IyMbmzwcOHAB27rTePViWISKysc6d\ngcWLgaeeAq5ds8492HInIpJg8mRgwADgf//XOtdncicikuSdd4B//AP47jvLX5tlGSIiSby8xKYe\nr7xi+Wuz5U5EJNGcOcBXXwFnzlj2umy5ExFJ1K0bMHcu8Nprlr2uTdeWyc1VMGCATW5HROQwLl8G\nvL2BvXsBP7/fv8a1ZYiIHFTXrsCzzwKvvmq5a9q05V5ZqcC92RXkiYhcT23rfc8ewN//+t87RMud\niZ2IqGG1rXdL1d5t2nK31b2IiBxRQ613h2i5ExFR47p2BZ57zjK1d7bciYjsSG3r3WgUOzdZreVu\nMBgi/Pz8zvj6+uYlJyfPa+icZ5555h1fX9+84ODgrBMnTtzS0iBcjdFolB2C3eCzuI7P4jpXfha1\nrfe21t6bTO7V1dVus2fPTjEYDBGnT58OWL9+fUxOTo5/3XPS09Mjz50755OXl+f7/vvvz3riiSfe\na1tIzs+Vf3Dr47O4js/iOld/FrNnA7t2AadOtf4aTSb3zMzM4T4+Pue8vLwK1Gp1ZXR09MdpaWmT\n656zdevWqEceeeQjAAgNDc0oLy/vUVxc7Nn6kIiIXFuXLm1vvTeZ3IuKijQ6nc5U+7lWqy0sKirS\nNHdOYWGhtvUhERHRU08Bu3e3/uubHHmuUqnM6gGtX+xv7OtUqhb3CTitpKQk2SHYDT6L6/gsruOz\naJsmk7tGoykymUy62s9NJpNOq9UWNnVOYWGhVqPRFNW/Vmt6e4mIqHWaLMuEhIQczcvL8y0oKPCq\nqKhov2HDhhlRUVFb654TFRW1ddWqVQ8DwOHDh0f06NGj3NPTs9iaQRMRUdOabLm7u7tXpaSkzA4P\nD/+yurraLS4ubrm/v39OampqAgAkJCSkRkZGpqenp0f6+Pic69y585WVK1fG2iZ0IiJqlKIoFj22\nb98eMXDgwDM+Pj55ixYtmtfQOU8//fQ7Pj4+eUFBQVnHjx+/xdIx2MvR3LNYs2bNA0FBQVmDBw/O\nvu222w5kZWUFyY5Z1rOoPTIzM4e5ublVbdq0aarsmGU9h927d+uHDBlyYtCgQd/ccccdRtkxy3oW\nFy5c8AgPDzcEBwd/PWjQoG9Wrlz5qOyYrXXExsauuPHGG4sDAwNPNnZOS/OmRQOsqqpy8/b2Ppef\nn+9VUVGhDg4O/vr06dP+dc/54osvIu+66650RVFw+PDh0NDQ0MOyH6w1DnOexcGDB8PKy8u7K4r4\nQXflZ1F73pgxY3bdfffdn2/cuHGa7LhlPIeysrIeAQEBp0wmk1ZRRIKTHbesZ7FgwYLEF1544Y3a\n59CrV6+LlZWV7rJjt8axd+/eUcePH7+lseTemrxp0bVlOC7+OnOeRVhY2KHu3bv/CIhn4axDSM15\nFgCwZMmSp6dPn76xT58+F2TEaW3mPId169bdP23atE21Axc8PDxK5ERrXeY8i759+/7n0qVL3QDg\n0qVL3Xr37n3R3d29Sk7E1jVq1Kh9PXv2LGvs9dbkTYsmd46Lv86cZ1HX8uXL4yIjI9NtE51tmftz\nkZaWNrl2hrO5w3AdiTnPIS8vz7e0tLTXmDFjdoeEhBxdvXr1Q7aP1PrMeRbx8fHLTp06Nahfv37f\nBwcHZy1evHiO7SO1D63JmxZdYd3S4+IdWUu+p927d49ZsWLFzAMHDoy0ZkyymPMs5s6d+/aiRYte\n+HWBOVX9nxFnYM5zqKysVB8/fnzozp07x129erVTWFjYoREjRhz29fXNs0WMtmLOs1i4cOGLQ4YM\n+dpoNOq//fZb7wkTJvwrKysruGvXrpdtEaO9aWnetGhyt+S4eEdnzrMAgOzs7KD4+PhlBoMhoqm3\nZY7MnGdx7NixW6Ojoz8GgJKSEo/t27ffpVarK+sPvXVk5jwHnU5n8vDwKOnYsePPHTt2/Hn06NF7\ns7Kygp0tuZvzLA4ePHjbSy+99DoAeHt7f3vzzTfn5+bmDgwJCTlq63hla1XetGSnQGVlpXv//v2/\nzc/P97p27Vr75jpUDx06NMJZOxHNeRb//ve/b/L29j536NChEbLjlf0s6h6PPvroSmccLWPOc8jJ\nyfEbN27cjqqqKrcrV650CgwMPHnq1KkA2bHLeBbPPvvs/yUmJi5QFAU//PCDp0ajKbx48WIv2bFb\n68jPz/cyp0PV3Lxp8QDT09PvGjBgQK63t/e5hQsX/k1RFCxdujRh6dKlCbXnPPXUUyne3t7ngoKC\nso4dOzZU9kO11tHcs4iLi/ugV69eF4cMGXJiyJAhJ4YNG5YpO2ZZz6Lu4azJ3dzn8NZbb/2/gICA\nU4GBgScXL178jOyYZT2LCxcueEycOHFbUFBQVmBg4Mm1a9feLztmax3R0dHr+/bt+71ara7QarWm\n5cuXz2xr3rTZZh1ERGQ73GaPiMgJMbkTETkhJnciIifE5E5E5ISY3ImInBCTOxGRE/r/EOWhpv26\nwocAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x38dd150>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the pdf to see how it looks ... \n",
    "x = np.linspace(0,1,50)\n",
    "plt.plot(x, betadist.pdf(x))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}