Random Variables.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "from scipy.stats import rv_continuous\n",
    "import random\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Triangle_rv(rv_continuous):\n",
    "    def _pdf(self, x):\n",
    "        if x <=0.0 or x >= 1.0:\n",
    "            return 0.0\n",
    "        elif x <= 0.5:\n",
    "            return 4.0*x\n",
    "        else: # 0.5 < x < 1.0\n",
    "            return 4 - 4.0*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "triangle = Triangle_rv(a=0.0,b=1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "triangle.cdf(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "triangle.pdf(0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check if scipy will catch an invalid pdf description\n",
    "class Unbounded_rv(rv_continuous):\n",
    "    def _pdf(self,x):\n",
    "        if x <= 0.5:\n",
    "            return 4.0*x\n",
    "        else:\n",
    "            return 4 - 4.0*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "bad_triangle = Unbounded_rv()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/patrick/anaconda3/lib/python3.6/site-packages/scipy/integrate/quadpack.py:385: IntegrationWarning: The algorithm does not converge.  Roundoff error is detected\n",
      "  in the extrapolation table.  It is assumed that the requested tolerance\n",
      "  cannot be achieved, and that the returned result (if full_output = 1) is \n",
      "  the best which can be obtained.\n",
      "  warnings.warn(msg, IntegrationWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-1.679743221855993e+20"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bad_triangle.cdf(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  This is unfortunate.  It appears that scipy does not check\n",
    "# the validity of continous random variable object pdfs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bad_triangle.pdf(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-8.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bad_triangle.pdf(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Horrible!!  It doesn't even catch negative values!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try creating a variable triangle \n",
    "class Concentrated_rv(rv_continuous):\n",
    "    \n",
    "    def __init__(self, a, b, m=0.5, seed=12345):\n",
    "        assert(a < m < b)\n",
    "        super().__init__(a=a, b=b)\n",
    "        self.m = m # midpoint (or other point between a and b)\n",
    "        \n",
    "        self.width = self.b - self.a\n",
    "        self.height = 2.0 / self.width\n",
    "        random.seed(seed)\n",
    "        \n",
    "    def _pdf(self,x):\n",
    "        if x <= self.a:\n",
    "            return 0.0\n",
    "        elif x >= self.b:\n",
    "            return 0.0\n",
    "        elif x <= self.m:\n",
    "            return (self.height / (self.m - self.a)) * (x - self.a)\n",
    "        else: # x > self.m\n",
    "            return self.height - (self.height / (self.b - self.m)) * (x - self.m)\n",
    "    \n",
    "    # vectorized pdf and cdf functions; these will compute the pdf, cdf\n",
    "    # function values for each x value in an array\n",
    "    #vpdf = numpy.vectorize(self.pdf)\n",
    "    #numpy.vectorize(self.cdf)\n",
    "    \n",
    "    def plotpdf(self):\n",
    "        x = numpy.linspace(0.0, 1.0, 100)\n",
    "        #y = self.pdf(x)\n",
    "        #y = numpy.apply_along_axis(self.pdf, 0, x)\n",
    "        vpdf = numpy.vectorize(self.pdf)\n",
    "        y = vpdf(x)\n",
    "        plt.plot(x,y)\n",
    "    \n",
    "    def plotcdf(self):\n",
    "        x = numpy.linspace(0.0, 1.0, 100)\n",
    "        #y = numpy.apply_along_axis(self.cdf, 0, x)\n",
    "        #y = self.vcdf(x)\n",
    "        vcdf = numpy.vectorize(self.cdf)\n",
    "        y = vcdf(x)\n",
    "        plt.plot(x,y)\n",
    "    \n",
    "    def sample(self):\n",
    "        return self.cdf(random.uniform(0,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "RV = Concentrated_rv(0.0,0.2,0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "0.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(RV._pdf(0.5))\n",
    "print(RV.pdf(0.5))\n",
    "RV.cdf(0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGY1JREFUeJzt3WtsZOdZB/D/MzO+zRnba3tOslevfdy06SYSNBhIKSrQtNCGquFDQakUCFJEaIFSLhJq1Q+tQAKEoFykCli1paWUthAqGlXlUnpRoGoCThPa7G5L43H25t3sGdtr78zYHs/Mw4eZ43Um9no8537m/5NWO7aP5zxn7f379XueeV9RVRARUfylwi6AiIi8wUAnIkoIBjoRUUIw0ImIEoKBTkSUEAx0IqKEYKATESUEA52IKCEY6ERECZEJ8mT5fF6npqaCPCURUew9/fTTRVU19zsu0ECfmprC3NxckKckIoo9ETnfyXGcciEiSggGOhFRQjDQiYgSYt9AF5GPicg1EXlux/vGReRLIvK91t9j/pZJRET76WSE/nEAb25733sBfFlV7wDw5dbbREQUon0DXVWfALDc9u4HAHyi9fgTAH7G47qIiOiAup1Dv11VrwBA6+/b9jpQRB4VkTkRmbNtu8vTERHRfny/Kaqqp1V1VlVnTXPfvvjI+dLZF3FppRJ2GURE++o20F8UkSMA0Pr7mnclRcd6tY5f/uQcTj9RCLsUIqJ9dRvojwN4uPX4YQCf96acaDl3dQ0NBebtUtilEBHtq5O2xU8D+AaAV4nIJRF5BMAfAniTiHwPwJtabyfOmcU1AMD8tXLIlRAR7W/ftVxU9R17fOg+j2uJnLOLqwCAq2sbKG/WYAwEuvQNEdGB8JWit3B2cQ3plAAAFoocpRNRtDHQ91CrN/CdqzfwulfkAXAenYiij4G+h3m7jM1aA/fffRgiQMHmCJ2Ioo2Bvoczrfnze06O4dihIRQ45UJEEcdA38OZxTUMZFKw8gYsM4cCp1yIKOIY6Hs4u7iGO4+MIJNuhvpCsQxVDbssIqI9MdB3oao4s7iKU0dGAAAzpoFKtY6raxshV0ZEtDcG+i4uraxjbaOGu442A90ycwB4Y5SIoo2BvgvnFaI3A90AAM6jE1GkMdB3cXZxFSkB7jzcDPTDI4PI9qcxzxE6EUUYA30XZ6+swTJzGOpPAwBEBJZpsHWRiCKNgb6LM4tr29MtDivP1kUiijYGepvVyhaurG7g1UfaAt00cPn6Oja26iFVRkR0awz0Ns+3RuGvaHW2OCwzB1XghSVOuxBRNDHQ2zjTKk5ni8PKO50uDHQiiiYGeptCsYxMSnBiPPuS97N1kYiijoHepmCXMDmRRV/6pf802f4MjowOcoRORJHFQG9TsMuw8rldP2aZBubZukhEEcVA36HeUJxfqmCmbf7c4bQucpEuIooiBvoOl1YqqNYbL7sh6rBMAzc2aiiWqgFXRkS0Pwb6Ds78uGXuNeXiLNLFG6NEFD0M9B2cfUOdFsV2262LnEcnoghioO9QKJYxOtSHcaN/148fOzSEgUyKI3QiiiQG+g4FuwTLNCAiu348lRJM5w22LhJRJDHQd1go7t2y6OCqi0QUVQz0ltJmDS+ube7Z4eKw8jlcWK6gWmsEVBkRUWcY6C0LrWmUvXrQHZZpoN5QXFiuBFEWEVHHGOgthaKzKNd+Uy5sXSSiaGKgt8zbZaQEODmRveVx24t0cR6diCKGgd5SsEs4PpbFQCZ9y+NGBvuQzw1whE5EkcNAbynY5X3nzx2WydZFIooeV4EuIr8pImdE5DkR+bSIDHpVWJAaDW22LO4zf+6YYesiEUVQ14EuIscA/DqAWVW9G0AawINeFRakq2sbWN+q79uy6LDyOSyXq7he4SJdRBQdbqdcMgCGRCQDIAtg0X1JwdtelGufFxU5nOCf57QLEUVI14GuqpcB/DGACwCuAFhV1X9vP05EHhWRORGZs227+0p95LQsdj6HztZFIooeN1MuYwAeADAN4CgAQ0Qeaj9OVU+r6qyqzpqm2X2lPirYZeQGMjCHBzo6/sTYEPrSwnl0IooUN1MubwSwoKq2qm4B+ByAH/GmrGDN77MoV7tMOoXJ8SxH6EQUKW4C/QKAe0UkK80kvA/AOW/KClZzH9HOplsclplj6yIRRYqbOfSnADwG4JsAvt16rtMe1RWY9Wodl6+vd9yy6LBMA+eXKqg3uL8oEUVDxs0nq+oHAHzAo1pCsVB0tp072Ah9Jp9Dtd7ApZUKTk4c7HOJiPzQ868U3V6Uq8OWRcf2mi6cdiGiiGCgtwJ5uos5dODmPqRERGFjoNslHDs0hKH+Wy/K1W7c6MehbB9bF4koMhjoxfKB588dVt5g6yIRRUZPB7qqdtWy6GDrIhFFSU8Hun1jE6XN2oFbFh2WaeDajU3c2NjyuDIiooPr6UB3Ftfqfsql+YNggfPoRBQBPR3one4jupcZti4SUYT0dqDbZQz2pXBkpLt9OSYnskinhDdGiSgSejzQS5jO55BKdbYoV7uBTBonxoYwzykXIoqA3g50Fy2LDna6EFFU9Gygb9bquLhcwUyXLYsOK29goVhCg4t0EVHIejbQLyxV0NDub4g6LDOHja0GrqxteFQZEVF3ejbQ3bYsOm4u0sUbo0QUrp4NdKdl8aCLcrXjqotEFBW9G+h2GbcND2B4sM/V85i5AQwPZDhCJ6LQ9XCgl1xPtwCAiMAyDa66SESh691AL5Zd3xB1sHWRiKKgJwN9uVzF9cpW16sstrPyBi5fX8d6te7J8xERdaMnA92Z757xcIQOcJEuIgpXjwa6Ny2Lju1OlyJvjBJReHoy0OeLJfSnUzg+lvXk+abzBkTYukhE4erJQC/YZZxsrZTohcG+NI6ODrF1kYhC1aOB7k3L4k6WaWy/+pSIKAw9F+hb9QbOL1U8a1l0zJg5FOwSVLlIFxGFo+cC/eJyBbWGetay6LBMA+VqHddubHr6vEREneq5QL/Z4eLtCN3ZX3Se8+hEFJLeC/Si04Pu/QgdYKcLEYWn9wLdLmPc6MehbL+nz3t4ZBBDfWkGOhGFpicD3ev5cwBIpQTTeYMvLiKi0LgKdBE5JCKPich3ROSciLzWq8L8Uih637LosEyDI3QiCo3bEfqfA/hXVb0TwPcBOOe+JP+srm+hWKp6fkPUYZk5XFqpYLPGRbqIKHhdB7qIjAB4PYCPAoCqVlX1uleF+cF5JacfUy5A80ZrQ4HzSxVfnp+I6FbcjNAtADaAvxGRZ0TkIyLiT1J6xK+WRYfTusglAIgoDG4CPQPgHgB/qaqvAVAG8N72g0TkURGZE5E527ZdnM69QrGEdEowOe7Nolztpltz81wCgIjC4CbQLwG4pKpPtd5+DM2AfwlVPa2qs6o6a5qmi9O5V7DLmBzPoj/jT3NPbiCD20cGeGOUiELRdbKp6lUAF0XkVa133QfgrCdV+cSvlsWdrHyOrYtEFAq3Q9V3A/iUiHwLwPcD+H33Jfmj3lAsLJV9a1l0OK2LXKSLiIKWcfPJqvosgFmPavHV4vV1VGsN326IOiwzh9X1LSyXq5jIDfh6LiKinXrmlaLzPrcsOm5uR8d5dCIKVs8Eut8ti44Zti4SUUh6J9CLJQwPZpDPebsoV7tjY0Poz6TY6UJEgeudQLfLsMwcRLzZR3Qv6ZRgaiLLXnQiClxPBbrXa6Dvha2LRBSGngj08mYNV9c2MOPz/LnDMg1cWKpgq94I5HxERECPBPpCq+PE7w4Xh2XmUGsoLi5zkS4iCk5PBPp2y2KAI3SA29ERUbB6ItALdhkiwMkJfxblarfdush5dCIKUG8EerGM42NDGOxLB3K+0WwfJox+jtCJKFC9Eeh2aXut8qBwOzoiClriA11VsVD0f1GudmxdJKKgJT7Qr65toFKtB3ZD1GGZBoqlKlbXtwI9LxH1rsQHujPtMRNQy6LD+QHCNV2IKCg9EOjBtiw62LpIREFLfKDP22UY/WncPhLs2uST41lkUsJ5dCIKTA8EegnTpuH7olzt+tIpTI5nOUInosAkPtCb+4gGO93iYOsiEQUp0YG+sVXH4up64C2LDsvMYWGpjHqD+4sSkf8SHegLxTJUg78h6rDyBqq1Bhavr4dyfiLqLYkO9O1t5wJuWXQ4P0jm2bpIRAFIeKA7LYthBTpbF4koOMkO9GIZR0YHke3PhHL+CaMfI4MZti4SUSCSHeh2KbTROQCICCwzxxE6EQUisYGuqqG2LDrYukhEQUlsoNulTdzYrIU6QgeAGTOHq2sbKG3WQq2DiJIvsYG+3eESUsuiw+mwWeAonYh8lvxAD6ll0bG96iJvjBKRzxIc6CUMZFI4dmgo1DpOTmSRkuYiYUREfkpuoBfLmM4bSKWCXZSr3WBfGsfHslwXnYh8l9xAD7llcSd2uhBREFwHuoikReQZEfmCFwV5oVpr4OLKeugtiw4rn8NCsYwGF+kiIh95MUJ/D4BzHjyPZy4sN1c4jNIIfX2rjqtrG2GXQkQJ5irQReQ4gJ8G8BFvyvHGfERaFh1c04WIguB2hP5nAH4HQMODWjzjBOd0yC2Ljhm2LhJRALoOdBF5K4Brqvr0Psc9KiJzIjJn23a3pzuQgl1CPjeA0aG+QM63n9uGB2D0pzlCJyJfuRmhvw7A20TkBQCfAfAGEfm79oNU9bSqzqrqrGmaLk7XuUKxHJn5c+DmIl1cF52I/NR1oKvq+1T1uKpOAXgQwFdU9SHPKnOhYJcwE6FAB9i6SET+S1wf+kq5ipXKVmRaFh1WPofL19exsVUPuxQiSihPAl1Vv6aqb/XiudxybjxGacoFuFnPQpGjdCLyR+JG6FFrWXSwdZGI/Ja4QC/YZfSlBSfGwl2Uq53TQsk1XYjILwkM9BImx7PIpKN1adn+DI6ODqLAKRci8km0Us8DzZbFaE23ONi6SER+SlSg1+oNnF+KVg/6Tk7roioX6SIi7yUq0C+trGOrrpiJWMuiw8obKG3WYN/YDLsUIkqgRAV6VFsWHc5UEHcvIiI/JCvQI9qy6NhuXeQiXUTkg0QF+rxdxqFsH8aN/rBL2dXR0SEM9qXYi05EvkhUoBfsEqyILJm7m1RKMDVhsBediHyRrECPcMuiY8bMsRediHyRmEC/sbEF+8ZmZG+IOizTwMXlCjZrXKSLiLyVmEDfviEa0ZZFh2UaaChwYakSdilElDDJCfRW50jU1kFv5/zAYesiEXktOYFul5ESYHIiG3Ypt8TWRSLyS6IC/cR4FgOZdNil3NLwYB/M4QG2LhKR5xIT6PMRb1ncycqzdZGIvJeIQG80FC8sRb9l0WGxdZGIfJCIQF9cXcfGViPyLYuOGdPA9coWlsvVsEshogRJRKDHpWXRcXM7Ok67EJF3EhLo8WhZdDg/eHhjlIi8lIxAL5YxPJCBOTwQdikdOT42hL60YJ6ti0TkoWQEut3cpUhEwi6lI5l0CicnDI7QichTCQn0Umw6XBxsXSQir8U+0CvVGhZXN2LTg+6wzBwuLFdQqzfCLoWIEiL2gb5QjPYuRXuxTANbdcXFlfWwSyGihIh9oN/cdi5eI/QZti4SkccSEegiwHTcplzYukhEHot/oBdLrb06o70oV7sxox9j2T6uukhEnol/oLdaFuPIMnNcF52IPBPrQFdVFOwSZmJ2Q9TRbF1koBORN7oOdBE5ISJfFZFzInJGRN7jZWGduHZjE+VqPdYj9GJpE2sbW2GXQkQJ4GaEXgPw26r6agD3AvhVETnlTVmdmW91iMRlUa52Nxfp4iidiNzrOtBV9YqqfrP1+AaAcwCOeVVYJ+Lasuhg6yIRecmTOXQRmQLwGgBPefF8nZq3SxjqS+PwyGCQp/XM5LiBdEo4QiciT7gOdBHJAfgnAL+hqmu7fPxREZkTkTnbtt2e7iUKdhnTeQOpVDwW5WrXn0nhxNgQWxeJyBOuAl1E+tAM80+p6ud2O0ZVT6vqrKrOmqbp5nQvUyiWYjvd4rDMHEfoROQJN10uAuCjAM6p6oe8K6kzG1t1XFpZj90aLu2svIGFYhmNhoZdChHFnJsR+usA/DyAN4jIs60/93tU177OL1WgGp9divZimTls1hq4fJ2LdBGRO5luP1FV/wtAaJPXhZi3LDq2WxeLZZwYz4ZcDRHFWWxfKVpoLZs7HfsROlsXicgbsQ30ebuE20cGkBvo+peMSDBzAxgeyPDGKBG5FttAL9jl2E+3AICIwDINti4SkWuxDHRnUa64tyw6Zti6SEQeiGWgL5WrWNuoxb5l0WGZBq6sbqBSrYVdChHFWCwDPe5ruLRzfjBxlE5EbsQ00JvzzTMJmEMHXtq6SETUrXgGerGM/kwKx8aGwi7FE1MTBkTYukhE7sQz0O0SpiaySMd0Ua52g31pHDs0xCkXInIlpoGejJbFnSwzx9ZFInIldoG+VW/gwnIlMTdEHVbewIJdhioX6SKi7sQu0C8sV1BraGJaFh0zpoFytY4X1zbDLoWIYip2gZ60lkXHzdZFTrsQUXdiGOjJall0OD+g5tm6SERdimGglzFh9GM02xd2KZ46PDKIbH+aI3Qi6lr8Aj0B287tRkQwnTfYukhEXYtfoCewZdHB1kUiciNWgb5a2cJSuZrIETrQbF28tLKOja162KUQUQzFKtDnW6PXpLUsOizTgCrwwhKnXYjo4GIV6EltWXTceXgEAPCdKzdCroSI4ihmgV5CJiWYTOhmyjOmgYFMCmcWV8MuhYhiKGaBXsbkeBZ96ViV3bFMOoU7Dw/jzOJa2KUQUQzFKhmT2rK406mjIzizuMY1XYjowGIT6PWG4oWlSmJviDpOHR3F6voWFlc3wi6FiGImNoF+eWUd1VoDVj7ZI/S7jjZvjJ65zHl0IjqY2AR60lsWHXceHoYIOI9ORAcWm0BPesuiI9ufgZU3GOhEdGAxCvQSRgYzmDD6wy7Fd3cdHcVZti4S0QHFKNDLsMwcRJKxj+it3HV0BIurG1gpV8MuhYhiJD6B3gMti45TrRujZ69w2oWIOheLQC9t1vDi2iZmEn5D1HHX0VEA4CtGiehAXAW6iLxZRL4rIs+LyHu9KqrdgnNDNOEti45xox9HRgd5Y5SIDqTrQBeRNIAPA3gLgFMA3iEip7wqbKdCj7Qs7nTqyAgDnYgOxM0I/YcAPK+qBVWtAvgMgAe8Keul5u0yRICTE8lclGs3dx0dQcEuYb3KtdGJqDMZF597DMDFHW9fAvDD7srZXcEu4fjYEAb70n48fSSdOjqKhgJv+fMnErsYGVEv+ejDP4hJnwelbgJ9t/7Bl60oJSKPAngUACYnJ7s60amjIziR0CVz9/Kjd+Txsz9wHOVqLexSiMgD/Rn/B2bS7ap+IvJaAB9U1Z9qvf0+AFDVP9jrc2ZnZ3Vubq6r8xER9SoReVpVZ/c7zs2PjP8BcIeITItIP4AHATzu4vmIiMiFrqdcVLUmIr8G4N8ApAF8TFXPeFYZEREdiJs5dKjqFwF80aNaiIjIBbZPEBElBAOdiCghGOhERAnBQCciSggGOhFRQnT9wqKuTiZiAzjf5afnARQ9LCcOeM29gdecfG6v96SqmvsdFGiguyEic528UipJeM29gdecfEFdL6dciIgSgoFORJQQcQr002EXEAJec2/gNSdfINcbmzl0IiK6tTiN0ImI6BYiF+j7bTwtIgMi8tnWx58Skangq/RWB9f8WyJyVkS+JSJfFpGTYdTppU43GBeRt4uIikisOyI6uV4R+bnW1/mMiPx90DV6rYPv60kR+aqIPNP63r4/jDq9JCIfE5FrIvLcHh8XEfmL1r/Jt0TkHk8LUNXI/EFzGd55ABaAfgD/C+BU2zG/AuCvWo8fBPDZsOsO4Jp/AkC29fhdvXDNreOGATwB4EkAs2HX7fPX+A4AzwAYa719W9h1B3DNpwG8q/X4FIAXwq7bg+t+PYB7ADy3x8fvB/AvaO74di+Ap7w8f9RG6J1sPP0AgE+0Hj8G4D4R2W07vLjY95pV9auqWmm9+SSA4wHX6LVONxj/PQB/BGAjyOJ80Mn1/hKAD6vqCgCo6rWAa/RaJ9esAEZaj0cBLAZYny9U9QkAy7c45AEAf6tNTwI4JCJHvDp/1AJ9t42nj+11jKrWAKwCmAikOn90cs07PYLmT/g42/eaReQ1AE6o6heCLMwnnXyNXwnglSLydRF5UkTeHFh1/ujkmj8I4CERuYTmvgrvDqa0UB30//uBuNrgwgedbDzd0ebUMdLx9YjIQwBmAfyYrxX575bXLCIpAH8K4BeDKshnnXyNM2hOu/w4mr+B/aeI3K2q132uzS+dXPM7AHxcVf+ktUfxJ1vX3PC/vND4ml9RG6FfAnBix9vH8fJfw7aPEZEMmr+q3epXnKjr5JohIm8E8H4Ab1PVzYBq88t+1zwM4G4AXxORF9Cca3w8xjdGO/2+/ryqbqnqAoDvohnwcdXJNT8C4B8AQFW/AWAQzTVPkqyj/+/dilqgd7Lx9OMAHm49fjuAr2jrbkNM7XvNremHv0YzzOM+twrsc82quqqqeVWdUtUpNO8bvE1V58Ip17VOvq//Gc2b3xCRPJpTMIVAq/RWJ9d8AcB9ACAir0Yz0O1Aqwze4wB+odXtci+AVVW94tmzh31XeI+7wP+H5h3y97fe97to/ocGml/0fwTwPID/BmCFXXMA1/wfAF4E8Gzrz+Nh1+z3Nbcd+zXEuMulw6+xAPgQgLMAvg3gwbBrDuCaTwH4OpodMM8C+Mmwa/bgmj8N4AqALTRH448AeCeAd+74On+49W/yba+/r/lKUSKihIjalAsREXWJgU5ElBAMdCKihGCgExElBAOdiCghGOhERAnBQCciSggGOhFRQvw/p6rhwuKa6UoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "RV.plotpdf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGBlJREFUeJzt3XtwVOd5x/HvoxviIglsSYDFRWCLi0Jw7CjYiTuOPXEz2E0hnUkz9kwml3pMk9TttMm0dScdN+P80TadNjOdum2YJk3jTuw4diahHlx3mho7cUJAxEaYi4wQYBYEktBdAnR7+ocWV5UFeySd3bN79veZ0cye3Zfd52XFz6/f857zmrsjIiLxUhB1ASIiEj6Fu4hIDCncRURiSOEuIhJDCncRkRhSuIuIxJDCXUQkhhTuIiIxpHAXEYmhoqg+uLKy0mtra6P6eBGRnHTgwIFOd69K1S6ycK+traWxsTGqjxcRyUlmdjpIO03LiIjEkMJdRCSGFO4iIjGkcBcRiSGFu4hIDKUMdzP7tpm1m9mb13jdzOzvzazFzJrM7PbwyxQRkZkIMnL/DrD1Oq/fD9Qlf3YA/zT3skREZC5SrnN391fNrPY6TbYD3/WJ/fr2mtliM1vu7m0h1ZjX+i+PcCjRy6GzvQxeGY26HBEJwUc2LuXWlYvT+hlhXMRUA5yZdJxIPveucDezHUyM7lm1alUIHx1fJzsH+cPvv0FTooer29yaRVuTiISjurw0J8J9usiZdtdtd98J7ARoaGjQztzXsLf1Ip//9wMUmPFH963j1pWL2VxTwZKFJVGXJiI5IoxwTwArJx2vAM6F8L556fkDCR77YROrbljAv352C6tuXBB1SSKSg8JYCrkL+HRy1cydQK/m22fnZ8c7+fIPDrJlzQ388It3KdhFZNZSjtzN7GngHqDSzBLAXwDFAO7+z8Bu4AGgBRgCPpeuYuNs4Moof/p8E2urFvKtz3yA0uLCqEsSkRwWZLXMQyled+D3QqsoT/3Vi0c513uJ5z7/IQW7iMyZrlDNAj9v6eTf977Nw3et4f2rl0RdjojEgMI9YkPDo/zJ802sqVzIlz+6PupyRCQmItusQyY8s+8Mie5LPP3Incwv0XSMiIRDI/cIjYyN862fnWRL7Q188OYboy5HRGJE4R6hF5rOcbbnEr/74bVRlyIiMaNwj4i7881XWqmrXsS966ujLkdEYkbhHpFX3urg2Pl+dty9loIC3TRGRMKlcI/IN19pZVl5KdvfVxN1KSISQwr3CLx5tpdftF7k4V9bQ0mRvgIRCZ+SJQLP/ypBSWEBn2xYmbqxiMgsKNwzbGzc+Y+Dbdy7oYqKBcVRlyMiMaVwz7Cfn+ikc+AKH9dcu4ikkcI9w370+jnK5hVx7wYtfxSR9FG4Z9DlkTFeOnyerZuW6c6PIpJWCvcM+snRdgaujPLx2zQlIyLppXDPoB+9cZbqsnncuVb3kRGR9FK4Z0jP0DB7mtvZdutNFOqKVBFJM4V7hvzXkQuMjDnb3ndT1KWISB5QuGfInuZ2lpbP4701FVGXIiJ5QOGeASNj4/z0eCf3rKvGTFMyIpJ+CvcM+NXpbvovj3LvhqqoSxGRPKFwz4CXmzsoKjDuuqUy6lJEJE8o3DNgT3M7H6i9gbJS3UtGRDJD4Z5mbb2XOHa+n3vWa0pGRDJH4Z5me5o7AHQvGRHJKIV7mu1pbqdm8XzqqhdFXYqI5BGFexoNj47zs+OdfHh9lZZAikhGKdzTqPFUF4PDY9y7XlMyIpJZCvc0eu1EJ0UFxodu1o3CRCSzFO5ptO9kF5tqKlg4ryjqUkQkzwQKdzPbambNZtZiZo9N8/oqM3vZzF43syYzeyD8UnPL5ZExDp7p5Y41N0RdiojkoZThbmaFwJPA/UA98JCZ1U9p9ufAs+5+G/Ag8I9hF5pr3jjTw/DYOFsU7iISgSAj9y1Ai7u3uvsw8AywfUobB8qTjyuAc+GVmJv2nezCDBpWK9xFJPOCTAbXAGcmHSeAO6a0+SrwX2b2+8BC4L5Qqsth+052sWFZORULdMsBEcm8ICP36RZo+5Tjh4DvuPsK4AHgKTN713ub2Q4zazSzxo6OjplXmyNGxsY5cLpb8+0iEpkg4Z4AVk46XsG7p10eBp4FcPdfAKXAu26B6O473b3B3RuqquJ7r5U3z/ZyaWRM8+0iEpkg4b4fqDOzNWZWwsQJ011T2rwNfATAzDYyEe7xHZqnsO9kFwAfqFW4i0g0Uoa7u48CjwIvAUeZWBVz2MyeMLNtyWZfBh4xs4PA08Bn3X3q1E3e2Heyi7WVC6kqmxd1KSKSpwJdXePuu4HdU557fNLjI8Bd4ZaWm8bHnf2nunjgvcujLkVE8piuUA1Z84V++i6Par5dRCKlcA/Z1fl2hbuIREnhHrLX3+6mumweK5YsiLoUEcljCveQNSV6uXXl4qjLEJE8p3APUe+lEVo7B3mfwl1EIqZwD9GhRC8Am1dURFyJiOQ7hXuIDiZ6ANhco5G7iERL4R6ig2d6WFO5UDcLE5HIKdxDdDDRw62akhGRLKBwD8n53stc6LvC5hWakhGR6CncQ3J1vl3LIEUkGyjcQ9KU6KGowHjPTeWpG4uIpJnCPSRNiV7WLS2jtLgw6lJERBTuYXB3Dp7p0ZSMiGQNhXsITl0cou/yqFbKiEjWULiH4OAZnUwVkeyicA9BU6KX0uIC6qoXRV2KiAigcA/FkbZeNiwrp6hQf50ikh2URnPk7hxt62fjci2BFJHsoXCfo7bey/ReGqF+eVnUpYiIvEPhPkdH2/oANHIXkayicJ+jq+G+QeEuIllE4T5HR9v6WXXDAhbNK4q6FBGRdyjc5+hoWx8bNd8uIllG4T4HQ8OjnLw4qPl2Eck6Cvc5aD7fj7tOpopI9lG4z8GR5MnUeoW7iGQZhfscHG3ro2xeESuWzI+6FBGR/0fhPgdH2/rZsLwMM4u6FBGR/0fhPkvj486xtj7Nt4tIVlK4z9KZ7iEGh8cU7iKSlQKFu5ltNbNmM2sxs8eu0eaTZnbEzA6b2ffCLTP76LYDIpLNUl5WaWaFwJPArwMJYL+Z7XL3I5Pa1AF/Btzl7t1mVp2ugrPFkbZ+CgzWL9UFTCKSfYKM3LcALe7e6u7DwDPA9iltHgGedPduAHdvD7fM7NN8vo/aGxcyv0QbYotI9gkS7jXAmUnHieRzk60D1pnZa2a218y2TvdGZrbDzBrNrLGjo2N2FWeJ4+0D1C3Vzksikp2ChPt06/x8ynERUAfcAzwE/IuZvWtDUXff6e4N7t5QVVU101qzxpXRMU5fHKKuWlMyIpKdgoR7Alg56XgFcG6aNj929xF3Pwk0MxH2sXSqc4ixcdfIXUSyVpBw3w/UmdkaMysBHgR2TWnzI+BeADOrZGKapjXMQrPJ8fZ+AG7RhtgikqVShru7jwKPAi8BR4Fn3f2wmT1hZtuSzV4CLprZEeBl4I/d/WK6io7a8QsDmMHNVQp3EclOgXaYcPfdwO4pzz0+6bEDX0r+xF5LxwCrblhAabFWyohIdtIVqrPQcmGAOk3JiEgWU7jP0OjYOK2dA9yilTIiksUU7jN0umuIkTHXyF1EsprCfYaOXxgA0DJIEclqCvcZakkug9RKGRHJZgr3GTrePkDN4vksnBdooZGISCQU7jN0/ILuKSMi2U/hPgNj486JDi2DFJHsp3CfgUT3EFdGx3XDMBHJegr3Gbi6UuYWTcuISJZTuM/A8fZkuGtaRkSynMJ9Bo6397O0fB7lpcVRlyIicl0K9xloaR/QqF1EcoLCPSB350T7gE6mikhOULgH1NZ7mcHhMW7WyF1EcoDCPaCWqydTddsBEckBCveAWrRSRkRyiMI9oJaOASrmF1O5qCTqUkREUlK4B3R1pYyZRV2KiEhKCveAWtp1TxkRyR0K9wC6BofpGhzWfLuI5AyFewBXT6ZqGaSI5AqFewBaBikiuUbhHkBL+wDziwupWTw/6lJERAJRuAfQ0jHA2qqFFBRopYyI5AaFewAndMMwEckxCvcUBq+McrbnkubbRSSnKNxTONGh2w6ISO5RuKdwdaVMnbbWE5EconBPoaV9gKICY/WNC6MuRUQksEDhbmZbzazZzFrM7LHrtPuEmbmZNYRXYrRa2gdYfeMCigv130ERyR0pE8vMCoEngfuBeuAhM6ufpl0Z8AfAL8MuMkrHtfuSiOSgIMPRLUCLu7e6+zDwDLB9mnZfA74OXA6xvkhdGh7j1MVBNixXuItIbgkS7jXAmUnHieRz7zCz24CV7v5CiLVF7q0L/bjDhmXlUZciIjIjQcJ9ussy/Z0XzQqAbwBfTvlGZjvMrNHMGjs6OoJXGZFj5/sA2KiRu4jkmCDhngBWTjpeAZybdFwGbAL2mNkp4E5g13QnVd19p7s3uHtDVVXV7KvOkKNt/SwoKWTlkgVRlyIiMiNBwn0/UGdma8ysBHgQ2HX1RXfvdfdKd69191pgL7DN3RvTUnEGHTvfx/plZbqnjIjknJTh7u6jwKPAS8BR4Fl3P2xmT5jZtnQXGBV359j5fjYs05SMiOSeoiCN3H03sHvKc49fo+09cy8rehf6rtAzNKKTqSKSk3RlzjUcTZ5M1chdRHKRwv0ajrX1A1oGKSK5SeF+Dc3n+7ipopSKBcVRlyIiMmMK92s4dr6fDcs1aheR3KRwn8bw6Dgt7QOabxeRnKVwn8aJjgFGx10jdxHJWQr3abxz2wGN3EUkRyncp3GsrZ+SwgLWVGqDDhHJTQr3aRw930/d0kUUaYMOEclRSq8p3J0j5/q0vl1EcprCfYoLfVfoHLjCe2sU7iKSuxTuUzQlegDYvHJxxJWIiMyewn2KQ2d7KSww6rUMUkRymMJ9iqZEL+uWllFaXBh1KSIis6Zwn8TdOXS2l801FVGXIiIyJwr3Sc72XKJrcJj3rlC4i0huU7hPcijRC8BmhbuI5DiF+yRNZ3spLjTW67YDIpLjFO6TNCV62LCsnHlFOpkqIrlN4Z7k7jQlejXfLiKxoHBPOn1xiP7Lo1opIyKxoHBPajo7cTJVI3cRiQOFe9KhRA8lRQWsW6qTqSKS+xTuSU2JXuqXl1Os2/yKSAwoyYCxcefNs71a3y4isaFwB4629TE4PMb7Vy+JuhQRkVAo3IHGU10AfKD2hogrEREJh8Id2H+6m5rF87lp8fyoSxERCUXeh7u703iqi4ZaTcmISHzkfbgnui9xoe8KDZqSEZEYyftw3//OfLtG7iISH4HC3cy2mlmzmbWY2WPTvP4lMztiZk1m9hMzWx1+qemx/1QXZaVFrKvWxUsiEh8pw93MCoEngfuBeuAhM6uf0ux1oMHdNwPPAV8Pu9B02X+qm4bVSygosKhLEREJTZCR+xagxd1b3X0YeAbYPrmBu7/s7kPJw73AinDLTI+uwWFa2gc03y4isRMk3GuAM5OOE8nnruVh4MXpXjCzHWbWaGaNHR0dwatMkwOnuwGtbxeR+AkS7tPNV/i0Dc0+BTQAfzPd6+6+090b3L2hqqoqeJVp0niqi5LCAt12QERipyhAmwSwctLxCuDc1EZmdh/wFeDD7n4lnPLSa/+pLjavqKC0WDsviUi8BBm57wfqzGyNmZUADwK7Jjcws9uAbwLb3L09/DLDNzQ8yqGzvbxfSyBFJIZShru7jwKPAi8BR4Fn3f2wmT1hZtuSzf4GWAT8wMzeMLNd13i7rPHL1i5Gxpxfu6Uy6lJEREIXZFoGd98N7J7y3OOTHt8Xcl1p98pbHZQWF+hkqojEUt5eofrq8Q7uWHOj5ttFJJbyMtwT3UO0dgxy97roV+yIiKRDXob7q291AnB3nebbRSSe8jTcO1heUcot1YuiLkVEJC3yLtxHx8Z57UQnd9dVYab7yYhIPOVduB9M9NB/eVTz7SISa3kX7q+81UmBofXtIhJreRfur77Vwa0rF1OxoDjqUkRE0iavwr1rcJimRA9312lKRkTiLa/C/T/fPM+4w0ffszTqUkRE0iqvwv2FpnOsrVxI/fLyqEsREUmrvAn3jv4r7G29yMc2L9cSSBGJvbwJ9/98s41xh9/YfFPUpYiIpF3ehPt/NLVRV72I9cvKoi5FRCTt8iLcL/RdZv+pLj6mUbuI5Im8CPfdh9pwh4/dujzqUkREMiIvwv2FpjY2Li/n5irdKExE8kPsw/1M1xAHTnfzsc0atYtI/oh9uD+19zSFBcZv3VYTdSkiIhkT63AfvDLK0/ve5v5Ny7hp8fyoyxERyZhYh/sPf5Wg//Ion7trTdSliIhkVGzDfXzc+dfXTnHrysXcvmpx1OWIiGRUbMP9leMdtHYO8jt31ep2AyKSd2Ib7t/+2Umqy+Zx/yatkhGR/BPLcD98rpefHu/k0x9cTUlRLLsoInJdsUu+8XHnL358mCULivnUnaujLkdEJBKxC/fnf5Wg8XQ3f/bARhYvKIm6HBGRSMQq3HuGhvnLF4/RsHoJn7h9RdTliIhEJlbh/vWXmum9NMLXPr6JggKtkBGR/BWbcN/T3M7T+97mcx+qZaO20RORPBco3M1sq5k1m1mLmT02zevzzOz7ydd/aWa1YRd6PS8fa2fHUwfYuKycP/z1dZn8aBGRrJQy3M2sEHgSuB+oBx4ys/opzR4Gut39FuAbwF+HXei1/PeRC/zuUwdYt3QR33vkDhbNK8rUR4uIZK0gSbgFaHH3VgAzewbYDhyZ1GY78NXk4+eAfzAzc3cPsdZ3jIyNs7f1IrsPnee5A2eoX17Od3/nDioWFKfj40REck6QcK8Bzkw6TgB3XKuNu4+aWS9wI9AZRpGTfX//2/zli8foGRphYUkhv7n5Jr66/T2UlyrYRUSuChLu0y07mToiD9IGM9sB7ABYtWpVgI9+t6Xlpdy7vpr7Ny3j7nVVlBYXzup9RETiLEi4J4CVk45XAOeu0SZhZkVABdA19Y3cfSewE6ChoWFWUzb3rK/mnvXVs/mjIiJ5I8hqmf1AnZmtMbMS4EFg15Q2u4DPJB9/AvifdM23i4hIailH7sk59EeBl4BC4NvuftjMngAa3X0X8C3gKTNrYWLE/mA6ixYRkesLtG7Q3XcDu6c89/ikx5eB3w63NBERma3YXKEqIiL/R+EuIhJDCncRkRhSuIuIxJDCXUQkhiyq5ehm1gGcnuUfryQNtzbIcupzflCf88Nc+rza3atSNYos3OfCzBrdvSHqOjJJfc4P6nN+yESfNS0jIhJDCncRkRjK1XDfGXUBEVCf84P6nB/S3uecnHMXEZHry9WRu4iIXEdWh3u2b8ydDgH6/CUzO2JmTWb2EzNbHUWdYUrV50ntPmFmbmY5v7IiSJ/N7JPJ7/qwmX0v0zWGLcDv9ioze9nMXk/+fj8QRZ1hMbNvm1m7mb15jdfNzP4++ffRZGa3h1qAu2flDxO3Fz4BrAVKgINA/ZQ2XwT+Ofn4QeD7UdedgT7fCyxIPv5CPvQ52a4MeBXYCzREXXcGvuc64HVgSfK4Ouq6M9DnncAXko/rgVNR1z3HPt8N3A68eY3XHwBeZGInuzuBX4b5+dk8cn9nY253Hwaubsw92Xbg35KPnwM+YmbTbfmXK1L22d1fdveh5OFeJnbGymVBvmeArwFfBy5nsrg0CdLnR4An3b0bwN3bM1xj2IL02YHy5OMK3r3jW05x91eZZke6SbYD3/UJe4HFZrY8rM/P5nCfbmPummu1cfdR4OrG3LkqSJ8ne5iJ//LnspR9NrPbgJXu/kImC0ujIN/zOmCdmb1mZnvNbGvGqkuPIH3+KvApM0swsX/E72emtMjM9N/7jATarCMioW3MnUMC98fMPgU0AB9Oa0Xpd90+m1kB8A3gs5kqKAOCfM9FTEzN3MPE/5391Mw2uXtPmmtLlyB9fgj4jrv/rZl9kInd3Ta5+3j6y4tEWvMrm0fuM9mYm+ttzJ1DgvQZM7sP+Aqwzd2vZKi2dEnV5zJgE7DHzE4xMTe5K8dPqgb93f6xu4+4+0mgmYmwz1VB+vww8CyAu/8CKGXiHixxFejf+2xlc7jn48bcKfucnKL4JhPBnuvzsJCiz+7e6+6V7l7r7rVMnGfY5u6N0ZQbiiC/2z9i4uQ5ZlbJxDRNa0arDFeQPr8NfATAzDYyEe4dGa0ys3YBn06umrkT6HX3ttDePeozyinONj8AvMXEWfavJJ97gol/3DDx5f8AaAH2AWujrjkDff5v4ALwRvJnV9Q1p7vPU9ruIcdXywT8ng34O+AIcAh4MOqaM9DneuA1JlbSvAF8NOqa59jfp4E2YISJUfrDwOeBz0/6jp9M/n0cCvv3WleoiojEUDZPy4iIyCwp3EVEYkjhLiISQwp3EZEYUriLiMSQwl1EJIYU7iIiMaRwFxGJof8FWXRdjJNr3xQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "RV.plotcdf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}