4K4K.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1 - Get Unit Weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shapes: (20, 1)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11310f1d0>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD8CAYAAABpcuN4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8XOV96P/Pd0b7LmuxrcWbZBYbDMYLGENuEgIhKbGbhARTEnACpVzictN7c1N45Ta/hm7ZaPtKQpsSQliyQEM20Zhih9BQbBYLbyDZRrLRZu2jfZdmvr8/5hx5kEfSSJpdz/v10sszZ845c2YszXee5/k+30dUFcMwDMOYjiPSF2AYhmFENxMoDMMwjBmZQGEYhmHMyAQKwzAMY0YmUBiGYRgzMoHCMAzDmJEJFIZhGMaMTKAwDMMwZmQChWEYhjGjhEhfQDDk5+frqlWrIn0ZhmEYMeXNN9/sVNWC2faLi0CxatUqKisrI30ZhmEYMUVE6gPZz3Q9GYZhGDMygcIwDMOYkQkUhmEYxoxMoDAMwzBmZAKFYRiGMSMTKAzDMIwZmUBhGIZhzMgEihjWOzzOr440RfoyDMOIcyZQxLCfvdHAXzxzjMauoUhfimEYccwEihj29tleABq7TaAwDCN0TKCIYdXNfQA094xE+EoMw4hnJlDEqMHRCd51DQJwtns4wldjGEY8M4EiRp1o6UPVe7u5xwQKwzBCxwSKGFVldTstz07hrAkUhmGEkAkUMaqquZfctESuWJFrWhSGYYSUCRQxqqq5j/VF2RTnpnK2Zxi1+6EMwzCCzASKGDQ24aGmbYD1RVkU56QyOuHBNTgW6csyDCNOmUARg2rbBxhze1hXlEVRTipgMp8MwwgdEyhiUFWzd6Ld+qJsinJSAJP5ZBhG6AQUKETkRhE5JSK1InK/n8d3i0iHiBy1fu6a8niWiDSJyPes+5k++x4VkU4R+edAzmV4xydSE52szk+nJCcNwGQ+GYYRMgmz7SAiTuBh4HqgCTgkIhWqWj1l12dUdc80p/kb4GX7jqr2A5f7PMebwC8DPNeiV93cx8XLM3E6hKzUBNKTnCZQGIYRMoG0KLYCtap6RlXHgKeBnYE+gYhsApYC+6Z5/AKgEPjvQM+5mHk8SnWLN+MJQES8mU9mjMIwjBAJJFAUA40+95usbVN9UkSOi8izIlIKICIO4CHgSzOcfxfeFoRvfud555pKRO4WkUoRqezo6AjgZcSHhq4hBkYnWFeUNbmtKCeV5l4TKAzDCI1gDWY/B6xS1Q3AfuAJa/u9wF5VnWnRhF3AzwI413uo6iOqullVNxcUFCz4BcQKe0b2ep9AUZyTagoDGoYRMrOOUQBnAd9v9SXWtkmq6vK5+yjwTev2NuBaEbkXyACSRGRAVe8HEJHLgARVfTOAcxl4M56cDuGCpZmT24pyUukaHGNobIK0pED+Sw3DMAIXSIviELBWRFaLSBLeFkCF7w4istzn7g7gBICq3qaqK1R1Fd7upyftIGG5lfe2JqY9l+FV3dLH2sIMUhKdk9tKcr1zKUyrwjCMUJj166eqTojIHuAFwAk8pqpVIvIgUKmqFcB9IrIDmAC6gN0BPv+ngY9O2Tbfcy0KVc19XLs2/z3bJifd9QxTXpgRicsyDCOOBdRPoap7gb1Ttn3V5/YDwAOznONx4PEp29b42W/Wcy1W7f0jdPSPTmY82Ypz7BaFGdA2DCP4zMzsGOJvIBugMDMZp0NMiqxhGCFhAkUMsZc+XTclUCQ4HSzLSjEtCsMwQsIEihhS1dzLiiVpZKUknvdYcU4qTSZQGIYRAiZQxJCq5j7WLc/y+1hxbqppURiGERImUMSIvpFx6l1D541P2IpyUmjtHcHtMQsYGYYRXCZQxIgT9kB28TQtipw0JjxKe7+ZS2EYRnCZQBEjqlvsjKdsv4/b61KYzCfDMILNBIoYUdXcR35GEoWZyX4ft2dnm3LjhmEEmwkUMaKquY91RdmIiN/HfWdnG4ZhBJMJFDFgdMJNTVv/tAPZAGlJCeSmJZrMJ8Mwgs4EihhQ0zbAhEdnDBTgbVWYMQrDMILNBIoYUNXcC0w/kG0rMutSGIYRAiZQxICq5j7Sk5ysXJI2437FOamc7RnmvYsFGoZhLIwJFDGgqrmPi5dn4XD4H8i2FeekMjA6Qd/IRJiuzDCMxcAEiijn8SgnWvpmHZ8AbxkPMHMpDMMILhMoolyda5ChMfes4xNwLkXWZD4ZhhFMJlBEuappSov7U7zI5lI0dg3xp09WMjBqutoMI5RMoIhyVc19JDqFC5ZmzrpvXnoSSQmORdOiePFEG/ur23irqTfSl2IYcS2gQCEiN4rIKRGpFZH7/Ty+W0Q6ROSo9XPXlMezRKRJRL7ns+2/rHPaxxRa25NF5BnruV4XkVULe4mxraq5l7WFmSQlzP5f5XDIolqXoqZ9ADBdbYYRarOumS0iTuBh4HqgCTgkIhWqWj1l12dUdc80p/kb4GU/229T1cop2+4EulW1XER2Ad8AbpntOuORqlLd3McHLyoM+JiinMWz0p0dKBZLV5thREogLYqtQK2qnlHVMeBpYGegTyAim4ClwL4AD9kJPGHdfha4TqYrcBTn2vpGcQ2OBZTxZCteRLOza+1AsUher2FESiCBohho9LnfZG2b6pMiclxEnhWRUgARcQAPAV+a5tw/srqd/sonGEw+n6pOAL1AXgDXGXfsGdnrAsh4shXlpNLeP8rohDtUlxUVXAOjdA2OAdDcawKFYYRSsAaznwNWqeoGYD/nWgT3AntVtcnPMbep6qXAtdbPZ+fyhCJyt4hUikhlR0fHAi49etkZTxcvn30g22ZnPrX2xncpD7vbKTs10bQoDCPEAgkUZ4FSn/sl1rZJqupS1VHr7qPAJuv2NmCPiNQB3wZuF5GvW8ectf7tB36Kt4vrPc8nIglANuCaelGq+oiqblbVzQUFBQG8jNhT3dzHqrw0MlMSAz5msaTI2oHimrX5pmyJYYRYIIHiELBWRFaLSBKwC6jw3UFElvvc3QGcAFDV21R1haquwtv99KSq3i8iCSKSbx2bCNwEvG0dXwHcYd2+Gfi9LtJPgaqW3oAm2vlaLLOza9v6SU9ysnllLqMTHlxWN5RhGME3a9aTqk6IyB7gBcAJPKaqVSLyIFCpqhXAfSKyA5gAuoDds5w2GXjBChJO4HfAD6zHfgg8JSK11rl2zf1lxb7e4XEau4bZtWXFnI5blu1dEjXeq8jWtA9QvjSTklxvocSz3cPkZ/hf/c8wjIWZNVAAqOpeYO+UbV/1uf0A8MAs53gceNy6Pci57qmp+40AnwrkuuJZdbO9RnbgGU8AyQlOCjOT4z5FtrZ9gPddUHBurfCeYS4rzYnwVRlGfDIzs6NUoGtQ+FNklRuPV71D47T3j7K2MIOSHG+LIt4Do2FEkgkUUaq6uY/CzGQKMufenVKcmxrXH5y1Hf0ArF2aQVZqAhnJCTTF+ZiMYUSSCRRRqqo5sNLi/sT7AkY1bd6Mp7WFmYgIRTkpcd2CMoxIM4EiCo2Mu6ntGAioYqw/xTmpcZ0JVNM+QEqiYzIVuDgnvltQhhFpJlBEoXfa+nF7dF7jE3BuXYp4TZGtaR+gvDBjcsW/4tz4HpMxjEgzgSIKVc0z48lWHOcLGNW29bO28Nxs9aKcVHqGxhk061IYRkiYQBGFqpp7yUxOoNSaIzBX8Tw7u39knObeEcoLMya3xXtgNIxIM4EiClU193FxUdZk18pcZaUmkJ7kjMtAcbpjEOA9gaLEmo2+WNbhMIxwM4Eiyrg9ysmW/nl3OwGIiLffPg7HKGrarNRYn0AR72MyhhFpJlBEmXc7Bxged897INtWlJMal+W3a9sHSHI6WLHkXLdcYWYKCQ4xXU+GESImUESZhQ5k2+J1AaPa9gHWFKST4Dz3q+t0CMvNXArDCBkTKKJMVXMfSU7He/rg56MoJ5XuoXGGxuIrE8hOjZ2qKDs+A6NhRAMTKKJMVXMvFyzLING5sP8ae4A3nrpjhsfcNHYPvSc11hbvZUsMI5JMoIgiqkp1cx/rly9sfAJ8BnjjqNz46Y4BVL01nqYqyUmltW+EcbcnAldmGPHNBIoo0tI7QvfQOOuLFzY+AT5zKeKoO6a23a7x5KfrKScVj8b/ErCGEQkmUESRYA1kAxRmJuOMs0ygmvZ+EhzCyrz08x4rjsOuNsOIFiZQRJGq5l5E4KJlCw8UCU4Hy7LiKxOopm2AVfnpJCWc/2sbz7PRDSPSTKCIIlXNfazOTyc9OaCFB2dVHGcLGNW2D1Be4D8bzEy6M4zQMYEiilQ39y14op2veJqdPTrhps416HcgGyAl0Ul+RlJcTjI0jEgLKFCIyI0ickpEakXkfj+P7xaRDhE5av3cNeXxLBFpEpHvWffTROS3InJSRKpE5OuBnitedQ+OcbZnOCjjE7ainBRa+0Zwe2J/AaN3OwfxKDPOLynOSTUr3RlGCMzaxyEiTuBh4HqgCTgkIhWqWj1l12dUdc80p/kb4OUp276tqi+JSBLwooh8RFWfD+Bccam6JXgD2bbinDTcHqWtb2SyayZW+a5qN52inFROWbWgDMMInkBaFFuBWlU9o6pjwNPAzkCfQEQ2AUuBffY2VR1S1Zes22PAYaBkLhceb6qtjKd1y4PbooD4yASqaR/AIbCm4PyMJ5u90l28LgFrGJESSKAoBhp97jdZ26b6pIgcF5FnRaQUQEQcwEPAl6Y7uYjkAB8DXpzpXPGuqrmXZVkp5GUkB+2c9uzseBjQPt0+wIolaaQkOqfdpzg3lZFxD11xugSsYURKsAaznwNWqeoGYD/whLX9XmCvqjb5O0hEEoCfAd9R1TOznGvqsXeLSKWIVHZ0dATpZUROVXNfULudwHd2duwHipr2fspn6HaC+Hq9hhFNAgkUZwHfb/Ul1rZJqupS1VHr7qPAJuv2NmCPiNQB3wZu9x24Bh4BalT1nwM413uo6iOqullVNxcUFATwMqLX8Jib0x0DQQ8UaUkJ5KYlxnzm07jbw7ud02c82cxKd4YRGoEEikPAWhFZbQ087wIqfHcQkeU+d3cAJwBU9TZVXaGqq/B2Pz2pqvdbx/wtkA18MZBzxbOTrX14FNYFMTXWVpQT+8Xy6l1DjLvVb+kOX5Mr3cV4YAzUyLibR14+zfCYO9KXYsS5WbOeVHVCRPYALwBO4DFVrRKRB4FKVa0A7hORHcAE0AXsnumcIlICfAU4CRwWEYDvqeqjcz1XPAhm6Y6pinNSqXMNBv284VTbbq9qN3PXU3ZqImlxugSsPxXHmvn7vSdJdDr43PbVkb4cI44FNAVYVfcCe6ds+6rP7QeAB2Y5x+PA49btJsDvgtCBnCveVDX3kZ2aOPmNOJiKclI5UNuJqmIF5Jhjp8aWFU6f8QTWErBx0IIK1L6qNgCefqOR3Vevitn/XyP6mZnZUaC6uZd1y7NC8odenJPK4JibvuHYXcCopn2A4pxU0pJm/15TnBtfZUumMzzm5pXaDgozkznV1s+xpt5IX5IRx0ygiLAJt4eTrf0h6XaCc1VVY/nDs6Z9YNaBbFtRnC4BO9XLNR2MjHt4cOclpCY6efqNhkhfkhHHTKCIsDOdg4xOeFgXokAR6ymjbo9yumNg1oFsW3GcLgE71f7qNrJSErju4kJu2rCcimPNDIzG92s2IscEigiravZ2GQSzGKCvWE8ZbewaYmzCM+tAti0el4CdasLt4cUTbVx38VISnQ52bS1laMzNb483R/rSjDhlAkWEVZ3tIznBQdkMpSkWIi89iaQER8x+cNZYq9qVz6HrCeI7RbayvpvuoXGuX7cUgCtW5LK2MIOnDzXOcqRhzI8JFBFW1dzHRcsySXCG5r/C4fBmAjXFaKCwlz+dqWqsr3MtqPhdEnV/dRtJCQ7ed4F3oqmIcMuWUo409HCq1RRFNILPBIoIUlWqmntDMtHOV1FOSgy3KPpZlpVCVkpiQPsvzUrB6RDO9gyF+MoiQ1XZV93K9rI8MnwWuPrEFSUkOoWnD5lBbSP4TKCIoKbuYfpGJkKW8WQrjuFMoNo5ZDwBOB3iXQI2Rl/vbE619dPYNcwN65e9Z/uS9CRuWL+MXx05y8i4maltBJcJFBEUyhnZvopyUmnvH2V0IrY+QDwe9S5/GmC3k604NzVuu572VbUhAtddXHjeY7duWUHP0Dj7qtsicGWxp3twjO++WEPv0HikLyXqmUARQXZpiguXBZbRM192v31rb2x9eDb3DjM05g4448lWEmdrhfvaV93KxtIcCjNTznvs6rI8SnJTzZyKAJxs7WPHw6/w0P53qDh2dvYDFjkTKCKozjXE0qzkgGYcL0RxjM6lsDOe5tL1BN4WVGvfCBNuTyguK2Kae4Z5+2zfed1ONodDuGVzKQdPu6iP8fpeofSfb7fwiX85yOi4h5REB6c7zHs1GxMoIqjeNcjKvNCkxfqanJ0dY/32tVaNp/KCuXc9uT1KW//o7DvHkP1Wl9INVlqsP5/aXIpD4N8rTarsVB6P8o/7TnHPjw9z4bJMnvvza7hgaSanOwYifWlRzwSKCKpzDbEqLy3kz7Ms214SNba6nmra+8nPSCI3PWlOx022oGIsMM5mX3UrZQXprJkhcC7LTuEDFxby88qmuGtRLUT/yDh3P/Um3/l9LZ/eXMLTd1/F0qwUygoyON1uAsVsTKCIkMHRCTr6R8PSokhOcFKYmRxzKaM18xjIBt+yJbH1emfSOzTO62e6pu128nXLllLa+0d56VTsr/wYDO92DvLxfznIS6fa+dqO9XzjkxtITvAuqVtWkE5z7wiDpvzJjEygiJB6l/dDbFUYAgXYCxjFTotCValtG5jzQDbE56S7l061M+HRGbudbB+4qJCCzGSeMXMq+MM7Hez83iu4BkZ56s6t3DGlHHuZ1Tp7t9OMU8zEBIoIsQcbV4ah6wlir/x2W98o/aMTcx7IBkhNcpKXnhRXZTz2V7dRmJnMZSU5s+6b6HTwqU0l/P5ke8xlugWLqvJvfzjN5370BkU5qVTsuYary/LP26/MarHWmu6nGZlAESF1VosibIHCShlV1bA830LVWKnD8+l6AqvceAwFxpmMjLv5r1PtfGjdUhyOwNYs+fTmUjwKvzjcFOKriz4j426++MxR/uH5k3zkkuX88t6rKV3i/+9sZV4aToeYAe1ZmEARIfWuQfIzksgMsDTFQhXnpDI24aFzYCwsz7dQ9je8+XQ9AXG10t2rp10MjrkniwAGYlV+OtvW5PH0oQY8ntj4chAMzT3DfOr7r1JxrJn/++EL+d6fbJwx/Tw5wcmKJWkmUMwioEAhIjeKyCkRqRWR+/08vltEOkTkqPVz15THs0SkSUS+57Ntk4i8ZZ3zO2J1HIrIEhHZLyI11r+5C32R0aguTKmxtqIYKzde0z5ATloi+Rlzy3iyFed6y5bESgtqJvuq20hPcnJ1Wd6cjtu1tZTGrmFePeMK0ZVFl0N1Xez43iu82znIo7dv5gsfKA9o1ciygnROt5sxipnMGihExAk8DHwEWAfcKiLr/Oz6jKpebv08OuWxvwFenrLtX4E/BdZaPzda2+8HXlTVtcCL1v24U+8aClu3E8TepDvvQHbGvJeHLcpJZXjcTXeMl2fweJTfnWjj/RcVTmbqBOrD65eRnZq4KMqP/+T1ev7kB6+RlZLIr7+wnesuDrz1VVaQwbudg7gXUctrrgJpUWwFalX1jKqOAU8DOwN9AhHZBCwF9vlsWw5kqepr6v3K9yTwx9bDO4EnrNtP+GyPGyPjblp6R1gdxhZFLC1gpKq8095P+Ty7nSC2Xu9Mjjb10NE/GlC201QpiU4+vrGYF95upXswNroc52pswsNXfvUWX/nV22wvz+dXX9g+53GtsoIMxtwemrrjJ5062AIJFMWA71eSJmvbVJ8UkeMi8qyIlAKIiAN4CPiSn3P6jrL5nnOpqrZYt1vxBpm40tBlDWTnhy9QZKUmkJ7kjIlMINfgGD1D4wEvf+qPvdJdLLzemeyraiPBIbz/wvOLAAbili2ljLk9/PJI8OsZdQ2O8c3/PElLb2Te497hcT7z6Ov85PUG/uf7y/jhHVvITp37mF9Zoffv0IxTTC9Yg9nPAatUdQOwn3MtgnuBvao6r9QLq7Xhtz0oIneLSKWIVHZ0xNbEojorZzscs7JtImJVVY3+D86atrktVuRPrK8VbttX3cq2srx5fQACXLw8i8tKc3jmUENQx2v6Rsa547E3+Jf/Os1tP3idzoHwlksZGpvg848f4khjN9+5dSN/eeNFOAPMCJvKnkthximmF0igOAuU+twvsbZNUlWXqtq/KY8Cm6zb24A9IlIHfBu4XUS+bh1fMs0526yuKbuLqt3fRanqI6q6WVU3FxQUBPAyokedPYdiSfhaFBA7KaN2Vd35zKGw5aYlkprojInAOJ3a9gHOdAzOKdvJn1u3lPJO2wBHGnuCcl3DY27uerySEy19fOmGC2jpHeGzP3wjbOW6R8bd3P3kmxxp6OY7uzay47KiBZ0vJy2J/IwkM5diBoEEikPAWhFZLSJJwC6gwncH+4PdsgM4AaCqt6nqClVdhbf76UlVvd/qWuoTkausbKfbgd9Yx1cAd1i37/DZHjfqXEPkpiWSnRae1FhbrKSM1rQPkJGcwLKs80tpB8puQcVyvSe7COCH5jAw689NlxWRluTkmTcWPqg9OuHmz378JpX1XfzzrsvZ88G1PHL7Jk63D/C5x98IeSmMcbeHP//ZEV6p7eSbN1/GRy5dPvtBAVhTkBGTXU9Vzb1hSX+eNVCo6gSwB3gBbwD4d1WtEpEHRWSHtdt9IlIlIseA+4DdATz3vXhbH7XAaeB5a/vXgetFpAb4kHU/roSrauxURTmpdA+NMzQW3XVtatq8NZ7mm/Fki5UW1HT2VbdyaXH2ZDfafGUkJ/CxDUU8d7yZgQV8kE+4PXzx6aO8/E4H//CJS7lpg/eb/LVrC/jOrRs51tTL3U9VhmyFPY9H+fKzx9lf3cbXdqzn5k0lsx8UoLIYDBRPv9HAzu8d4NFXzoT8uQIao1DVvap6gaqWqerfWdu+qqoV1u0HVHW9ql6mqh9Q1ZN+zvG4qu7xuV+pqpdY59xjjUfY3VjXqepaVf2QqnYF56VGj7rO8FSNncoe4I32VkVN+8CCBrJtsdKC8qe9b4SjjT3zynby55atpQyNuXnuWPO8jvd4lPt/+RbPv93KX920jlu2rHjP4zdesoxv3byBA7Uu9vz0CONBrlyrqny14m1+deQs//fDF3LH1auCev6ygnS6h8bpioHsMLdH+bvfVnP/L99iW1keu7aumP2gBTIzs8NsdMJNc+9wxFoUEN2ZQD1DY3QOjC5ofMJWkpuKa3CM4bHYWgIW4Hcn2lEloGqxgdhYmsOFSzPnNadCVXnwP6p59s0m/uJDF3DnNav97veJK0r4m53r+d2JNr7082NB7RL5xn+e4sevNXDP/yjj3veXBe28NrvmU7S3KgZGJ/izpyr5wX+/yx3bVvKj3VvICkN1BxMowqyxaxhVWJUf/hZFLFRVXWjpDl9FOd4xjljsftpf3cqKJWlcEISACd4xm1u2lHKssYcTLX1zOvYf97/D4wfruOua1dx3XfmM+3522yq+fOOF/OZoM//vN28HJdPq4Zdq+f4fTvOZq1bwlzdeuOAuSX/KJzOfojdQnO0Z5uZ/PchLpzp4cOd6vrbzEhKc4fkIN4EizM5VjQ1/i6IwMxmnQ6J6nQZ7+dOFpMbainO8wTjWup8GRic4UOvihnVLg/qh+PGNxSQ5HTwzh1bF9/9wmu/+vpZdW0r5yh9dHND13Pv+cv7n+8v46esNfP35kwsKFk8crONbL5zi4xuLeXDHJSEJEuBtbScnOKK2RXG4oZud3zvA2e5hfrR7C7dvWxXW5zeBIszqwrwOha8Ep4NlWSlR3aKoaRsgNdE52fpZiMklYGMsUPzhVAdjbk/Qup1suelJfPiSZfzycFNAA84/fq2erz9/kps2LOfvPn7pnD6kv/zhC/nsVSv5t5fP8PBLtfO63mffbOL/q6jihnVL+dbNGwKunDsfToewOj89KlNkK441s+uR10hLcvKrL1zN+y4I/3QAEyjCrN41SGZKArlhTo21FedEd8poTXs/5YUZQflQWGq3oKL49fqzv7qVJelJbFoZ/HqYt24ppW9kgheqWmfc79dHzvJXv3mb6y4q5J9uuXzOk9lEhK/tWM8nNhbz7X3v8KMD787p+OffauHLzx7jmvJ8vvsnG8PSxVJemMHpjuiZdKeq/OP+d7jvZ0e4vCSHX39h+4LK2iyECRRh5l0nOz1kTejZRPsCRrVByngC3xZU9L7eqcbdHn5/sp3rLiqc90zjmVy1Jo8VS9J4eoY5FfuqWvk/Pz/GVavzePi2K0ic54e0wyF88+YNfHj9Ur72XDX/XhlYl9cf3ungvqePsHFFLo/cvmnOxRDnq6wgg8buoZCl987FyLibP//ZEb7zYg03byrhqbu2smSOa8cHkwkUYeadQxH+gWxbUU4KrX0jTAQ5fTEY+kfGaekdmcxACYbinFSaYihQvPFuF30jEwuejT0dh8M7qP3qGddkKRlfr9R0suenR7i0OJsf3LGZlMSFfUgnOB1859aNXLs2n/t/cZy9b7XMuP+hui7+7KlK1hZm8tjuLTOuJRFsZYUZqJ6rnBAp7X0j3PLIa/z2rRbu/8hFfOvmDWELltMxgSKMxt0emrqHIzI+YSvOScPtUdr7w1ubJxDnMp6CFyiKclJiqutpX1UrKYkOrl0bun7omzeV4HQIz0z5hv9mfRd/+mQlawrSefxzW8hIDs6HdHKCk3/77CauWJHL/3r6CC+d8luVh7eaevn8jw5RlJPKk3dunXd9q/kqK7CKA0aw5lNVcy87Hz7AO639fP8zm7jnf5RFrPfBlwkUYXS2exi3RyPeooDozASyM57WLg1eP2xxbiqtfSMxsdaAqrK/uo33rS0gNSl03yCXZqXwgQsLefbNpsmJcVXNvez+0SGWZiXz5J1byUkLbjdHWlICP9y9hQtw6ufIAAAgAElEQVSWZnLPU2/y+pTFlGra+rn9sdfJSk3kJ3ddSX5GclCfPxBr8iM7l2J/dRuf+v6rAPz8nm18OMjJDAthAkUY2U3aVWEsLz5VSRRnAtW2D5CU4KA0d+EZTza7BdXWF72ZXraq5j6ae0dC1u3ka9eWUjr6R3npZDunOwa4/YdvkJmcwI/vupLCzPnX2JpJdmoiT35+KyW5qdz5RCXHm7xFChtcQ3zmh6+T4HTwk7uuZHl28P7/5yI1yZttF+5Aoar82x9Oc/dTlawtzOA3X9jOJcXZYb2G2ZhAEUb1VmpsZFsU0Rsoatr6WZOfHtQMl1iadLevqhWHMKfV2ebr/RcWUJiZ7J3I9ujriMCP77qSktzQ/m7mZSTz47uuJCctkdsfe4OX3+ngth++xuiEhx/feWVEv0SBd5winCmyYxMe/vIXx/mH50/y0UuX8/Td2yhcQDHMUDGBIozqXIOkJTkpiECz2paW5E3NjcZ++5r2gaB2O0Hs1LcC79rYW1YtCUt2S4LTwac2l3C4oYfB0QmeuvNK1hQEb2xoJsuzU/nJXVeS5HRw+2Nv0D04zhOf28qFyyKT+umrvCCDMx2DYanIOuH2cMdjb/DvlU3cd91avrtrY0i7HBfCBIow8q6THbnUWFtRFBbLGxqboKl7OKgD2RAb9a3A2/1ysrU/LN1Ots9etYr3X1jA45/fysXLs8L2vOCtTPDju65ke3keP7xjM5eV5oT1+adTVpjO8LibljB0VR5r6uXVMy7+3x9dzP++/oKQTihcqPDlnhnUuQa5MMjfmOejOCc14imAU52xJjoFO1BMtqCiLDBOta/aOwHuhnXhG8Bclp3C45/bGrbnm+qCpZn85K6rIvb8/pT51HwKRnWAmRxp6AZY8MJL4WBaFGHi9iiNXUMRqfE0VZE1OzuYS2MuVE0QVrWbTiwsAbu/uo2LlmWyIoLjV4ZPoAjDgPbhhm6Kc1KjckxiKhMowqS5Z5hxt0ZkHYqpinNSGRxz0zccPQsY1bQNkOCQkATSaC9b0jU4xqG6rqCtPWHMX35GElkpCeEJFPU9XBGCMi2hYAJFmJzLeIp8iyIai+XVtA+wOj993uUiZmKvdBdNLShfL55owxPEtSeM+RMRygozQj7prqV3mNa+Ea5YER1jM7MxgSJM7DGB1RFO/4PoTJGtbR8ISmlxf4pzUhkac9M7PB6S8y/U/uo2irJTWF8U3gFlw79wLIt6uN47h+SKFaZFYfiodw2SkuigMDNyqbG2cwsYRUegGBl3U+8aDPpAts1OkY3GzKfhMTcv13RwfZDXnjDmr7wwg/b+UfpGQvfF4nBDN8kJjrBnm81XQIFCRG4UkVMiUisi9/t5fLeIdIjIUevnLmv7ShE5bG2rEpF7rO2ZPvseFZFOEfnnmc4V6+pcQ6xckh4VKXB56UkkJTiipkXxbucgHoXyEGWERWMLyvZKbScj48Ffe8KYv7IwrHZ3uKGbS4uzSUqIje/qs6bHiogTeBi4HmgCDolIhapWT9n1GVXdM2VbC7BNVUdFJAN42zq2Gbjc5zneBH45y7liWr1rMKLFAH05HOId4I2SD86aEBQD9BVtLShf+6payUpJYOvqJZG+FMMyWRywY5CNIegaGp1wU3W2j93bVwX93KESSDjbCtSq6hlVHQOeBnYGcnJVHVNVu0xpsr/nE5ELgELgvwO75Njj8Sj1rqGIlyfwFU1VVWvb+nFI6MZvlqQnkZLoiJrXa3N7lBdPtvPBiwpDMohvzE/pkjQSnRKycYqq5j7G3J6YGciGwAJFMeBbj7jJ2jbVJ0XkuIg8KyKl9kYRKRWR49Y5vmG1JnztwtuC0NnOFata+0YYnfBEtMbTVMVRNDu7pn2AlXnpC177YDoiMpn5FE0O1XXRNTjG9WGcZGfMLtHpYGVeesi6ng7XeyfahaK1EirB+hrzHLBKVTcA+4En7AdUtdHaXg7cISJTk8V3AT8L5Fy+RORuEakUkcqOjo4gvYzQmKwaGyVdT+Dtt2/vH2V0IvKredWEMOPJFk2B0fbcsWZSE5184KLwr4FszKysID1kLYojDT0U56SyNAYm2tkCCRRnAd9v9SXWtkmq6vLpYnoU2DT1JFZL4m3gWnubiFwGJKjqm3M5l7XfI6q6WVU3FxRE9x9aNFSNncrut2/tjWz57XG3h7rO0GU82UqibAnYcbeHvW+18KF1S8O6ipsRmLKCDOpdQ5PrdQTT4YZuNsZQtxMEFigOAWtFZLWIJOFtAVT47iAiy33u7gBOWNtLRCTVup0LXAOc8tn3Vt7bmpj2XLGszjVIktMRsTr7/tiBItL99vWuQSY8GpLSHb6KslPpHBiLivWQwZvt1D00zs4YqPOzGJUXZjBhjS0GU0vvMC29IzEzf8I261cZVZ0QkT3AC4ATeExVq0TkQaBSVSuA+0RkBzABdAG7rcMvBh4SEQUE+LaqvuVz+k8DH53ylNOdK2bVdw5RuiQVZxSkxtqiZXZ2TZud8RTaYonFPuXGw1VOeyYVR5vJTk3kfRdEd2t4sfKt+RTMbtEjDdZEuxgp3WELqM2rqnuBvVO2fdXn9gPAA36O2w9smOG8a/xs83uuWFYXRamxtmXZ0bGgT037ACLn/jBDpdhnLkWkA8XwmJt9Va187LKimMmjX2zWTKbIBnec4nB9N0kJDtbFyEQ7m/ktDTFVnVyHIpokJzgpzEyO+ABvjVXOOdQLthRFSVcbwO9PtjM45o6J8tKLVWZKIkuzkoNe8ynWJtrZYutqY1BH/yjD425W5UfPQLYtGlJGa9r6Qz6QDd4WlEOiY9JdxbGzFGYmc+WavEhfijGDYNd8Gp1w83ZzX0zNn7CZQBFidVFUNXYq7zoNkct6mnB7ONM5GPTlT/1JdDpYlpVCU4QDRe/wOC+d6uCmDUVRNWZlnM8OFMGqOlzd3MfYhCfmBrLBBIqQOzeHIvpaFHYZj3CsD+xPY/cwYxOekM+hsBVFwboUL1S1MjbhYcflptsp2pUVpNM/MkHHwOjsOwfgcIwOZIMJFCFX7xokwaqtFG2Kc1IZm/AE7Q9hrk61ele1C1egKM5Npbk3soHiuWPNrMxL47KS7IhehzG7cisTrzZIM7QPN3RTlJ0SUxPtbCZQhFida4iS3FQSorCWzyXF3g+ryrruiDz/G+92kRzGDJDinFRaekZwR6gF1dE/yoHaTj62ociUFI8BZYXnigMGw5H6bjbGYGsCTKAIuXrXYFSOTwBcVpJNRnICB053RuT5D57uZPOq3JDVeJqqKCeVCY/S3h+ZcZm9b7XgUdhpup1iwrKsFNKSnEGp+dTaO0JzDE60s5lAEUKqSn3nUFSOTwAkOB1ctWYJB2rDHyg6+kc52drP9vL8sD2n76S7SKg41sxFyzLDMnhvLJyIBC3z6UiDt9UeixlPYAJFSHUNjtE/OhG1LQqAq8vyqXcN0dgV3FIFszlotWK2l4UvUJTkRG6lu8auId6s7zaD2DGmrCCdM0Hoejrc4J1ot74oNsemTKAIITs1NhrnUNiuWev9oD4Y5u6nA7WdZKUkTI6ThEMkV7p77ri3uv7HNphAEUvKCjI42zPM0NjEgs5zuKGHS4qyYm6inS02rzpG1FupsdHcolhbmEFBZjIHal1he05V5UCti21leWGdS5CenEBOWmJEup4qjjZzxYocSpdE75cG43xlVkbeQloVYxMe3jrbG7PjE2ACRUjVuYZwiLfEdbQSEbaX5XHwdGfQJhbNpqFriLM9w1wTxvEJW3EE5lLUtPVzsrWfnZf7W+/LiGZ26vZCximqW6yJdjGa8QQmUIRUvWuQopxUkhPCk9UzX1eX59M5MMaptv6wPN8r1uD51REIFJEoW1JxrBmHwEcvXT77zkZUWZmXhkNYUOaTvaKdaVEYftW5hkK2DnQw2ZlH4ep+OljrYllWCmsi8N7YLYpwtZ5Uld8cbWZ7eT4FmclheU4jeJITnKxYkraguRSHG7pZnp0yWbE5FplAEULeORTR3yddnJPK6vz0sKTJejzKwdOdbC/Pj8iks5LcVAbH3PQNL2xwMlDHmnpp6BriY6ZSbMxaaIrskYaemG5NgAkUIdMzNEbP0HjUrUMxnavL8nj9jCskSz/6qm7po3tonO3lkamcamc+NfWEJx244mgzSU4HH16/LCzPZwRfWWEGZzoH5zWjv71vhLM9wzG39OlUJlCESH0UV43155ryfAbH3Bxr7Anp89itlnBOtPNl19wKR9Vct0f5j+PNfOCiArJTE0P+fEZolBWkMzbhmVcSxGF7ol0MD2SDCRQhE81VY/3ZVpaHSOjHKQ6cdlFemBGxwmiTS8B2h75F8fq7Ltr7R9lxmcl2imW+y6LO1eGGHpKcDtYXxdaKdlOZQBEi9a4hRIiZvPmctCQuKcoO6TjF6ISbQ+92RSQt1paXnkRygiMsmU8VR5tJT3Jy3cWFIX8uI3QWFCjqu1lfnBX1mY+zCShQiMiNInJKRGpF5H4/j+8WkQ4ROWr93GVtXykih61tVSJyj88x/2Wd0z6m0NqeLCLPWM/1uoisCs5LDa+6zkGWZ6WEreBdMFxdnseRxu4Fz0KdzpGGHobH3VxdFrmV3US8Jd9D3fU0NuHh+bdbuWH9spj6HTDOl5ueRF560pzLjcfDRDvbrIFCRJzAw8BHgHXArSKyzs+uz6jq5dbPo9a2FmCbql4OXAncLyK+6R+3+RzTbm27E+hW1XLgn4BvzO+lRVZdFFeNnc415fmMu5U33u0KyfkP1nbiELgqgoECvN1PoV7p7uV3OugdHjfrYseJ+WQ+nWjpYzRGV7SbKpAWxVagVlXPqOoY8DSwM5CTq+qYqtqr4iQH+Hw7gSes288C10kMFu+vdw1FdY0nfzavXEKS0xGy7qdXajvZUJJDVkpkB3aLskM/O7viWDO5aYmTtbSM2FZWmD7nuRTnBrJjO+MJAvvgLgYafe43Wdum+qSIHBeRZ0Wk1N4oIqUictw6xzdUtdnnmB9Z3U5/5RMMJp9PVSeAXuC8r6AicreIVIpIZUdHRwAvI3z6RsZxDY7FXIsiNcnJFStzQjKg3T8yzrGm3oilxfoqzk2lc2CUkXF3SM4/NDbB/uo2PnrpchKjcMEqY+7KCjLoGhyja3As4GMON/SwPDuF5dnRW8InUMH6LX4OWKWqG4D9nGsRoKqN1vZy4A4RWWo9dJuqXgpca/18di5PqKqPqOpmVd1cUFAQlBcRLA121dgYyXjydU15PtUtfbiCvDzqG+924fZoxNJifdkpsi29oRmn2F/dxvC423Q7xRF7QPvMHLqfjjR0x/z8CVsggeIsUOpzv8TaNklVXT5dTI8Cm6aexGpJvI03KKCqZ61/+4Gf4u3ies/ziUgCkA2Er7RpENTFQNXY6dj1l149E9y3/JXaTpITHFHRXztZbjxE3U/PHWtmeXYKW1YtCcn5jfCba+ZTe/8ITd3DUfH7HgyBBIpDwFoRWS0iScAuoMJ3BxHxrXa2AzhhbS8RkVTrdi5wDXBKRBJEJN/angjchDeIYJ37Duv2zcDvNVyFeYLk3GS72GtRbCjOJjM5IejjFAdrXWxdvSQqMoBKQrjSXc/QGH94p4ObNizHEcYS6kZoFeemkpTgCHic4nC9d+LqxjgJFAmz7aCqEyKyB3gBcAKPqWqViDwIVKpqBXCfiOwAJoAuYLd1+MXAQyKigADfVtW3RCQdeMEKEk7gd8APrGN+CDwlIrXWuXYF6bWGTV3nIIWZyaQlzfr2Rp0Ep4Mr1+QFdZyivX+EU239/PHG6Jh4tiw7BRFCkvn0n2+3Mu5WU1I8zjgdwpr89ICryB5p6CbJ6eCS4tieaGcL6JNMVfcCe6ds+6rP7QeAB/wctx/Y4Gf7IH66p6zHRoBPBXJd0areNRQzNZ78uaY8j9+daKOxaygoEwZfPe2yzhv58QmARKeDpZkpIel6qjjWzJr89JifiWucr6wwg7fP9ga07+GGbtYVxf5EO5tJyQiBuhipGjudc2XHg9P99EpNJ9mpiayLog/P4tzUoHc9tfWN8OoZFx+7rCgilXGN0CoryKCxa2jWbLlxt4fjTfEx0c5mAkWQDY1N0N4/yqoYWIdiOuWFGRRmJnPg9MK7n1SVg6ddbFsT3mVPZ1McggWM/uN4C6qw43KT7RSPygrS8ei5McjpTE60i4P5EzYTKIIslgeybSLC9vJ8DtZ24plHaWVf9S7vsqfbo2ziWVFOKi29wwt+fb4qjjWzvihrMkPGiC+BZj7Fw4p2U5lAEWT1k1VjY7dFAd71KVyDC18e1V72dHuEy3ZMVZybyrhb6QjSfJF61yDHGnvYaVoTcWtNgfdverYB7cMNPSzLSplMw44HJlAEWZ3VolgRwy0KCN44xcHTnRRlp0TdkrAl9gJGQRrQfu6Yt+DATRtMoIhXaUkJFOekztqiONIYPxPtbCZQBFm9a5C89KSI1zNaqKKcVNYscHlU77KnLq6O0LKnM5mcdBeEcQp7Xeytq5bE1bdI43xlhRkzzqXo6B+lsSt+JtrZTKAIsrrOoZgen/C1vTyf19/tYmxifsujVrf00TM0HjVpsb6Kgzjp7mRrPzXtA3zMdDvFvbKCdE53DEw7thVPhQB9mUARZPWuwZgfn7BtL89jaMzNsab5LY9qt0Yiuf7EdDKSE8hOTQzKXIqKY804HcJHLzHrYse7soIMhsbctPb5rxN2uKGbRKewvig7zFcWWiZQBNHIuJvm3pGYrPHkz7Y1+dbyqPPrfnqltpMLlmZQGKFlT2dTFIQUWVWl4mgz167NJy8jOUhXZkSr2TKfjjT0sK4oOypK1QSTCRRB1NhlVY2NsXUoppOdlsilxfNbHnV0ws2hui6uLou+biebd6W7hQWKww3dnO0ZNpViF4mywukzn7wT7Xq4Is4GssEEiqCqmywvHh8tCoCry/I50tDD4Ojclkc9XN/DyLgnKsqKT6ckd+ELGFUcbSY5wcEN602302JQkJFMZkqC3wHtky39jIzHx4p2U8Ve1booFi9zKHxdU57P9/9wmjfe7eIDFxUGfNzB0504HcKVa6K31HZRTgr9oxM0dg2RkTz3PwWPKr99q4XrLi6c1/FG7BGRaZdFPTeQbQKFMYM61yA5aYlkp8V2aqyvzatySUrwLo86l0DhXfY0O6rThFdYBQ+v/eZLCzqP6XZaXMoLM/jvmvNX1Tzc0M3SrGSKsqNzTG4hTKAIonrXUNwMZNtSEp1sXpk7OcM6EP0j4xxv6uXe95eF8MoW7gMXFfL1T1y6oCVR05MTuH6d6XZaTMoKMnj2zSb6Rsbf80XoSEMPG0tzo27OUDCYQBFEda7BuOyf3F6ez7deOEXnwCj5AWT2vH7Gu+xpNA9kAyQnONm1dUWkL8OIMWVWKY8zHYNcXuoduO4cGKWha4jPXBWfv09mMDtIxiY8nO0ejrsWBZwr5/FqgNVkX6ntJCXREXeTjgwDvLOz4b2ZT/FYCNCXCRRB0tg9hEdhVZzMyvZ1aXE2mSmBL4968HQnW1YtiZtFWwzD14olaSQ45D0D2ocbekh0CpcUx9dEO5sJFEFiZzzFY4vC6RCuWpPHgdOzB4r2vhHeaRuI6rRYw1iIRKeDlXlp7wkURxq6Wbc8K+4m2tkCChQicqOInBKRWhG538/ju0WkQ0SOWj93WdtXishha1uViNxjbU8Tkd+KyElr+9dnO1e0q+u051DEX4sCvGmyjV3DNMyyaMvBKFv21DBCwZsi6/1yOGGtaLcxTrudIIDBbBFxAg8D1wNNwCERqVDV6im7PqOqe6ZsawG2qeqoiGQAb4tIBdADfFtVXxKRJOBFEfmIqj4/w7miWr1rkMzkBJakJ0X6UkJie7m3XtOB052syJt+wO5AbSc5aYmsWx49y54aRrCVF2bw0ql2xt0eTrX2Mzzujsv5E7ZAWhRbgVpVPaOqY8DTwM5ATq6qY6pqrwyTbD+fqg6p6kv2PsBhoGSuFx9N6lxDrMxPi8vUOPB+g1qalTxjmqyqcqC2k6vL8nBE0bKnhhFsZQUZjLuVxq6hcxPt4rB0hy2QQFEMNPrcb7K2TfVJETkuIs+KSKm9UURKReS4dY5vqGqz70EikgN8DHhxtnMF2+mOAf5+74mgLIdZ7xqMy/EJm4iwvSyfV0+7pn2/6lxDNPeORH1arGEslJ35VNs+wOH6bgozkymO47VIgjWY/RywSlU3APuBJ+wHVLXR2l4O3CEiS+3HRCQB+BnwHVU9M9u5fInI3SJSKSKVHR3nz5IMxEsn23nk5TP8w/Mn5nW8bdztoal7OG7HJ2zby/PpGhzjRGuf38ft1oYZnzDi3eSyqB2DHGnsYeOKnLjtTYDAAsVZwPdbfYm1bZKquny6mB4FNk09idWSeBu41mfzI0CNqv7zXM5l7feIqm5W1c0FBQUBvIzz3XnNanZfvYof/Pe7PPLy6XmdA7yL30x4NK5bFHBuPsXBWv/zKQ7WdlKckxo3CzcZxnSyUhIpzEzmjXdd1LuG4nb+hC2QQHEIWCsiq62B511Ahe8OIrLc5+4O4IS1vUREUq3bucA1wCnr/t8C2cAXAzlXKIgIX71pHTdtWM7f7z3JL95smtd54rFqrD/LslMoK0j3mybrtpY93V6eF9ffrAzDVlaQwcs13r+FeB7IhgCynlR1QkT2AC8ATuAxVa0SkQeBSlWtAO4TkR3ABNAF7LYOvxh4SEQUELyZTm+JSAnwFeAkcNj6YPmeqj46w7lCwuEQHvr0ZXQPjfHlXxxnSXrSnIrfgW/V2Pj/Jr29PJ+fVzYxNuEhKeHc94zq5j56h8fN/Alj0SgrTOfVMy4SHMKlcTrRzhZQrSdV3QvsnbLtqz63HwAe8HPcfmCDn+1NeAOHv+fye65QSk5w8v3PbOLWH7zGvT85zE/+9Mo5NSXrOodITXRSkBn/K5xdXZbPk6/Wc7Sxh62rz5UQt1sZ26Jw2VPDCIVya7W7dUXxO9HOZmZmWzJTEvnR7q0UZiXz+ccPUetnBavpeDOe4jc11te2NXk4hPPSZA/UdnLh0kwKM+OvxLJh+GNnPsX7+ASYQPEeBZnJPPn5rSQ4HNzx2Bu09vpfQH2qOtdg3I9P2OzlUQ/6BIqRcWvZ03LTmjAWj3XLs8hMTuCDc+yqjkUmUEyxMi+dxz+3hd7hce547A16h8Zn3N/tURq7hlkZJ+tkB2J7eT5HG3sYsJZHPdzQzci4x6TFGotKXkYyx//6Bt53wfyyLmOJCRR+XFKczSOf3cS7nYPc9eShGRe2aekdZsztWTQtCvAGigmP8sa73jTZg7UunA55z5iFYSwGi6G7GUygmNbV5fn80y2XU1nfzZ6fHmHC7fG7X72VGruY5g5sWplLcoKDA9Z8igOnO7m8NIfMKF721DCM+TOBYgZ/tGE5X9uxnt+daOP//fptVM8vXVE3mRq7eFoUKYlONq/K5UBtJ30j4xxr7GG7yXYyjLhlAsUsbt+2ij//YDlPH2rkH/e/c97j9a4hkhMcLMtaXNk+V5flc7K1n98eb8Gj3haYYRjxyQSKAPzv6y9g15ZSvvv7Wp44WPeex+o6vamxi61aqj1w/Z0Xa0hNdLIxjitnGsZiF9CEu8VORPjbP74E1+AYf/1cFfkZyfzRBm+lkXrXUNzXePLnkuJsslISaOkd4X0XFJhlTw0jjpkWRYASnA6+e+tGNq/M5S+eOcrB2k48HqW+a3BRlO6YyumQyVnY15j5E4YR10ygmIOURCeP3r6F1fnp3P3Um7x0qp2Rcc+ibFEAXLvWmz9+TXn855EbxmJmAsUcZacl8sTnt5Kdmsg9P34TWFwZT75u2VLKz+/Zxrois+ypYcQzEyjmYVl2Ck98fisZyd4hnsU0h8JXotPBllVmkp1hxDszmD1P5YUZPHXnlfz2rRZKcuN3CUTDMAwTKBbgkuJsLonzOvSGYRim68kwDMOYkQkUhmEYxoxMoDAMwzBmZAKFYRiGMaOAAoWI3Cgip0SkVkTu9/P4bhHpEJGj1s9d1vaVInLY2lYlIvf4HLNJRN6yzvkdsQq7i8gSEdkvIjXWv/G/zqBhGEYUmzVQiIgTeBj4CLAOuFVE1vnZ9RlVvdz6edTa1gJsU9XLgSuB+0WkyHrsX4E/BdZaPzda2+8HXlTVtcCL1n3DMAwjQgJpUWwFalX1jKqOAU8DOwM5uaqOqeqodTfZfj4RWQ5kqepr6l3k4Ungj639dgJPWLef8NluGIZhREAggaIYaPS532Rtm+qTInJcRJ4VkVJ7o4iUishx6xzfUNVm6/imac65VFVbrNutwFJ/FyUid4tIpYhUdnR0BPAyDMMwjPkI1oS754CfqeqoiPwZ3pbABwFUtRHYYHU5/VpEng30pKqqInL+snLexx4BHgGwxkfq53nt+UDnPI8NB3N9C2Oub+Gi/RrN9c3fykB2CiRQnAVKfe6XWNsmqarL5+6jwDennkRVm0XkbeBa4IB1Hn/nbBOR5araYnVRtc92gao67/KlIlKpqpvne3yometbGHN9Cxft12iuL/QC6Xo6BKwVkdUikgTsAip8d7A+0G07gBPW9hIRSbVu5wLXAKesrqU+EbnKyna6HfiNdXwFcId1+w6f7YZhGEYEzNqiUNUJEdkDvAA4gcdUtUpEHgQqVbUCuE9EdgATQBew2zr8YuAhq/tIgG+r6lvWY/cCjwOpwPPWD8DXgX8XkTuBeuDTC36VhmEYxrwFNEahqnuBvVO2fdXn9gPAA36O2w9smOaclcAlfra7gOsCua4geSSMzzUf5voWxlzfwkX7NZrrCzHxZqcahmEYhn+mhIdhGIYxo0UTKAIoQ5IsIs9Yj78uIqvCeG2lIvKSiFRbpU7+l5993i8ivT5lUr7q71whvMY6q+TKURGp9PO4WKVYaq35NFeE8dou9HlfjopIn4MQp2YAAAQ2SURBVIh8cco+YX//ROQxEWm3sv3sbQGVqBGRO6x9akTkDn/7hODaviUiJ63/v1+JSM40x874uxDia/xrETnr8//40WmOnfHvPYTX94zPtdWJyNFpjg3Lexg0qhr3P3gH4U8Da4Ak4Biwbso+9wLft27vwluSJFzXtxy4wrqdCbzj5/reD/xHBN/DOiB/hsc/ijchQYCrgNcj+H/dCqyM9PsHvA+4AnjbZ9s3gfut2/fjnYQ69bglwBnr31zrdm4Yru0GIMG6/Q1/1xbI70KIr/GvgS8F8Dsw4997qK5vyuMPAV+N5HsYrJ/F0qIIpAyJb+mQZ4HrrNTdkFPVFlU9bN3ux5te7G/2ezTbCTypXq8BOVPSpsPlOuC0qs53AmbQqOrLeLMAfQVSoubDwH5V7VLVbmA/52qhhezaVHWfqk5Yd1/jvXOdwm6a9y8Q8y47NBczXZ/12fFp4GfBft5IWCyBIpAyJJP7WH8svUBeWK7Oh9XltRF43c/D20TkmIg8LyLrw3phoMA+EXlTRO7283igpV5CbRfT/3FG8v2zBVKiJhrey89zLmV9qtl+F0Jtj9U99tg0XXfR8P5dC7Spas00j0f6PZyTxRIoYoKIZAC/AL6oqn1THj6MtzvlMuC7wK/DfHnXqOoVeKsIf0FE3hfm55+VNSF0B/BzPw9H+v07j3r7IKIu7VBEvoJ3TtRPptklkr8L/wqUAZfjrU79UBifey5uZebWRNT/PflaLIFi1jIkvvuISAKQDbgIExFJxBskfqKqv5z6uKr2qeqAdXsvkCgi+eG6PlU9a/3bDvwKb/PeVyDvcah9BDisqm1TH4j0++ejze6Sk+lL1ETsvRSR3cBNwG1WIDtPAL8LIaOqbarqVlUP8INpnjuiv4vW58cngGem2yeS7+F8LJZAMWsZEt5bOuRm4PfT/aEEm9Wf+UPghKr+4zT7LLPHTERkK97/u7AEMhFJF5FM+zbeQc+3p+xWAdxuZT9dBfT6dLGEy7Tf4iL5/k0RSImaF4AbRCTX6lq5wdoWUiJyI/BlYIeqDk2zTyC/C6G8Rt9xr49P89yB/L2H0oeAk6ra5O/BSL+H8xLp0fRw/eDNynkHbzbEV6xtD+L9owBIwdtlUQu8AawJ47Vdg7cL4jhw1Pr5KHAPcI+1zx6gCm8Gx2vA1WG8vjXW8x6zrsF+/3yvT/AucHUaeAvYHOb/33S8H/zZPtsi+v7hDVotwDjefvI78Y57vQjUAL8Dllj7bgYe9Tn289bvYi3wuTBdWy3evn37d9DOAiwC9s70uxDG9+8p6/frON4P/+VTr9G6f97feziuz9r+uP1757NvRN7DYP2YmdmGYRjGjBZL15NhGIYxTyZQGIZhGDMygcIwDMOYkQkUhmEYxoxMoDAMwzBmZAKFYRiGMSMTKAzDMIwZmUBhGIZhzOj/B4PhsR5XgbURAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "weights = np.array([[41, 43, 46, 45, 44, 47, 44, 30, 32, 32,\n",
    "                     43, 45, 35, 27, 21, 27, 46, 25, 35, 34,\n",
    "                    ]]\n",
    "                  ).transpose()\n",
    "numbers = np.array([[77, 81, 85, 82, 82, 88.5, 82, 55, 60, 60, \n",
    "                     80, 83, 65, 50, 39, 50, 85, 47, 65, 63,\n",
    "                    ]]\n",
    "                  ).transpose()\n",
    "print(\"shapes:\", weights.shape)\n",
    "\n",
    "fish_weights = weights/numbers\n",
    "plt.plot(fish_weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([4., 2., 0., 3., 4., 4., 1., 0., 1., 1.]),\n",
       " array([0.5308642 , 0.53265583, 0.53444746, 0.53623908, 0.53803071,\n",
       "        0.53982234, 0.54161397, 0.5434056 , 0.54519723, 0.54698886,\n",
       "        0.54878049]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEYxJREFUeJzt3X+MZWV9x/H3h2VBrT9QdlTK7jgaaBp/KxOKMW2ohgTBQhsxwbQKRrLxB1ETmwZsipE0KdjWHxRTsgErWqtYNGaVJZZWiPoH6C5dVgGtK6FhKS3roihRMWu//WPOwvU6M/fMzJ29M0/fr+Rkz3nOc5/z3TPPfvbMub9SVUiS2nLEpAuQJI2f4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lq0JGTOvCmTZtqZmZmUoeXpHVp165dP6iqqVH9JhbuMzMz7Ny5c1KHl6R1Kcl/9unnbRlJapDhLkkNMtwlqUGGuyQ1yHCXpAb1DvckG5L8e5IvzbPv6CTXJdmb5LYkM+MsUpK0NEu5cn8XcPcC+94C/LCqTgA+BFy+0sIkScvXK9yTbAbOBK5eoMvZwLXd+vXAq5Nk5eVJkpaj75X7h4E/A/53gf3HA/cBVNVB4GHg2BVXJ0lalpHvUE3yWuDBqtqV5NSVHCzJVmArwPT09LLHmbnohpWUsSL3XnbmxI79/80kf86T4vzSuPS5cn8lcFaSe4HPAK9K8o9Dfe4HtgAkORJ4GnBgeKCq2lZVs1U1OzU18qMRJEnLNDLcq+riqtpcVTPAucBXqupPhrptB87r1s/p+tRYK5Uk9bbsDw5Lcimws6q2A9cAn0yyF3iIuf8EJEkTsqRwr6pbgFu69UsG2n8OvH6chUmSls93qEpSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDRoZ7kick+UaSO5LcmeT98/Q5P8n+JLu75YLVKVeS1Eefr9l7FHhVVT2SZCPw9SQ3VtWtQ/2uq6oLx1+iJGmpRoZ7VRXwSLe5sVtqNYuSJK1Mr3vuSTYk2Q08CNxUVbfN0+11SfYkuT7JlrFWKUlakl7hXlW/rKqXApuBk5O8cKjLF4GZqnoxcBNw7XzjJNmaZGeSnfv3719J3ZKkRSzp1TJV9SPgZuD0ofYDVfVot3k1cNICj99WVbNVNTs1NbWceiVJPfR5tcxUkmO69ScCpwHfGepz3MDmWcDd4yxSkrQ0fV4tcxxwbZINzP1n8Nmq+lKSS4GdVbUdeGeSs4CDwEPA+atVsCRptD6vltkDvGye9ksG1i8GLh5vaZKk5fIdqpLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktSgPt+h+oQk30hyR5I7k7x/nj5HJ7kuyd4ktyWZWY1iJUn99LlyfxR4VVW9BHgpcHqSU4b6vAX4YVWdAHwIuHy8ZUqSlmJkuNecR7rNjd1SQ93OBq7t1q8HXp0kY6tSkrQkI78gGyDJBmAXcALw0aq6bajL8cB9AFV1MMnDwLHAD4bG2QpsBZienl5Z5TpsZi66YdIlSFqiXk+oVtUvq+qlwGbg5CQvXM7BqmpbVc1W1ezU1NRyhpAk9bCkV8tU1Y+Am4HTh3bdD2wBSHIk8DTgwDgKlCQtXZ9Xy0wlOaZbfyJwGvCdoW7bgfO69XOAr1TV8H15SdJh0uee+3HAtd199yOAz1bVl5JcCuysqu3ANcAnk+wFHgLOXbWKJUkjjQz3qtoDvGye9ksG1n8OvH68pUmSlst3qEpSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KD+nyH6pYkNye5K8mdSd41T59TkzycZHe3XDLfWJKkw6PPd6geBN5TVbcneQqwK8lNVXXXUL+vVdVrx1+iJGmpRl65V9UDVXV7t/4T4G7g+NUuTJK0fEu6555khrkvy75tnt2vSHJHkhuTvGCBx29NsjPJzv379y+5WElSP73DPcmTgc8B766qHw/tvh14TlW9BPg74AvzjVFV26pqtqpmp6amlluzJGmEXuGeZCNzwf6pqvr88P6q+nFVPdKt7wA2Jtk01kolSb31ebVMgGuAu6vqgwv0eXbXjyQnd+MeGGehkqT++rxa5pXAG4FvJdndtb0XmAaoqquAc4C3JTkI/Aw4t6pqFeqVJPUwMtyr6utARvS5ErhyXEVJklbGd6hKUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg/p8h+qWJDcnuSvJnUneNU+fJLkiyd4ke5K8fHXKlST10ec7VA8C76mq25M8BdiV5Kaqumugz2uAE7vld4C/7/6UJE3AyCv3qnqgqm7v1n8C3A0cP9TtbOATNedW4Jgkx429WklSL32u3B+TZAZ4GXDb0K7jgfsGtvd1bQ8MPX4rsBVgenp6aZWuETMX3TCR49572ZkTOa6k9an3E6pJngx8Dnh3Vf14OQerqm1VNVtVs1NTU8sZQpLUQ69wT7KRuWD/VFV9fp4u9wNbBrY3d22SpAno82qZANcAd1fVBxfoth14U/eqmVOAh6vqgQX6SpJWWZ977q8E3gh8K8nuru29wDRAVV0F7ADOAPYCPwXePP5SJUl9jQz3qvo6kBF9CnjHuIqSJK2M71CVpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDWoz9fsfSzJg0m+vcD+U5M8nGR3t1wy/jIlSUvR52v2Pg5cCXxikT5fq6rXjqUiSdKKjbxyr6qvAg8dhlokSWMyrnvur0hyR5Ibk7xgTGNKkpapz22ZUW4HnlNVjyQ5A/gCcOJ8HZNsBbYCTE9Pj+HQkqT5rPjKvap+XFWPdOs7gI1JNi3Qd1tVzVbV7NTU1EoPLUlawIrDPcmzk6RbP7kb88BKx5UkLd/I2zJJPg2cCmxKsg94H7ARoKquAs4B3pbkIPAz4NyqqlWrWJI00shwr6o3jNh/JXMvlZQkrRG+Q1WSGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaNDLck3wsyYNJvr3A/iS5IsneJHuSvHz8ZUqSlqLPlfvHgdMX2f8a4MRu2Qr8/crLkiStxMhwr6qvAg8t0uVs4BM151bgmCTHjatASdLSjeOe+/HAfQPb+7o2SdKEHHk4D5ZkK3O3bpienj6ch5bWhZmLbpjIce+97MyJHBcm93eepMNxvsdx5X4/sGVge3PX9muqaltVzVbV7NTU1BgOLUmazzjCfTvwpu5VM6cAD1fVA2MYV5K0TCNvyyT5NHAqsCnJPuB9wEaAqroK2AGcAewFfgq8ebWKlST1MzLcq+oNI/YX8I6xVSRJWjHfoSpJDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkN6hXuSU5P8t0ke5NcNM/+85PsT7K7Wy4Yf6mSpL76fIfqBuCjwGnAPuCbSbZX1V1DXa+rqgtXoUZJ0hL1uXI/GdhbVfdU1S+AzwBnr25ZkqSV6BPuxwP3DWzv69qGvS7JniTXJ9kyluokScsyridUvwjMVNWLgZuAa+frlGRrkp1Jdu7fv39Mh5YkDesT7vcDg1fim7u2x1TVgap6tNu8GjhpvoGqaltVzVbV7NTU1HLqlST10CfcvwmcmOS5SY4CzgW2D3ZIctzA5lnA3eMrUZK0VCNfLVNVB5NcCHwZ2AB8rKruTHIpsLOqtgPvTHIWcBB4CDh/FWuWJI0wMtwBqmoHsGOo7ZKB9YuBi8dbmiRpuXyHqiQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDWoV7gnOT3Jd5PsTXLRPPuPTnJdt/+2JDPjLlSS1N/IcE+yAfgo8Brg+cAbkjx/qNtbgB9W1QnAh4DLx12oJKm/PlfuJwN7q+qeqvoF8Bng7KE+ZwPXduvXA69OkvGVKUlaij7hfjxw38D2vq5t3j5VdRB4GDh2HAVKkpbuyMN5sCRbga3d5iNJDgA/OJw1LNMmJlxn+t/omnitPa2XOmH91LrsOpcwv8al+XO6mBWe7+f06dQn3O8Htgxsb+7a5uuzL8mRwNOAA8MDVdU2YNuh7SQ7q2q2T6GTtF7qhPVT63qpE9ZPreulTlg/ta6XOufT57bMN4ETkzw3yVHAucD2oT7bgfO69XOAr1RVja9MSdJSjLxyr6qDSS4EvgxsAD5WVXcmuRTYWVXbgWuATybZCzzE3H8AkqQJ6XXPvap2ADuG2i4ZWP858PplHH/b6C5rwnqpE9ZPreulTlg/ta6XOmH91Lpe6vw18e6JJLXHjx+QpAatKNx7fCzB+Un2J9ndLRd07c9JcnvXdmeStw485qQk3+rGvOLQm6GSPCPJTUm+1/359EnVmeRJSW5I8p2u/bJRY034nN7SjXnoMc/s2pf9sRGrcE6fMtB3d5IfJPnwJM/pwP6nJtmX5MqBtjUzTxeqcy3O04Vq7drWzDxdqM7VnKdjV1XLWph7cvX7wPOAo4A7gOcP9TkfuHKexx4FHN2tPxm4F/jNbvsbwClAgBuB13TtHwAu6tYvAi6fVJ3Ak4DfH+jztYE65x1rwuf0FmB2nse8HbiqWz8XuG6SdQ712wX83iTP6cD+jwD/NNhnLc3Thepci/N0xDldM/N0sTpXY56uxrKSK/c+H0swr6r6RVU92m0eTfcbRJLjgKdW1a01d7Y+Afxh12/wIw6uHWg/7HVW1U+r6uZDfYDbmXv9/0qNvdYRlvuxEataZ5LfAp7JXBit1LJr7Wo5CXgW8C8DbWtqni5U51qcpwvVOsJhn6d96hzzPB27lYR7n48lAHhdkj1Jrk/y2JuhkmxJsqcb4/Kq+q/u8fsWGPNZVfVAt/7fzJ30SdXJwP5jgD8A/m3UWBOu9R+6XxX/YuAfxnI/NmJVzymPX53VqLFWs9YkRwB/C/zpPGOumXm6SJ2PWSvztEeta2Ke9jmnjHeejt1qP6H6RWCmql4M3MTj//tSVfd17ScA5yXp+4+A7mSO82U+y6ozc+/G/TRwRVXdM2qsCdb6x1X1IuB3u+WNY65pXHUeci5z53XkWKtc69uBHVW1b8FHLuIwztNF61xj83SxWtfSPO3zsz/c83RJVhLuIz+WoKoODPwKfjVw0vAg3VXbt5n7Yd7Pr/7aODjm/3S/Dh/6tfjBCdZ5yDbge1X14aWMdbhrrar7uz9/wtz9w5OHj5dFPjbicNXZ1fES4Miq2rWUsVap1lcAFya5F/gb4E2Ze1Jyrc3Theo8ZC3N0wVrXWPzdNFzugrzdPyWepP+0MLcG6DuAZ7L409WvGCoz3ED638E3Nqtbwae2K0/HfgP4EXd9vATVWd07X/Nrz5R9YEJ1/mXwOeAI/qMNalauzE3de0bmbtn+dZu+x386hNVn53kOe3aLgPevxbO6VCf81n8CdWJzdMRda6pebpQrWttni52Tldjnq7GsrIHwxndP87vA3/etV0KnNWt/xVwZ3dSbwZ+u2s/DdjTte8Btg6MOcvc1dz3gSt5/I1WxzJ3v/B7wL8Cz5hUncwFVAF3A7u75YLFxppgrb/B3DP6e7rHfQTY0O17AvDPwF7mwup5k/zZd/vvGT5nkzqni/0DX0vzdKE61+I8XaTWNTVPF/vZr9Y8HffiO1QlqUG+Q1WSGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoP8D7XmAY+8DsicAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(fish_weights, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean: 0.5378740636259465\n",
      "std: 0.004662838260396349\n",
      "min/max: 0.5308641975308642 0.5487804878048781\n"
     ]
    }
   ],
   "source": [
    "mu = fish_weights.mean()\n",
    "std = fish_weights.std()\n",
    "\n",
    "print(\"mean:\", mu)\n",
    "print(\"std:\", std)\n",
    "print(\"min/max:\", fish_weights.min(), fish_weights.max())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "expected weight: 0.5454545454545454\n",
      "50% interval: (0.5252198529486629, 0.5656892379604279)\n",
      "sample size needed: 138.2976\n"
     ]
    }
   ],
   "source": [
    "expected_weight = 30/55  # according to the nutrition guide\n",
    "confidence_interval = stats.norm.interval(0.50, loc=expected_weight,\n",
    "                                          scale=0.03)  # just some guesses\n",
    "\n",
    "print(\"expected weight:\", expected_weight)\n",
    "print(\"50% interval:\", confidence_interval)\n",
    "\n",
    "n_needed = (2*1.96*0.03)**2 / 0.01**2  # 95% confidence on .03 std, < -.01 error\n",
    "print(\"sample size needed:\", n_needed)  # don't have time pho dis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2 - Get Density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.42982101 0.48456228 0.46221063]\n",
      " [0.44935833 0.47446724 0.46221063]\n",
      " [0.44935833 0.4879273  0.46221063]\n",
      " [0.43958967 0.49129231 0.45962845]\n",
      " [0.42982101 0.48456228 0.45962845]\n",
      " [0.45912699 0.47110222 0.45962845]\n",
      " [0.44935833 0.48119727 0.45962845]]\n",
      "container means: [0.44377624 0.4821587  0.4607351 ]\n"
     ]
    }
   ],
   "source": [
    "WEIGHT = mu\n",
    "#WEIGHT = expected_weight\n",
    "volumes = np.array([\n",
    "    [10.4 * 6.1 * 3],\n",
    "    [13 * 8.5 * 5], # [13 * 8.5 * 5.3],\n",
    "    [16 * 9 * 5]\n",
    "]).transpose()\n",
    "weights = np.array([\n",
    "    [44, 46, 46, 45, 44, 47, 46],\n",
    "    [144, 141, 145, 146, 144, 140, 143],\n",
    "    [179, 179, 179, 178, 178, 178, 178]\n",
    "]).transpose()\n",
    "cnumbers = weights / WEIGHT\n",
    "\n",
    "d = cnumbers / volumes\n",
    "container_means = d.mean(axis=0)\n",
    "print(d)\n",
    "print(\"container means:\", container_means)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean: 0.462223347183501\n",
      "std: 0.01723460541842234\n",
      "min/max: 0.42982101403693374 0.49129231456850464\n",
      "interval: (0.45059878248010127, 0.47384791188690073)\n"
     ]
    }
   ],
   "source": [
    "print(\"mean:\", d.mean())\n",
    "print(\"std:\", d.std())\n",
    "print(\"min/max:\", d.min(), d.max())\n",
    "\n",
    "interval = stats.norm.interval(0.5, loc=d.mean(), scale=d.std())\n",
    "print(\"interval:\", interval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3 - Estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "per-container estimates: [4051.18893345 4401.57855743 4206.00465143]\n",
      "overall estimate: 4219.5907141034595\n",
      "guess interval: [4107.82010582 4331.65289798]\n",
      "GUESSES: [4100, 4132, 4164, 4196, 4228, 4260, 4292, 4324]\n"
     ]
    }
   ],
   "source": [
    "tank_volume = 21.5 * 22 * 19.3\n",
    "\n",
    "#tank_interval = np.array(\n",
    "#    stats.norm.interval(.80, loc=tank_volume, scale=200))\n",
    "tank_interval = np.array(\n",
    "    stats.norm.interval(.05, loc=tank_volume, scale=200))\n",
    "\n",
    "pce = container_means * tank_volume\n",
    "print(\"per-container estimates:\", pce)\n",
    "print(\"overall estimate:\", pce.mean())\n",
    "\n",
    "gint = interval * tank_interval\n",
    "print(\"guess interval:\", gint);\n",
    "\n",
    "\n",
    "gint = (4100, 4330)\n",
    "step = (gint[1] - gint[0]) / 7\n",
    "guesses = range(int(gint[0]), int(gint[1])+1, int(step))\n",
    "print(\"GUESSES:\", [i for i in guesses])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4 - Pick Up the TV\n",
    "\n",
    "... or not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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
}