Bootstrap.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Demonstration of Boostrap \n",
    "# Contact: Michael Pyrcz, University of Texas at Austin, Geostatistics Course\n",
    "#\n",
    "# Steps:\n",
    "# 1. Build an initial sample set with $ndata$ samples.\n",
    "# 2. Draw from this initial sample set, with replacement, $ndata$ times to build a new realization of the sample.  \n",
    "#    Repeat this $nreal$ times to make realizations of the sample.\n",
    "# 3. Calculate the statistic of interest for each realization. This demonstration considers with mean and variance.  \n",
    "#    We could have considered any sstatistic including median, 13th percentile, skew etc. \n",
    "# 4. - 6. Quantify and visualize uncertainty with histograms and summary statistics.\n",
    "#\n",
    "# Efron, 1982, The jackknife, the bootstrap, and other resampling plans, Society of Industrial and Applied Math, \n",
    "# CBMS-NSF Monographs, 38."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Imports \n",
    "import scipy\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import math\n",
    "import random as rand\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The initial sample set:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 0.78905746, -0.74536691, -0.3672344 , -1.7200198 , -1.45637476,\n",
       "        0.596527  ,  1.26733395,  0.00919964,  1.65619443, -0.97839249])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 1. Make Dataset with User Specified Number of Data (assume Gaussian N[0,1] for this example, but could be any distribution)\n",
    "nreal = 100                                # number of bootstrap realizations\n",
    "ndata = 10                                 # number of data samples available\n",
    "data = np.zeros((ndata))                    \n",
    "for idata in range(0, ndata):\n",
    "    data[idata] = np.random.normal()\n",
    "print('The initial sample set:')\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The realizations of the sample set:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 1.26733395,  0.596527  ,  0.00919964,  0.00919964,  0.00919964,\n",
       "        -0.74536691, -0.3672344 ,  0.78905746, -1.7200198 , -1.7200198 ,\n",
       "        -0.97839249,  0.596527  , -0.97839249,  1.26733395,  0.596527  ,\n",
       "         0.78905746,  1.65619443, -1.7200198 , -1.45637476, -1.7200198 ,\n",
       "         0.00919964,  1.26733395, -0.74536691, -0.3672344 ,  0.00919964,\n",
       "        -0.74536691,  0.596527  ,  1.65619443, -1.7200198 ,  0.596527  ,\n",
       "        -0.3672344 ,  0.78905746,  0.00919964,  1.26733395, -1.45637476,\n",
       "         0.00919964,  0.00919964,  0.596527  ,  0.00919964, -0.74536691,\n",
       "        -1.45637476,  1.26733395,  0.596527  , -1.45637476,  0.596527  ,\n",
       "         0.596527  ,  1.26733395,  1.26733395,  1.65619443, -0.74536691,\n",
       "         1.65619443,  0.596527  , -1.7200198 , -1.45637476, -0.97839249,\n",
       "        -0.3672344 , -0.97839249, -0.74536691, -1.7200198 , -0.3672344 ,\n",
       "         0.596527  ,  0.596527  ,  1.26733395,  1.26733395,  0.596527  ,\n",
       "         1.26733395,  0.00919964,  1.65619443, -0.74536691,  1.26733395,\n",
       "        -1.7200198 ,  1.65619443, -1.7200198 ,  1.65619443, -1.7200198 ,\n",
       "        -0.74536691, -0.74536691, -1.7200198 , -1.45637476, -0.3672344 ,\n",
       "        -0.74536691,  1.65619443, -1.45637476, -0.97839249, -0.3672344 ,\n",
       "         0.78905746,  0.596527  , -1.7200198 , -0.74536691,  0.596527  ,\n",
       "        -0.74536691, -1.45637476,  1.65619443, -0.74536691,  1.26733395,\n",
       "         1.26733395,  1.26733395, -1.45637476,  1.65619443, -0.97839249],\n",
       "       [ 0.596527  ,  0.596527  , -1.7200198 ,  1.65619443,  0.78905746,\n",
       "         0.78905746, -1.7200198 ,  1.26733395,  1.26733395,  0.78905746,\n",
       "         1.65619443, -0.74536691,  1.65619443, -1.7200198 , -0.74536691,\n",
       "        -0.3672344 , -1.45637476, -1.45637476, -0.74536691, -0.74536691,\n",
       "        -1.7200198 ,  0.596527  , -0.3672344 , -0.74536691,  0.596527  ,\n",
       "        -0.74536691, -0.3672344 ,  0.596527  , -0.3672344 , -0.74536691,\n",
       "        -0.3672344 , -0.3672344 ,  1.65619443,  0.596527  , -1.45637476,\n",
       "         0.78905746,  0.00919964, -0.3672344 ,  0.596527  , -1.7200198 ,\n",
       "         0.00919964, -1.7200198 ,  1.26733395, -1.45637476,  0.78905746,\n",
       "        -0.74536691, -0.3672344 ,  0.596527  ,  0.596527  ,  1.65619443,\n",
       "         1.65619443, -0.74536691,  0.00919964, -1.45637476, -0.74536691,\n",
       "        -1.7200198 , -1.45637476, -0.3672344 , -1.45637476,  1.65619443,\n",
       "         0.596527  ,  0.78905746, -0.74536691, -1.45637476, -0.97839249,\n",
       "         0.78905746,  1.65619443, -1.7200198 ,  1.26733395, -1.45637476,\n",
       "        -1.7200198 , -0.3672344 , -0.97839249, -0.97839249,  1.26733395,\n",
       "         0.78905746, -0.3672344 ,  1.65619443,  0.596527  , -0.74536691,\n",
       "        -0.97839249,  0.78905746, -1.45637476, -1.45637476,  0.78905746,\n",
       "         0.596527  , -0.97839249, -0.97839249, -1.45637476,  0.00919964,\n",
       "        -0.3672344 , -0.74536691, -0.97839249,  0.78905746,  1.65619443,\n",
       "        -0.97839249,  1.65619443,  1.26733395,  0.00919964, -0.74536691],\n",
       "       [-0.3672344 , -0.97839249,  0.596527  , -1.7200198 , -1.7200198 ,\n",
       "         0.78905746,  0.00919964,  0.596527  ,  1.26733395, -0.74536691,\n",
       "         1.26733395, -1.7200198 , -1.45637476,  0.596527  , -0.97839249,\n",
       "         0.78905746, -0.3672344 , -1.45637476, -0.3672344 , -0.3672344 ,\n",
       "         0.596527  ,  0.596527  ,  1.26733395, -0.74536691, -0.74536691,\n",
       "         0.00919964,  0.00919964,  0.00919964,  1.65619443, -0.3672344 ,\n",
       "         0.596527  , -0.3672344 ,  1.26733395, -0.74536691, -1.45637476,\n",
       "        -1.7200198 , -1.7200198 , -1.45637476, -1.45637476,  1.65619443,\n",
       "        -0.97839249,  0.596527  ,  0.78905746, -0.97839249, -0.3672344 ,\n",
       "        -0.3672344 ,  1.65619443, -0.74536691,  1.26733395,  0.00919964,\n",
       "         0.78905746,  0.00919964, -1.7200198 ,  1.65619443, -1.7200198 ,\n",
       "        -0.3672344 ,  1.26733395, -0.97839249,  0.596527  ,  0.596527  ,\n",
       "         0.78905746, -0.3672344 , -0.97839249,  0.596527  , -0.74536691,\n",
       "         0.00919964,  1.26733395,  1.65619443,  1.65619443, -0.74536691,\n",
       "        -0.74536691, -1.7200198 ,  0.596527  ,  0.596527  , -1.7200198 ,\n",
       "        -0.97839249, -1.45637476,  1.26733395, -1.45637476, -0.97839249,\n",
       "         1.26733395,  0.00919964,  0.78905746, -1.7200198 ,  0.596527  ,\n",
       "         1.65619443,  0.596527  ,  1.65619443,  0.78905746,  1.26733395,\n",
       "         1.26733395, -0.74536691, -1.7200198 ,  0.596527  , -0.3672344 ,\n",
       "        -0.97839249,  0.00919964, -1.7200198 ,  1.65619443,  0.78905746],\n",
       "       [ 1.65619443,  0.596527  ,  1.26733395, -0.74536691, -0.97839249,\n",
       "        -0.74536691,  0.00919964,  0.00919964, -0.97839249,  0.596527  ,\n",
       "        -1.7200198 ,  0.78905746,  0.78905746,  0.596527  ,  0.00919964,\n",
       "         1.26733395,  0.00919964, -0.97839249,  1.26733395,  0.596527  ,\n",
       "         0.596527  , -0.74536691,  1.65619443,  0.596527  ,  1.26733395,\n",
       "        -1.45637476, -0.3672344 ,  0.78905746, -0.74536691, -1.7200198 ,\n",
       "         1.26733395, -0.3672344 ,  0.00919964, -0.74536691, -0.97839249,\n",
       "        -0.74536691, -1.45637476, -0.3672344 ,  0.78905746,  0.78905746,\n",
       "        -1.7200198 ,  0.596527  , -1.7200198 ,  1.65619443,  1.26733395,\n",
       "        -1.7200198 ,  0.596527  ,  0.00919964,  1.65619443, -1.45637476,\n",
       "        -1.45637476, -0.97839249, -1.7200198 ,  1.65619443, -0.3672344 ,\n",
       "        -0.3672344 ,  0.78905746,  0.596527  ,  0.596527  ,  0.78905746,\n",
       "        -1.45637476,  1.65619443, -0.97839249, -0.74536691, -0.74536691,\n",
       "         1.65619443, -0.74536691,  0.78905746,  1.26733395, -1.7200198 ,\n",
       "         0.00919964, -0.74536691, -0.74536691,  0.00919964,  0.00919964,\n",
       "        -0.3672344 ,  0.00919964,  1.65619443, -0.74536691,  0.596527  ,\n",
       "         1.26733395,  0.78905746,  0.00919964,  0.00919964,  0.00919964,\n",
       "        -1.7200198 ,  0.596527  ,  0.78905746, -0.3672344 , -1.45637476,\n",
       "         0.78905746, -1.7200198 , -0.74536691,  1.26733395, -0.3672344 ,\n",
       "         1.26733395, -1.7200198 , -1.45637476, -0.74536691,  1.65619443],\n",
       "       [-0.74536691, -0.3672344 , -0.3672344 , -0.74536691, -1.45637476,\n",
       "        -1.45637476,  0.596527  ,  0.00919964, -0.3672344 ,  0.00919964,\n",
       "         1.26733395, -1.45637476, -0.74536691, -1.7200198 ,  0.00919964,\n",
       "        -1.45637476,  0.78905746,  0.596527  , -1.7200198 ,  1.65619443,\n",
       "         0.78905746, -1.45637476,  0.596527  ,  0.596527  ,  0.596527  ,\n",
       "         1.26733395,  0.00919964,  1.26733395, -0.74536691,  1.65619443,\n",
       "        -0.3672344 ,  1.65619443, -0.3672344 , -1.7200198 , -0.3672344 ,\n",
       "        -0.3672344 , -0.74536691,  0.596527  ,  0.00919964,  0.596527  ,\n",
       "         0.00919964, -1.45637476,  0.596527  , -1.45637476, -0.3672344 ,\n",
       "        -0.74536691, -0.3672344 ,  1.65619443, -1.45637476,  0.78905746,\n",
       "         0.78905746,  0.00919964,  1.65619443,  0.78905746, -0.97839249,\n",
       "        -1.45637476,  0.00919964, -0.97839249,  0.78905746,  0.596527  ,\n",
       "        -0.3672344 , -1.7200198 , -0.74536691,  1.26733395, -0.74536691,\n",
       "         0.00919964,  0.78905746,  0.00919964, -1.7200198 , -0.74536691,\n",
       "        -1.45637476, -0.3672344 ,  0.78905746, -0.74536691, -0.74536691,\n",
       "        -1.45637476,  0.78905746, -1.7200198 ,  0.78905746,  1.65619443,\n",
       "         1.65619443,  1.65619443, -1.7200198 , -0.3672344 ,  0.596527  ,\n",
       "         0.78905746,  0.78905746, -0.97839249,  0.00919964,  0.00919964,\n",
       "         0.78905746, -1.45637476, -0.3672344 ,  0.00919964,  1.65619443,\n",
       "        -1.45637476, -1.7200198 ,  0.00919964, -0.3672344 ,  0.596527  ],\n",
       "       [-0.97839249,  1.65619443, -1.7200198 ,  0.78905746,  1.26733395,\n",
       "        -1.7200198 ,  1.65619443, -0.74536691, -1.7200198 ,  0.00919964,\n",
       "        -0.3672344 , -0.97839249, -0.97839249,  1.26733395, -0.3672344 ,\n",
       "         0.596527  , -1.45637476, -0.3672344 , -1.45637476, -1.45637476,\n",
       "         1.26733395, -0.74536691,  0.00919964,  1.26733395,  0.596527  ,\n",
       "        -0.3672344 ,  0.78905746, -1.45637476, -0.97839249, -0.97839249,\n",
       "         0.00919964, -0.3672344 , -1.7200198 ,  0.596527  ,  0.78905746,\n",
       "        -0.74536691, -0.3672344 ,  1.65619443,  0.00919964,  0.78905746,\n",
       "         0.00919964,  1.65619443, -1.45637476,  0.00919964,  0.596527  ,\n",
       "         0.00919964,  1.65619443, -0.74536691,  1.65619443, -0.3672344 ,\n",
       "        -0.74536691, -0.74536691,  0.00919964, -0.3672344 ,  1.26733395,\n",
       "        -0.74536691,  1.26733395,  0.00919964,  0.596527  , -0.97839249,\n",
       "        -0.97839249,  1.26733395,  1.26733395, -0.74536691, -0.97839249,\n",
       "        -1.7200198 , -1.7200198 , -1.45637476,  1.26733395, -0.74536691,\n",
       "        -0.97839249,  0.00919964, -1.7200198 , -0.74536691,  0.00919964,\n",
       "        -1.45637476, -0.3672344 ,  0.00919964,  0.596527  ,  1.65619443,\n",
       "         0.78905746,  1.26733395,  1.65619443,  0.78905746,  1.65619443,\n",
       "        -0.3672344 ,  1.65619443,  1.65619443,  1.65619443, -0.3672344 ,\n",
       "        -1.45637476, -0.74536691,  1.26733395, -1.45637476, -0.97839249,\n",
       "         1.65619443,  0.596527  ,  1.26733395, -1.7200198 ,  0.596527  ],\n",
       "       [-1.7200198 ,  0.596527  ,  0.00919964,  1.26733395,  1.26733395,\n",
       "        -1.7200198 , -0.97839249, -0.3672344 ,  0.00919964,  0.00919964,\n",
       "        -0.3672344 , -1.7200198 , -0.74536691,  0.596527  , -1.45637476,\n",
       "         0.00919964,  0.00919964,  0.78905746,  0.596527  ,  1.65619443,\n",
       "         0.78905746, -0.97839249,  0.596527  , -1.45637476,  0.00919964,\n",
       "         0.596527  ,  0.596527  ,  0.78905746, -0.74536691, -0.74536691,\n",
       "         1.65619443,  0.00919964, -0.3672344 ,  1.26733395,  0.00919964,\n",
       "        -0.74536691,  1.65619443,  0.78905746, -1.7200198 , -1.45637476,\n",
       "        -1.7200198 , -1.7200198 ,  0.78905746, -1.45637476, -1.45637476,\n",
       "         0.00919964,  1.65619443, -0.3672344 , -0.3672344 ,  0.78905746,\n",
       "        -0.97839249,  1.65619443,  1.26733395, -0.74536691, -1.45637476,\n",
       "        -1.7200198 , -1.7200198 ,  1.26733395,  0.78905746,  0.596527  ,\n",
       "        -0.3672344 , -1.45637476, -0.97839249, -1.45637476, -1.7200198 ,\n",
       "         1.65619443,  0.596527  , -0.3672344 ,  0.00919964,  0.00919964,\n",
       "         1.26733395,  0.596527  ,  0.78905746,  1.65619443, -0.74536691,\n",
       "         0.78905746,  1.65619443,  1.65619443,  0.78905746, -0.74536691,\n",
       "        -0.3672344 ,  1.65619443, -1.7200198 , -1.45637476,  1.65619443,\n",
       "        -0.3672344 , -0.74536691,  0.00919964,  0.78905746,  0.596527  ,\n",
       "         1.65619443,  0.78905746, -0.97839249, -0.97839249, -0.74536691,\n",
       "         1.65619443, -1.7200198 , -0.3672344 , -1.45637476,  0.596527  ],\n",
       "       [-0.74536691,  0.596527  ,  0.596527  ,  1.26733395, -0.3672344 ,\n",
       "        -0.74536691,  1.65619443,  0.78905746, -0.3672344 , -0.74536691,\n",
       "         1.65619443,  1.26733395, -0.97839249,  0.78905746,  1.65619443,\n",
       "        -0.97839249,  1.65619443,  0.00919964, -0.97839249, -1.45637476,\n",
       "         0.00919964,  1.65619443,  0.78905746,  1.65619443, -0.3672344 ,\n",
       "        -0.3672344 , -0.74536691,  0.00919964, -1.45637476,  1.65619443,\n",
       "        -0.74536691, -1.7200198 , -0.97839249, -0.74536691, -0.3672344 ,\n",
       "         1.65619443, -0.3672344 ,  0.78905746,  1.65619443,  1.65619443,\n",
       "         0.78905746, -1.7200198 ,  1.26733395,  1.26733395,  1.26733395,\n",
       "        -1.45637476, -0.74536691, -1.7200198 ,  1.65619443,  0.00919964,\n",
       "         1.26733395,  0.596527  , -0.3672344 ,  0.00919964, -0.97839249,\n",
       "         0.596527  , -0.97839249,  0.78905746,  0.00919964,  1.26733395,\n",
       "        -0.97839249, -0.74536691, -1.45637476, -0.74536691, -0.3672344 ,\n",
       "         0.78905746,  1.26733395,  1.65619443, -0.3672344 , -1.7200198 ,\n",
       "         1.26733395,  1.65619443,  0.78905746,  0.78905746, -0.74536691,\n",
       "         0.596527  , -0.3672344 ,  0.596527  , -1.45637476, -0.97839249,\n",
       "         0.00919964, -0.74536691,  0.596527  , -1.45637476, -0.3672344 ,\n",
       "         0.78905746, -1.7200198 ,  1.26733395, -1.45637476,  0.596527  ,\n",
       "        -1.45637476,  0.00919964, -1.7200198 ,  1.26733395,  1.26733395,\n",
       "         1.65619443,  0.00919964,  0.78905746, -1.45637476, -1.45637476],\n",
       "       [-1.45637476, -0.74536691, -0.97839249, -0.74536691, -0.74536691,\n",
       "        -0.74536691, -1.45637476,  1.65619443,  0.00919964, -1.45637476,\n",
       "         0.596527  , -0.74536691,  0.00919964, -1.7200198 , -0.74536691,\n",
       "         0.00919964,  1.65619443,  0.596527  ,  0.78905746,  0.596527  ,\n",
       "         1.26733395,  0.78905746,  0.596527  , -1.45637476,  1.26733395,\n",
       "         0.00919964, -1.7200198 , -1.7200198 , -1.45637476,  1.26733395,\n",
       "        -0.74536691, -0.74536691,  1.26733395,  0.78905746,  0.78905746,\n",
       "        -0.97839249, -1.45637476, -0.74536691, -0.3672344 , -1.7200198 ,\n",
       "        -0.74536691, -1.45637476, -1.7200198 , -1.45637476, -0.97839249,\n",
       "         0.00919964, -1.45637476, -1.45637476,  0.596527  ,  0.78905746,\n",
       "         0.596527  , -1.45637476,  1.65619443,  0.78905746, -0.97839249,\n",
       "        -1.7200198 , -0.3672344 ,  0.596527  ,  1.26733395,  0.78905746,\n",
       "        -1.7200198 , -0.97839249,  1.26733395,  1.65619443, -0.74536691,\n",
       "        -1.7200198 ,  1.26733395,  1.65619443,  1.26733395,  0.00919964,\n",
       "        -0.97839249, -1.7200198 , -1.45637476,  0.596527  ,  1.65619443,\n",
       "         0.78905746, -1.7200198 ,  0.78905746, -1.7200198 , -1.7200198 ,\n",
       "        -0.97839249, -0.3672344 ,  1.26733395,  1.65619443, -1.45637476,\n",
       "         0.596527  ,  0.596527  ,  1.65619443,  0.596527  , -1.45637476,\n",
       "        -0.74536691, -1.45637476,  0.78905746, -1.7200198 , -1.7200198 ,\n",
       "         0.78905746, -0.3672344 , -0.3672344 ,  0.00919964,  0.596527  ],\n",
       "       [ 0.78905746, -0.74536691,  0.78905746, -0.74536691,  1.26733395,\n",
       "        -0.3672344 ,  0.00919964, -0.74536691, -1.7200198 ,  0.596527  ,\n",
       "         0.596527  , -0.3672344 , -0.3672344 ,  0.00919964, -1.45637476,\n",
       "        -0.74536691, -1.45637476, -1.7200198 , -0.74536691,  0.78905746,\n",
       "        -0.74536691, -0.97839249, -0.74536691,  0.78905746,  1.26733395,\n",
       "        -1.7200198 ,  0.00919964, -0.74536691, -0.97839249, -0.74536691,\n",
       "        -0.97839249,  1.65619443, -0.74536691,  0.78905746,  0.596527  ,\n",
       "         1.65619443,  1.26733395, -0.3672344 ,  1.26733395,  1.65619443,\n",
       "        -1.45637476, -0.74536691, -1.45637476,  1.65619443, -0.97839249,\n",
       "        -1.7200198 ,  0.78905746,  0.00919964,  1.65619443, -1.45637476,\n",
       "        -0.97839249,  1.65619443, -0.3672344 ,  0.596527  ,  1.26733395,\n",
       "        -1.7200198 , -1.45637476,  0.596527  , -0.97839249, -0.97839249,\n",
       "        -1.45637476, -1.7200198 , -1.7200198 ,  0.78905746,  0.00919964,\n",
       "         1.26733395,  1.65619443,  0.78905746, -1.45637476,  0.596527  ,\n",
       "        -0.74536691, -1.45637476, -0.74536691,  0.00919964,  0.596527  ,\n",
       "         0.78905746, -0.74536691,  1.26733395, -0.97839249, -0.74536691,\n",
       "        -1.7200198 ,  1.26733395,  1.26733395,  1.26733395, -1.7200198 ,\n",
       "        -1.45637476,  0.596527  , -0.3672344 ,  0.596527  , -0.97839249,\n",
       "         1.26733395, -0.97839249, -1.7200198 ,  0.78905746,  1.65619443,\n",
       "         0.00919964, -0.74536691, -0.3672344 ,  0.596527  , -1.45637476]])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 2. Perform ndata random draws with replacement, nreal times.  Here we aquire the nreal realizations of the distribution of \n",
    "# ndata, samples.\n",
    "draw = np.zeros((ndata,nreal)) \n",
    "for ireal in range(0, nreal):\n",
    "    for isample in range(0, ndata):\n",
    "        draw[isample,ireal] = rand.choice(data)\n",
    "print('The realizations of the sample set:')\n",
    "draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The statistics of interest for each realization of the sample set:\n",
      "[ -1.70364243e-01   1.80246874e-01  -1.51782181e-01   2.87631993e-02\n",
      "  -6.67129427e-02  -6.66700148e-01  -5.85506667e-02   3.25860136e-01\n",
      "  -4.31985353e-01  -2.65741800e-01   3.60722967e-01  -5.07985667e-01\n",
      "  -3.79506890e-01  -3.75534229e-03  -3.47798951e-01  -8.69934163e-03\n",
      "   1.03968134e-01  -5.70710492e-01  -4.81621163e-01  -4.50870316e-02\n",
      "   2.85884938e-01   1.74628855e-04   3.65339827e-01   1.34922096e-02\n",
      "   4.49738081e-01  -3.51933695e-01  -1.19014513e-01   1.19480810e-01\n",
      "  -7.53669499e-01  -1.25497597e-02  -4.15744790e-03   1.76321649e-02\n",
      "   3.10136079e-03   1.34971628e-01  -3.89814402e-01  -1.19110145e-01\n",
      "  -3.17067738e-01   1.12391847e-01   7.93082788e-02   1.50144393e-01\n",
      "  -7.25989215e-01  -4.70159347e-01  -1.04695233e-01  -3.67134385e-01\n",
      "   3.69150812e-02  -6.13025666e-01   4.68529123e-01  -1.49590812e-01\n",
      "   8.91775094e-01   1.64152470e-03   2.59583808e-01   5.98341087e-02\n",
      "  -1.29640612e-01   1.47087958e-01  -5.66789793e-01  -9.58699709e-01\n",
      "  -3.62386372e-01   7.85785767e-02   4.89442446e-02   3.96720492e-01\n",
      "  -5.34191165e-01  -2.67829533e-01  -3.80030400e-01   4.27596535e-02\n",
      "  -6.41978016e-01   4.00353134e-01   6.04378809e-01   4.66846331e-01\n",
      "   2.44573398e-01  -5.25025486e-01  -5.80006562e-01  -2.45813457e-01\n",
      "  -4.40184130e-01   2.84377331e-01  -2.13768566e-01  -1.25098650e-01\n",
      "  -3.31438005e-01   5.45799568e-01  -5.04173457e-01  -2.37122403e-01\n",
      "   1.99713337e-02   7.97796443e-01  -7.67142716e-02  -3.71298551e-01\n",
      "   1.39283659e-01   1.30555743e-01   1.98410770e-01   2.99013516e-01\n",
      "   4.11212152e-02  -1.18306218e-01   9.98259495e-02  -8.50538020e-01\n",
      "  -4.51685986e-01  -1.81644516e-02   3.32500319e-01   4.88834855e-01\n",
      "  -2.73420606e-01  -2.40154755e-01  -1.81805549e-01   1.94850973e-02] [ 1.23200839  0.64005696  0.97192761  1.17504292  1.21237765  0.71954446\n",
      "  1.19793851  0.61345089  1.14487049  0.68190067  1.22941512  1.02453174\n",
      "  0.81642093  1.37777637  0.83346641  0.69753278  1.52415707  0.93283847\n",
      "  0.96687817  1.46908881  0.78660902  1.08300904  0.59506179  1.11696552\n",
      "  0.45538982  0.71839372  0.49766721  1.14426632  0.79794062  1.31636614\n",
      "  0.71778504  1.01590317  1.06963378  0.96121931  0.76716562  1.16840806\n",
      "  1.13697822  0.76547556  1.0479284   1.81981099  0.7024077   1.65233372\n",
      "  1.52185632  1.73066925  0.88381974  0.58820621  1.15646575  1.08068553\n",
      "  1.03359846  0.95434452  1.26104286  1.03514399  1.60693243  1.18032516\n",
      "  0.95986999  0.60798482  1.1771616   0.58395557  0.99317687  0.71373654\n",
      "  0.78654475  1.42958934  1.24365697  1.28241933  0.34646916  1.42577332\n",
      "  1.11442092  1.52521297  1.4339373   0.89103106  1.09118902  1.41009043\n",
      "  1.0392418   0.80801207  1.18318108  0.85795341  0.88623873  1.5357716\n",
      "  1.02493777  1.19540825  1.21373632  0.70121599  1.70339318  1.40068489\n",
      "  1.19755636  1.05012737  0.92224375  1.45715137  0.96155552  0.78680422\n",
      "  1.25888814  0.53340649  1.44118657  1.14024429  1.507666    1.3709843\n",
      "  1.35233781  1.08633302  1.34023971  1.05400514]\n"
     ]
    }
   ],
   "source": [
    "# 3. Calculate summary statistic - average and variance\n",
    "mean = np.zeros((nreal))\n",
    "var = np.zeros((nreal))\n",
    "for ireal in range(0, nreal):\n",
    "    s = 0.0\n",
    "    ss = 0.0\n",
    "    count = 0.0\n",
    "    for isample in range(0, ndata):\n",
    "        s = s + draw[isample,ireal]\n",
    "        ss = ss + draw[isample,ireal] * draw[isample,ireal]\n",
    "        count = count + 1.0\n",
    "    mean[ireal] = s / count \n",
    "    var[ireal] = ss / count - mean[ireal] * mean[ireal]\n",
    "print('The statistics of interest for each realization of the sample set:')\n",
    "print(mean,var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHRtJREFUeJzt3XucHGWd7/HPlxBEIAQhQ7iEEFQOS5SLOAbloIQ9wknA\nGDnqmhxEQTCiZhX3yDGuvoBd1z24rusF0BAxi0EBXSUQIYDgCrjyQhIQgXDREMIm4ZJwDQGWGPid\nP+oZLTo9M0/PdHX3ZL7v16tfU/Vcun5d1TO/qaeqn1ZEYGZm1p+t2h2AmZkNDU4YZmaWxQnDzMyy\nOGGYmVkWJwwzM8vihGFmZlmcMMyGMEknSvqPFm/zOEmrJG2Q9KZWbtvaywnDkLRS0gvpD8BTkq6S\ntFcTnneCpJC0dWb7kPT6wW43l6TJkl5Or/tZSfdLOqlV2x/C/hmYHRE7RMRvayvTcVxbPu6SRqYy\nf/BrCHPCsB7TImIHYHfgMeCcNsezmdzE06CH0+veEfgM8F1J+1WwnS3J3sCyfto8BUwtrU9NZTaE\nOWHYK0TEfwE/ASb2lEkaLWmBpHWSHpL0RUlbpbqt0vpD6T/IBZJGp643pZ9Pp//i3ybp9ZJulPSM\npMcl/Sg9T0/b36W2H0hnAKslfU7So8C/SnqNpCtTLE+l5XGlWG+Q9P8k3SppvaQrJO2c8bojIhYD\nTwIHlp7vLyRdJ+nJdAbyV6W6YyTdk85O1kj6bCrPifEfJN2cXuvPJO0i6Ycp5iWSJpTah6RPSVqR\n9tlXe/Z/rYHEW+c56h5TSa+StAEYkY7TA33s0ouAD5XWPwQsqNnOaEnfk/RIiucfJI1Ida+T9O+S\nnkiv+YeSdir1XSnps5LuTO+lH0nato94rBkiwo9h/gBWAu9My9sB3wcWlOoXAFcAo4AJwO+Bk1Pd\nR4DlwGuBHYDLgItS3QQggK1Lz3UJ8AWKf1a2BQ4v1QXw+tL6ZGAT8BXgVcCrgV2A96Y4RwH/Blxe\n6nMDsAZ4I7A98FPgB7287snA6rS8FfBu4GXgTalse2AVcBKwNfAm4HFgYqp/BHh7Wn4NcEhazolx\nOfA6YDRwT9qn70zbWQD8a81++SWwMzA+tT0l1Z0I/Mdg4q2zX3o9pvWOU53+kfb/Y8BOaVuPpbIo\ntVsInJ/i3hW4FfhYqns9cFQ67l0U/3x8o+Y9eyuwR9ov9wKntvt3aUt/tD0AP9r/SL98G4CngT8C\nDwMHpLoRwMaePzqp7GPADWn5F8AnSnX7pefYmvoJYwEwDxhXJ456CWMjsG0fsR8MPFVavwE4u7Q+\nMT3HiDp9J1MkiKeBF4GXgNNK9R8AflXT53zgzLT8n2lf7NjP/q0X4xdK618Dri6tTwPuqNkvU0rr\nnwB+kZZP5M8Jo1nx9npM6x2n3o4jcEHa3qnAd1NZpDZj0z5/danfTOCXvTzne4Df1rxnP1ha/ydg\nbrt/l7b0h4ekrMd7ImIniv/6ZwM3StoNGAOMBB4qtX0I2DMt71GnbmuKPwj1/F9AwK2Slkn6SD9x\nrYtimAwASdtJOj8Nl6yn+M9zp56hjGRVTTwj0+uo5+H0uncEvgX8Zalub+BQSU/3PIDjgd1S/XuB\nY4CH0jDb2xqI8bHS8gt11neoibP2Ne1R57UMKN46Gj2mvVlAMRS12XBUinUk8Egp1vMpzjSQNFbS\npWmoaj3wAzY/ho+Wlp9n831mTeaEYa8QES9FxGUU/20fTjGk8UeKX/Ae4ymGfaA4G6mt20TxB3Cz\nO2Ii4tGI+GhE7EHx3+e31fedUbXP8X8o/uM9NCJ2BN6RylVqU77Da3yK//E+tkFEvAh8DjhA0ntS\n8SrgxojYqfTYISI+nvosiYjpFH/kLgd+3ECMjap9TQ/XaTPQeGv1dUwb8SuKmyjGArW3/q6iOMMY\nU4p1x4h4Q6r/R4pjf0Dahx9kcPvPmsAJw15BhekU4873RsRLFH9YvixplKS9gb+h+I8PimsSn5G0\nj6QdKH7RfxQRm4B1FEM+ry09//tLF4Cfovij8HJaf6zcthejKP4DfzpdzD6zTpsPSpooaTvg74Gf\npNfRp4jYSDE8dEYquhL4b5JOUHFb6EhJb5G0v6RtJB0vaXRE/BFYX3odOTE26vR0MX0v4NPAj+q0\nGWi8tfo6ptmiGCuaBrw7LZfrHgF+DnxN0o7pQvvrJB2RmoyiGCZ9RtKewOmNbNuq4YRhPX6W7oBZ\nD3wZ+HBE9Nw6+dfAc8AKiv8ULwbmp7r5FHfE3AQ8CPxXak9EPJ+e69dp2OGtwFuA36RtLQI+HREr\n0nOdBXw/tf3T3T01vkFx8ftx4BbgmjptLgIupBiy2Bb4VAP7YT4wXtK0iHgWOBqYQfFf96P8+QI8\nwAnAyjRkcirF8E9ujI26ArgNuAO4CvhebYNBxFur12PaqIhYVnof1foQsA3FRf+nKO7O2z3V/R1w\nCPAMxeu9bCDbt+ZSTeI3G9Ik3UBxV9QF7Y6lWVR82G3fiFje7lhsePMZhpmZZXHCMDOzLB6SMjOz\nLD7DMDOzLFVM5tY2Y8aMiQkTJrQ7DDOzIeO22257PCK6ctpuUQljwoQJLF26tN1hmJkNGZIe6r9V\nwUNSZmaWxQnDzMyyOGGYmVkWJwwzM8vihGFmZlmcMMzMLIsThpmZZXHCMDOzLE4YZmaWZYv6pLdZ\np5ow56oB91159rFNjMRs4HyGYWZmWZwwzMwsixOGmZllccIwM7MsThhmZpalsrukJM0H3gWsjYg3\n1qk/HTi+FMf+QFdEPClpJfAs8BKwKSK6q4rTzMzyVHmGcSEwpbfKiPhqRBwcEQcDnwdujIgnS02O\nTPVOFmZmHaCyhBERNwFP9tuwMBO4pKpYzMxs8Np+DUPSdhRnIj8tFQdwvaTbJM1qT2RmZlbWCZ/0\nngb8umY46vCIWCNpV+A6SfelM5bNpIQyC2D8+PHVR2tmNky1/QwDmEHNcFRErEk/1wILgUm9dY6I\neRHRHRHdXV1dlQZqZjactTVhSBoNHAFcUSrbXtKonmXgaODu9kRoZmY9qryt9hJgMjBG0mrgTGAk\nQETMTc2OA34eEc+Vuo4FFkrqie/iiLimqjjNzCxPZQkjImZmtLmQ4vbbctkK4KBqojIzs4HqhGsY\nZmY2BDhhmJlZFicMMzPL4oRhZmZZnDDMzCyLE4aZmWVxwjAzsyxOGGZmlsUJw8zMsjhhmJlZFicM\nMzPL4oRhZmZZnDDMzCyLE4aZmWVxwjAzsyxOGGZmlsUJw8zMsjhhmJlZFicMMzPL4oRhZmZZKksY\nkuZLWivp7l7qJ0t6RtId6XFGqW6KpPslLZc0p6oYzcwsX5VnGBcCU/pp86uIODg9/h5A0gjgPGAq\nMBGYKWlihXGamVmGyhJGRNwEPDmArpOA5RGxIiI2ApcC05sanJmZNazd1zAOk3SnpKslvSGV7Qms\nKrVZncrqkjRL0lJJS9etW1dlrGZmw1o7E8btwPiIOBA4B7h8IE8SEfMiojsiuru6upoaoJmZ/Vnb\nEkZErI+IDWl5MTBS0hhgDbBXqem4VGZmZm3UtoQhaTdJSsuTUixPAEuAfSXtI2kbYAawqF1xmplZ\nYeuqnljSJcBkYIyk1cCZwEiAiJgLvA/4uKRNwAvAjIgIYJOk2cC1wAhgfkQsqypOMzPLU1nCiIiZ\n/dSfC5zbS91iYHEVcZmZ2cC0+y4pMzMbIpwwzMwsixOGmZllccIwM7MsThhmZpbFCcPMzLI4YZiZ\nWRYnDDMzy+KEYWZmWZwwzMwsixOGmZllccIwM7MsThhmZpbFCcPMzLI4YZiZWRYnDDMzy+KEYWZm\nWSr7xj0za44Jc64acN+VZx/bxEhsuPMZhpmZZaksYUiaL2mtpLt7qT9e0p2S7pJ0s6SDSnUrU/kd\nkpZWFaOZmeWr8gzjQmBKH/UPAkdExAHAl4B5NfVHRsTBEdFdUXxmZtaAyq5hRMRNkib0UX9zafUW\nYFxVsZiZ2eB1yjWMk4GrS+sBXC/pNkmz+uooaZakpZKWrlu3rtIgzcyGs7bfJSXpSIqEcXip+PCI\nWCNpV+A6SfdFxE31+kfEPNJwVnd3d1QesJnZMNXWMwxJBwIXANMj4ome8ohYk36uBRYCk9oToZmZ\n9WhbwpA0HrgMOCEifl8q317SqJ5l4Gig7p1WZmbWOpUNSUm6BJgMjJG0GjgTGAkQEXOBM4BdgG9L\nAtiU7ogaCyxMZVsDF0fENVXFaWZmeaq8S2pmP/WnAKfUKV8BHLR5DzMza6dOuUvKzMw6nBOGmZll\nyUoYkg6oOhAzM+tsuWcY35Z0q6RPSBpdaURmZtaRshJGRLwdOB7YC7hN0sWSjqo0MjMz6yjZ1zAi\n4g/AF4HPAUcA35J0n6T/VVVwZmbWOXKvYRwo6evAvcBfAtMiYv+0/PUK4zMzsw6R+zmMcyim8Pjb\niHihpzAiHpb0xUoiMzOzjpKbMI4FXoiIlwAkbQVsGxHPR8RFlUVnZmYdI/caxvXAq0vr26UyMzMb\nJnITxrYRsaFnJS1vV01IZmbWiXITxnOSDulZkfRm4IU+2puZ2RYm9xrGacC/SXoYELAb8IHKojIz\ns46TlTAiYomkvwD2S0X3R8QfqwvLtmQT5lw14L4rzz62iZGYWSMamd78LcCE1OcQSUTEgkqiMjOz\njpOVMCRdBLwOuAN4KRUH4IRhZjZM5J5hdAMTIyKqDMbMzDpX7l1Sd1Nc6DYzs2Eq9wxjDHCPpFuB\nF3sKI+LdlURlZmYdJzdhnNXoE0uaD7wLWBsRb6xTL+CbwDHA88CJEXF7qpuS6kYAF0TE2Y1u38zM\nmiv3+zBuBFYCI9PyEuD2frpdCEzpo34qsG96zAK+AyBpBHBeqp8IzJQ0MSdOMzOrTu705h8FfgKc\nn4r2BC7vq09E3AQ82UeT6cCCKNwC7CRpd2ASsDwiVkTERuDS1NbMzNoo96L3J4H/DqyHP32Z0q6D\n3PaewKrS+upU1lt5XZJmSVoqaem6desGGZKZmfUmN2G8mP7bB0DS1hSfw2i7iJgXEd0R0d3V1dXu\ncMzMtli5F71vlPS3wKvTd3l/AvjZILe9huI7wnuMS2Ujeyk3M7M2yk0Yc4CTgbuAjwGLKb6BbzAW\nAbMlXQocCjwTEY9IWgfsK2kfikQxA/jfg9yW2aANZg4ssy1B7uSDLwPfTY8ski4BJgNjJK0GzqQ4\neyAi5lIknWOA5RS31Z6U6jZJmg1cS3Fb7fyIWJa7XTMzq0buXFIPUueaRUS8trc+ETGzr+dM04x8\nspe6xRQJxczMOkQjc0n12BZ4P7Bz88MxM7NOlfvBvSdKjzUR8Q3AX0xgZjaM5A5JHVJa3YrijKOR\n79IwM7MhLveP/tdKy5sopgn5q6ZHY2ZmHSv3Lqkjqw7EzMw6W+6Q1N/0VR8R/9KccMzMrFM1cpfU\nWyg+bAcwDbgV+EMVQZmZWefJTRjjgEMi4lkASWcBV0XEB6sKzMzMOkvu5INjgY2l9Y2pzMzMhonc\nM4wFwK2SFqb19wDfryYkMzPrRLl3SX1Z0tXA21PRSRHx2+rCMjOzTpM7JAWwHbA+Ir4JrE6zyZqZ\n2TCR+xWtZwKfAz6fikYCP6gqKDMz6zy5ZxjHAe8GngOIiIeBUVUFZWZmnSc3YWxM05EHgKTtqwvJ\nzMw6UW7C+LGk84GdJH0UuJ4GvkzJzMyGvty7pP45fZf3emA/4IyIuK7SyMzMrKP0mzAkjQCuTxMQ\nOkmYmQ1T/Q5JRcRLwMuSRrcgHjMz61C5n/TeANwl6TrSnVIAEfGpvjpJmgJ8ExgBXBARZ9fUnw4c\nX4plf6ArIp6UtBJ4FngJ2BQR5a+JNTOzFstNGJelR7Y0lHUecBSwGlgiaVFE3NPTJiK+Cnw1tZ8G\nfCYiniw9zZER8Xgj2zUzs2r0mTAkjY+I/4yIgcwbNQlYHhEr0nNdCkwH7uml/UzgkgFsx8zMWqC/\nM4zLgUMAJP00It7bwHPvCawqra8GDq3XUNJ2wBRgdqk4gOslvQScHxHzeuk7C5gFMH78+AbCswlz\nrhpw35VnH9vESKxT+T1iZf1d9FZp+bUVxjEN+HXNcNThEXEwMBX4pKR31OsYEfMiojsiuru6uioM\n0cxseOsvYUQvyznWAHuV1selsnpmUDMcFRFr0s+1wEKKIS4zM2uT/hLGQZLWS3oWODAtr5f0rKT1\n/fRdAuwraR9J21AkhUW1jdLtukcAV5TKtpc0qmcZOBq4O/9lmZlZs/V5DSMiRgz0iSNik6TZwLUU\nt9XOj4hlkk5N9XNT0+OAn0fEc6XuY4GFknpivDgirhloLGZmNni5t9UOSEQsBhbXlM2tWb8QuLCm\nbAVwUJWxmZlZYxr5AiUzMxvGnDDMzCyLE4aZmWVxwjAzsyxOGGZmlsUJw8zMslR6W61Zsw1mbiMz\nGxyfYZiZWRYnDDMzy+KEYWZmWZwwzMwsixOGmZllccIwM7MsThhmZpbFCcPMzLI4YZiZWRYnDDMz\ny+KEYWZmWTyXlA2I53QyG34qPcOQNEXS/ZKWS5pTp36ypGck3ZEeZ+T2NTOz1qrsDEPSCOA84Chg\nNbBE0qKIuKem6a8i4l0D7GtmZi1S5RnGJGB5RKyIiI3ApcD0FvQ1M7MKVJkw9gRWldZXp7Jah0m6\nU9LVkt7QYF8kzZK0VNLSdevWNSNuMzOro913Sd0OjI+IA4FzgMsbfYKImBcR3RHR3dXV1fQAzcys\nUGXCWAPsVVofl8r+JCLWR8SGtLwYGClpTE5fMzNrrSoTxhJgX0n7SNoGmAEsKjeQtJskpeVJKZ4n\ncvqamVlrVXaXVERskjQbuBYYAcyPiGWSTk31c4H3AR+XtAl4AZgREQHU7VtVrGZm1r9KP7iXhpkW\n15TNLS2fC5yb29fMzNqn3Re9zcxsiHDCMDOzLJ5Lyswq0a75xlaefWxbtjsc+AzDzMyyOGGYmVkW\nJwwzM8vihGFmZlmcMMzMLIsThpmZZXHCMDOzLE4YZmaWxQnDzMyyOGGYmVkWJwwzM8viuaSGsHbN\n1WNmw5PPMMzMLIsThpmZZXHCMDOzLE4YZmaWpdKEIWmKpPslLZc0p0798ZLulHSXpJslHVSqW5nK\n75C0tMo4zcysf5XdJSVpBHAecBSwGlgiaVFE3FNq9iBwREQ8JWkqMA84tFR/ZEQ8XlWMZmaWr8oz\njEnA8ohYEREbgUuB6eUGEXFzRDyVVm8BxlUYj5mZDUKVCWNPYFVpfXUq683JwNWl9QCul3SbpFkV\nxGdmZg3oiA/uSTqSImEcXio+PCLWSNoVuE7SfRFxU52+s4BZAOPHj29JvGZmw1GVZxhrgL1K6+NS\n2StIOhC4AJgeEU/0lEfEmvRzLbCQYohrMxExLyK6I6K7q6urieGbmVlZlQljCbCvpH0kbQPMABaV\nG0gaD1wGnBARvy+Vby9pVM8ycDRwd4WxmplZPyobkoqITZJmA9cCI4D5EbFM0qmpfi5wBrAL8G1J\nAJsiohsYCyxMZVsDF0fENVXF2k6eD8qqNBzfX4N9zSvPPrZJkWx5Kr2GERGLgcU1ZXNLy6cAp9Tp\ntwI4qLbczMzax5/0NjOzLE4YZmaWxQnDzMyyOGGYmVkWJwwzM8vihGFmZlmcMMzMLIsThpmZZXHC\nMDOzLE4YZmaWxQnDzMyydMT3YXSC4ThJm5k1V7v+jrRqwkSfYZiZWRYnDDMzy+KEYWZmWZwwzMws\nixOGmZllccIwM7MsThhmZpbFCcPMzLJUmjAkTZF0v6TlkubUqZekb6X6OyUdktvXzMxaq7KEIWkE\ncB4wFZgIzJQ0sabZVGDf9JgFfKeBvmZm1kJVnmFMApZHxIqI2AhcCkyvaTMdWBCFW4CdJO2e2dfM\nzFqoyrmk9gRWldZXA4dmtNkzsy8AkmZRnJ0AbJB0fx8xjQEe7zfy9nKMzeEYm2coxNm0GPWVZjxL\nXZXtx0HGvHduwyE/+WBEzAPm5bSVtDQiuisOaVAcY3M4xuYZCnE6xtaoMmGsAfYqrY9LZTltRmb0\nNTOzFqryGsYSYF9J+0jaBpgBLKppswj4ULpb6q3AMxHxSGZfMzNrocrOMCJik6TZwLXACGB+RCyT\ndGqqnwssBo4BlgPPAyf11bcJYWUNXbWZY2wOx9g8QyFOx9gCioh2x2BmZkOAP+ltZmZZnDDMzCzL\nFpcwJL1f0jJJL0vq9Ra23qYekbSzpOsk/SH9fE0FMfa7DUn7Sbqj9Fgv6bRUd5akNaW6Y9oRY2q3\nUtJdKY6ljfavOkZJe0n6paR70vvi06W6yvbjUJgWJyPG41Nsd0m6WdJBpbq6x70NMU6W9EzpGJ6R\n27eFMZ5eiu9uSS9J2jnVtWQ/Nk1EbFEPYH9gP+AGoLuXNiOAB4DXAtsAvwMmprp/Auak5TnAVyqI\nsaFtpHgfBfZO62cBn614P2bFCKwExgz2NVYVI7A7cEhaHgX8vnSsK9mPfb2/Sm2OAa4GBLwV+E1u\n3xbGeBjwmrQ8tSfGvo57G2KcDFw5kL6tirGm/TTg31u5H5v52OLOMCLi3ojo69Pe0PfUI9OB76fl\n7wPvqSDMRrfxP4AHIuKhCmLpzWD3Q0fsx4h4JCJuT8vPAvdSzCRQpaEwLU6/24mImyPiqbR6C8Xn\noVppMPuiY/ZjjZnAJRXE0RJbXMLI1NuUJABjo/gsCBT/1Y+tYPuNbmMGm7/J/joNF8yvYriH/BgD\nuF7SbSqmaWm0fytiBEDSBOBNwG9KxVXsx77eX/21yenbqhjLTqY4I+rR23FvptwYD0vH8GpJb2iw\nb6tiRNJ2wBTgp6XiVuzHphmSU4NIuh7YrU7VFyLiimZtJyJC0oDuO+4rxka2oeKDi+8GPl8q/g7w\nJYo325eArwEfaVOMh0fEGkm7AtdJui8ibmqgfytiRNIOFL+op0XE+lTclP24pZN0JEXCOLxU3O9x\nb5HbgfERsSFdg7qcYvbrTjQN+HVEPFkq65T9mGVIJoyIeOcgn6KvaUsek7R7RDyShgjWNjtGSY1s\nYypwe0Q8VnruPy1L+i5wZbtijIg16edaSQspTtFvooP2o6SRFMnihxFxWem5m7If6xgK0+LkxIik\nA4ELgKkR8URPeR/HvaUxlpI/EbFY0rcljcnp26oYSzYbKWjRfmya4Tok1dfUI4uAD6flDwNNO2Mp\naWQbm415pj+OPY4D7m5qdIV+Y5S0vaRRPcvA0aVYOmI/ShLwPeDeiPiXmrqq9uNQmBan3+1IGg9c\nBpwQEb8vlfd13Fsd427pGCNpEsXftCdy+rYqxhTbaOAISu/RFu7H5mn3VfdmPyh+8VcDLwKPAdem\n8j2AxaV2x1DcMfMAxVBWT/kuwC+APwDXAztXEGPdbdSJcXuKN//omv4XAXcBd1K8OXdvR4wUd4b8\nLj2WdeJ+pBhGibSv7kiPY6rej/XeX8CpwKlpWRRfEvZAiqG7r74V/a70F+MFwFOl/ba0v+Pehhhn\npxh+R3Fh/rBO249p/UTg0pp+LduPzXp4ahAzM8syXIekzMysQU4YZmaWxQnDzMyyOGGYmVkWJwwz\nM8vihGHWABUz3/7PmrLTJH2njz4bqo/MrHpOGGaNuYTiw1ll9eb6MtviOGGYNeYnwLHpU709Exru\nAfxW0i8k3Z6+32CzGUtVfHfDlaX1cyWdmJbfLOnGNAndtTWfQjfrCE4YZg2IYuK4Wynm+ILi7OLH\nwAvAcRFxCHAk8LWeKSv6k+a6Ogd4X0S8GZgPfLnZsZsN1pCcfNCszXqGpa5IP0+mmOrjHyW9A3iZ\nYorrsRTTrvdnP+CNFLOVQvGlPI/02cOsDZwwzBp3BfB1FV+rul1E3JaGlrqAN0fEHyWtBLat6beJ\nV57V99QLWBYRb6s2bLPB8ZCUWYMiYgPwS4qho56L3aOBtSlZHAnsXafrQ8BESa+StBPFNykC3A90\nSXobFENUpS8CMusYPsMwG5hLgIX8+Y6pHwI/k3QXsBS4r7ZDRKyS9GOKKawfBH6byjdKeh/wrTQN\n9tbANyhmMDXrGJ6t1szMsnhIyszMsjhhmJlZFicMMzPL4oRhZmZZnDDMzCyLE4aZmWVxwjAzsyz/\nH2bgaUMkwe4BAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1b1bcd9acf8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHU9JREFUeJzt3XuYHVWZ7/Hvj4RbIBAgLZeQ0KAMGpVLaC4qaMDLkIAi\nM4zAIDc5BhAUZjxnyKgP+MwcZvAZUQ+gYsDIRQwgd0yAAYbrKEKCARIgEkOQJECaawgwQsN7/qjV\nUNns3b2603tXJ/l9nmc/XVVrVdW7a3ftd9eqqlWKCMzMzHqzVtUBmJnZqsEJw8zMsjhhmJlZFicM\nMzPL4oRhZmZZnDDMzCyLE4bZICPpaEn3tHidB0l6StJySbu0YH17S5rX7PXYwHLCWE1JWijp9fQF\n8KKk6ZJGD8By2yWFpKGZ9UPSB1Z2vbkkjZf0dnrfr0iaJ+mYVq1/FfZ94KSI2DAi/lAukPSYpK/U\nziDpZEkz+7OyiLg7InboZ6xWESeM1dvnI2JDYEvgWeCciuN5j9zE00dL0vveCPgH4HxJ/nLq2TbA\n3AZlFwFH1pl+RCrrkyZ95tYCThhrgIj4H+BKYGz3NEkbS7pYUqekJyV9R9JaqWytNP6kpKWp3sZp\n1rvS35fSr/iPSfqApDslvSzpOUmXp+V0130w1T0kHQEsknSqpGeAX0jaRNJvUiwvpuGtS7HeIenf\nJd0naZmk6yRtmvG+IyJmAC8AO5aW90FJt0h6IR2BfKlUNlHSI+noZLGk/52m58T4fyX9Nr3XGyRt\nJunSFPP9ktpL9UPSNyQtSNvsP7q3f63+xFtnGXU/U0nrSloODEmf05/qzH4JsJekbUrLG5u26bQ0\nfoykR1McCyQdV6pb7zMfL2lRqc5kSX9K8z8i6aBS2dGS7pH0/bTtn5A0oVS+qaRfSFqSyq8tlR0g\nabakl9Jn887/gfVDRPi1Gr6AhcBn0vAwil+CF5fKLwauA4YD7cAfgWNT2VeA+cB2wIbA1cAlqawd\nCGBoaVnTgG9T/ABZD9irVBbAB0rj44Eu4HvAusD6wGbA36Y4hwO/Bq4tzXMHsBj4CLABcBXwywbv\nezywKA2vBXwBeBvYJU3bAHgKOAYYCuwCPAeMTeVPA3un4U2AcWk4J8b5wPuBjYFH0jb9TFrPxcAv\narbL7cCmwJhU93+lsqOBe1Ym3jrbpeFnWu9zqjP/LcB3SuP/XvP+90/vXcCngNdK267eZ/7O55Tq\n/B2wVfrMDgFeBbYsbY83ga9SJLYTgCWAUvl04PL0/tcGPpWm7wIsBfZI8x1FsV+sW/X+uaq+Kg/A\nryZ9sMWOsRx4Ke1sS4CPprIhwBvdXzpp2nHAHWn4NuBrpbId0jKGUj9hXAxMAbauE0e9hPEGsF4P\nse8MvFgavwM4szQ+Ni1jSJ15x1MkiJeAvwBvAaeUyg8B7q6Z52fA6Wn4z2lbbNTL9q0X47dL42cB\nN5bGPw/Mrtku+5XGvwbcloaP5t2EMVDxNvxM631Odeb/MjAvDa+V1ntQD/WvBU5u9JlTkzDqzD8b\nOLC0PeaXyoaleLegaG59G9ikzjJ+CvxrzbR5pITiV99fbpJavX0xIkZQ/Oo/CbhT0hbASIpfYk+W\n6j4JjErDW9UpGwps3mA9/0Txy/I+SXNV5wRpjc4omskAkDRM0s9Sc8kyimavEZKGlOZ5qiaetdP7\nqGdJet8bAWcD+5bKtgH2SE0UL0l6CTic4ssHiqOIicCTqZntY32I8dnS8Ot1xjesibP2PW1V5730\nK946+vqZ1roa2FLSnhRf9sMoftkDIGmCpHtTs9lLKaby57PCZ15L0pGlpqOXKI4my/M/0z0QEa+l\nwQ2B0cALEfFincVuA3yzZtuNpv52tgxOGGuAiHgrIq6m+LW9F0WTxpsUO1S3MRTNPlAcjdSWdVF8\nAb6ne+OIeCYivhoRW1H82v2Jer4yqnYZ36T4xbtHRGwEfDJNV6lO+QqvMSn+53pYBxHxF+BU4KOS\nvpgmPwXcGREjSq8NI+KENM/9EXEg8D6KX8lX9CHGvqp9T0vq1OlvvLV6+kx7lb6kr6Q4+X0EcFlE\nvAEgaV2KZsLvA5unZD2DFbdNw26x07mR8yl+1GyW5p9D3rZ9CthU0ogGZWfUbLthETEtY7lWhxPG\nGkCFAynaeB+NiLcovljOkDQ87bD/CPwyzTIN+AdJ20raEPg34PKI6AI6KZoAtist/+/07gngFym+\nHN5O48+W6zYwnOIX+EsqTmafXqfOlyWNlTQM+BfgyvQ+epS+1M4CTkuTfgP8laQjJK2dXrtJ+pCk\ndSQdLmnjiHgTWFZ6Hzkx9tX/UXEyfTRwMkU7fK3+xlurp88010UUTWR/y4pXR61DcW6iE+hKJ6Q/\n14flbkDxP9MJxQl0iiOMXkXE08CNFD9SNknbpzuZnw8cL2mPtA9sIGl/ScP7EJuVOGGs3m5QcQXM\nMuAM4KiI6L508usUJxYXAPcAvwKmprKpFFfG3AU8AfxPqt/9S/MM4L/TYf6ewG7A79O6rqdou16Q\nlvVd4KJU952re2r8iOJE6HPAvcBNdepcAlxI0TSxHvCNPmyHqcAYSZ+PiFcovswOpfjV/QzvnoyF\n4tfzwtTsdDxF809ujH11HTCLor1+OvDz2gorEW+thp9pH9wFvExx7uH+mhi/QfEj5EXg7yn+D7JE\nxCMUSf13FD8wPgr8dx/iOoLiiPMxipPcp6TlzqQ4UX5uims+xfkQ66fuqwzMBi1Jd1BcFXVB1bEM\nFEkBbB8R86uOxSyXjzDMzCyLE4aZmWVxk5SZmWXxEYaZmWVZrToBGzlyZLS3t1cdhpnZKmPWrFnP\nRURbTt3VKmG0t7czc2a/els2M1sjSXqy91oFN0mZmVkWJwwzM8vihGFmZlmcMMzMLIsThpmZZXHC\nMDOzLE1LGJJGS7pdxfN550o6OU3fVMXziR9PfzdpMP9+Kp5fPF/S5GbFaWZmeZp5hNEFfDMixgJ7\nAieqeHD8ZIpHUW5P8djI9ySD9BSzHwMTKB7HeVia18zMKtK0hBERT0fEA2n4FeBRikeAHsi7D1+5\nCPhindl3p3iG74L0AJzL0nxmZlaRltzpLakd2AX4PcUjHJ9ORc9Q/5nCo1jxeceLgD0aLHsSMAlg\nzJgxAxOw2QBrnzy990pNsPDM/StZr62emn7SOz0O8irglIhYVi6LoqvcleouNyKmRERHRHS0tWV1\nh2JmZv3Q1IQhaW2KZHFpRFydJj8ractUviXFIxVrLQZGl8a3TtPMzKwizbxKShTPKH40In5QKroe\nOCoNH0XxXONa9wPbpwfWr0PxPOPsZwSbmdnAa+YRxicoHs6+r6TZ6TUROBP4rKTHgc+kcSRtJWkG\nQER0AScBN1OcLL8iIuY2MVYzM+tF0056R8Q9gBoUf7pO/SXAxNL4DGBGc6IzM7O+8p3eZmaWxQnD\nzMyyOGGYmVkWJwwzM8uyWj3T26w3Vd1xbbY68BGGmZllccIwM7MsThhmZpbFCcPMzLI4YZiZWRYn\nDDMzy+KEYWZmWZwwzMwsixOGmZllccIwM7MsThhmZpalaX1JSZoKHAAsjYiPpGmXAzukKiOAlyJi\n5zrzLgReAd4CuiKio1lxmplZnmZ2PnghcC5wcfeEiDike1jSWcDLPcy/T0Q817TozMysT5r5iNa7\nJLXXK5Mk4EvAvs1av5mZDayqzmHsDTwbEY83KA/gVkmzJE1qYVxmZtZAVc/DOAyY1kP5XhGxWNL7\ngFskPRYRd9WrmBLKJIAxY8YMfKRmZgZUcIQhaSjwN8DljepExOL0dylwDbB7D3WnRERHRHS0tbUN\ndLhmZpZU0ST1GeCxiFhUr1DSBpKGdw8DnwPmtDA+MzOro2kJQ9I04HfADpIWSTo2FR1KTXOUpK0k\nzUijmwP3SHoQuA+YHhE3NStOMzPL08yrpA5rMP3oOtOWABPT8AJgp2bFZWZm/eM7vc3MLIsThpmZ\nZXHCMDOzLE4YZmaWxQnDzMyyOGGYmVkWJwwzM8vihGFmZlmcMMzMLIsThpmZZXHCMDOzLE4YZmaW\nxQnDzMyyOGGYmVkWJwwzM8vihGFmZlmcMMzMLEszH9E6VdJSSXNK074rabGk2ek1scG8+0maJ2m+\npMnNitHMzPI18wjjQmC/OtN/GBE7p9eM2kJJQ4AfAxOAscBhksY2MU4zM8vQtIQREXcBL/Rj1t2B\n+RGxICLeAC4DDhzQ4MzMrM+qOIfxdUkPpSarTeqUjwKeKo0vStPqkjRJ0kxJMzs7Owc6VjMzS1qd\nMH4KbAfsDDwNnLWyC4yIKRHREREdbW1tK7s4MzNroKUJIyKejYi3IuJt4HyK5qdai4HRpfGt0zQz\nM6tQSxOGpC1LowcBc+pUux/YXtK2ktYBDgWub0V8ZmbW2NBmLVjSNGA8MFLSIuB0YLyknYEAFgLH\npbpbARdExMSI6JJ0EnAzMASYGhFzmxWnmZnlaVrCiIjD6kz+eYO6S4CJpfEZwHsuuTUzs+r4Tm8z\nM8vihGFmZlmcMMzMLIsThpmZZXHCMDOzLE4YZmaWxQnDzMyyOGGYmVkWJwwzM8vihGFmZlmcMMzM\nLIsThpmZZXHCMDOzLE4YZmaWxQnDzMyyNO15GDb4tU+eXsl6F565fyXrNbOV4yMMMzPLkpUwJH20\nrwuWNFXSUklzStP+Q9Jjkh6SdI2kEQ3mXSjpYUmzJc3s67rNzGzg5R5h/ETSfZK+JmnjzHkuBPar\nmXYL8JGI2BH4I/DPPcy/T0TsHBEdmeszM7MmykoYEbE3cDgwGpgl6VeSPtvLPHcBL9RM+8+I6Eqj\n9wJb9z1kMzOrQvY5jIh4HPgOcCrwKeDs1Lz0N/1c91eAGxutDrhV0ixJk3paiKRJkmZKmtnZ2dnP\nUMzMrDe55zB2lPRD4FFgX+DzEfGhNPzDvq5U0reBLuDSBlX2ioidgQnAiZI+2WhZETElIjoioqOt\nra2voZiZWabcI4xzgAeAnSLixIh4ACAillAcdWSTdDRwAHB4RES9OhGxOP1dClwD7N6XdZiZ2cDL\nvQ9jf+D1iHgLQNJawHoR8VpEXJK7Mkn7Af8EfCoiXmtQZwNgrYh4JQ1/DviX3HWYmVlz5B5h3Aqs\nXxoflqY1JGka8DtgB0mLJB0LnAsMB25Jl8yel+puJWlGmnVz4B5JDwL3AdMj4qbsd2RmZk2Re4Sx\nXkQs7x6JiOWShvU0Q0QcVmfyzxvUXQJMTMMLgJ0y4zIzsxbJPcJ4VdK47hFJuwKvNyckMzMbjHKP\nME4Bfi1pCSBgC+CQpkVlZmaDTlbCiIj7JX0Q2CFNmhcRbzYvLDMzG2z60lvtbkB7mmecJCLi4qZE\nZWZmg05WwpB0CfB+YDbwVpocgBOGmdkaIvcIowMY2+hGOzMzW/3lXiU1h+JEt5mZraFyjzBGAo9I\nug/4S/fEiPhCU6IyM7NBJzdhfLeZQZiZ2eCXe1ntnZK2AbaPiFvTXd5DmhuamZkNJrlXSX0VmARs\nSnG11CjgPODTzQvNzFZl7ZOnV7LehWfuX8l61wS5J71PBD4BLIN3Hqb0vmYFZWZmg09uwvhLRLzR\nPSJpKMV9GGZmtobITRh3SvoWsH56lvevgRuaF5aZmQ02uQljMtAJPAwcB8ygj0/aMzOzVVvuVVJv\nA+enl5mZrYFyr5J6gjrnLCJiuwGPyMzMBqXcJqkOit5qdwP2Bs4GftnTDJKmSloqaU5p2qaSbpH0\nePq7SYN595M0T9J8SZMzYzQzsybKShgR8XzptTgifgT0drHzhcB+NdMmA7dFxPbAbWl8BZKGAD8G\nJgBjgcMkjc2J08zMmie3SWpcaXQtiiOOHueNiLsktddMPhAYn4YvAu4ATq2pszswPz3bG0mXpfke\nyYnVzMyaI7cvqbNKw13AQuBL/Vjf5hHxdBp+Bti8Tp1RwFOl8UXAHo0WKGkSxV3ojBkzph8hWatV\ndQfwmsjb2gZS7lVS+wz0iiMiJK30zX8RMQWYAtDR0eGbCc3MmiS3SeofeyqPiB9kru9ZSVtGxNOS\ntgSW1qmzGBhdGt86TTMzswr15SqpEyiai0YBxwPjgOHplet64Kg0fBRwXZ069wPbS9pW0jrAoWk+\nMzOrUO45jK2BcRHxCoCk7wLTI+LLjWaQNI3iBPdISYuA04EzgSskHQs8SToPImkr4IKImBgRXZJO\nAm6m6EJ9akTM7c+bMzOzgZObMDYH3iiNv0H9E9bviIjDGhS9p0v0iFgCTCyNz6DofsTMzAaJ3IRx\nMXCfpGvS+BcpLos1M7M1RO5VUmdIupHiLm+AYyLiD80Ly8zMBpvck94Aw4BlEfH/gEWStm1STGZm\nNghlJQxJp1Pckf3PadLa9NKXlJmZrV5yjzAOAr4AvArvnKTuy+W0Zma2istNGG9ERJC6OJe0QfNC\nMjOzwSg3YVwh6WfACElfBW7FD1MyM1uj5F4l9f30LO9lwA7AaRFxS1MjMzOzQaXXhJGeT3Fr6oDQ\nScLMbA3Va5NURLwFvC1p4xbEY2Zmg1Tund7LgYcl3UK6UgogIr7RlKjMzGzQyU0YV6eXmZmtoXpM\nGJLGRMSfI8L9RpmZreF6O4dxbfeApKuaHIuZmQ1ivSUMlYa3a2YgZmY2uPV2DiMaDNsAaZ88veoQ\nzMyy9JYwdpK0jOJIY/00TBqPiNioqdGZmdmg0WOTVEQMiYiNImJ4RAxNw93j/UoWknaQNLv0Wibp\nlJo64yW9XKpzWn/WZWZmAyf3stoBExHzgJ3hnbvIFwPX1Kl6d0Qc0MrYzMyssb48QKkZPg38KSKe\nrDgOMzPrRdUJ41BgWoOyj0t6SNKNkj7caAGSJkmaKWlmZ2dnc6I0M7PqEoakdSgeyvTrOsUPAGMi\nYkfgHEr3g9SKiCkR0RERHW1tbc0J1szMKj3CmAA8EBHP1hZExLKIWJ6GZwBrSxrZ6gDNzOxdVSaM\nw2jQHCVpC0lKw7tTxPl8C2MzM7MaLb9KCt55xOtngeNK044HiIjzgIOBEyR1Aa8Dh6ZHxJqZWUUq\nSRgR8SqwWc2080rD5wLntjouMzNrrOqrpMzMbBXhhGFmZlmcMMzMLIsThpmZZXHCMDOzLE4YZmaW\nxQnDzMyyOGGYmVkWJwwzM8vihGFmZlmcMMzMLIsThpmZZXHCMDOzLE4YZmaWxQnDzMyyOGGYmVkW\nJwwzM8tSScKQtFDSw5JmS5pZp1ySzpY0X9JDksZVEaeZmb2rkke0JvtExHMNyiYA26fXHsBP018z\nM6vIYG2SOhC4OAr3AiMkbVl1UGZma7KqEkYAt0qaJWlSnfJRwFOl8UVp2ntImiRppqSZnZ2dTQjV\nzMyguoSxV0TsTNH0dKKkT/Z3QRExJSI6IqKjra1t4CI0M7MVVJIwImJx+rsUuAbYvabKYmB0aXzr\nNM3MzCrS8oQhaQNJw7uHgc8Bc2qqXQ8cma6W2hN4OSKebnGoZmZWUsVVUpsD10jqXv+vIuImSccD\nRMR5wAxgIjAfeA04poI4zcyspOUJIyIWADvVmX5eaTiAE1sZl5mZ9WywXlZrZmaDjBOGmZllccIw\nM7MsThhmZpalyr6kzMwGXPvk6ZWte+GZ+1e27lbwEYaZmWVxwjAzsyxOGGZmlsUJw8zMsjhhmJlZ\nFicMMzPL4oRhZmZZnDDMzCyLE4aZmWVxwjAzsyzuGsTMbIBU1S1Jq7ok8RGGmZllqeKZ3qMl3S7p\nEUlzJZ1cp854SS9Lmp1ep7U6TjMzW1EVTVJdwDcj4gFJw4FZkm6JiEdq6t0dEQdUEJ+ZmdXR8iOM\niHg6Ih5Iw68AjwKjWh2HmZn1TaXnMCS1A7sAv69T/HFJD0m6UdKHe1jGJEkzJc3s7OxsUqRmZlZZ\nwpC0IXAVcEpELKspfgAYExE7AucA1zZaTkRMiYiOiOhoa2trXsBmZmu4ShKGpLUpksWlEXF1bXlE\nLIuI5Wl4BrC2pJEtDtPMzEqquEpKwM+BRyPiBw3qbJHqIWl3ijifb12UZmZWq4qrpD4BHAE8LGl2\nmvYtYAxARJwHHAycIKkLeB04NCKigljNzCxpecKIiHsA9VLnXODc1kRUqPLB8WZmqwLf6W1mZlmc\nMMzMLIsThpmZZXHCMDOzLE4YZmaWxQnDzMyyOGGYmVkWJwwzM8vihGFmZlmcMMzMLIsThpmZZXHC\nMDOzLE4YZmaWxQnDzMyyOGGYmVkWJwwzM8vihGFmZlkqSRiS9pM0T9J8SZPrlEvS2an8IUnjqojT\nzMze1fKEIWkI8GNgAjAWOEzS2JpqE4Dt02sS8NOWBmlmZu9RxRHG7sD8iFgQEW8AlwEH1tQ5ELg4\nCvcCIyRt2epAzczsXUMrWOco4KnS+CJgj4w6o4CnaxcmaRLFUQjAcknzBi7UhkYCz7VgPSvDMQ6c\nVSFOxzhwVoU4V4hR31upZW2TW7GKhDGgImIKMKWV65Q0MyI6WrnOvnKMA2dViNMxDpxVIc6qYqyi\nSWoxMLo0vnWa1tc6ZmbWQlUkjPuB7SVtK2kd4FDg+po61wNHpqul9gRejoj3NEeZmVnrtLxJKiK6\nJJ0E3AwMAaZGxFxJx6fy84AZwERgPvAacEyr4+xFS5vA+skxDpxVIU7HOHBWhTgriVERUcV6zcxs\nFeM7vc3MLIsThpmZZXHCaKC37ktK9XaT1CXp4FbGV1p/r3FKGi9ptqS5ku4cbDFK2ljSDZIeTDG2\n/JyVpKmSlkqa06C88u5qMmI8PMX2sKTfStqp1TGmOHqMs1Svsn0nJ8ZBsN/09nm3fr+JCL9qXhQn\n4/8EbAesAzwIjG1Q778oTtIfPBjjBEYAjwBj0vj7BmGM3wK+l4bbgBeAdVoc5yeBccCcBuUTgRsB\nAXsCv6/g8+4txo8Dm6ThCVXEmBNn6f+iyn2nt21Z6X6TGWPL9xsfYdSX030JwNeBq4ClrQyuJCfO\nvweujog/A0REq2PNiTGA4ZIEbEjxj9/VyiAj4q603kYq766mtxgj4rcR8WIavZfi/qWWy9iWUPG+\nkxFj1ftNTowt32+cMOpr1DXJOySNAg6i2o4Re40T+CtgE0l3SJol6ciWRVfIifFc4EPAEuBh4OSI\neLs14WXLeR+DybEUR0SDziDZd3pT9X6To+X7zSrfNUiFfgScGhFvFwl+0BoK7Ap8Glgf+J2keyPi\nj9WGtYK/BmYD+wLvB26RdHdELKs2rFWTpH0oEsZeVcfSwKqw73i/qcMJo76crkk6gMvSP/xIYKKk\nroi4tjUhAnlxLgKej4hXgVcl3QXsBLTqHz8nxmOAM6NojJ0v6Qngg8B9rQkxyyrRXY2kHYELgAkR\n8XzV8TQwGPad3lS93+Ro+X7jJqn6eu2+JCK2jYj2iGgHrgS+VsE/fE43K9cBe0kaKmkYRc/Ajw6y\nGP9M8UsOSZsDOwALWhhjjkHfXY2kMcDVwBGD7JfwCgbJvtObqvebHC3fb3yEUUfkdV9SuZw4I+JR\nSTcBDwFvAxdERI+XO7Y6RuBfgQslPUxxFdKpEdHS7qUlTQPGAyMlLQJOB9YuxVh5dzUZMZ4GbAb8\nJP1674oKejTNiLNyvcVY9X6TEyMV7DfuGsTMzLK4ScrMzLI4YZiZWRYnDDMzy+KEYWZmWZwwzMws\nixOGWR9Iul3SX9dMO0VSw24uJC1vfmRmzeeEYdY30yhuPiw7NE03W605YZj1zZXA/umudSS1A1sB\nf5B0m6QH0vMo3tO7cXq+wm9K4+dKOjoN7yrpztTR3c2t7gnXLIcThlkfRMQLFH31TEiTDgWuAF4H\nDoqIccA+wFnK7FlP0trAORTPhdgVmAqcMdCxm60sdw1i1nfdzVLXpb/HUnTN8G+SPknRlcQoYHPg\nmYzl7QB8hKK3USi6UBlU/VSZgROGWX9cB/wwPaZ1WETMSk1LbcCuEfGmpIXAejXzdbHiUX13uYC5\nEfGx5oZttnLcJGXWRxGxHLidoumo+2T3xsDSlCz2AbapM+uTwFhJ60oaQeppFJgHtEn6GBRNVJI+\n3NQ3YdYPPsIw659pwDW8e8XUpcANqefQmcBjtTNExFOSrgDmAE8Af0jT35B0MHC2pI0p9ssfAXOb\n/i7M+sC91ZqZWRY3SZmZWRYnDDMzy+KEYWZmWZwwzMwsixOGmZllccIwM7MsThhmZpbl/wOK68lq\nevhjSAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1b1bcd9aef0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 4. Visualize the summary statistic, calculate the variance or any other measure of spread to represent uncertainty in\n",
    "#    the summary statistic.\n",
    "plt.hist(mean,bins=20,normed=True)\n",
    "plt.title(\"Bootstrap Resamples of Mean\")\n",
    "plt.xlabel(\"Value\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.show()\n",
    "\n",
    "plt.hist(var)\n",
    "plt.title(\"Bootstrap Resamples of Variance\")\n",
    "plt.xlabel(\"Value\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>100.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>-0.080045</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.359072</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>-0.958700</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>-0.354547</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>-0.031626</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>0.131660</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>0.891775</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             mean\n",
       "count  100.000000\n",
       "mean    -0.080045\n",
       "std      0.359072\n",
       "min     -0.958700\n",
       "25%     -0.354547\n",
       "50%     -0.031626\n",
       "75%      0.131660\n",
       "max      0.891775"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 5. Summary statistics to describe uncertainty in the mean from bootstrap\n",
    "columns = ['mean']\n",
    "df_mean = pd.DataFrame(mean,columns=columns)\n",
    "df_mean.describe()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>var</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>100.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>1.068578</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.313462</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>0.346469</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>0.814319</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>1.075160</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>1.259427</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>1.819811</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              var\n",
       "count  100.000000\n",
       "mean     1.068578\n",
       "std      0.313462\n",
       "min      0.346469\n",
       "25%      0.814319\n",
       "50%      1.075160\n",
       "75%      1.259427\n",
       "max      1.819811"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 6. Summary statistics to describe uncertainty in the variance from bootstrap\n",
    "columns = ['var']\n",
    "df_var = pd.DataFrame(var,columns=columns)\n",
    "df_var.describe()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# What did we learn?\n",
    "# 1. Uncertainty decreases rapidly as the number of samples increases.\n",
    "# 2. Bootstrap result for uncertainty in the mean is the same as standard error.\n",
    "# 3. Bootstrap provides uncertainty in any statistic.\n",
    "# 4. Bootstrap does not account for spatial context, location of data, local nonstationarity, spatial correlation etc."
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}