Spatial_Bootstrap.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GeostatsPy: Basic Univariate Statistics and Distribution Plotting for Subsurface Data Analytics in Python \n",
    "\n",
    "\n",
    "### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "\n",
    "#### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Demonstration of Spatial 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.  Specify locations.  Sort the data in a data only array for \n",
    "   an equal weighted cumulative density function (CDF) for backtransform\n",
    "2. Calculate the covariance matrix (between the data) for LU simulation.  This provides a fast and easy method\n",
    "   to calculate spaitally correlated, unconditional, Gaussian realizations at the data locations.\n",
    "3. Perform LU (lower upper) decomposition of the covariance matrix.  Realizations only require L x random vector.\n",
    "4. 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.  Use LU simulation to ensure spaital correlation \n",
    "   between samples.  Backtransform from Gaussian to original data values (use percentiles to sample).\n",
    "5. Calculate the statistic of interest for each realization. This demonstration considers the mean only.  \n",
    "   We could have considered any statistic including median, 13th percentile, skew etc. \n",
    "6. (and 7.) Quantify and visualize uncertainty with histograms and summary statistics.\n",
    "\n",
    "References:\n",
    "\n",
    "Efron, 1982, The jackknife, the bootstrap, and other resampling plans, Society of Industrial and Applied Math, \n",
    "CBMS-NSF Monographs, 38.\n",
    "\n",
    "Journel, A. G. Resampling from stochastic simulations. Environmental and Ecological Statistics, 1:63–84, 1994.\n",
    "\n",
    "#### Import Packages\n",
    "\n",
    "We will also need some standard packages. These should have been installed with Anaconda 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports \n",
    "import pprint\n",
    "import math\n",
    "import scipy\n",
    "import scipy.linalg                # SciPy Linear Algebra Library\n",
    "from scipy.stats import norm\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set the working directory\n",
    "\n",
    "I always like to do this so I don't lose files and to simplify subsequent read and writes (avoid including the full address each time). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'os' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-2-8e1ada8ab58b>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mos\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mchdir\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"c:/PGE383\"\u001b[0m\u001b[1;33m)\u001b[0m             \u001b[1;31m# set the working directory\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'os' is not defined"
     ]
    }
   ],
   "source": [
    "os.chdir(\"c:/PGE383\")             # set the working directory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Assume Data Values and Locations\n",
    "\n",
    "For demonstration I assume 4 data values in 2D and make a CDF by sorting for the back transform from Gaussian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[25.  50.   2.3]\n",
      " [75.  80.   5.5]\n",
      " [10.  25.   7. ]\n",
      " [95.  15.   8.7]]\n",
      "[2.3 5.5 7.  8.7]\n"
     ]
    }
   ],
   "source": [
    "# 1. Make Spatial Dataset and Define Number of Realizations (just need a sufficient number to summarize uncertainty)\n",
    "nreal = 100                                               # User Specified\n",
    "# Assign Data Values - assume 4 data locations\n",
    "ndata = 4                                                 # User Specified\n",
    "output = np.zeros((ndata,nreal))  # x, y, value \n",
    "data = np.zeros((ndata,3))  # x, y, value \n",
    "data[0,0] = 25.0;  data[0,1] = 50.0;  data[0,2] = 2.3     # User Specified\n",
    "data[1,0] = 75.0;  data[1,1] = 80.0;  data[1,2] = 5.5     # User Specified\n",
    "data[2,0] = 10.0;  data[2,1] = 25.0;  data[2,2] = 7.0     # User Specified\n",
    "data[3,0] = 95.0;  data[3,1] = 15.0;  data[3,2] = 8.7     # User Specified\n",
    "# Set up backtransform lookup table\n",
    "value = np.zeros(ndata) \n",
    "for i in range(0, ndata):\n",
    "    value[i] = data[i,2]\n",
    "print(data) \n",
    "value.sort()\n",
    "print(value)\n",
    "\n",
    "# data in ascending "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate the Covariance Matrix\n",
    "\n",
    "The covariance matrix is calculated from the distance matrix assuming a isotropic spherical variogram with range set as $var_range$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.2244834 , 0.57506938, 0.06574285],\n",
       "       [0.2244834 , 1.        , 0.03145365, 0.13715671],\n",
       "       [0.57506938, 0.03145365, 1.        , 0.0296663 ],\n",
       "       [0.06574285, 0.13715671, 0.0296663 , 1.        ]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 2. Calculate Symmetric Covariance Array - assuming variogram with spherical structure with range specified, var_range\n",
    "cov = np.zeros((ndata,ndata))\n",
    "var_range = 100.0                                          # User Specified\n",
    "for i in range(0, ndata):\n",
    "    for j in range(0, ndata):\n",
    "        distance = math.sqrt(math.pow((data[i,0]-data[j,0]),2) + math.pow((data[i,1]-data[j,1]),2))\n",
    "        cova = 0.0\n",
    "        if(distance < var_range):\n",
    "            hr = distance / var_range\n",
    "            cova = 1.0 - hr*(1.5 - 0.5*hr*hr)  # spherical structure, no nugget\n",
    "            cov[i,j] = cova\n",
    "cov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### LU Decomposition\n",
    "\n",
    "Peform LU decomposition to calculate the needed lower matrix $L$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.          0.          0.          0.        ]\n",
      " [ 0.2244834   1.          0.          0.        ]\n",
      " [ 0.57506938 -0.10282133  1.          0.        ]\n",
      " [ 0.06574285  0.12889386  0.00674212  1.        ]]\n",
      "[[ 1.          0.2244834   0.57506938  0.06574285]\n",
      " [ 0.          0.9496072  -0.09763988  0.12239854]\n",
      " [ 0.          0.          0.65925575  0.00444478]\n",
      " [ 0.          0.          0.          0.97987149]]\n"
     ]
    }
   ],
   "source": [
    "# 3. Perform LU Decomposition using scipy (used tutorial at www.quantstart.com) \n",
    "P, L, U = scipy.linalg.lu(cov)\n",
    "print(L);  print(U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test the LU Decomposition (available for testing only)\n",
    "#Test = cov - np.matmul(L,U)\n",
    "#Test  # should be zeros"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate the Realizations\n",
    "\n",
    "The following steps:\n",
    "\n",
    "1. Matrix multiplication of the nxn lower matrix with a 1D matrix of random standard normal values to calculate a realization with spatial correlation at data locations $i = 0,\\ldots,n-1$.\n",
    "\n",
    "2. Convert to p-values with the standard normal CDF\n",
    "\n",
    "3. Transform to discrete values in the data CDF.\n",
    "\n",
    "Realizations stored in output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 4. Calculate Unconditional Realization at the Specified Locations from the Data Distribution with Replacement\n",
    "#    Realizations in Gaussian space are converted to p-values and then applied to draw from the original distribution\n",
    "rand = np.zeros(ndata); pvalue = np.zeros(ndata)\n",
    "for ireal in range(0, nreal):    \n",
    "    for i in range(0, ndata):\n",
    "        rand[i] = np.random.normal()\n",
    "    realization = np.matmul(L,rand)\n",
    "# Convert to p-values    \n",
    "    for i in range(0, ndata):\n",
    "        pvalue[i] = norm.cdf(realization[i])\n",
    "# Complete backtransform by selecting percentiles from original cdf (drawing with replacement from provided data values)       \n",
    "        index = int(round( ndata * pvalue[i] - 0.5))\n",
    "        output[i,ireal] = value[index]\n",
    "idata = 0; ireal = 0\n",
    "output[idata,ireal]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate the Target Statistic for Each Realization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5.875, 4.7  , 3.1  , 5.875, 6.675, 6.25 , 6.675, 6.25 , 6.25 ,\n",
       "       5.875, 7.1  , 6.675, 5.875, 7.9  , 5.075, 6.675, 4.65 , 5.875,\n",
       "       7.85 , 5.45 , 8.275, 7.05 , 5.875, 3.9  , 5.075, 3.1  , 8.7  ,\n",
       "       7.85 , 6.675, 6.25 , 8.275, 5.45 , 8.7  , 7.05 , 8.275, 7.475,\n",
       "       5.5  , 6.675, 7.05 , 4.7  , 7.475, 5.075, 5.075, 5.875, 6.3  ,\n",
       "       3.9  , 5.875, 7.05 , 5.075, 8.275, 7.475, 4.65 , 6.3  , 5.45 ,\n",
       "       8.275, 5.45 , 7.475, 8.275, 7.1  , 7.85 , 5.875, 7.05 , 5.075,\n",
       "       7.05 , 6.625, 6.3  , 6.3  , 3.1  , 4.275, 6.3  , 3.9  , 6.3  ,\n",
       "       5.875, 5.875, 8.275, 5.45 , 7.475, 6.25 , 5.875, 6.25 , 4.7  ,\n",
       "       5.825, 4.275, 7.475, 8.275, 3.9  , 7.475, 5.5  , 4.275, 4.7  ,\n",
       "       5.875, 6.25 , 5.075, 3.1  , 3.475, 5.075, 8.275, 6.3  , 3.475,\n",
       "       4.275])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 5. Calculate the summary statistic of interest (mean only for this example)\n",
    "# Any summary statistic may be calculated.\n",
    "summary = np.zeros(nreal) \n",
    "for ireal in range(0, nreal): \n",
    "    mean = 0\n",
    "    for i in range(0, ndata):\n",
    "        mean = mean + output[i,ireal]    \n",
    "    mean = mean / ndata\n",
    "    summary[ireal] = mean\n",
    "summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\pm27995\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\matplotlib\\axes\\_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
      "  warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGmJJREFUeJzt3Xu8ZXVd//HXmwHkIoIKqdwcNKQmw8ARK28k2ANChyz9ieUFy7B+Epn2MLRf/JDSB5mlVlSi+fOWIFDamOMP89GPSgucESm5iE2IMozKQCIXERj8/P5Y37PcHM5ln5mzZs8cXs/HYz/OXmt991qftffZ673XPVWFJEkAO026AEnS9sNQkCT1DAVJUs9QkCT1DAVJUs9QkCT1DIUHqSRvTPKeMdu+L8nvD12Ttr0k1yc5dhtOb/ckH0/y7SQXbqvpanyGwnYkydOT/Gv7wvx3ks8mecoijPfoJBtG+1XVW6rqlYsw7pOT3Jfkjva4Lsmvbe1427jHDqMkZyb50GJMd1xJLkny3TbfNyf52ySP2ZY17IBeADwKeGRVvXD6wPY5VpLTpvV/Tet/5jaq80HLUNhOJHkY8PfAnwKPAA4A3gTcPcm6xvRvVfXQqnoo3Zf+rUmOmHRRo9IZ4v/91DbfPwg8FHjbANNYSh4LfLmqNs/R5svAy6f1e1nrr4EZCtuPJwBU1XlVdV9V3VVVn6qq/4D+F/lnk/xpW5P4UpJjpl6c5BVJrklye/u1/qrWf0/gk8D+I7/m95/+yzrJhUm+0cb9z0l+ZEtmoqouB64Bfnhk3KuSXJXk1vbrenTYD7d+t7Y2q1r/U4BfBF7fav546//bSW5s83ltkmOSHAe8EXhRa/vvre0lSd6c5LPAd4DHzfY+tfZHJ9nQNq3d3Dat/OKY830r8DHgx0bGt1OS05P8V5JbklyQ5BFt2G5JPtT635pkbZJHtWHj1Pj6JDcl+XqSn03yM0m+3NYw3zjS/swkFyX5SBvf5UmeNNM8bGm9M4xnts/0TcAZI5/TL8/ydq4F9pj6H2x/d2/9R6fz3CRXtOn8a5LDR4ZNzcftSa5O8vyRYScn+UyStyX5VpKvJDl+lloefKrKx3bwAB4G3AK8HzgeePi04ScDm4HfBHYBXgR8G3hEG34C8HggwLPoFoJHtmFHAxumje9M4EMj3b8E7AU8BHgHcMXIsPcBvz9L3ScDnxnpfgpwK/CE1v0E4E7gOa3u1wPrgV1b93q6BfquwLOB24HDZpoucBhwA7B/614OPH6m+Wn9LgG+BvwIsHOb3nzv02bgj9v78KxW+2GzzPslwCvb80cCnwb+bmT4a4BLgQPb+N4FnNeGvQr4OLAHsAx4MvCwMT/LzXQL112AXwE2AR9un9+PAN8FHjfyvtxLtwa3C/BbwFeAXdrw64Fjt6beae/JfJ/pAz6nmf4v2+v/oPV7K/CG1v/M1u9I4Cbgqa2el7d5eUgb/kJgf7ofvi9qn+NjRv5n723v3TLg14CNQCa9HNgeHhMvwMfIh9H9un4fsKF98VcDj2rDTp7+jwt8DnjpLOP6GPAb7fnRzBMK04btAxSwd+t+H3OHwma6ILijve5Pp+oEfhe4YKT9TsCNraZnAN8AdhoZft7IF/9+06XbRHMTcCxtoTbX/NAttM+a5z2f/j5tBvYcGX4B8LuzvPYSugX2t9t8XwEcPDL8GuCYke7HtIXRznQh/K/A4WP8X0yv8S5gWeveq037qSPtPw/87Mj7cum09//rwDNa9/V8PxS2ut4xPtNZ/+9GhwMH0wX6Lu3vQdw/FP4C+L1pr70WeNYs470COHHkf3b9yLA92nv46KG/4zvCw81H25GquqaqTq6qA4En0v3SecdIkxur/Rc3X21tSHJ8kkvb5oNbgZ8B9h1nukmWJTm7rW7fRregYNzX0y109qlu2/qj6X6tvqUN27/VOTWP36P7tX9AG3ZD6zc6TwfMNJGqWk/3a/ZM4KYk5yfZf57abhjtGON9+lZV3TmtnrmmcVpV7Q0cDjyc7lf2lMcCH22bN26lW+jeR7ej9YPAxcD5STYmeWuSXcas8Zaquq89v6v9/ebI8Lvo9m884D1o7/WGWeZpi+qdZkGf6Wyq6mt0axxvAf6zqm6Y1uSxwOumam31HsT3vw8vG9m0dCvd92n0PfzGyLS+056OvmcPWobCdqqqvkT3S/mJI70PSJKR7oOBjUkeAvwN3U7OR1XVPsAaus0P0P0KmssvACfS/QLfm26zDCOvX0jd32y1PK/12kj3Be5G2NV/EN3awkbgoNx/B/DBbdiMdVfVh6vq6W2cBfzBbG2n9x/jfQJ4eLr9MKP1bJxtfkfq+iLw+8A5I5/RDcDxLTCnHrtV1Y1VdW9VvamqVgA/CTwXeNmYNS7UQVNP2nt94CzztOB6ZxjHfJ/pQnwAeF37O1Otb55W6x5VdV6SxwLvBk6lO8ppH+BKtu49fNAwFLYTSX4oyeuSHNi6DwJeTLeNd8oPAKcl2SXJC+k2N62h23b7ELpty5vbTrOfHnndN4FHJtl7lsnvRXeU0y10q9JvmaXdOPPxSOD5wFWt1wXACel2CO9C9yW/m25TxGV023pf3+bpaLowOX+k7seNjPuwJM9uC87v0v0ivm+k7fLMfYTRfO/TlDcl2TXJM+gWfuMeT/9+us9oVev+S+DNbSFFkv2SnNie/1SSH02yDLiNbjPNfQuocSGenOTnkuxMt6Z1N/f/v5qyJfVON99nuhAfoZv3C2YY9m7gV5M8NZ09k5yQZC9gT7ofA5ta7a/g/j+uNAdDYftxO91Os8uS3En3pb2SbiE65TLgUOBm4M3AC6rqlqq6HTiN7svzLbpf/qunXtTWOs4Drmur09M3HXyAbhX/RuBqZl5gzOUn2tEkd9BtctgE/Hqb9rXAS+j2M9xMt4B4XlXdU1X30C1Aj2/D/hx4WasX4K+AFa3mj9EtLM9ubb9BtwCeOtJmasF9S5LLZypyvvep+UYbthH4a+BXR+qZU5ufP6HbjwLwzjb+TyW5ne59fWob9mjgIroF7DXAP9Ftax+nxoX6O7qdrd8CXgr8XFXdO0O7Bdc7fQRjfKZjq+4IvE9X1V0zDFtHt6P4z9p8rafbV0BVXQ38EfBvdD8WfhT47EKn/2A1tTNQ27kkJ9Md6fL0SdeyVLVftR9q+3SWhHQne/1gVb1k0rVox+CagiSpZyhIknpuPpIk9VxTkCT1dp50AQu177771vLlyyddhiTtUD7/+c/fXFX7zdduhwuF5cuXs27dukmXIUk7lCRfnb+Vm48kSSMMBUlSz1CQJPUMBUlSz1CQJPUMBUlSb9BQSHJcuvvork9y+gzDT06yqd0M44okrxyyHknS3AY7T6Fdd/0cunvzbgDWJlndLms76iNVdepQdUiSxjfkmsJRdPdBva5dY/18urt7SZK2U0Oe0XwA978/7ga+f8OOUT+f5JnAl4HfnOFerCQ5BTgF4OCDDx6gVKmz/PRPbPFrrz/7hB122tKUIdcUZrof6vRLsn4cWF5VhwOfprud4QNfVHVuVa2sqpX77TfvpTskSVtoyFDYwMgNw5nhZuHtVpJ3t853A08esB5J0jyGDIW1wKFJDkmyK3AS0+41m+QxI52r6O79KkmakMH2KVTV5iSnAhcDy4D3VtVVSc4C1lXVauC0JKuAzcB/0268LUmajEEvnV1Va4A10/qdMfL8DcAbhqxBkjQ+z2iWJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSb9BQSHJckmuTrE9y+hztXpCkkqwcsh5J0twGC4Uky4BzgOOBFcCLk6yYod1ewGnAZUPVIkkaz5BrCkcB66vquqq6BzgfOHGGdr8HvBX47oC1SJLGMGQoHADcMNK9ofXrJTkCOKiq/n6uESU5Jcm6JOs2bdq0+JVKkoBhQyEz9Kt+YLIT8HbgdfONqKrOraqVVbVyv/32W8QSJUmjhgyFDcBBI90HAhtHuvcCnghckuR64MeB1e5slqTJGTIU1gKHJjkkya7AScDqqYFV9e2q2reqllfVcuBSYFVVrRuwJknSHAYLharaDJwKXAxcA1xQVVclOSvJqqGmK0nacjsPOfKqWgOsmdbvjFnaHj1kLZKk+XlGsySpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpZyhIknqGgiSpt/M4jZI8saquHLoYSdreLT/9E1v82uvPPmERKxnGuGsKf5nkc0n+Z5J9Bq1IkjQxY4VCVT0d+EXgIGBdkg8nec6glUmStrmx9ylU1X8C/wv4beBZwJ8k+VKSnxuqOEnStjVWKCQ5PMnbgWuAZwPPq6ofbs/fPmB9kqRtaNw1hT8DLgeeVFWvrqrLAapqI93aw4ySHJfk2iTrk5w+w/BfTfLFJFck+UySFVsyE5KkxTHW0UfAzwB3VdV9AEl2Anarqu9U1QdnekGSZcA5wHOADcDaJKur6uqRZh+uqr9s7VcBfwwct2WzIknaWuOuKXwa2H2ke4/Wby5HAeur6rqqugc4HzhxtEFV3TbSuSdQY9YjSRrAuGsKu1XVHVMdVXVHkj3mec0BwA0j3RuAp05vlOTVwGuBXen2UUiSJmTcNYU7kxw51ZHkycBd87wmM/R7wJpAVZ1TVY+nO6ppxv0TSU5Jsi7Juk2bNo1ZsiRpocZdU3gNcGGSja37McCL5nnNBrrzGqYcCGycpS10m5f+YqYBVXUucC7AypUr3cQkSQMZKxSqam2SHwIOo1sD+FJV3TvPy9YChyY5BLgROAn4hdEGSQ5t5z8AnAD8J5KkiRl3TQHgKcDy9pojklBVH5itcVVtTnIqcDGwDHhvVV2V5CxgXVWtBk5NcixwL/At4OVbOB+SpEUw7gXxPgg8HrgCuK/1LmDWUACoqjXAmmn9zhh5/hsLKVaSNKxx1xRWAiuqyu352iaW+pUol4qt+ZzAz2p7NO7RR1cCjx6yEEnS5I27prAvcHWSzwF3T/WsqlWDVCVJmohxQ+HMIYuQJG0fxj0k9Z+SPBY4tKo+3c5mXjZsaZKkbW3cS2f/CnAR8K7W6wDgY0MVJUmajHF3NL8aeBpwG/Q33PmBoYqSJE3GuKFwd7vSKQBJdsYrmkrSkjNuKPxTkjcCu7d7M18IfHy4siRJkzBuKJwObAK+CLyK7izlWe+4JknaMY179NH3gHe3hyRpiRr32kdfYeZ7ITxu0SuSJE3MQq59NGU34IXAIxa/HEnSJI21T6Gqbhl53FhV78BbZ0rSkjPu5qMjRzp3oltz2GuQiiRJEzPu5qM/Gnm+Gbge+B+LXo20CLb2cs47oh31UuM7at1L2bhHH/3U0IVIkiZv3M1Hr51reFX98eKUI0mapIUcffQUYHXrfh7wz8ANQxQlSZqMhdxk58iquh0gyZnAhVX1yqEKkyRte+Ne5uJg4J6R7nuA5YtejSRposZdU/gg8LkkH6U7s/n5wAcGq0qSNBHjHn305iSfBJ7Rer2iqr4wXFmSpEkYd/MRwB7AbVX1TmBDkkMGqkmSNCHj3o7zfwO/Dbyh9doF+NBQRUmSJmPcNYXnA6uAOwGqaiNe5kKSlpxxQ+Geqira5bOT7DlcSZKkSRk3FC5I8i5gnyS/Anwab7gjSUvOuEcfva3dm/k24DDgjKr6h0ErkyRtc/OGQpJlwMVVdSxgEEjSEjZvKFTVfUm+k2Tvqvr2tihK24etvQS1lzaWdjzjntH8XeCLSf6BdgQSQFWdNkhVkqSJGDcUPtEekqQlbM5QSHJwVX2tqt6/JSNPchzwTmAZ8J6qOnva8NcCr6S7m9sm4Jeq6qtbMi1J0tab75DUj009SfI3Cxlx20F9DnA8sAJ4cZIV05p9AVhZVYcDFwFvXcg0JEmLa75QyMjzxy1w3EcB66vquqq6BzgfOHG0QVX9v6r6Tuu8FDhwgdOQJC2i+UKhZnk+jgO4/53ZNrR+s/ll4JMzDUhySpJ1SdZt2rRpgWVIksY1347mJyW5jW6NYff2nNZdVfWwOV6bGfrNGCxJXkJ3y89nzTS8qs4FzgVYuXLlQsNJkjSmOUOhqpZtxbg3AAeNdB8IbJzeKMmxwO8Az6qqu7diepKkrbSQ+yks1Frg0CSHJNkVOAlYPdogyRHAu4BVVXXTgLVIksYwWChU1WbgVOBi4Brggqq6KslZSVa1Zn8IPBS4MMkVSVbPMjpJ0jYw7slrW6Sq1gBrpvU7Y+T5sUNOX5K0MENuPpIk7WAMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSz1CQJPUMBUlSb9DbcW5vlp/+ia16/fVnnzCRaU9qupIefFxTkCT1DAVJUs9QkCT1DAVJUs9QkCT1DAVJUs9QkCT1DAVJUs9QkCT1DAVJUs9QkCT1DAVJUs9QkCT1Br1KapLjgHcCy4D3VNXZ04Y/E3gHcDhwUlVdNGQ9kpaOSV15eKkbbE0hyTLgHOB4YAXw4iQrpjX7GnAy8OGh6pAkjW/INYWjgPVVdR1AkvOBE4GrpxpU1fVt2PcGrEOSNKYh9ykcANww0r2h9VuwJKckWZdk3aZNmxalOEnSAw0ZCpmhX23JiKrq3KpaWVUr99tvv60sS5I0myFDYQNw0Ej3gcDGAacnSdpKQ4bCWuDQJIck2RU4CVg94PQkSVtpsFCoqs3AqcDFwDXABVV1VZKzkqwCSPKUJBuAFwLvSnLVUPVIkuY36HkKVbUGWDOt3xkjz9fSbVaSJG0HPKNZktQzFCRJPUNBktQzFCRJPUNBktQzFCRJvUEPSV1qtuZSvZK0tcuQbXHJb9cUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1DMUJEk9Q0GS1PPS2ZIedLwM/uxcU5Ak9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QwFSVLPUJAk9QYNhSTHJbk2yfokp88w/CFJPtKGX5Zk+ZD1SJLmNlgoJFkGnAMcD6wAXpxkxbRmvwx8q6p+EHg78AdD1SNJmt+QawpHAeur6rqqugc4HzhxWpsTgfe35xcBxyTJgDVJkuaQqhpmxMkLgOOq6pWt+6XAU6vq1JE2V7Y2G1r3f7U2N08b1ynAKa3zMODaMUrYF7h53lY7nqU6X7B058352vEsxXl7bFXtN1+jIe+nMNMv/ukJNE4bqupc4NwFTTxZV1UrF/KaHcFSnS9YuvPmfO14lvK8zWfIzUcbgINGug8ENs7WJsnOwN7Afw9YkyRpDkOGwlrg0CSHJNkVOAlYPa3NauDl7fkLgH+sobZnSZLmNdjmo6ranORU4GJgGfDeqroqyVnAuqpaDfwV8MEk6+nWEE5axBIWtLlpB7JU5wuW7rw5XzuepTxvcxpsR7MkacfjGc2SpJ6hIEnqLalQSLJbks8l+fckVyV506RrWmxJliX5QpK/n3QtiyXJ9Um+mOSKJOsmXc9iSbJPkouSfCnJNUl+YtI1LYYkh7XPaupxW5LXTLquxZDkN9uy48ok5yXZbdI1bWtLap9COxt6z6q6I8kuwGeA36iqSydc2qJJ8lpgJfCwqnrupOtZDEmuB1ZOP2lxR5fk/cC/VNV72hF4e1TVrZOuazG1y9ncSHfS6VcnXc/WSHIA3TJjRVXdleQCYE1VvW+ylW1bS2pNoTp3tM5d2mPJpF6SA4ETgPdMuhbNLcnDgGfSHWFHVd2z1AKhOQb4rx09EEbsDOzezpvagweeW7XkLalQgH7zyhXATcA/VNVlk65pEb0DeD3wvUkXssgK+FSSz7dLmiwFjwM2Af+nbe57T5I9J13UAE4Czpt0EYuhqm4E3gZ8Dfg68O2q+tRkq9r2llwoVNV9VfVjdGdQH5XkiZOuaTEkeS5wU1V9ftK1DOBpVXUk3RV1X53kmZMuaBHsDBwJ/EVVHQHcCTzg8vE7srZJbBVw4aRrWQxJHk53kc5DgP2BPZO8ZLJVbXtLLhSmtFX1S4DjJlzKYnkasKptfz8feHaSD022pMVRVRvb35uAj9JdYXdHtwHYMLKmehFdSCwlxwOXV9U3J13IIjkW+EpVbaqqe4G/BX5ywjVtc0sqFJLsl2Sf9nx3ug/5S5OtanFU1Ruq6sCqWk63yv6PVbXD/4pJsmeSvaaeAz8NXDnZqrZeVX0DuCHJYa3XMcDVEyxpCC9miWw6ar4G/HiSPdpBK8cA10y4pm1uyKukTsJjgPe3IyJ2Ai6oqiVz6OYS9Sjgo+02GjsDH66q/zvZkhbNrwN/3TazXAe8YsL1LJokewDPAV416VoWS1VdluQi4HJgM/AFHoSXu1hSh6RKkrbOktp8JEnaOoaCJKlnKEiSeoaCJKlnKEiSeoaCNIskleSDI907J9m0lK5QK01nKEizuxN4YjsRErrj8m+cYD3S4AwFaW6fpLsyLUw7g7edjf3eJGvbRe9ObP2XJ/mXJJe3x0+2/kcnuWTkHgt/3c6clbYbhoI0t/OBk9rNVg4HRq+6+zt0lxt5CvBTwB+2S3XcBDynXeTvRcCfjLzmCOA1wAq6K6k+bfhZkMa31C5zIS2qqvqPJMvp1hLWTBv803QXKfyt1r0bcDDdNfj/LMmPAfcBTxh5zeeqagNAu8T7crobu0jbBUNBmt9quuvsHw08cqR/gJ+vqmtHGyc5E/gm8CS6tfHvjgy+e+T5ffgd1HbGzUfS/N4LnFVVX5zW/2Lg16f2CyQ5ovXfG/h6VX0PeCmwbJtVKm0lQ0GaR1VtqKp3zjDo9+hu+fofSa5s3QB/Drw8yaV0m47u3DaVSlvPq6RKknquKUiSeoaCJKlnKEiSeoaCJKlnKEiSeoaCJKlnKEiSev8fpPF0Sk76sWEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 6. Visualize the summary statistic, calculate the variance or any other measure of spread to represent uncertainty in\n",
    "#    the summary statistic.\n",
    "plt.hist(summary,bins=20,normed=True)\n",
    "plt.title(\"Spatial Bootstrap Resamples of Mean\")\n",
    "plt.xlabel(\"Mean\")\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\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.00000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>6.07875</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>1.41530</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>3.10000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>5.07500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>6.25000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>7.05000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>8.70000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            mean\n",
       "count  100.00000\n",
       "mean     6.07875\n",
       "std      1.41530\n",
       "min      3.10000\n",
       "25%      5.07500\n",
       "50%      6.25000\n",
       "75%      7.05000\n",
       "max      8.70000"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 7. Summary statistics to describe uncertainty in the mean from bootstrap\n",
    "columns = ['mean']\n",
    "df_mean = pd.DataFrame(summary,columns=columns)\n",
    "df_mean.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summary of Exercise\n",
    "\n",
    "What did we learn?\n",
    "1. Uncertainty decreases rapidly as the number of data samples increases.\n",
    "2. Bootstrap result for uncertainty in the mean is the same as standard error if uncorrelated and more if spatially correlated.\n",
    "3. Spatial bootstrap provides uncertainty in any statistic.\n",
    "4. Spatial bootstrap does account for spatial context, location of data and spatial correlation.  Resampling from realizations\n",
    "   taught by Journel is available to account for spatial trends, conditioning data, and secondary information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This was a basic demonstration of spatial bootstrap. Much more could be done, I have other demonstrations on basics of working with DataFrames, ndarrays and many other workflows availble at https://github.com/GeostatsGuy/PythonNumericalDemos and https://github.com/GeostatsGuy/GeostatsPy.\n",
    "\n",
    "I hope this was helpful,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "Michael Pyrcz, Ph.D., P.Eng. Associate Professor The Hildebrand Department of Petroleum and Geosystems Engineering, Bureau of Economic Geology, The Jackson School of Geosciences, The University of Texas at Austin\n",
    "\n",
    "#### More Resources Available at: [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig)  | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)"
   ]
  },
  {
   "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
}