{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using the dcnorm module\n",
    "\n",
    "* This notebook teaches some of the basic functionality of the `dcnorm` module.\n",
    "\n",
    "* The `NCGR` package (https://adirkson.github.io/sea-ice-timing/) should be installed before running this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a first step, we'll import the modules used in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from NCGR.dcnorm import dcnorm_gen\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First define variables for the **minimum** ($a$) and **maximum** ($b$) values that the DCNORM distribution takes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=120 # minimum\n",
    "b=273 # maximum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now instantiate the `dcnorm_gen` class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "dcnorm = dcnorm_gen(a=a, b=b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now make a DCNORM distribution object with the fixed `a` and `b` values, and some arbitrary parameter values $\\mu$ and $\\sigma$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mu variable for the DCNORM distribution\n",
    "m=132.\n",
    "# sigma variable for the DCNORM distribution\n",
    "s=20.\n",
    "# freeze a dcnorm distribution object\n",
    "rv = dcnorm(m,s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`rv` is now a frozen object representing a random variable described by the DCNORM distribution with the given parameters (try some others, too!); it has several methods (see documentation for `dcnorm_gen`. We'll go through some main ones now.\n",
    "\n",
    "For instance, its PDF can be calculated and plotted simply with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD4CAYAAAD2FnFTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAd90lEQVR4nO3de5gU1ZnH8e/LDIyCK67DaFxuwwom4l0nmBjjmrAqXtEVFR+S4EoWY8C7Ioq6CnEjeEPXK4oJCBGJ0RWTddGo0cREZFC8IKIjIgwQBEHEG4i8+8epkaatmemZ6Znqy+/zPPV0d/Xp6rcLih9Vp6qOuTsiIiLp2iVdgIiI5CYFhIiIxFJAiIhILAWEiIjEUkCIiEis0qQLaIouXbp4ZWVl0mWIiOSVefPmrXH3iqZ+Lq8CorKykurq6qTLEBHJK2b2XnM+p0NMIiISSwEhIiKxFBAiIhJLASEiIrEUECIiEqvwA2L6dKishHbtwuP06UlXJCKSF/LqNNcmmz4dhg+HTz8Nr997L7wGGDIkubpERPJAYe9BjBmzNRzqfPppmC8iIg0q7IBYurRp80VE5CuFHRA9ejRtvoiIfKWwA+Laa6Fjx23ndewY5ouISIMKOyCGDIFJk6BnTzALj5MmqYNaRCQDhX0WE4QwUCCIiDRZYe9BiIhIsykgREQklgJCRERiKSBERCSWAkJERGIpIEREJJYCQkREYikgREQklgJCRERiKSBERCRWRgFhZgPMbJGZ1ZjZ6Jj3y8zswej9OWZWGc0/wszmmdlr0eMPUz7zp2iZ86Npl2z9KBERablG78VkZiXA7cARQC0w18xmufsbKc2GAevcvbeZDQbGA6cBa4Dj3X2Fme0NzAa6pnxuiLtXZ+m3iIhIFmWyB9EPqHH3xe6+CZgBDExrMxCYEj1/COhvZubuL7v7imj+AmA7MyvLRuEiItK6MgmIrsCylNe1bLsXsE0bd98MrAfK09qcDLzs7htT5v0qOrx0pZlZ3Jeb2XAzqzaz6tWrV2dQroiIZEMmARH3D7c3pY2Z7UU47HRWyvtD3H0f4PvR9OO4L3f3Se5e5e5VFRUVGZQrIiLZkElA1ALdU153A1bU18bMSoHOwNrodTfgEeAn7v5O3QfcfXn0uAH4DeFQloiI5IhMAmIu0MfMeplZB2AwMCutzSxgaPR8EPC0u7uZ7QT8AbjM3Z+va2xmpWbWJXreHjgOeL1lP0VERLKp0YCI+hRGEs5AWgjMdPcFZjbWzE6Imk0Gys2sBrgQqDsVdiTQG7gy7XTWMmC2mb0KzAeWA/dk84eJiEjLmHt6d0Luqqqq8upqnRUrItIUZjbP3aua+jldSS0iIrEUECIiEksBISIisRQQIiISSwEhIiKxFBAiIhJLASEiIrEUECIiEksBISIisRQQIiISSwEhIiKxFBAiIhJLASEiIrEUECIiEksBISIisRQQIiISSwEhIiKxFBAiIhJLASEiIrEUECIiEksBISIisRQQIiISSwEhIiKxFBAiIhJLASEiIrEUECIiEksBISIisTIKCDMbYGaLzKzGzEbHvF9mZg9G788xs8po/hFmNs/MXosef5jymYOi+TVmdquZWbZ+lIiItFyjAWFmJcDtwNFAX+B0M+ub1mwYsM7dewM3A+Oj+WuA4919H2AocH/KZ+4EhgN9omlAC36HiIhkWWkGbfoBNe6+GMDMZgADgTdS2gwEro6ePwTcZmbm7i+ntFkAbGdmZcDOwI7u/rdomVOBE4HHW/Bb2taSJfDss1BTA2vXQrt20L077LsvHHYYdOyYdIUiIi2SSUB0BZalvK4FDq6vjbtvNrP1QDlhD6LOycDL7r7RzLpGy0ldZte4Lzez4YQ9DXr06JFBua3o889h+nS4/XZ4Ocq+khLYaSfYvBnWrw/zyspg0CA491zo1y+5ekVEWiCTPoi4vgFvShsz24tw2OmsJiwzzHSf5O5V7l5VUVGRQbmtwB0eegi+9S346U9DGNx0E7z2Gnz2GaxZAx9+COvWwezZMGwYPPYYHHwwnHZa2NsQEckzmQRELdA95XU3YEV9bcysFOgMrI1edwMeAX7i7u+ktO/WyDJzw4cfwqmnwimnwI47whNPwCuvwAUXwN57Q/v2W9vutBMceWTYw1i+HP7zP0NQ7Ltv2PMQEckjmQTEXKCPmfUysw7AYGBWWptZhE5ogEHA0+7uZrYT8AfgMnd/vq6xu68ENpjZd6Kzl34CPNrC35J9CxbAgQfC//wPXHcdvPQSHHEEZHLC1Q47wNVXw8KFsN9+8KMfwVlnhb0PEZE80GhAuPtmYCQwG1gIzHT3BWY21sxOiJpNBsrNrAa4EKg7FXYk0Bu40szmR9Mu0XtnA/cCNcA75FoH9Z//DIceGvodnnsOLr0USjPpsknTsyc88wyMHg2TJsHxx8OGDdmvV0Qky8w99tB/TqqqqvLq6urW/6Knn4Zjjw3/uM+eHR6z4Z574Oyz4ZBD4PHHoVOn7CxXRKQBZjbP3aua+jldSZ3ub3+DE06A3XeHv/wle+EA8B//AdOmwfPPh+/4/PPsLVtEJMsUEKkWLoRjjoHddoMnn4QuXbL/HYMHw69/HfZS/v3fwxlSIiI5qBkH1QvU2rWhf6BDhxAOu+3Wet/14x/DihWhX2KPPeCaa1rvu0REmkkBAeHMolNPhWXLQodyZWXrf+eoUfDWWzB2LOy1V/h+EZEcokNMEK5XeOopuPvu0IHcFszgzjvhu98NF9+9/XbbfK+ISIYUEE8/Db/8ZfhH+owz2va7O3SAGTPCxXannqpOaxHJKcUdEGvWhAvYvvlNmDgxmRp69IApU2D+/NAnISKSI4o3INzDaacffAAPPJDsNQnHHQcjR8Itt4Q7xIqI5IDiDYjf/jbcQuMXv4D990+6mnArj913D6e+fvxx0tWIiBRpQHzwAZxzDhx0ULjpXi7o1ClcH7FkSTjDSUQkYcUZEBdeGK57mDy5efdXai2HHgrnnx/ObnrhhaSrEZEiV3wB8cQTMHVquPnefvslXc3XjR0LXbuGezbpzq8ikqDiCohNm8Khpd694Yorkq4m3g47hDOq5s8PexIiIgkproD47/8OVy9PnAjbbZd0NfU7+eQw8NAVV8DKlUlXIyJFqngCYtWqcPjm6KPDrbxzmRncdlu4cE4d1iKSkOIJiMsvD+NH33xz0pVkpk+f0Jk+bRrMm5d0NSJShIojIObOhfvug/POC1dN54vLLoOKCrjoIt0WXETaXHEExMyZsOuucOWVSVfSNDvuGMa1fvZZeOyxpKsRkSJTHEOOusPy5dCtW/aLam1ffAH77gtbtsDrr4cb+4mINIGGHG2IWX6GA4RAmDAhnH11991JVyMiRaQ4AiLfHXccHH44jBsHn3ySdDUiUiQUEPnADK69Ft5/P5z+KiLSBhQQ+eKQQ8I1HBMmwEcfJV2NiBQBBUQ+GTcu3GQwqcGNRKSoKCDyyUEHwUknwY03hqAQEWlFCoh8c801sGED3HBD0pWISIFTQOSbffaB004Lw5OuXp10NSJSwDIKCDMbYGaLzKzGzEbHvF9mZg9G788xs8pofrmZPWNmH5vZbWmf+VO0zPnRtEs2flBRuOqqcF8p9UWISCtqNCDMrAS4HTga6AucbmZ905oNA9a5e2/gZmB8NP9z4Erg4noWP8Td94+m95vzA4rSnnvCoEHh9uXr1iVdjYgUqEz2IPoBNe6+2N03ATOAgWltBgJToucPAf3NzNz9E3f/CyEoJJvGjAl9EbouQkRaSSYB0RVYlvK6NpoX28bdNwPrgfIMlv2r6PDSlWZmcQ3MbLiZVZtZ9Wodc99qv/3g+OPDYaYNG5KuRkQKUCYBEfcPd/od/jJpk26Iu+8DfD+afhzXyN0nuXuVu1dVVFQ0WmxRGTMmnO56111JVyIiBSiTgKgFuqe87gasqK+NmZUCnYEGT9R39+XR4wbgN4RDWdIUBx8MRxwRrov47LOkqxGRApNJQMwF+phZLzPrAAwGZqW1mQUMjZ4PAp72Bu4jbmalZtYlet4eOA54vanFC2Hc6lWrYPLkpCsRkQLTaEBEfQojgdnAQmCmuy8ws7FmdkLUbDJQbmY1wIXAV6fCmtkS4CbgDDOrjc6AKgNmm9mrwHxgOXBP9n5WETnsMPj+92H8eNi0KelqRKSAFMeAQYVu9mwYMADuvReGDUu6GhHJMRowqJgdeSQccABcf30YeU5EJAsUEIXADEaNgkWLNHa1iGSNAqJQDBoElZVhvAgRkSxQQBSK0lK46CL461/h+eeTrkZECoACopCceSaUl2svQkSyQgFRSDp2hHPOgVmz4I03kq5GRPKcAqLQjBgB22+vAYVEpMUUEIWmS5dwLcS0abB8edLViEgeU0AUoosuCtdD3HJL0pWISB5TQBSiyko49dRwl9f165OuRkTylAKiUF1ySRgnQrcCF5FmUkAUqgMOCLcCnzgRNm5MuhoRyUMKiEI2ahT8/e9w//1JVyIieUgBUcj694cDDwynvOomfiLSRAqIQpZ6E79Z6WM8iYg0TAFR6E4+GXr1CgMK5dHYHyKSPAVEoSsthYsvhhde0E38RKRJFBDF4IwzwhXW48cnXYmI5BEFRDGou4nf738PCxYkXY2I5AkFRLEYMSIExfXXJ12JiOQJBUSxKC+Hn/4Upk+H2tqkqxGRPKCAKCYXXBDOZJo4MelKRCQPKCCKSWUlDB4Md98N69YlXY2I5DgFRLG55BL4+GPdxE9EGqWAKDb77QdHHRXGivj886SrEZEcpoAoRpdeCqtWwdSpSVciIjlMAVGMDj8cqqrCTfy+/DLpakQkR2UUEGY2wMwWmVmNmY2Oeb/MzB6M3p9jZpXR/HIze8bMPjaz29I+c5CZvRZ95lYzs2z8IMlA3U383n4bHn006WpEJEc1GhBmVgLcDhwN9AVON7O+ac2GAevcvTdwM1B3T4fPgSuBi2MWfScwHOgTTQOa8wOkmf7t32D33XUTPxGpVyZ7EP2AGndf7O6bgBnAwLQ2A4Ep0fOHgP5mZu7+ibv/hRAUXzGz3YAd3f1v7u7AVODElvwQaaKSknATvxdfhOeeS7oaEclBmQREV2BZyuvaaF5sG3ffDKwHyhtZZurlvHHLBMDMhptZtZlVr169OoNyJWNDh8Iuu8CECUlXIiI5KJOAiOsbSD8mkUmbZrV390nuXuXuVRUVFQ0sUpps++3h3HPhf/8XXnst6WpEJMdkEhC1QPeU192AFfW1MbNSoDOwtpFldmtkmdIWzj4bOnXSXoSIfE0mATEX6GNmvcysAzAYSB+/chYwNHo+CHg66luI5e4rgQ1m9p3o7KWfADqdJgk77wzDh8MDD8DixUlXIyI5pNGAiPoURgKzgYXATHdfYGZjzeyEqNlkoNzMaoALga9OhTWzJcBNwBlmVptyBtTZwL1ADfAO8Hh2fpI02cUXh5Hnrrsu6UpEJIdYA//RzzlVVVVeXV2ddBmFacQIuOceqKmBHj2SrkZEssjM5rl7VVM/pyupJbj00vCoYUlFJKKAkKBHjzB29eTJsELnC4iIAkJSjR4NmzdrWFIRARQQkuqf/xl+9KMwoNCqVUlXIyIJU0DIti6/HDZuhBtvTLoSEUmYAkK2tccecNppcMcdsGZN0tWISIIUEPJ1V1wBn36qvQiRIqeAkK/r2xcGD4Zbb4W//z3pakQkIQoIiXfNNaEv4pe/TLoSEUmIAkLi9ekTrou46y5YujTpakQkAQoIqd9VV4XHceOSrUNEEqGAkPr16AE/+xn86ldh/GoRKSoKCGnY5ZdDWRlcfXXSlYhIG1NASMN23TWMOvfAAxp1TqTIKCCkcZdcAp07b73jq4gUBQWENG7nncPFc48/Dk8+mXQ1ItJGFBCSmZEjoVevMPrcl18mXY2ItAEFhGSmrCwMSfrqqzBlStLViEgbUEBI5k45BQ4+OBxu+vjjpKsRkVamgJDMmcFNN8HKlbqRn0gRUEBI0xxyCAwaBBMmQG1t0tWISCtSQEjTTZgAW7bARRclXYmItCIFhDRdr15w2WUwcyb88Y9JVyMirUQBIc0zalQYw/qcc2DTpqSrEZFWoICQ5tluuzCg0JtvwsSJSVcjIq1AASHNd+yxMHAgjB0Ly5YlXY2IZJkCQlpm4sTQYT1iBLgnXY2IZFFGAWFmA8xskZnVmNnomPfLzOzB6P05ZlaZ8t5l0fxFZnZUyvwlZvaamc03s+ps/BhJQGUl/OIX8NhjMGNG0tWISBY1GhBmVgLcDhwN9AVON7O+ac2GAevcvTdwMzA++mxfYDCwFzAAuCNaXp0fuPv+7l7V4l8iyTnvvHCF9bnnwurVSVcjIlmSyR5EP6DG3Re7+yZgBjAwrc1AoO4GPQ8B/c3Movkz3H2ju78L1ETLk0JSUgKTJ8NHH4WQEJGCkElAdAVSeyBro3mxbdx9M7AeKG/ksw48YWbzzGx4fV9uZsPNrNrMqlfrf6e5a6+94Morw2GmRx9NuhoRyYJMAsJi5qX3RtbXpqHPfs/dDyQcuhphZofFfbm7T3L3KnevqqioyKBcScyll8J++8FZZ8H77yddjYi0UCYBUQt0T3ndDVhRXxszKwU6A2sb+qy71z2+DzyCDj3lv/btYdo0+PBDGDZMZzWJ5LlMAmIu0MfMeplZB0Kn86y0NrOAodHzQcDT7u7R/MHRWU69gD7Ai2bWycz+AcDMOgFHAq+3/OdI4vbeG66/Hn7/e7jzzqSrEZEWKG2sgbtvNrORwGygBLjP3ReY2Vig2t1nAZOB+82shrDnMDj67AIzmwm8AWwGRrj7l2a2K/BI6MemFPiNu/9fK/w+ScLIkWF40osugsMPh77pJ72JSD4wz6PDAFVVVV5drUsm8sKqVbDPPvCNb8ALL0DHjklXJFK0zGxecy4n0JXU0jp23TX0R7z+eui0zqP/iIhIoICQ1nPkkeE+TdOmwR13JF2NiDSRAkJa1+WXw3HHwfnnw1//mnQ1ItIECghpXe3awf33Q8+ecPLJsHRp0hWJSIYUENL6dtopXF396adhb+Kjj5KuSEQyoICQtrHXXvC738HChXDKKfDFF0lXJCKNUEBI2/nXf4W774YnnoCf/1xnNonkuEYvlBPJqjPPhHffDWNI7Lgj3HADWNwtu0QkaQoIaXtjx8L69XDTTdCpU3gtIjlHASFtzywMVfrZZzBuHGy3XTgdVkRyigJCktGuHdx1VwiJMWPCGU7jxulwk0gOUUBIckpKYMqUcJ+ma6+FtWvhtttCeIhI4hQQkqySknBm0847w/jxISR+/etw2ElEEqWAkOSZwXXXQXk5jBoFS5bAI4/AbrslXZlIUdO+vOSOSy6Bhx8Od4CtqgLd2l0kUQoIyS0nnRRu6te+PRx6aOiT0AV1IolQQEju2XdfmDsX+veHc86BgQNhzZqkqxIpOgoIyU0VFWFc64kTYfbsMDrdww9rb0KkDSkgJHeZwXnnwZw5YYS6k0+GE0+EZcuSrkykKCggJPftv3845DRhAjz5JOy5Z7io7pNPkq5MpKApICQ/tG8fznJasACOOgquugr69IFJk2DTpqSrEylICgjJL716hXElnn8eKivhrLOgd+/QV/Hxx0lXJ1JQFBCSnw45JITEH/4QQuOCC8KwppdeCm+/nXR1IgVBASH5ywyOOQaefTZcO/Ev/wI33gh77AE/+AFMnQoffph0lSJ5SwEhheG73w2nwS5dCv/1X+Fx6FDYZRc4+mi4915YuTLpKkXyinkenVdeVVXl1br9gmRiy5Zw5tNDD4U+i3ffDfP33DNcgNe/P3zve+F6C5ECZ2bz3L2qyZ9TQEjBc4dXXgljYT/1FPz5z2EcCgj9Ft/+drj30/77w7e+Bd2765bjUlBaNSDMbABwC1AC3Ovu16W9XwZMBQ4CPgBOc/cl0XuXAcOAL4Fz3X12JsuMo4CQrNi4MVx89+KLYS+juhoWL976/vbbwze/GaaePUNg1E3dukGXLgoQaRvTp4cBtZYuhR49wrgpQ4Y0eTGtFhBmVgK8BRwB1AJzgdPd/Y2UNj8H9nX3n5nZYOAkdz/NzPoCDwD9gH8C/gjsEX2swWXGUUBIq/ngg3CNxZtvhmnRojAtW/b16yzatQvjV3TpEm5RXvfYuTPssEP81LEjlJVBhw4NT6WlGlVPgunTYfjwMNpinY4dw7U/TQyJ5gZEJuNB9ANq3H1x9EUzgIFA6j/mA4Gro+cPAbeZmUXzZ7j7RuBdM6uJlkcGyxRpO+XlcNhhYUq1ZQusXh2Com5avTrcPPCDD8Lju++GPZGPPgpXd7f0sG1JSQih1Mns6/Pqm8ziQ6a+4GnK/NZq29D8YvXWW7B587bzPv007FE0Yy+iOTIJiK5A6s1vaoGD62vj7pvNbD1QHs1/Ie2zXaPnjS0TADMbDgwH6NGjRwblimRRu3bhPlC77hr6KRqzZUvo3/jkk3DhXuq0aVOYvvhi6/O46csvw3LiJvf630ud0tUXWk2Z31ptG5pfzN6o5//LS5e2WQmZBERcrKf/adbXpr75cQdwY/+GuPskYBKEQ0z1lymSA9q1g06dwrTLLklXI/msshLee+/r89vwP8qZ9LTVAt1TXncDVtTXxsxKgc7A2gY+m8kyRUSK17XXhj6HVB07hvltJJOAmAv0MbNeZtYBGAzMSmszCxgaPR8EPO2h93sWMNjMysysF9AHeDHDZYqIFK8hQ0KHdM+eoX+mZ89mdVC3RKOHmKI+hZHAbMIpqfe5+wIzGwtUu/ssYDJwf9QJvZbwDz5Ru5mEzufNwAh3/xIgbpnZ/3kiInlsyJA2DYR0ulBORKTANfc0V13tIyIisRQQIiISSwEhIiKxFBAiIhIrrzqpzWw1EHPlSEa6AGuyWE625XJ9uVwbqL6WUn3Nl8u1wdb6erp7k+9tn1cB0RJmVt2cXvy2ksv15XJtoPpaSvU1Xy7XBi2vT4eYREQklgJCRERiFVNATEq6gEbkcn25XBuovpZSfc2Xy7VBC+srmj4IERFpmmLagxARkSZQQIiISKyCCQgzu8/M3jez11PmXW9mb5rZq2b2iJntlPLeZWZWY2aLzOyotq4t5b2LzczNrEv02szs1qi2V83swNasraH6zOycaP0sMLMJKfPbbN3VV5+Z7W9mL5jZfDOrNrN+0fw2XX9m1t3MnjGzhdF6Oi+av7OZPWlmb0eP/5hj9eXKthFbX8r7iW4fDdWXC9tHA3++2dk+3L0gJuAw4EDg9ZR5RwKl0fPxwPjoeV/gFaAM6AW8A5S0ZW3R/O6EW56/B3SJ5h0DPE4Yje87wJyE1t0PgD8CZdHrXZJYdw3U9wRwdMo6+1MS6w/YDTgwev4PwFvROpoAjI7mj075u5cr9eXKthFbX/Q68e2jgfWXE9tHA/VlZfsomD0Id3+OMBZF6rwn3L1u1O8XCCPXAQwEZrj7Rnd/F6gB+rVlbZGbgVFsO9zqQGCqBy8AO5nZbq1VWwP1nQ1c5+4bozbvp9TXZuuugfoc2DF63pmtIxK26fpz95Xu/lL0fAOwkDDu+kBgStRsCnBiLtWXQ9tGfesPcmD7aKC+nNg+GqgvK9tHwQREBs4kJCeEFbgs5b1atv6lbBNmdgKw3N1fSXsr8doiewDfN7M5ZvasmX07mp8r9Z0PXG9my4AbgMui+YnVZ2aVwAHAHGBXd18JYSMG6gaozpX6UuXEtpFaXy5uH2nrL+e2j7T6srJ9FEVAmNkYwoh20+tmxTRrs/N9zawjMAa4Ku7tmHlJnItcCvwjYTf0EmCmmRm5U9/ZwAXu3h24gDCqISRUn5ntAPwOON/dP2qoacy8xOrLlW0jtb6onpzaPmLWX05tHzH1ZWX7KPiAMLOhwHHAEI8OwhFSs3tKs25s3QVrC7sTjk++YmZLou9/ycy+kQO11akFHo52RV8EthBu/JUr9Q0FHo6e/5atu/FtXp+ZtSdsnNPdva6mVXW77tFj3SGIXKkvZ7aNmPpyavuoZ/3lzPZRT33Z2T5aq/MkiQmoZNuOzAGE8bAr0trtxbYdSYtp/Y7WbWpLe28JWzvhjmXbTqQXE1p3PwPGRs/3IOyWWhLrrp76FgKHR8/7A/OSWH/R90wFJqbNv55tO6kn5Fh9ObFt1FdfWpvEto8G1l9ObB8N1JeV7aPVVmxbT8ADwErgC0JKDiN0EC0D5kfTXSntxxDOMFhE1NvflrWlvZ+6ARhwe1Tba0BVQuuuAzANeB14CfhhEuuugfoOBeZFG+Mc4KAk1l9UhwOvpvw9OwYoB54C3o4ed86x+nJl24itL1e2jwbWX05sHw3Ul5XtQ7faEBGRWAXfByEiIs2jgBARkVgKCBERiaWAEBGRWAoIERGJpYAQEZFYCggREYn1/8qHy6Q0S0VtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(a, b, 1000) # discretize the range from a to b\n",
    "x_sub = x[(x!=a)&(x!=b)] # extract from those values where x is not a or b\n",
    "plt.figure()\n",
    "plt.plot(x_sub, rv.pdf(x_sub), color='r') # plot for a<x<b\n",
    "plt.plot(a, rv.pdf(a)*1e-1, 'o', color='r') # point mass at a (re-scale by 1/10 for plotting)\n",
    "plt.plot(b, rv.pdf(b)*1e-1, 'o', color='r') # point mass at b (re-scale by 1/10 for plotting)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the maginute of the circles have been reduced by a factor of 1e-1 to make it easier to see the shape of the PDF. It's not actually necessary to plot the different components of the PDF seperately; the reason for doing so was purely cosmetic. We could have also typed \n",
    "```python \n",
    "plt.figure()\n",
    "plt.plot(x,rv.pdf(x), color='r')\n",
    "plt.show()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A random sample can be drawn from the distribution with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = rv.rvs(size=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll now plot the this sample that was generated, along with the true PDF:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAeG0lEQVR4nO3de5xVdb3/8dcbRi6igOFoBeigoieM1BiQ+uWVMsCTk+dIolT2Uw+paRf1mB6LH1JWWCesh5RxUiMi0bwcp0Tx5z3PUWTwAqIRI6FM+NNREEVEHPn8/liLw2a7h9kwe9ibWe/n47Efe63v+q7Zn72YNW/WXRGBmZllV5dyF2BmZuXlIDAzyzgHgZlZxjkIzMwyzkFgZpZxDgIzs4xzEJiZZZyDwDJP0gpJb0taJ+llSTdI2kPSg5I2SHpT0huSFkq6VFL3nHknS3o3nXfz65Jyfh+z7eUgMEt8LiL2AD4ODAe+k7afHxF7Ah8CLgLGA3MlKWfemyJij5zXVTu1crN2chCY5YiIvwN3AR/Na38rIh4ETgI+AZy486sz6xgOArMckgYCY4EnC02PiBeBBuConVmXWUdyEJgl/lPS68AjwEPAD7bRdxXwgZzxL0h6Pef14Y4s1KzUqspdgFmF+HxE3JvbsPVhgK30B/47Z/zmiPhiRxVm1tG8RWC2HdJdR8OAP5e7FrNScRCYFUHS7pKOAe4AHgfmlrkks5JxEJht2zWS3gReBq4GbgVGR8Sm8pZlVjryg2nMzLLNWwRmZhnnIDAzyzgHgZlZxjkIzMwyruIuKNt7772jpqam3GWYme1SFi5c+GpEVO/IvBUXBDU1NTQ0NJS7DDOzXYqkF3Z0Xu8aMjPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllnIPAzCzjHARmZhnnIDAzy7iKu7K4vSZPLm0/M7POzlsEZmYZ5yAwM8s4B4GZWcY5CMzMMs5BYGaWcUUFgaTRkpZKapR0aYHp3SXdlE6fL6kmbd9N0kxJiyU9J+my0pZvZmbt1WYQSOoKTAfGAEOA0yQNyet2FrAmIg4CpgFT0/ZxQPeIGAoMA766OSTMzKwyFLNFMAJojIjlEbERmAPU5fWpA2amw7cAoyQJCKCXpCqgJ7AReKMklZuZWUkUEwT9gZU5401pW8E+EdECrAX6kYTCW8BLwIvATyJidf4HSJooqUFSQ3Nz83Z/CTMz23HFBIEKtEWRfUYA7wEfBgYBF0k64H0dI2ZERG1E1FZX79Czl83MbAcVEwRNwMCc8QHAqtb6pLuB+gCrgdOBuyPi3Yh4BfgvoLa9RZuZWekUEwQLgMGSBknqBowH6vP61ANnpMOnAPdHRJDsDjpeiV7ASOAvpSndzMxKoc0gSPf5nw/MA54Dbo6IJZKmSDop7XYd0E9SI3AhsPkU0+nAHsAzJIFyQ0QsKvF3MDOzdijq7qMRMReYm9c2KWd4A8mpovnzrSvUbmZmlcNXFpuZZZyDwMws4xwEZmYZ5yAwM8s4B4GZWcY5CMzMMs5BYGaWcQ4CM7OMcxCYmWWcg8DMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjCsqCCSNlrRUUqOkSwtM7y7ppnT6fEk1afsESU/lvDZJOry0X8HMzNqjzSCQ1JXkkZNjgCHAaZKG5HU7C1gTEQcB04CpABExOyIOj4jDgS8BKyLiqVJ+ATMza59itghGAI0RsTwiNgJzgLq8PnXAzHT4FmCUJOX1OQ24sT3FmplZ6RUTBP2BlTnjTWlbwT7pw+7XAv3y+pxKK0EgaaKkBkkNzc3NxdRtZmYlUkwQ5P/PHiC2p4+kI4H1EfFMoQ+IiBkRURsRtdXV1UWUZGZmpVJMEDQBA3PGBwCrWusjqQroA6zOmT4e7xYyM6tIxQTBAmCwpEGSupH8Ua/P61MPnJEOnwLcHxEBIKkLMI7k2IKZmVWYqrY6RESLpPOBeUBX4PqIWCJpCtAQEfXAdcAsSY0kWwLjc37E0UBTRCwvfflmZtZebQYBQETMBebmtU3KGd5A8r/+QvM+CIzc8RLNzKwj+cpiM7OMcxCYmWWcg8DMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllnIPAzCzjHARmZhnnIDAzyzgHgZlZxjkIzMwyrqggkDRa0lJJjZIuLTC9u6Sb0unzJdXkTPuYpEclLZG0WFKP0pVvZmbt1WYQSOoKTAfGAEOA0yQNyet2FrAmIg4CpgFT03mrgN8B50TEocCxwLslq97MzNqtmC2CEUBjRCyPiI0kzx6uy+tTB8xMh28BRkkScAKwKCKeBoiI1yLivdKUbmZmpVBMEPQHVuaMN6VtBftERAuwFugHHAyEpHmSnpB0SaEPkDRRUoOkhubm5u39DmZm1g7FBIEKtEWRfaqATwET0veTJY16X8eIGRFRGxG11dXVRZRkZmalUkwQNAEDc8YHAKta65MeF+gDrE7bH4qIVyNiPTAX+Hh7izYzs9IpJggWAIMlDZLUDRgP1Of1qQfOSIdPAe6PiADmAR+TtHsaEMcAz5amdDMzK4WqtjpERIuk80n+qHcFro+IJZKmAA0RUQ9cB8yS1EiyJTA+nXeNpJ+ShEkAcyPizg76LmZmtgPaDAKAiJhLslsnt21SzvAGYFwr8/6O5BRSMzOrQL6y2Mws4xwEZmYZ5yAwM8s4B4GZWcY5CMzMMs5BYGaWcQ4CM7OMcxCYmWWcg8DMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllXFFBIGm0pKWSGiVdWmB6d0k3pdPnS6pJ22skvS3pqfR1bWnLNzOz9mrzCWWSugLTgc+QPIx+gaT6iMh99vBZwJqIOEjSeGAqcGo67fmIOLzEdZuZWYkUs0UwAmiMiOURsRGYA9Tl9akDZqbDtwCjJKl0ZZqZWUcpJgj6AytzxpvStoJ9IqIFWAv0S6cNkvSkpIckHVXoAyRNlNQgqaG5uXm7voCZmbVPMUFQ6H/2UWSfl4D9IuII4ELg95J6v69jxIyIqI2I2urq6iJKMjOzUikmCJqAgTnjA4BVrfWRVAX0AVZHxDsR8RpARCwEngcObm/RZmZWOsUEwQJgsKRBkroB44H6vD71wBnp8CnA/RERkqrTg81IOgAYDCwvTelmZlYKbZ41FBEtks4H5gFdgesjYomkKUBDRNQD1wGzJDUCq0nCAuBoYIqkFuA94JyIWN0RX8TMzHZMm0EAEBFzgbl5bZNyhjcA4wrMdytwaztrNDOzDuQri83MMs5BYGaWcQ4CM7OMcxCYmWWcg8DMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllnIPAzCzjHARmZhnnIDAzyzgHgZlZxhUVBJJGS1oqqVHSpQWmd5d0Uzp9vqSavOn7SVon6eLSlG1mZqXSZhCkzxyeDowBhgCnSRqS1+0sYE1EHARMA6bmTZ8G3NX+cs3MrNSK2SIYATRGxPKI2AjMAery+tQBM9PhW4BRkgQg6fMkD6xfUpqSzcyslIoJgv7AypzxprStYJ+IaAHWAv0k9QK+DVyxrQ+QNFFSg6SG5ubmYms3M7MSKCYIVKAtiuxzBTAtItZt6wMiYkZE1EZEbXV1dRElmZlZqVQV0acJGJgzPgBY1UqfJklVQB9gNXAkcIqkq4C+wCZJGyLimnZXbmZmJVFMECwABksaBPwdGA+cntenHjgDeBQ4Bbg/IgI4anMHSZOBdQ4BM7PK0mYQRESLpPOBeUBX4PqIWCJpCtAQEfXAdcAsSY0kWwLjO7JoMzMrnWK2CIiIucDcvLZJOcMbgHFt/IzJO1CfmZl1MF9ZbGaWcQ4CM7OMcxCYmWWcg8DMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnEOAjOzjHMQmJllnIPAzCzjHARmZhnnIDAzyzgHgZlZxhUVBJJGS1oqqVHSpQWmd5d0Uzp9vqSatH2EpKfS19OSTi5t+WZm1l5tBoGkrsB0YAwwBDhN0pC8bmcBayLiIGAaMDVtfwaojYjDgdHAr9JnGpuZWYUoZotgBNAYEcsjYiMwB6jL61MHzEyHbwFGSVJErI+IlrS9BxClKNrMzEqnmCDoD6zMGW9K2wr2Sf/wrwX6AUg6UtISYDFwTk4wmJlZBSgmCFSgLf9/9q32iYj5EXEoMBy4TFKP932ANFFSg6SG5ubmIkoyM7NSKSYImoCBOeMDgFWt9UmPAfQBVud2iIjngLeAj+Z/QETMiIjaiKitrq4uvnozM2u3YoJgATBY0iBJ3YDxQH1en3rgjHT4FOD+iIh0nioASfsDhwArSlK5mZmVRJtBkO7TPx+YBzwH3BwRSyRNkXRS2u06oJ+kRuBCYPMppp8Cnpb0FHA7cF5EvFrqLwHA7NlQU8P/uaIL37y6hqGLZ3fIx5iZdTZFncoZEXOBuXltk3KGNwDjCsw3C5jVzhrbNns2TJwI69cjoO/aF/jcHycCsHjohA7/eDOzXVnnOKf/8sth/fqtmrq9u55R911esUEweXLH9DUz216d4xYTL75YsLnP2sLtZma2RecIgv32K9i8tk/hdjMz26JzBMGVV8Luu2/VtHG33blv1JVlKsjMbNfROYJgwgSYMQP2359AvN5nf/74uRkVe3zAzKySdI6DxZCEwYQJXDG53IWYme1aOscWgZmZ7TAHgZlZxjkIzMwyzkFgZpZxDgIzs4xzEJiZZZyDwMws4xwEZmYZ5yAwM8s4B4GZWcYVFQSSRktaKqlR0qUFpneXdFM6fb6kmrT9M5IWSlqcvh9f2vLNzKy92gwCSV2B6cAYYAhwmqQhed3OAtZExEHANGBq2v4q8LmIGEryTOOOf1qZmZltl2JuOjcCaIyI5QCS5gB1wLM5feqAyenwLcA1khQRT+b0WQL0kNQ9It5pd+U7w7vvwksvwauvwjvvJK9u3aB3b+jbFz74QajqPPftM7NsKuavWH9gZc54E3Bka30iokXSWqAfyRbBZv8MPFkoBCRNBCYC7NfKQ2Y63Ouvw4MPwoIF8MQTsHhxEgKbNrU+z267waBBcMghcOSR8IlPJO+9eu20ss3M2quYIFCBttiePpIOJdlddEKhD4iIGcAMgNra2vyf3XFWrIAbb4Q//Qnmz4f33oOuXeHQQ+H445M/8gMGwD77QI8eydbAxo2wdi2sWZPMv2wZLFkCf/xj8jO7dUvmPekkOPnkZKvBzKyCFRMETcDAnPEBwKpW+jRJqgL6AKsBJA0Abge+HBHPt7vidqp6922GLv49HPUbeOSRpLG2Fi67DE44AYYPT/7ob681a+Cxx+Dee+GOO+C88+CCC+DEE+HMM2Hs2GQLwsyswhQTBAuAwZIGAX8HxgOn5/WpJzkY/ChwCnB/RISkvsCdwGUR8V+lK3v77fnG3xm+4BfULvwVu7/9GnzkI8kjLk8/HWpq2v8Be+0FY8Ykr5/8BJ59FmbNgpkzob4eBg6Eb30Lzj4b9tyz/Z9nZlYiimh7T4ykscDVQFfg+oi4UtIUoCEi6iX1IDkj6AiSLYHxEbFc0neAy4BlOT/uhIh4pbXPqq2tjYaGhh3+QpMnbz2+x5svcdSff8CwJ2bQZVMLSw+p47GR3+SF/Y4CFdqj1T75n09LC9x5J0ybBg89lBxk/vrX+eHGi3ine++O//ydrNjPL3edZp2NpIURUbsj8xZ1yktEzAXm5rVNyhneAIwrMN/3ge/vSGHt1fPt1Rz18JUMb/gFXTa18NTh/5s/f+oyXt9r0M4tpKoK6uqS1/z5MHUqTJnCN3pO5+GjLqdh+Lm0VO3ArigzsxLpdOc+dtnUwrCFMzjuge/SY8PrLPrYl3jomEms2euAcpeWnFF0223Q0MBLp/4bo++5kJHzr2beZ6fx3D+c3CFbKGZmbelcQfDgg0yc8Q0++PIi/lZzHHeN/hmv7Du03FW9X20ts750D4OW38dn77mQU2/+ZxoP/Cx3jfk5r/U7uNzVmVnGdJ57Da1fD+PG0WPD69w07hZmfvm+ygyBHH87YBQzJi7krtE/Y0DTo5z7y6Ec98AkurbsGtfbmVnn0Hm2CHbfHe6+m2tuG0LLbj3LXU3RNnWpYv6RX2fJoV/ghHsu5piHv8c//OV27qi7gVUf3qHjPmZm26XzbBEADBu2S4VArnV7fJDb/ul3zD7tT/R8ezVn/3okx993ubcOzKzDda4g6ASWHXwivzhvCU8f9mWOfuQHfHXGMPZ9eVG5yzKzTsxBUIE29OjLHXXX87vT59Lz7df4l/8YzpGP/QyKuObDzGx7OQgqWOPgMfzynEU8f+BnGTPvm0z4/Vh6rXu53GWZWSfjIKhw63tVc+P4O7hz7HRqVjzIeb8cyuBlc9ue0cysSA6CXYHEguHnMeNfGnhzzw8x4fcnMvqub/hAspmVhINgF9K8z6H8+uz5PHbkNxj5+M85+7qR9Ht1abnLMrNdnINgF9NS1YO7R1/N70/7I33WruSrM4Zx2FMzfSDZzHaYg2AX9deD/5FfnvM0qz5cy8l3fIV/uv1LdHvnzXKXZWa7IAfBLuzN3v2Z+eX7uP/YKXz0mRs551dHQDtu4W1m2eQg2MVFl648fMx3+c1XHqLrexvhk5+Ef//3bT9r2cwsh4Ogk3hxv09x7TlPJY/GvPji5P2VVp//Y2b2P4oKAkmjJS2V1Cjp0gLTu0u6KZ0+X1JN2t5P0gOS1km6prSlW763e34ged7B9OnwwANw2GFw333lLsvMKlybQSCpKzAdGAMMAU6TNCSv21nAmog4CJgGTE3bNwDfBS4uWcW2bRKcdx48/njyHOXPfAb+7d/g3XfLXZmZVahitghGAI0RsTwiNgJzgLq8PnXAzHT4FmCUJEXEWxHxCEkg2M70sY/BggVw5pnwwx/C0UfDihXlrsrMKlAxQdAfWJkz3pS2FewTES3AWqBfsUVImiipQVJDc3NzsbNZW3r1gl//GubMgWefhcMPhz/8odxVmVmFKSYICj1IN//qpWL6tCoiZkREbUTUVldXFzubFevUU+HJJ+GQQ+ALX4CvfjV5opuZGcUFQRMwMGd8ALCqtT6SqoA+wOpSFGglcsAB8MgjcMklMGMGjBgBzzxT7qrMrAIUEwQLgMGSBknqBowH6vP61ANnpMOnAPdH+J4HFWe33WDqVJg3D5qbYfhwuPZa357CLOPaDIJ0n//5wDzgOeDmiFgiaYqkk9Ju1wH9JDUCFwL/c4qppBXAT4GvSGoqcMaR7WwnnACLFiUHkM89Fz7/eViVv5FnZllR1MPrI2IuMDevbVLO8AZgXCvz1rSjPuso++4Ld90F06bBd74DQ4YkVySfeWZyCqqZZYavLM6yLl3goouSrYPDDoOzz4ZPfxqWLy93ZWa2EzkIDAYPTq5Evvba5NqDj34UrroKNm4sd2VmthM4CCzRpUtyWumzzyZXI3/72zB0KNx9d7krM7MO5iCwrQ0YAHfcAXPnJmcTjRmTHEz27iKzTstBYIWNGQOLF8OPfgT33pscTL7oInj11XJXZmYl5iCw1nXvnuwiWroUTj8drr46uTDte9+DdevKXZ2ZlYiDwNrWvz9cf32yhTBqFEyaBAceCD/+Mbzpx2Oa7eocBFa8IUPg9tvh0UeTA8mXXAL77Qff/a53GZntwhwEtv1GjkyOG8yfD8cdB9//fhIIX/sa+7zi+xeZ7WocBLbjRoxInoj27LPJHU6vu47zfjmUr/zmGA595qbkGcpmVvEcBNZ+H/kI3HADNDVxz6evos/alYy7dTwX/nQAY+ZeQP+m+b6xnVkFK+peQ2ZF2Xtv/vt//SuPfvIiDmq8m8OfuoFhT/wHRy64htV7HciioRP4y0dO5v/texiFH2FhZuXgILCSC3Vh2eCxLBs8lu4b1jLkuVsZung2xzz8PY59eAprew+EVz8HJ52U3AG1Z89yl2yWaQ4C61Dv9OjDk0ecyZNHnEmvdS9z8LI7OWRpPX1+8xv4xS+gWzf4xCfg2GOTA88jRybXL5jZTuNjBLbTvLXHvjx5xJnMGf+fyemmd94JF1yQXIswZUoSBr17Jwehv/a15LjD4sXQ0lLu0s06NW8RWHn07AljxyYvgNdfh4cfTh6n2dAAs2YlWwyQbCEcfHByUHrz65BDoKYG+vYt21cw6yyKCgJJo4GfAV2BX0fEj/Kmdwd+CwwDXgNOjYgV6bTLgLOA94CvR8S8klVvnUffvskxg5PSh95t2gTLliW3xX76aXjuuSQg/vCHrc9A6t0b9t9/y2vAANhnny2v6urkvVev8nwvs2LMng2XXw4vvphck3PllTBhwk77+DaDQFJXYDrwGZKH1C+QVB8Rz+Z0OwtYExEHSRoPTAVOTR9LOR44FPgwcK+kgyPivVJ/EetkunRJ/td/yCHwxS9uaX/7bfjrX5P7H734IrzwwpbXI48kWxaF9OyZhE3v3smrT5+t3/fYA3r0aPu1225QVQVdu27fu7TlyW+bh7c1btkxezZMnAjr1yfjL7yQjMNOC4NitghGAI0RsRxA0hygDsgNgjpgcjp8C3CNJKXtcyLiHeBv6TONRwCPlqZ8y5yePZOnqR12WOHpb70Fzc3J65VXtryam2HtWnjjjS2vl17aMvzmm8lWSCXZnuDIHe/MOuN3XL/+/dfZrF+fbCFUUBD0B1bmjDcBR7bWJyJaJK0F+qXtj+XN2z//AyRNBNIIZJ2kpUVVX9jeQNlufHPFFduc3OG1tfH5bdlpy24H6yzrv20RSlvf5j8OpbkYL1vLrvQ6rL5hyS7193vhBRZKC4v8MXsD++9oDcUEQaEIzv/NbK1PMfMSETOAGUXU0iZJDRFRW4qfVWqVXBu4vvaq5PoquTZwfe2V1lezo/MXc/poEzAwZ3wAsKq1PpKqgD7A6iLnNTOzMiomCBYAgyUNktSN5OBvfV6feuCMdPgU4P6IiLR9vKTukgYBg4HHS1O6mZmVQpu7htJ9/ucD80hOH70+IpZImgI0REQ9cB0wKz0YvJokLEj73UxyYLkF+NpOOGOoJLuYOkgl1waur70qub5Krg1cX3u1qz6F7wppZpZpvsWEmVnGOQjMzDJulwoCSddLekXSMzltP5b0F0mLJN0uqW/OtMskNUpaKumz5agvZ9rFkkLS3um4JP08rW+RpI+Xqz5JF6TLaImkq3Lay778JB0u6TFJT0lqkDQibd+py0/SQEkPSHouXU7fSNs/IOn/SlqWvu9VYfVVxPrRWn0508u2fmyrtkpYN7bxb1u6dSMidpkXcDTwceCZnLYTgKp0eCowNR0eAjwNdAcGAc8DXXd2fWn7QJKD7S8Ae6dtY4G7SK61GAnML9PyOw64F+ieju9TScsPuAcYk7PMHizH8gM+BHw8Hd4T+Gu6jK4CLk3bL835/auU+ipi/WitvnS8rOvHNpZdRawb26ivZOvGLrVFEBEPk5yVlNt2T0Rsvk/xYyTXKkDO7S0i4m/A5ttb7NT6UtOAS9j6Yro64LeReAzoK+lDZajvXOBHkdwGhIh4Jae+Slh+AfROh/uw5TqUnbr8IuKliHgiHX4TeI7kKvk6YGbabSbw+Uqqr1LWj20sPyjz+rGN2ipi3dhGfSVbN3apICjCmSRJCIVvjfG+21t0NEknAX+PiKfzJlVEfcDBwFGS5kt6SNLwtL1S6vsm8GNJK4GfAJel7WWrT1INcAQwH9g3Il6CZIUF9qmw+nJVxPqRW1+lrR95y67i1o28+kq2bnSaIJB0Ocm1CrM3NxXotlPPlZW0O3A5MKnQ5AJt5TiXtwrYi2QT8l+BmyWJyqnvXOBbETEQ+BbJNStQpvok7QHcCnwzIt7YVtcCbWWrr1LWj9z60noqZv0osOwqat0oUF/J1o1OEQSSzgD+EZgQ6U4yKuP2FgeS7EN8WtKKtIYnJH2wQuojreO2dDPycWATyQ2sKqW+M4Db0uE/sGUTfKfXJ2k3khVxdkRsrunlzZvd6fvm3QeVUl/FrB8F6quY9aOVZVcx60Yr9ZVu3eioAxwd9QJq2Ppg4miSK5er8/odytYHdJbTwQc7C9WXN20FWw6GncjWB3QeL9PyOweYkg4fTLJJqUpZfiT7Q49Nh0cBC8ux/NLP+S1wdV77j9n6YPFVFVZfRawfrdWX16cs68c2ll1FrBvbqK9k60aH/WJ20AK5EXgJeJck9c4iOVCzEngqfV2b0/9ykiP6S0mPru/s+vKm5/6ii+SBP88Di4HaMi2/bsDvgGeAJ4DjK2n5AZ8CFqYr3nxgWDmWX1pHAItyftfGktxu/T5gWfr+gQqrryLWj9bqq4T1YxvLriLWjW3UV7J1w7eYMDPLuE5xjMDMzHacg8DMLOMcBGZmGecgMDPLOAeBmVnGOQjMzDLOQWBmlnH/H0MErhTbsSZiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.hist(X, density=True, color='b', alpha=0.5)\n",
    "plt.plot(x_sub, rv.pdf(x_sub), color='r') # plot for a<x<b\n",
    "plt.plot(a, rv.pdf(a)*1e-1, 'o', color='r') # point mass at a (re-scale by 1/10 for plotting)\n",
    "plt.plot(b, rv.pdf(b)*1e-1, 'o', color='r') # point mass at b (re-scale by 1/10 for plotting)\n",
    "plt.title('PDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll use the `dcnorm.fit` function to fit the sample of data to a DCNORM distribution using Maximum Likelihood estimation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_fit, s_fit = dcnorm.fit(X) # fit parameters to data\n",
    "rv_fit = dcnorm(m_fit, s_fit) # create new distribution object with fitted parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and recreate the previous plot, but also include the fitted distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXwV1fnH8c9DwhJAATEiq4mKKIsEDeAGuItLQatWkKoIikrd6tJiqRS1WhdaN7T+qLj8LBYVt6ggtgIiVYGgICCigIBBfhhW2aJAnt8fM+DlckMu5Ca55H7fr9d93ZkzZ+48d2Duk5kzc465OyIikrqqVXYAIiJSuZQIRERSnBKBiEiKUyIQEUlxSgQiIilOiUBEJMUpEYiIpDglAkl5ZrbYzDab2QYzW2Fmz5pZXTObZGZFZrbezH4wsxlmNsjMakasO9TMtoTrbn/9rjK/j8ieUiIQCfzC3esCxwAdgT+G5de7+35AY+BWoBcw1swsYt2X3L1uxOvBCo1cpIyUCEQiuPsyYBzQNqp8o7tPAnoAxwPnVnx0IuVDiUAkgpk1B84BPou13N2XAvlAl4qMS6Q8KRGIBN4ws7XAFOAD4L7d1P0OOCBi/ldmtjbi1aQ8AxVJtPTKDkAkSZzv7v+JLNi5GWAnTYGPIuZfdvdfl1dgIuVNZwQieyC8dHQs8GFlxyKSKEoEInEws9pm1g14E5gGjK3kkEQSRolAZPeGm9l6YAXwCPAq0N3diys3LJHEMQ1MIyKS2nRGICKS4pQIRERSnBKBiEiKUyIQEUlxSfdA2YEHHuhZWVmVHYaIyD5lxowZK909c2/WTbpEkJWVRX5+fmWHISKyTzGzJXu7ri4NiYikOCUCEZEUp0QgIpLikq6NQESqhi1btlBQUEBRUVFlh1Kl1KpVi2bNmlG9evWEfaYSgYiUi4KCAvbbbz+ysrJ216W37AF3Z9WqVRQUFJCdnZ2wz9WlIREpF0VFRTRs2FBJIIHMjIYNGyb8LEuJQETKjZJA4pXHPlUiEBFJcUoEIlIlrV27lieffLLCtpeVlcXKlSsBOOGEE3Zb9777djckdsWrco3FQ4cmtp6I7Ju2J4KBAwfusmzbtm2kpaWV27Y/+uij3S6/7777+MMf/lBu299TOiMQkSpp0KBBLFy4kJycHG6//XYmTZrEKaecwqWXXkq7du1YvHgxbdu23VF/2LBhDA3/Qly4cCHdu3fn2GOPpUuXLnz55Ze7fP6qVas488wz6dChA9dccw2Rg3zVrVsXgOXLl9O1a1dycnJo27YtH374IYMGDWLz5s3k5OTQp0+f8t0JcapyZwQikoRuvhlmzkzsZ+bkwCOPlLj4/vvvZ86cOcwMtztp0iSmTZvGnDlzyM7OZvHixSWuO2DAAJ566ilatmzJ1KlTGThwIBMmTNipzl133cVJJ53EkCFDeOeddxgxYsQun/Piiy9y1llnMXjwYLZt28amTZvo0qULw4cP3xFXMlAiEJGU0alTp1Lvv9+wYQMfffQRF1988Y6yH3/8cZd6kydP5rXXXgPg3HPPpUGDBrvU6dixI/369WPLli2cf/755OTklPEblA8lAhEpf7v5y70i1alTZ8d0eno6xcXFO+a335tfXFxM/fr14/qLvbRbObt27crkyZN55513uOyyy7j99tu5/PLL9zL68qM2AhGpkvbbbz/Wr19f4vJGjRrx/fffs2rVKn788UfefvttAPbff3+ys7N55ZVXgOBp3lmzZu2yfteuXRk1ahQA48aNY82aNbvUWbJkCQcddBBXX301/fv359NPPwWgevXqbNmypczfMVHiSgRm1t3M5pvZAjMbFGN5TTN7KVw+1cyywvLqZva8mc02s3lmdkdiwxcRia1hw4aceOKJtG3blttvv32X5dWrV2fIkCF07tyZ8847jyOPPHLHslGjRjFy5Ejat29PmzZtePPNN3dZ/09/+hOTJ0/mmGOO4b333qNFixa71Jk0aRI5OTl06NCBV199lZtuugkI2iCOPvroHY3F55xzDt99912ivvoes8iW7pgVzNKAr4AzgAJgOtDb3b+IqDMQONrdrzWzXsAF7n6JmV0K9HD3XmZWG/gCONndF5e0vdzcXC/LwDS6fVQkOcybN4+jjjqqssOokmLtWzOb4e65e/N58ZwRdAIWuPsid/8JGA30jKrTE3g+nB4DnGbBxTMH6phZOpAB/AT8sDeBiohI+YgnETQFvo2YLwjLYtZx963AOqAhQVLYCCwHlgLD3H119AbMbICZ5ZtZfmFh4R5/CRER2XvxJIJYzeLR15NKqtMJ2AY0AbKBW83s0F0quo9w91x3z83M3Kuxl0VEZC/FkwgKgOYR882A6FaNHXXCy0D1gNXApcC77r7F3b8H/gvs1TUsEREpH/EkgulASzPLNrMaQC8gL6pOHnBFOH0RMMGDVuilwKkWqAMcB+z6rLaIiFSaUhNBeM3/emA8MA942d3nmtndZtYjrDYSaGhmC4BbgO23mD4B1AXmECSUZ9398wR/BxERKYO4niNw97HufoS7H+bu94ZlQ9w9L5wucveL3f1wd+/k7ovC8g1heRt3b+3uD5XfVxER2dljjz3GUUcdRZ8+fcjLy+P+++8H4I033uCLL3bcAc9zzz23x/fxR3daV5K+ffsyZswYAK666qqdthttb+JIBHUxISJV1pNPPsm4ceN29C/Uo0dwEeONN97gvPPOo3Xr1kDwA9y2bVuaNGlSrvE8/fTTu11eUXFEUyIQkSrp2muvZdGiRfTo0YN+/frRoEED8vPzufTSS8nLy+ODDz7gz3/+M7179yY/P58+ffqQkZHBxx9/zBdffMEtt9zChg0bOPDAA3nuuedo3LgxM2bMoF+/ftSuXZuTTjop5nbdnRtuuIEJEyaQnZ29U/fUJ598MsOGDaNDhw7079+f/Px8zIx+/frRvHnzXeLIyMiokH2lRCAi5a4SeqHmqaee4t1332XixIk7fswhGD2sR48enHfeeVx00UVA0FfQsGHDyM3NZcuWLdxwww28+eabZGZm8tJLLzF48GCeeeYZrrzySh5//HG6desWs9sKgNdff5358+cze/ZsVqxYQevWrenXr99OdWbOnMmyZcuYM2cOEAyiU79+fYYPH74jjoqkRCAiEmH+/PnMmTOHM844AwhGM2vcuDHr1q1j7dq1dOvWDYDLLruMcePG7bL+5MmT6d27N2lpaTRp0oRTTz11lzqHHnooixYt4oYbbuDcc8/lzDPPLN8vVQolAhEpd0nSC3Vc3J02bdrw8ccf71S+du3aUrud3q60eg0aNGDWrFmMHz+eJ554gpdffplnnnlmr2MuK3VDLSIpJ7qL6sj5Vq1aUVhYuCMRbNmyhblz51K/fn3q1avHlClTAHZ0QR2ta9eujB49mm3btrF8+XImTpy4S52VK1dSXFzMhRdeyD333LOje+rSus4uL0oEIpJyevXqxUMPPUSHDh1YuHAhffv25dprryUnJ4dt27YxZswYfv/739O+fXtycnJ2DEb/7LPP8pvf/Ibjjz++xIbcCy64gJYtW9KuXTuuu+66HZeSIi1btoyTTz6ZnJwc+vbty1/+8heAneLYvHkzQ4YMIS8v+vndxCu1G+qKpm6oRaoGdUNdfiqjG2oREanClAhERFKcEoGISIpTIhARSXFKBCIiKU6JQEQkxenJYhGpEIm+ZXtPP2/o0KHUrVuX2267LebyN954gyOOOGJHj6SpRGcEIiLsOkZBKokrEZhZdzObb2YLzGxQjOU1zeylcPlUM8sKy/uY2cyIV7GZ5ST2K4iIxHbvvffSqlUrTj/9dObPnw/AP/7xDzp27Ej79u258MIL2bRpEx999BF5eXncfvvt5OTksHDhwpj1qqpSE4GZpREMOXk20BrobWbR5079gTXufjjwMPAAgLuPcvccd88BLgMWu3uCO6MVEdnVjBkzGD16NJ999hmvvfYa06dPB+CXv/wl06dPZ9asWRx11FGMHDlyR9fUDz30EDNnzuSwww6LWa+qiqeNoBOwYPvwk2Y2GugJRJ5D9QSGhtNjgOFmZr5z/xW9gX+VOWIRkTh8+OGHXHDBBdSuXRv4eXSyOXPm8Mc//pG1a9eyYcMGzjrrrJjrx1uvKojn0lBT4NuI+YKwLGadcLD7dUDDqDqXUEIiMLMBZpZvZvmFhYXxxC0iUqpY3UH37duX4cOHM3v2bP70pz9RVFQUc91461UF8SSCWB1rR/dUt9s6ZtYZ2OTuc2JtwN1HuHuuu+dmZmbGEZKIyO517dqV119/nc2bN7N+/XreeustANavX0/jxo3ZsmXLTl1JR3cBXVK9qiieS0MFQPOI+WbAdyXUKTCzdKAesDpieS90WUgkpVV0j7/HHHMMl1xyCTk5ORxyyCF06dIFgHvuuYfOnTtzyCGH0K5dux0//r169eLqq6/mscceY8yYMSXWq4pK7YY6/GH/CjgNWAZMBy5197kRdX4DtHP3a82sF/BLd/9VuKwasBTour2dYXfUDbVI1aBuqMtPoruhLvWMwN23mtn1wHggDXjG3eea2d1AvrvnASOBF8xsAcGZQK+Ij+gKFMSTBEREpOLF9WSxu48FxkaVDYmYLgIuLmHdScBxex+iiIiUJz1ZLCLlJtlGQKwKymOfKhGISLmoVasWq1atUjJIIHdn1apV1KpVK6Gfq07nRKRcNGvWjIKCAvRsUGLVqlWLZs2aJfQzlQhEpFxUr16d7Ozsyg5D4qBLQyIiKU6JQEQkxSkRiIikOCUCEZEUp0QgIpLilAhERFKcEoGISIpTIhARSXFKBCIiKU6JQEQkxSkRiIikuLgSgZl1N7P5ZrbAzAbFWF7TzF4Kl081s6yIZUeb2cdmNtfMZptZYrvNExGRMik1EZhZGvAEcDbQGuhtZq2jqvUH1rj74cDDwAPhuunAP4Fr3b0NcDKwJWHRi4hImcVzRtAJWODui9z9J2A00DOqTk/g+XB6DHCamRlwJvC5u88CcPdV7r4tMaGLiEgixJMImgLfRswXhGUx67j7VmAd0BA4AnAzG29mn5rZ72JtwMwGmFm+meWr73IRkYoVTyKwGGXRQw6VVCcdOAnoE75fYGan7VLRfYS757p7bmZmZhwhiYhIosSTCAqA5hHzzYDvSqoTtgvUA1aH5R+4+0p33wSMBY4pa9AiIpI48SSC6UBLM8s2sxpALyAvqk4ecEU4fREwwYOBSscDR5tZ7TBBdAO+SEzoIiKSCKUOVenuW83seoIf9TTgGXefa2Z3A/nungeMBF4wswUEZwK9wnXXmNnfCJKJA2Pd/Z1y+i4iIrIX4hqz2N3HElzWiSwbEjFdBFxcwrr/JLiFVEREkpCeLBYRSXFKBCIiKU6JQEQkxSkRiIikOCUCEZEUp0QgIpLilAhERFKcEoGISIpTIhARSXFKBCIiKU6JQEQkxSkRiIikOCUCEZEUp0QgIpLilAhERFKcEoGISIqLKxGYWXczm29mC8xsUIzlNc3spXD5VDPLCsuzzGyzmc0MX08lNnwRESmrUkcoM7M04AngDILB6KebWZ67R4493B9Y4+6Hm1kv4AHgknDZQnfPSXDcIiKSIPGcEXQCFrj7Inf/CRgN9Iyq0xN4PpweA5xmZpa4MEVEpLzEkwiaAt9GzBeEZTHruPtWYB3QMFyWbWafmdkHZtYl1gbMbICZ5ZtZfmFh4R59ARERKZt4EkGsv+w9zjrLgRbu3gG4BXjRzPbfpaL7CHfPdffczMzMOEISEZFEiScRFADNI+abAd+VVMfM0oF6wGp3/9HdVwG4+wxgIXBEWYMWEZHEiScRTAdamlm2mdUAegF5UXXygCvC6YuACe7uZpYZNjZjZocCLYFFiQldREQSodS7htx9q5ldD4wH0oBn3H2umd0N5Lt7HjASeMHMFgCrCZIFQFfgbjPbCmwDrnX31eXxRUREZO+UmggA3H0sMDaqbEjEdBFwcYz1XgVeLWOMIiJSjvRksYhIilMiEBFJcUoEIiIpTolARCTFKRGIiKQ4JQIRkRSnRCAikuKUCEREUpwSgYhIilMiEBFJcUoEIiIpTolARCTFKRGIiKQ4JQIRkRSnRCAikuKUCEREUlxcicDMupvZfDNbYGaDYiyvaWYvhcunmllW1PIWZrbBzG5LTNgiIpIopSaCcMzhJ4CzgdZAbzNrHVWtP7DG3Q8HHgYeiFr+MDCu7OGKiEiixXNG0AlY4O6L3P0nYDTQM6pOT+D5cHoMcJqZGYCZnU8wYP3cxIQc26hRkJUFd90FjzwCs2eX59ZERKqOeBJBU+DbiPmCsCxmHXffCqwDGppZHeD3wF2724CZDTCzfDPLLywsjDf2HUaNggEDYMmSYH7dOnjrLSUDEZF4xJMILEaZx1nnLuBhd9+wuw24+wh3z3X33MzMzDhC2tngwbBp085lW7bA++/v8UeJiKSc9DjqFADNI+abAd+VUKfAzNKBesBqoDNwkZk9CNQHis2syN2HlznyCEuXxi5fty6RWxERqZriSQTTgZZmlg0sA3oBl0bVyQOuAD4GLgImuLsDXbZXMLOhwIZEJwGAFi1+viwUqV69RG9JRKTqKfXSUHjN/3pgPDAPeNnd55rZ3WbWI6w2kqBNYAFwC7DLLabl6d57oXbtncuqV4fTTqvIKERE9k3xnBHg7mOBsVFlQyKmi4CLS/mMoXsRX1z69AneBw8Ozgzq1QuSQLt25bVFEZGqI65EsC/o0yd4DR1a2ZGIiOxb1MWEiEiKUyIQEUlxSgQiIilOiUBEJMUpEYiIpDglAhGRFKdEICKS4pQIRERSnBKBiEiKUyIQEUlxSgQiIilOiUBEJMUpEYiIpDglAhGRFBdXIjCz7mY238wWmNkug86YWU0zeylcPtXMssLyTmY2M3zNMrMLEhu+iIiUVamJwMzSgCeAs4HWQG8zax1VrT+wxt0PBx4GHgjL5wC57p4DdAf+JxzTWEREkkQ8ZwSdgAXuvsjdfwJGAz2j6vQEng+nxwCnmZm5+6ZwqEuAWoAnImgREUmceBJBU+DbiPmCsCxmnfCHfx3QEMDMOpvZXGA2cG1EYhARkSQQTyKwGGXRf9mXWMfdp7p7G6AjcIeZ1dplA2YDzCzfzPILCwvjCElERBIlnkRQADSPmG8GfFdSnbANoB6wOrKCu88DNgJtozfg7iPcPdfdczMzM+OPXkREyiyeRDAdaGlm2WZWA+gF5EXVyQOuCKcvAia4u4frpAOY2SFAK2BxQiIXEZGEKDURhNf0rwfGA/OAl919rpndbWY9wmojgYZmtgC4Bdh+i+lJwCwzmwm8Dgx095WJ/hIAjBoFWVn86a5q3PxIFu1mjyqXzYiIVDVx3crp7mOBsVFlQyKmi4CLY6z3AvBCGWMs3ahRMGAAbNqEAfXXLeEXbw0AYHa7PuW+eRGRfVnVuKd/8GDYtGmnohpbNnHa+4OTNhEMHVo+dUVE9lTV6GJi6dKYxfXWxS4XEZGfVY1E0KJFzOJ19WKXi4jIz6pGIrj3Xqhde6ein6rX5v3T7q2kgERE9h1VIxH06QMjRsAhh+AYa+sdwlu/GJG07QMiIsmkajQWQ5AM+vThrqGVHYiIyL6lapwRiIjIXlMiEBFJcUoEIiIpTolARCTFKRGIiKQ4JQIRkRSnRCAikuKUCEREUpwSgYhIilMiEBFJcXElAjPrbmbzzWyBmQ2Ksbymmb0ULp9qZllh+RlmNsPMZofvpyY2fBERKatSE4GZpQFPAGcDrYHeZtY6qlp/YI27Hw48DDwQlq8EfuHu7QjGNC7/0cpERGSPxNPpXCdggbsvAjCz0UBP4IuIOj2BoeH0GGC4mZm7fxZRZy5Qy8xquvuPZY68ImzZAsuXw8qV8OOPwatGDdh/f6hfHw4+GNKrTr99IpKa4vkVawp8GzFfAHQuqY67bzWzdUBDgjOC7S4EPouVBMxsADAAoEUJg8yUu7VrYdIkmD4dPv0UZs8OkkBxccnrVK8O2dnQqhV07gzHHx+816lTYWGLiJRVPInAYpT5ntQxszYEl4vOjLUBdx8BjADIzc2N/uzys3gx/Otf8PbbMHUqbNsGaWnQpg2cemrwI9+sGRx0ENSqFZwN/PQTrFsHa9YE63/9NcydC2+9FXxmjRrBuj16wAUXBGcNIiJJLJ5EUAA0j5hvBnxXQp0CM0sH6gGrAcysGfA6cLm7LyxzxGWUvmUz7Wa/CF2egylTgsLcXLjjDjjzTOjYMfjR31Nr1sAnn8B//gNvvgkDB8INN8C550K/fnDOOcEZhIhIkoknEUwHWppZNrAM6AVcGlUnj6Ax+GPgImCCu7uZ1QfeAe5w9/8mLuw9t98Py+g4/UlyZ/wPtTevgqOOCoa4vPRSyMoq+wYaNICzzw5ew4bBF1/ACy/A889DXh40bw6//S1cdRXst1/ZtycikiDmXvqVGDM7B3gESAOecfd7zexuIN/d88ysFsEdQR0IzgR6ufsiM/sjcAfwdcTHnenu35e0rdzcXM/Pz9/rLzR06M7zddcvp8uH93HspyOoVryV+a168slxN7OkRRewWFe0yiZ6+2zdCu+8Aw8/DB98EDQy33gjf/npVn6suX/5b7+Cxbv9yo5TpKoxsxnunrs368Z1y4u7jwXGRpUNiZguAi6Osd6fgT/vTWBllbF5NV0m30vH/CepVryVmTlX8uFJd7C2QXbFBpKeDj17Bq+pU+GBB+Duu7kp4wkmdxlMfsfr2Jq+F5eiREQSpMrd+1iteCvHzhjBKRPvpFbRWj4/+jI+6DaENQ0OrezQgjuKXnsN8vNZfskf6P7eLRw39RHGn/Uw8468oFzOUERESlO1EsGkSQwYcRMHr/icb7JOYVz3R/m+UbvKjmpXubm8cNl7ZC96n7Peu4VLXr6QBYedxbizH2NVwyMqOzoRSTFVp6+hTZvg4oupVbSWly4ew/OXv5+cSSDCN4eexogBMxjX/VGaFXzMdX9vxykTh5C2dd943k5Eqoaqc0ZQuza8+y7DX2vN1uoZbNkChYWwenVw2//atbB+PWzeHDwgXFQUPDgcqXr14M7R7a/994d69YL23QYNIDMz8XeAFldLZ2rnG5nb5lec+d5tdJt8D0d++Tpv9nyW75rsVbuPiMgeqTKJwB3ue/dY/vUGrFgRJIDIG6K2/7BnZAQ/7LVqBe242y/Luwc3+BQVBa/tz4v9GPHHuRk0bAiNGkHjxtCiBTRpEjyDVlYb6h7Ma7/8J7Pb9uYXbw/gqqePY8qJv+eDbkPYll6z7BsQESlBlUkEZvDss0ECaNQI2rYNHghu2DD44a+5l7+lRUXBGcWqVUGC+f57WLYseJgYgmTStGnwEHLLlkGPFNXKcMHt6yPO5cmBczlr/C10nXIfR85/k1cvfJEVjY7e+w8VEdmNKpMIIHiG6777EvuZ2y8TNWoErSP6XN24EZYuhSVLgvdJk4LX228HDxH36AFnnbV3CaioVn3e7PkMc1tfTM+8flz9j478+/QHmdr5Rt1ZJCIJV6USQY0aFbetOnWCh5OPOiqY37gx6HaoRg145RUYOTJoX7jwQujdG045Zc8vIS1oeTZ/v/Zzeub15+zxN3P4wnd5o+dzbKzbKPFfSERSVtW5a6iS1akDOTnw8stBr9XjxgXPkL3yCpxxRnD56NZb4csv9+xzN9XJ5F+93uSdc54ga/EkBv69HS2/Hlv6iiIicVIiKAfVq0P37kE3QytWwJgxcMIJ8NhjwRlEt27w+edB43RczJjecSAjrs5n/X6N6fPiuXQfd5NuMxWRhFAiKGcZGcHloddeg4ICuP/+oLH59dfhr3+F994LGqPjUXhQG56+aiqfdL6J46Y9xlUjj6Phyvnl+wVEpMpTIqhAjRrB738PX30Fl18Ohx4a9Fz96KPBWUNBQemfsTW9Fu92f4QXe79FvXXfcs2IY2k/8/md75UVEdkDVaqxeF9RrVpwu2l2dnA2MHVqMCja3LnBODjHHRdcQtrdbahfHXEef792Fhe+1ocL3uzLYYv+zdvn/h1QF9cismeUCCpZvXrBeDjdusHMmUFSGDMGDjgATjwRjj665GGR1+/flOcvf58uH97HyR8MpVnBJ3Dh6GCgHRGROOnSUJKoWTPonPT66+FXvwrm33oraGD+5JNghMxYvFoak7vdyXN9PyBt209Bq/Rf/7r7sZZFRCIoESSZatWCy0JXXw2//nVwZjB+PDzyCEyeHDzpHMvSFifx1LUzg6Exb7steP++xPF/RER2iCsRmFl3M5tvZgvMbFCM5TXN7KVw+VQzywrLG5rZRDPbYGbDExt61WYGhx0GffvClVcGbQcTJwYDnf3nP8EDbNE2ZxwQ3J70xBNB5fbt4f33Kzx2Edm3lJoIzCwNeAI4G2gN9Daz1lHV+gNr3P1w4GHggbC8CLgTuC1hEaegFi2CoZWvuSboz+i//w3OEN59F374IaqyGQwcCNOmBV2mnnEG/OEPu3a1KiISiueMoBOwwN0XuftPwGigZ1SdnsDz4fQY4DQzM3ff6O5TCBKClNHBB8NFF8FvfhN0qjdtWtCG8PbbQW+pOzn6aJg+Hfr1g7/8Bbp2DbpTFRGJEk8iaAp8GzFfEJbFrOPuW4F1QMN4gzCzAWaWb2b5hYWF8a6Wsg48MOi+4oYbgm4tZs6Exx+HN96A+ZHPl9WpA08/DaNHBz3y5eQEfV6IiESIJxHE6u4y+umleOqUyN1HuHuuu+dmZmbGu1rKa9AAzjsPbrwxuONo7tygofmSS4IuLHa45BL47DNo1Sq4Jemaa4IR3UREiC8RFADNI+abAd+VVMfM0oF6wOpEBCil23//oMvrm2+GQYOCDu/atw/OGqZNCysdeihMmQK/+x2MGAGdOsGcOZUat4gkh3gSwXSgpZllm1kNoBeQF1UnD7ginL4ImOCuPg8qWp06wXgMS5bA3XcHv/udOwcPrE2eTNAb3gMPBPejFhZCx47w1FPqnkIkxZWaCMJr/tcD44F5wMvuPtfM7jazHmG1kUBDM1sA3ALsuMXUzBYDfwP6mllBjDuOJMEaNIA77wzahh98EGbNCp5c7to16OTOz4GOmMQAAAs2SURBVDgzuHbUtStcdx2cfz58F32SJyKpIq7nCNx9rLsf4e6Hufu9YdkQd88Lp4vc/WJ3P9zdO7n7ooh1s9z9AHev6+7N3P2L8vkqEm2//eD224OE8Pjj8M03wSWkTp3gzU8aUfzOOBg2LMgOrVsHo+no7EAk5ejJ4hSQkRF0XbFwIfzjH8G4zuefDznHVOOlZrey7bPPg0aFq66C00+HRYtK/1ARqTKUCFJIjRrBb/38+fDCC8HAOL16wRHntuTRCybyw8Mjg2cP2rYNrimV1MGRiFQpSgQpKD096Mdozpygp9PGjeHm31aj2ZB+3PyrZSw84bJg4IR27YLHl0WkSlMiSGHVqgWjp02ZEtxm2qMHPPH8frSc8D/0PO7/mLixE3722cF1JF0uEqmylAgECO4k/ec/g1tPBw+GjxY04tRlL9A6cyV/G3ckK4/qArfeCitXVnaoIpJgSgSykyZN4J57YOlSePZZaHB4Q2796X6abl1Mr7914v0WV1J81z2wYUNlhyoiCaJEIDFlZARdYH/0EcyeDddeX5339r+I0ze/xRFDe3PPwcNZOOgfsH59ZYcqImWkRCClatsWHn0Ulv1fGi+8AM06HMSQjYM4/IGrOa7Blzze/W2+n7eqssMUkb2kRCBxy8gI7jaa9On+LF0KD1y/lKK6B3Lj+PNo0roe3Vt8wdeTlumqkcg+RoPXy15p3hx+93gLfvc4zM1byIt3fsG/Pm/LN982xT4o5oiGqzgkpwGtWqdzwAGVHa2I7I4SgZRZmx6HcW+Pw/hz4UpG9HiWJTPXMH5VN957P5P33oeD622mRataHHaYkZUVPNgmIslDiUASxjIPZPlZV1LzzGJuWfAuB027jjkLM/j3utOZNO0Upk3LIK1aMc2aGenpRteuwW2rGRmVHblIalMikIRzq8bXLc/h65bnULNoHdfMe5W/zrqA75Zs4d/FZzCu4FzuvLMdANWrO8cea5x0Epx4IpxwAhx0UCV/AZEUo0Qg5erHWvX4rEM/PuvQjzobVnD81+/Qd/6dZC6dwUebc5iyrRv//eIsHpt+FMOGBf8dmzaFY46BDh1+fm/eHCzWOHgiUmZKBFJhNtZttCMpDP39Zn4xcSK/mDABJval6NMvmMExfJJ2Ep9tPplPp+TwztuNKPbgxrYDDgiG4TzyyJ1fWVlB30kisvd0CEnlyMiAc84JXkCttWs5cfJkTpwyBfIfgvx8Nvo2PudoPkvvyEzvyvyvWvPWzBaM3Fh3x8fUqAGHHBK8srJ2nm7RAg4+GGrWrJyvKLKviCsRmFl34FEgDXja3e+PWl4T+F/gWGAVcIm7Lw6X3QH0B7YBN7r7+IRFL1VH/fpBr3c9wkHvioup8/XXHD99OsfPmgXznoN582DlN6ymPvNpxXxa8SXt+eb7o1hcmM1b/23Mis37x/ho5+CDjUaNoFGjIDk0agQNGwajudWvv+u7zjKkIo0aOIXBI7JYuq0JLdK+494Bi+nz5EkVtv1S/7ubWRrwBHAGwSD1080sL2qksf7AGnc/3Mx6AQ8Al4TDUvYC2gBNgP+Y2RHuvi3RX0SqmGrVoFWr4PXrX/9cvnkzB3z1FcfPn8/xS5fCkkWwZGLQW97SpWzeXMRSWrCEQ1hKC/6Pg1mxthH/90NTVixqwmd+MCuKD+SHbXVL3jZQt+ZP1Ku9hTq1iqmTUUyd2sXUzoA6tZ3atT14zwjGic6oDTVrGjVqGdVrVKNGLaNGzWrUqBW8qtesRo2aRo0awRlMjZpGejqkpRvVqkG1NAum04y0tGC+WjWC6Yj3ksq2t51Ev5dUJsll1MApDPh7BzZRB4Al25ox4O8NgCkVlgzi+bunE7Bg+/CTZjYa6AlEJoKewNBwegww3MwsLB/t7j8C34RjGncCPk5M+JJyMjKC0dTat4+9eONGWhUW0qqwEL7//udX4Uewbh388AP88AOb1xSxeo2x9odqrFmfztpNNVjj9VhLfdbQgLU/1mfdj/XYSB02UodN1GYddVhO7Z3KNlKHYtIqeCeUnVEcMe27fd/bOhKfTZyAR3XysIk6DB6RRZ8nKyaGeBJBU+DbiPkCoHNJddx9q5mtAxqG5Z9Erds0egNmNgAYEM5uMLP5cUUf24FApfWVfNddu11c7rGVsv3SVNi+28s4K/XfNg7JHN9OsUX+VCfJz3Yy7zso1/iOPTZW6ZJtYDZjRpwfciBwyN5GEE8iiHVSGf1/p6Q68ayLu48ARsQRS6nMLN/dcxPxWYmWzLGB4iurZI4vmWMDxVdWYXxZe7t+PJ3OFQDNI+abAd+VVMfM0oF6wOo41xURkUoUTyKYDrQ0s2wzq0HQ+JsXVScPuCKcvgiY4O4elvcys5pmlg20BKYlJnQREUmEUi8Nhdf8rwfGE9w++oy7zzWzu4F8d88DRgIvhI3BqwmSBWG9lwkalrcCv6mAO4YScompnCRzbKD4yiqZ40vm2EDxlVWZ4rPgD3cREUlVGphGRCTFKRGIiKS4fSoRmNkzZva9mc2JKHvIzL40s8/N7HUzqx+x7A4zW2Bm883srMqIL2LZbWbmZnZgOG9m9lgY3+dmdkxlxWdmN4T7aK6ZPRhRXun7z8xyzOwTM5tpZvlm1iksr9D9Z2bNzWyimc0L99NNYfkBZvZvM/s6fG+QZPElxfFRUnwRyyvt+NhdbMlwbOzm3zZxx4a77zMvoCtwDDAnouxMID2cfgB4IJxuDcwCagLZwEIgraLjC8ubEzS2LwEODMvOAcYRPGtxHDC1kvbfKcB/gJrh/EHJtP+A94CzI/bZpMrYf0Bj4Jhwej/gq3AfPQgMCssHRfz/S5b4kuL4KCm+cL5Sj4/d7LukODZ2E1/Cjo196ozA3ScT3JUUWfaeu28NZz8heFYBIrq3cPdvgO3dW1RofKGHgd+x88N0PYH/9cAnQH0za1wJ8V0H3O9BNyC4+/cR8SXD/nNge09y9fj5OZQK3X/uvtzdPw2n1wPzCJ6S7wk8H1Z7Hjg/meJLluNjN/sPKvn42E1sSXFs7Ca+hB0b+1QiiEM/gkwIsbvG2KV7i/JmZj2AZe4+K2pRUsQHHAF0MbOpZvaBmXUMy5MlvpuBh8zsW2AYcEdYXmnxmVkW0AGYCjRy9+UQHLDA9vHVkiW+SElxfETGl2zHR9S+S7pjIyq+hB0bVSYRmNlggmcVRm0vilGtQu+VNbPawGBgSKzFMcoq417edKABwSnk7cDLZmYkT3zXAb919+bAbwmeWYFKis/M6gKvAje7+w+7qxqjrNLiS5bjIzK+MJ6kOT5i7LukOjZixJewY6NKJAIzuwI4D+jj4UUykqN7i8MIriHOMrPFYQyfmtnBSRIfYRyvhaeR04Bigg6skiW+K4DXwulX+PkUvMLjM7PqBAfiKHffHtOK7afd4fv2ywfJEl/SHB8x4kua46OEfZc0x0YJ8SXu2CivBo7yegFZ7NyY2J3gyeXMqHpt2LlBZxHl3NgZK76oZYv5uTHsXHZu0JlWSfvvWuDucPoIglNKS5b9R3A99ORw+jRgRmXsv3A7/ws8ElX+EDs3Fj+YZPElxfFRUnxRdSrl+NjNvkuKY2M38SXs2Ci3/5jltEP+BSwHthBkvf4EDTXfAjPD11MR9QcTtOjPJ2xdr+j4opZH/kc3ggF/FgKzgdxK2n81gH8Cc4BPgVOTaf8BJwEzwgNvKnBsZey/MA4HPo/4v3YOQXfr7wNfh+8HJFl8SXF8lBRfMhwfu9l3SXFs7Ca+hB0b6mJCRCTFVYk2AhER2XtKBCIiKU6JQEQkxSkRiIikOCUCEZEUp0QgIpLilAhERFLc/wP6RfnWoxDgugAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.hist(X, density=True, color='b', alpha=0.5, label='data')\n",
    "plt.plot(x_sub, rv.pdf(x_sub), color='r', label='true dist.') # plot for a<x<b\n",
    "plt.plot(a, rv.pdf(a)*1e-1, 'o', color='r') # point mass at a (re-scale by 1/10 for plotting)\n",
    "plt.plot(b, rv.pdf(b)*1e-1, 'o', color='r') # point mass at b (re-scale by 1/10 for plotting)\n",
    "\n",
    "plt.plot(x_sub, rv_fit.pdf(x_sub), color='b', label='fitted dist.') # plot for a<x<b\n",
    "plt.plot(a, rv_fit.pdf(a)*1e-1, 'o', color='b') # point mass at a (re-scale by 1/10 for plotting)\n",
    "plt.plot(b, rv_fit.pdf(b)*1e-1, 'o', color='b') # point mass at b (re-scale by 1/10 for plotting)\n",
    "\n",
    "plt.legend()\n",
    "plt.title('PDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can make an analogous plot for the CDF's as well, making used of the `dcnorm.ecdf` function to plot the empirical CDF for the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3hUxfrA8e9LCCVIkyYQSlRUEDVoAAsiFhQEUfBeBZErgmABRKx4QaTotcBPvShFkBIURQTBqHDFAmIFAgYMIApYCKCG0GuS3fn9MRuyJJtkCdv3/TzPPrvnnDm7bw6cd2fnzJkRYwxKKaUiX5lgB6CUUiowNOErpVSU0ISvlFJRQhO+UkpFCU34SikVJTThK6VUlNCEr5RSUUITvooqInKHiKSKyEER2Skii0WkjYiMFJEcETngevwsIq+JSF23fduJiNO1b97jw2D+PUqdDE34KmqIyMPAK8B/gDpAQ2AicLOryLvGmMrA6UBX4AxgtXvSB3YYY05ze9wUuL9AqVOjCV9FBRGpCowGBhhj3jfGHDLG5BhjPjTGPOZe1rV+PXA7kAk8EoSQlfI5TfgqWlwGVAAWeLuDMcYBfABc6a+glAokTfgqWtQAdhljck9yvx3YJp489URkr9vjNt+FqJR/lQ12AEoFSBZQU0TKnmTSrw/sdlveYYyJ921oSgWG1vBVtPgOOArc4u0OIlIGuAn4yl9BKRVIWsNXUcEYs09ERgATRCQXWALkANcBVwOH88qKSCxwNjAS21PnpYAHrJQfaA1fRQ1jzEvAw8BwbO+bbcBAYKGryO0ichDYC6Rgm4EuMcbsCEK4Svmc6AQoSikVHbSGr5RSUUITvlJKRQlN+EopFSU04SulVJQIWrfMmjVrmsaNGwfr45VSKiytXr16lzGmVmn2DVrCb9y4MampqcH6eKWUCksi8ntp99UmHaWUihKa8JVSKkpowldKqSihCV8ppaKEJnyllIoSJSZ8EZkuIn+LSHoR20VExovIZhFZJyIX+z5MpZRSp8qbGv5MoEMx2zsCTVyP/sCkUw9LKaWUr5XYD98Ys1xEGhdT5GZglrHDbn4vItVEpK4xZqePYlRhau9eSEuDoA7IagySk03MkYOUyTlGmdxsJCfb9XyMMsdf5z+LIwdxOsE4EacT43SSnS0czRZycoRcB+TmgsMh5OSC02HIdQi5uUKuU+w2U4YcZwy5zjI4nYITwWkEjMFpXK8RnAacRjBGMAhOJxhwK5P/2iAYV3kPf6aHv93j4fB8mDyu9e5zjKdynj7DQ9zRqtPdtblhSLOAf64vbryqjx1XPE+Ga12hhC8i/bG/AmjYsKEPPlqFstWr4auvQHxwnpdx5FB5/3YqH9hB3OFd9nEk6/jriod3Uf7YAcplH6Rctn0ue+wQB7LLs8vUIJNa7KUa+6jKXqodf+Qt76cKhzmdI1TkMHEnPB8h7tT/AAWA4Ax2CCGhzidfc8OQwH+uLxK+p9PZY4XBGDMFmAKQlJSkA/FHuNxcKFcO/v1vL3fYtQt++gk2bbLPv/4K27bZx59/nlC9zKEs22jAb7HnsKnS+fxeriW/0Zg/nbX4K6cGf+VUIzO7Mg4TU+THVa6YQ7XTcqlW2UHlSk4qxxlqV4S4ilCxorHPFY4RF3eMihUhLs5QvoIQGyvElst/lI0VYsuXOXFdWYiNtY+YGIgpK4hAmTKGMmXyXtuHlJH81wJlYiT/dRm7XLD8CV+iroXj69w2evqyLeoL2NuyJ/OehWk/EattUD7VFwk/A2jgthwP6AxBCofDJrtCjIGtWyE1Nf/x44+QlZVfpnx5SEhg9xnN2NDiFjaUac7Gowls2H0GP+2oTMbf5XA6xU5SuNcmwnr17KNhHUiqDXXq5D9q14Zq1fIfVapATEwsEBugo6FU8Pki4acAA0VkDtAa2Kft9wrA6XQlfGNg40b44gtYuhSWLYPdu22h8uUhMRG6dWN3w0RSpSWpe85i1ebqrEoVtv+U/35xcdC0KVx5LZx1FjRubB+NGkF8vP01oZQqWokJX0TeAdoBNUUkA3gaV7XIGDMZWATcCGzGTgR9t7+CVWEkO5vTU5fS/tOF8HIK7HD96GvUCLp0gSuuYNdZrVn2V1OWflWWZctgw9T83Zs0gbZt4eKLoVkz+2jY0NbklVKl400vnR4lbDfAAJ9FpMKXMbB8OcycCQsWcPm+fWSXqwQ3dYCOHXG2u4bUrARSUuCj12DtWrtbpUrQpg307AmtWsEll0D16kH9S5SKSEEbHllFkD//hKlTbaLfuhUqV4Zbb+Xr2t1YW+s6zk2syNy58OFwWzQmxib4Z56Bq6+Gli3txU2llH9pwlell54OL70Es2dDdjZccw2MGgXdupH2cxzPD4evv4Z9++x3QMeOtjWnY0c4/fRgB69U9NGEr07eunUwfDh8+CFUrAh9+8KQIRyJb8KcOTDxKtvxJiYGmjeHYcOgc2dbVCkVPJrwlfe2bIERI+Cdd2y/xlGjYMAAth2uwX//C9Onw5499gLrq6/aC6zly8M//xnswJVSoAlfeePQIXj2WRg3zja2Dx0Kjz3GL7uq8/zj8Oabtgtm164wYABcdZW9EWfmzGAHrpRypwlfFe+DD+DBB+GPP+Cuu+C559hyuC5PDYB337V93++9Fx591Pa4dOdwaN94pUKJJnzl2Z49MGiQvSDbvDksX07meVcyZgxMmmQT+aOPwsMP2ztZPXE4tN+8UqFEE74q7LPPoHdv24dy5EiyH/03r0yI5ZlOcPgw3HmnvTGqYkX7fVCUgwfhnHMCFrVSqgSa8FU+hwNGj7aP886DhQtZfjiJ+1vBhg22p82LL9pm/LfegjPPtMMdFOfCCwMTulKqZJrwlZWVZW91/eQT6N2b3WMm8MhTccycadvmU1Lgppts0U2b7HO7dlC3brACVkqdLE34Ctavt9X3HTvg9df5X4N+9GklZGbaDjnDh9vhD/I4HPbZ40iYSqmQpZfUot2yZXDFFXD0KIeWfMMDaf3peKNQvTqsWAHPPXdisgdN+EqFK0340eydd+CGG6BePX6cuZqL+ycxebLtebN6tR2p0hNN+EqFJ23SiVYTJx6/S+qtHh/Tv2slqla1HXSuuab4XfMSvna5VCq86CkbjcaPhwEDOHZjVx4493N63VeJli1hzZqSkz1oDV+pcKUJP9q8/DIMHkzmjXdx7d55TJoSw2OPweefe9/jRhO+UuFJm3SiycSJ8PDD/HT9g3T66RV27BDefRduu+3k3sbptM+a8JUKL5rwo8XcuTBwIEsv+zfdVj5DuXLC0qVw6aUn/1Zaw1cqPHnVpCMiHURkk4hsFpGhHrY3EpHPRWSdiCwTkXjfh6pK7bPP4M47efecp7gh9Rnq1hW+/750yR5swhexD6VU+Cgx4YtIDDAB6Ag0A3qISLMCxcYBs4wxFwKjged8HagqpR9+gFtuYVrtJ+nx80guvVT49ltISCj9WzoctnavCV+p8OJNk04rYLMxZiuAiMwBbgY2uJVpBgxxvV4KLPRlkMo6fBh27vS+fEzmn8R37cL4mEd5ZPtI2rSxzfhZWfZRWrt3a3OOUuHIm4RfH9jmtpwBtC5QZi1wK/BfoCtQWURqGGNOSCsi0h/oD9CwYcPSxhy1Fi2y08h6o2zuUe5K7spzf/dhhHMkTZvaiUnmzfNNLNWq+eZ9lFKB403C9/TD3RRYfhR4TUR6A8uB7UBuoZ2MmQJMAUhKSir4HqoER45ArVr5g5gVyRhOf+ReJmZcyQhG0bWrnayqrA8v0WvCVyr8eJMCMoAGbsvxwA73AsaYHUA3ABE5DbjVGLPPV0Eqy+m0wxGX+ONo/Ku8Or8KT/Ai3bvboYy1CUYp5U0vnVVAExFJEJFyQHcgxb2AiNQUkbz3ehKY7tswFXg5g9TKlbzx8Hoe5FVuucUwa5Yme6WUVWLCN8bkAgOBT4CNwFxjzHoRGS0iXVzF2gGbRORnoA7wrJ/ijWp5vWOKtGcPczq9SX/HJDpcm82cOUJsbMDCU0qFOK9adY0xi4BFBdaNcHs9D/DR5UBVlGITvjEs7fx//GvX/3Fl4kHe/7AK5csHNDylVIjTsXTCSHEJf/3wd+j67aM0qXOAhV9UoWLFwMamlAp9mvDDSFEJf8fyzXR87koqlney6PvTqV498LEppUKfjqUTRjwl/AN7cunU0cEeU53lKYdp1Fhvf1VKeaY1/DBSMOE7nXDXFZtZd/hs3ntyDS2urx284JRSIU8TfhhxOk9M+M8M2MmCjecxrsXbdPhP2+AFppQKC5rww4h7P/yF8x08PbkuvSq8x0OfdgpuYEqpsKBt+GEkr0ln/Xro1dNBEmt4/Y0YpMbpwQ5NKRUGNOGHkEOH7ND1OTmet+fmwtGj0PWmHCod28OCaydQ8Y4ZgQ1SKRW2NOGHkN9/t8PXV6vmuftlzZrw+mTDll/LsLTCv4ifPlUHpVdKeU0TfgjJmzrwzjttci/o9dfh40XCszxF2+c6ejGKmlJK5dOEH0KKmyt27VoYPNhwffkvGXr+pzDo+8AGp5QKe5rwQ4jTaZ8L3Vx1AP75T6gRu583D95GmQkpOgSmUuqkacIPIXk1/IJDIA8aBFu2GJaW6UbtXh1KP/u4UiqqacIPIZ6adObPh+RkeOrc92ibsQKe/zk4wSmlwp7eeBVCCib8nTvh3nvhkib7eWpTTxg2DOrVC16ASqmwpgk/hLgnfGOgb184dMjwlulJbEIDGDIkuAEqpcKaNumEEPc2/ClTYPFiGN/9O86b85Ft26lQIbgBKqXCmtbwQ0je0AmbN8PDD0P7axwM+PI2uOwy6No12OEppcKc1vBDiNNpb5zt2xfKlYMZrSdT5ovtMOdtvaNWKXXKvKrhi0gHEdkkIptFZKiH7Q1FZKmI/CAi60TkRt+HGvkcDli9Gr76Cl4afZD6E4dBp07QVoc+VkqduhITvojEABOAjkAzoIeINCtQbDgw1xjTAugOTPR1oJHEGDsQWsHHjh223f6666D3H6Nh/3547rlgh6uUihDeNOm0AjYbY7YCiMgc4GZgg1sZA1Rxva4K7PBlkJFm8WJYufLEdcbAO+/Y5ymjdiLXvgq9esEFFwQnSKVUxPEm4dcHtrktZwCtC5QZCSwRkUFAJeA6T28kIv2B/gANo3jgr127oGpVaNkyf91nn8Evv8Dw4ZAwa5Rt0B81KnhBKqUijjcJ39PVQlNguQcw0xjzfyJyGfCmiDQ3xjhP2MmYKcAUgKSkpILvETUcDqheHdq0scuZmbYTTuvWMPLu3+G86dCvHzRuHNQ4lVKRxZuLthlAA7fleAo32fQF5gIYY74DKgAeBvhVUHgy8scfh3374I03IOZFV5v90ELXxpVS6pR4k/BXAU1EJEFEymEvyqYUKPMHcC2AiDTFJvxMXwYaSdwnI//mG5g5Ex55BJpX+QOmT4d77oEGDYp9D6WUOlklJnxjTC4wEPgE2IjtjbNeREaLSBdXsUeAfiKyFngH6G2Midomm5Lk1fBzc+GBB2xuHz6c/B45Tz4Z1PiUUpHJqxuvjDGLgEUF1o1we70BuMK3oUUuh8MOnzBhAqxbZ0dNqJT1B0ybprV7pZTf6J22QeBw2C72Tz0FHTq4Rk14QNvulVL+pWPpBIHDAW++CceOwfjxINszbO2+b1+dp1Yp5Tea8IPg55/hyy/hiSegSRPgpZfsldwnngh2aEqpCKYJP8Byc2HhQqhTx3VtdvduOxZy9+7a714p5Vea8ANs2jT46y97X1XFisDEiXDokO2Mr5RSfqQJP4D27bMXahs1gquuAo4csY34HTvChRcGOzylVITThB9Azz5rx9G54QbXjVczZthxFbTtXikVAJrwA2TLFvjvf20XzLp1oSy5MG6cHUBHx7tXSgWAJvwAefxxiI2FAQPs8kWb58Gvv9ravc5mpZQKAE34AfDll/D++/aeqtNPB4wh7rUX4dxz4eabgx2eUipKaML3M6fTTkjeoIEdIM3hgMa/LaPsuh/g0UftGAtKKRUAOrSCn82ZA2vWwFtv2W6YTidcuuIVTI2ayJ13Bjs8pVQU0eqlH2Vn21EwExOhRw+7ruzvWzh304c4+t0HFSoEN0ClVFTRGr4fvf66vS67eHF+y02dua/iLBODue/+4AanlIo6WsP3kwMHYMwYaNfO9rsHYP9+an88nfXn306Z+HrBDE8pFYU04fvJSy/Ze6peeMGt1+WMGZQ9fIDvWw/Wa7VKqYDTtOMHf/9t76m69VZo1cq10uGAV19lb7PL+athS+16r5QKOE34fvDMM3aYnGefdVu5aBFs2cLvtww+YQJzpZQKFK8Svoh0EJFNIrJZRApNySQiL4tImuvxs4js9X2o4WHrVpg8Gfr0sfdVHffKKxAfz85Lu2rCV0oFRYm9dEQkBpgAtAcygFUikuKaxxYAY8wQt/KDgBZ+iDUsjBxpB0Z7+mm3lT/+CF98Ac8/T67EasJXSgWFNzX8VsBmY8xWY0w2MAcobjyAHsA7vggu3Pz0E8yebcfLqV/frluyBDYMnEhubAWSY+9hwwa9uVYpFRze9MOvD2xzW84AWnsqKCKNgATgiyK29wf6AzSMwLlbR4+2d9O6j3a8ceUB2q14iz8uvZ3YM2oQjx0PXymlAs2bhO+pP4kpomx3YJ4xxuFpozFmCjAFICkpqaj3CEvr19thFJ54AmrVyl9/3uq3KHfsIGeOvZ8zPX5NKqVUYHjTuJABNHBbjgd2FFG2O1HanDNqFFSqZMdDO84YWnw3if1ntXDrn6mUUsHhTcJfBTQRkQQRKYdN6ikFC4nIuUB14Dvfhhj61q2D996Dhx6CGjXcNnz7LbX/+pGMLvfrmPdKqaArMeEbY3KBgcAnwEZgrjFmvYiMFpEubkV7AHOMMRHVVOONkSOhalU7DLI7M2kSR8tXIfO6O4ISl1JKufNq8DRjzCJgUYF1Iwosj/RdWOFjzRpYsMA26VSv7rYhMxPee4+1F/VHTqsUtPiUUiqPdhA8RU8/bRP94MEFNsyYgWRnk5p0n/a7V0qFBE34p2DlSvjoI3uhtmpVtw1OJ7z+Os42bcmsfb72u1dKhQRNRadg5EioWRMGDSqwYckS2LqV7L52zHut4SulQoEm/FJavdpObPLww1C5coGNkyZB7dpkd+4GaMJXSoUGTfil9OyzUK2aHUbhBNu323aePn1wxJQDNOErpUKDJvxSSE+3PXMGD4YqVQpsnDnTtuH37YvTaVdpwldKhQKd09YL+/fb8e3zjBhh76rt3h3++sutoNNJjSnTcF5xNXsqn83uv+1qTfhKqVCgCb8Ehw/boezzautZWbBwIVx+uR07x13Cr8u4649fmX/JGH6clL++fPnAxauUUkXRhF+CI0dssr/0UmjYEIYNswn8uecKDKMA1H/sDRxVqtFseDeaVrDrYmMhISHwcSulVEGa8EvgcI372aABxMXBhx/CwIFw5ZUFCu7eDZ+/D/360fTiigGPUymlSqIXbUuQl/BjYuCFF+zzY495KDh7Nhw7BvfcE9D4lFLKW5rwS5DXdp+ZCdOn27lq82azOs4YmDoVLrkELroo4DEqpZQ3NOGXIK+GP2OGfe0+m9Vxqal23lqt3SulQpgm/BI4HHDoELzzDvTqBY0beyg0bZqd27BHj0CHp5RSXtOEXwKHA77/3jbPDx3qocChQ/D223DbbQVGUFNKqdCiCb8E+/bBqlVw441w7rkeCrz3Hhw4AH37Bjw2pZQ6GZrwSzB3rq3dFxozJ88bb8A550CbNgGNSymlTpYm/GLk5MCsWdCoEbRs6aHATz/BN9/Yi7U6Z61SKsR5lfBFpIOIbBKRzSLiqSUbEblNRDaIyHoRedu3YQbH3Lnw5592GAWP4+HMmGE3/OtfAY9NKaVOVol32opIDDABaA9kAKtEJMUYs8GtTBPgSeAKY8weEantr4ADxRgYO9YOi9CkiYeE73DAW2/Zxv06dYISo1JKnQxvavitgM3GmK3GmGxgDnBzgTL9gAnGmD0Axpi/fRtm4H32GaxdC1dfbVtrCk1T+PnnsGOH1u6VUmHDm4RfH9jmtpzhWufuHOAcEflGRL4XkQ6e3khE+otIqoikZmZmli7iABk71k5OXq8e1K5tu9mfIDnZFrjppqDEp5RSJ8ubhO/paqQpsFwWaAK0A3oAb4hItUI7GTPFGJNkjEmqVavWycYaMGlp8Omn0KED1K0LDzxQoEln/347A0r37jr2sVIqbHiT8DOABm7L8cAOD2U+MMbkGGN+BTZhvwDC0rhxcNpptjnH48XaefPsuMnanKOUCiPeJPxVQBMRSRCRckB3IKVAmYXA1QAiUhPbxLPVl4EGyh9/2IlN+vWDChWKSPjJybbvfevWAY9PKaVKq8SEb4zJBQYCnwAbgbnGmPUiMlpEuriKfQJkicgGYCnwmDEmy19B+9Mrr9jnhx6yI2UWSvi//grLl8Ndd2nfe6VUWPFqAhRjzCJgUYF1I9xeG+Bh1yNs7d1rRznu3t3ObuVwQNmCR2jWLJvo77wzKDEqpVRp6Z22biZPhoMH8yc4cTgK1PCNsQn/6qvtN4JSSoURTfgux47Bf/8L7dvnz2HicBTof//NN7B1q23OUUqpMKNz2rrMnm2HUZg1K39doRp+cjJUqgTdugU8PqWUOlVaw8denB03DhIT4brr8tefkPCPHLGD69x6q+2zqZRSYUZr+MCiRbBxox0ax73jzQkJ/4MP7A1X2pyjlApTWsMHXnzRXoO97bYT15+Q8JOToUEDaNcu0OEppZRPRH3CX7ECvvoKhgyB2NgTtx3vh79jByxZYie1LTSKmlJKhYeoa9J5+217N22e2bPtHbX79sHzz59Y9uhRV8J/+22b/XUoBaVUGIu6hP/HH1CtGjRubCvu69fbphxPoySIwMUtDDyaDJdeWsSktkopFR6iLuE7HHDmmXD99XYUzNhYePllOyqmR2t+gPR0mDQpoHEqpZSvRV2DdN6F2MxMO0Nhr17FJHuwHfPLlSt8RVcppcJMVCV8Y/IvxE6YYNvoH3mkmB1ycmz7fZcucPrpAYtTKaX8IaoSvsNhn3Ny4LXX7GRVTZsWs8PixfangPa9V0pFgKhK+E6nff7kE8jKyh8krUjJyVCrFtxwg99jU0opf4uqhO9w2KQ/Z47tldOmTTGFd++GDz+Enj0Ld9BXSqkwFFW9dBwO+Okn2L7djoxZ7Pwlc+bYth9tzlFKRYioquHn5sK339oREm65pYTCyclw4YV2RDWllIoAUZXwv/7a1u579Spirto8P/0EK1fqnbVKqYjiVcIXkQ4isklENovIUA/be4tIpoikuR73+D7UU/fqqxAXB127llBw1iz7jdCzZ0DiUkqpQCixDV9EYoAJQHsgA1glIinGmA0Fir5rjBnohxh9YsMGO/7ZVVfZOUyK5HDAm2/anjlnnBGw+JRSyt+8qeG3AjYbY7YaY7KBOcDN/g3Ld4yBLVvgqaegfHlo1aqE5pylSyEjQ5tzlFIRx5uEXx/Y5rac4VpX0K0isk5E5olIA09vJCL9RSRVRFIzMzNLEe7J++MPmDjRzl9ywQW2SScurpgdZs2CqlXh5rD5TlNKKa94k/A9dV40BZY/BBobYy4EPgOSPb2RMWaKMSbJGJNUq1atk4u0lI4ds2PeG2Pvrn3wQYiPL6LwgQMwfz7cfrsdM1kppSKINwk/A3CvsccDO9wLGGOyjDHHXItTgUt8E96p278fUlPhxhvtzVbFDokzfz4cPqzNOUqpiORNwl8FNBGRBBEpB3QHUtwLiIj7eJNdgI2+C/HUvPuureUP9OZy8qxZcPbZcPnlfo9LKaUCrcReOsaYXBEZCHwCxADTjTHrRWQ0kGqMSQEeFJEuQC6wG+jtx5i9lpNjJyZv1Ahatiyh8O+/2wu2o0eXcAuuUkqFJ6+GVjDGLAIWFVg3wu31k8CTvg3t1L37Lvz5J/ToUULPHLBdMcHelaWUUhEoYu+0NQbGjoWEBGjSpISEb4xtzrnqKjv3oVJKRaCITfiffgrr1tmbZUVKSPjffw+//KIDpSmlIlrEJvyxY+3Uhddea5fLFPeXJidDxYrwj38EJDallAqGiEz4P/wAn30Ggwfn1+yLrOEfPWob+7t1g8qVAxajUkoFWkQm/HHj4LTT4N5786c1LDLhp6TA3r3Qu3egwlNKqaCIuIT/+++2wt6/P1SrZhN+mTLF9LRMTra33l59dUDjVEqpQIu4hP/KKza5P/SQXXY4iqnd79wJ//ufFwPkK6VU+AvbKQ7377e5Ojc3f92hQzBpkh1C4auv7Lq//y4ml8+ebSe51d45SqkoELYJf9s2O8Z9rVr5c4wvXmyHUbj6ajh40K6Li4NzzvHwBsbY5pxLL4Vzzw1Y3EopFSxhm/DzLsZ27w41atjONiNGQIcOMGaMF2/www+Qnm5/EiilVBQI2zb8gr1vZs2Cv/6Cxx/38g2Sk+2MKLff7pf4lFIq1EREwnc4bFfMpCRo186LnbOzbft9ly5Qvbo/w1RKqZAR9k06MTF2NqtffoG5c70c6HLRIsjK0ou1SqmoEvYJv0wZeOEFOPNMe7OsV5KToU4dO1G5UqrUcnJyyMjI4OjRo8EOJeJUqFCB+Ph4YvN6pfhA2Cf8776DlSvtvLVedaXPzISPPrLjLpQN2z9fqZCQkZFB5cqVady4MaLzSPiMMYasrCwyMjJISEjw2fuGfRv+uHG2a6bXIyO8847tvK/NOUqdsqNHj1KjRg1N9j4mItSoUcPnv5zCNuE7nbayvmgRDBpkB7v0SnIytGgBF1zg1/iUihaa7P3DH8c1bBO+w2Gbc+Li4IEHvNwpPR3WrNHavVIqKoVtwv/zT1i7Fu65x9545ZXkZNtuf8cdfo1NKRU8I0eOZNy4cUVuX7hwIRs2bAhgRKHDq4QvIh1EZJOIbBaRocWU+4eIGBFJ8l2Ini1YYJ+HDPFyh9xcO6N5p0620V8pFZWiOeGX2E1FRGKACUB7IANYJSIpxpgNBcpVBh4EVvgj0Dw5OXZAtMWL4aKLTmIK2iVL7M8Cbc5Ryi/+9z97ivnSGR7N+2kAABAiSURBVGfY4VJK8uyzzzJr1iwaNGhArVq1uOSSS5g6dSpTpkwhOzubs88+mzfffJO0tDRSUlL48ssveeaZZ5g/fz5ffPFFoXJxcXG+/UNChDc1/FbAZmPMVmNMNjAHuNlDuTHAi4BfO+SuXGlz9tGjcN11J7FjcrJt++nUyW+xKaUCb/Xq1cyZM4cffviB999/n1WrVgHQrVs3Vq1axdq1a2natCnTpk3j8ssvp0uXLowdO5a0tDTOOussj+UilTcd0esD29yWM4DW7gVEpAXQwBjzkYg8WtQbiUh/oD9Aw4YNTz5aoF49SEuzQyAPGuTlTllZsHAh3HcflCtXqs9VShXPm5q4P3z11Vd07dr1eK28S5cuAKSnpzN8+HD27t3LwYMHuaGIGy29LRcJvEn4nvoGmeMbRcoALwO9S3ojY8wUYApAUlKSKaG4R59/bvP3u+9CgwZe7jR7th0/p0+f0nykUirEeerC2Lt3bxYuXMhFF13EzJkzWbZsmcd9vS0XCbxp0skA3FNrPLDDbbky0BxYJiK/AZcCKf66cHvuuXD//XDNNV7uYAxMmwaXXGIb/ZVSEaVt27YsWLCAI0eOcODAAT788EMADhw4QN26dcnJyWH27NnHy1euXJkDBw4cXy6qXCTypoa/CmgiIgnAdqA7cLxfozFmH1Azb1lElgGPGmNSfRuqdeWV9uG1NWtg3To79oJSKuJcfPHF3H777SQmJtKoUSOudCWIMWPG0Lp1axo1asQFF1xwPMl3796dfv36MX78eObNm1dkuUgkxpTcsiIiNwKvADHAdGPMsyIyGkg1xqQUKLsMLxJ+UlKSSU31y3fCiR54AGbMsPPXVqvm/89TKops3LiRpk2bBjuMiOXp+IrIamNMqVpQvBo9zBizCFhUYN2IIsq2K00gfnHkCLz9Ntx6qyZ7pVTUC9s7bb0yfz7s2wd9+wY7EqWUCrrITvjTp9uB8q+6KtiRKKVU0EVuwt+yBZYuhbvvtrOkKKVUlIvcTDhjhk30Xg+Ur5RSkS0yE77DATNn2ikM4+ODHY1SSoWEyEz4S5bA9u16Z61SEW7v3r1MDOA9No0bN2bXrl0AXH755cWW/c9//hOIkE5KZCb8adOgZk1wjamhlIpMxSV8R948qH7y7bffFrs9FBN+5M3i/ddfkJICAwfqQGlKBdJDD9mRDX0pMRFeeaXIzUOHDmXLli0kJibSvn17OnXqxKhRo6hbty5paWksWrSIzp07k56eDsC4ceM4ePAgI0eOZMuWLQwYMIDMzEzi4uKYOnUq55133gnvn5WVRY8ePcjMzKRVq1a436h62mmncfDgQXbu3Mntt9/O/v37yc3NZdKkSXz88cccOXKExMREzj///JAZsiHyEv6MGXbQ/P79gx2JUsrPnn/+edLT00lzfdEsW7aMlStXkp6eTkJCAr/99luR+/bv35/JkyfTpEkTVqxYwQMPPMAXX3xxQplRo0bRpk0bRowYwccff8yUKVMKvc/bb7/NDTfcwLBhw3A4HBw+fJgrr7yS11577XhcoSKyEr7TCVOm2H73Bb6plVJ+VkxNPJBatWpFQkJCsWUOHjzIt99+yz//+c/j644dO1ao3PLly3n//fcB6NSpE9WrVy9UpmXLlvTp04ecnBxuueUWEhMTT/Ev8J/IasP/9FP49Vc77r1SKipVqlTp+OuyZcvidDqPLx89audncjqdVKtWjbS0tOOPjRs3enw/T0Mvu2vbti3Lly+nfv369OrVi1mzZvngr/CPyEr4r79uL9Z27RrsSJRSAVBwqOOC6tSpw99//01WVhbHjh3jo48+AqBKlSokJCTw3nvvAWCMYe3atYX2b9u27fH298WLF7Nnz55CZX7//Xdq165Nv3796Nu3L2vWrAEgNjaWnJycU/4bfSlyEv6OHfZibZ8+UL58sKNRSgVAjRo1uOKKK2jevDmPPfZYoe2xsbGMGDGC1q1b07lz5xMuys6ePZtp06Zx0UUXcf755/PBBx8U2v/pp59m+fLlXHzxxSxZssTjTH3Lli0jMTGRFi1aMH/+fAYPHgzYawQXXnghPXv2BODGG29kx44dhfYPJK+GR/YHnw+PPGYMjBgBv/wCZ5/tu/dVShVJh0f2L18PjxwZNXyHA6ZOhfbtNdkrpVQRIiPhL14M27bBvfcGOxKllApZkZHwX38dzjhD76xVSqlihH/C/+MPWLTITnISGxvsaJRSKmR5lfBFpIOIbBKRzSIy1MP2+0TkRxFJE5GvRaSZ70MtwtSpYAz06xewj1RKqXBUYsIXkRhgAtARaAb08JDQ3zbGXGCMSQReBF7yeaSeHDtm76zt3BkaNQrIRyqlVLjypobfCthsjNlqjMkG5gA3uxcwxux3W6wEBKav59y58Pff8OCDAfk4pVToGT9+PE2bNqVnz56kpKTw/PPPA7Bw4UI2bNhwvNzMmTNPuh/8b7/9RvPmzUss17t3b+bNmwfAPffcc8LnFlSaOHzFm7F06gPb3JYzgNYFC4nIAOBhoBxwjac3EpH+QH/A4w0MJ8UYGD8emjaFa689tfdSSoWtiRMnsnjx4uPj53Rxdd5YuHAhnTt3plkz2yAxc+ZMmjdvTr169fwazxtvvFHs9kDF4Yk3Cd/TQBKFavDGmAnABBG5AxgO3OWhzBRgCtgbr04u1AJWrIDUVJgwAUoY60Ip5X9BGB2Z++67j61bt9KlSxf69OlD9erVSU1N5Y477iAlJYUvv/ySZ555hh49epCamkrPnj2pWLEi3333HRs2bODhhx/m4MGD1KxZk5kzZ1K3bl1Wr15Nnz59iIuLo02bNh4/1xjDoEGD+OKLL0hISDhh2OR27doxbtw4WrRoQd++fUlNTUVE6NOnDw0aNCgUR8WKFX170IrhTZNOBtDAbTkeKO73yBzgllMJyiuvvgpVqsC//uX3j1JKhabJkydTr149li5dypAhQ46vv/zyy+nSpQtjx44lLS2NJ554gqSkJGbPnk1aWhply5Zl0KBBzJs373iCHzZsGAB3330348eP57vvvivycxcsWMCmTZv48ccfmTp1qsfJUNLS0ti+fTvp6en8+OOP3H333fzjH/84IY5AJnvwroa/CmgiIgnAdqA7cId7ARFpYoz5xbXYCfgFf/rzT3jvPXjgATjtNL9+lFLKOyEyOrJXNm3aRHp6Ou3btwfs7Fh169Zl37597N27l6uuugqAXr16sXjx4kL7L1++nB49ehATE0O9evW45prCrdhnnnkmW7duZdCgQXTq1Inrr7/ev3+UF0pM+MaYXBEZCHwCxADTjTHrRWQ0kGqMSQEGish1QA6wBw/NOT71+ut2kpMBA/z6MUqpyGSM4fzzzy9Ui9+7d2+JwyHnKalc9erVWbt2LZ988gkTJkxg7ty5TJ8+vdQx+4JX/fCNMYuMMecYY84yxjzrWjfClewxxgw2xpxvjEk0xlxtjFnvt4izs2HyZOjYEZo08dvHKKXCW8Ghk92Xzz33XDIzM48n/JycHNavX0+1atWoWrUqX3/9NUCRUxO2bduWOXPm4HA42LlzJ0uXLi1UZteuXTidTm699VbGjBlzfNjkkoZ09qfwu9N2/nzbpDNoULAjUUqFsO7duzN27FhatGjBli1b6N27N/fddx+JiYk4HA7mzZvHE088wUUXXURiYuLxdvgZM2YwYMAALrvssiLb2Lt27UqTJk244IILuP/++483Abnbvn077dq1IzExkd69e/Pcc88BnBDHkSNHGDFiBCkpKf47EG7Cb3jkjz6CadNs4i8Tft9XSkUSHR7Zv3w9PHL4zWnbubN9KKWUOilaRVZKqSihCV8pdUqC1Swc6fxxXDXhK6VKrUKFCmRlZWnS9zFjDFlZWVSoUMGn7xt+bfhKqZARHx9PRkYGmZmZwQ4l4lSoUIH4+HifvqcmfKVUqcXGxh4ftEyFPm3SUUqpKKEJXymlooQmfKWUihJBu9NWRDKB30u5e01glw/D8TWN79SEcnyhHBtofKcqHOKrZIypVZqdg5bwT4WIpJb21uJA0PhOTSjHF8qxgcZ3qiI9Pm3SUUqpKKEJXymlokS4JvwpwQ6gBBrfqQnl+EI5NtD4TlVExxeWbfhKKaVOXrjW8JVSSp0kTfhKKRUlQjLhi8h0EflbRNLd1o0VkZ9EZJ2ILBCRam7bnhSRzSKySURuCEZ8btseFREjIjVdyyIi413xrRORi4MRm4gMch2f9SLyotv6oB87EUkUke9FJE1EUkWklWt9QI+d6zMbiMhSEdnoOlaDXetPF5FPReQX13P1QMdYTGwhcW4UFZ/b9mCfG0XGFwrnRzH/vr47P4wxIfcA2gIXA+lu664HyrpevwC84HrdDFgLlAcSgC1ATKDjc61vAHyCvaGspmvdjcBiQIBLgRVBOHZXA58B5V3LtUPp2AFLgI5ux2tZMI6d6zPrAhe7XlcGfnYdpxeBoa71Q93+/wUsxmJiC4lzo6j4XMuhcG4UdfxC4vwoJj6fnR8hWcM3xiwHdhdYt8QYk+ta/B7IGzf0ZmCOMeaYMeZXYDPQKtDxubwMPA64Xwm/GZhlrO+BaiJSN8Cx3Q88b4w55irzt1tsoXDsDFDF9boqsMMtvoAdO1d8O40xa1yvDwAbgfquWJJdxZKBWwIdY1Gxhcq5Ucyxg9A4N4qKLyTOj2Li89n5EZIJ3wt9sN9sYA/INrdtGeT/JwsYEekCbDfGrC2wKRTiOwe4UkRWiMiXItIyhGIDeAgYKyLbgHHAk671QY1PRBoDLYAVQB1jzE6wJyZQO5gxFojNXUicG+7xheK5UeD4hdz5USA+n50fYZfwRWQYkAvMzlvloVhA+5qKSBwwDBjhabOHdYHuC1sWqI792fcYMFdEhNCIDWwNa4gxpgEwBJjmWh+0+ETkNGA+8JAxZn9xRT2s82uMRcUWKueGe3yueELq3PBw/ELq/PAQn8/Oj7BK+CJyF9AZ6GlcjVjYb7UGbsXiyf/JEyhnYdv41orIb64Y1ojIGSESXwbwvuun30rAiR2EKRRiA7gLeN/1+j3yfzYHJT4RicWecLONMXlx/ZX3c9n1nPezP6AxFhFbyJwbHuILqXOjiOMXMudHEfH57vzw1wWIU30AjTnxwl4HYANQq0C58znxwspW/Hzh0VN8Bbb9Rv6FqU6ceGFlZRCO3X3AaNfrc7A/AyVUjh22rbKd6/W1wOogHjsBZgGvFFg/lhMv2r4Y6BiLiS0kzo2i4itQJmjnRjHHLyTOj2Li89n54beDe4p/+DvATiAH+y3WF3vBZBuQ5npMdis/DHsFfROuq9mBjq/Advf/1AJMcMX3I5AUhGNXDngLSAfWANeE0rED2gCrXSfXCuCSYBw712e2wf4sXuf2f+1GoAbwOfCL6/n0QMdYTGwhcW4UFV8InRtFHb+QOD+Kic9n54cOraCUUlEirNrwlVJKlZ4mfKWUihKa8JVSKkpowldKqSihCV8ppaKEJnyllIoSmvCVUipK/D82NxArp9klgAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(x, dcnorm.ecdf(x,X), color='b', alpha=0.5, label='data')\n",
    "plt.plot(x, rv.cdf(x), color='r', label='true dist.') # plot for a<x<b\n",
    "plt.plot(x, rv_fit.cdf(x), color='b', label='fitted dist.') # plot for a<x<b\n",
    "plt.legend()\n",
    "plt.title('CDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we'll compute the statistical moments of the distribution (mean, variance, skewness, kurtosis), noting that the mean and variance are calculalated with analytic expressions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean, variance, skewness, kurtosis (array(135.37345464), array(238.43710092), array(0.9142011), array(0.24932382))\n"
     ]
    }
   ],
   "source": [
    "m, v, s, k = rv.stats('mvsk')\n",
    "print(\"mean, variance, skewness, kurtosis\", (m, v, s, k))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}