Interactive_Monte_Carlo_simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<p align=\"center\">\n",
    "    <img src=\"https://github.com/GeostatsGuy/GeostatsPy/blob/master/TCG_color_logo.png?raw=true\" width=\"220\" height=\"240\" />\n",
    "\n",
    "</p>\n",
    "\n",
    "## Data Analytics \n",
    "\n",
    "### Monte Carlo Simulation 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)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Monte Carlo Simulation in Python \n",
    "\n",
    "Here's a simple workflow, demonstration of Monte Carlo simulation for subsurface uncertainty modeling workflows. This should help you get started with building subsurface models that integrate uncertainty sources.  \n",
    "\n",
    "#### Monte Carlo Simulation\n",
    "\n",
    "Definition: random sampling from a distribution\n",
    "\n",
    "Procedure: \n",
    "\n",
    "1. Model the representative distribution (CDF)\n",
    "2. Draw a random value from a uniform [0,1] distribution (p-value)\n",
    "3. Apply the inverse of the CDF to calculate the associated realization\n",
    "\n",
    "In practice, Monte Carlo simulation refers to the workflow with multiple realizations drawn to buld an uncertainty model. \n",
    "\n",
    "\\begin{equation}\n",
    "X^\\ell = F_x(p^\\ell),  \\, \\forall \\, \\ell = 1,\\ldots, L\n",
    "\\end{equation}\n",
    "\n",
    "where $X^\\ell$ is the realization of the variable $X$ drawn from its CDF, $F_x$, with cumulative probability, p-value, $p^\\ell$.  \n",
    "\n",
    "#### Monter Carlo Simulation for Uncertainty Modeling\n",
    "\n",
    "It would be trivial to apply Monte Carlo simulation to a single random variable, after many realizations one would get back the original distribution. The general approach is to:\n",
    "\n",
    "1. Model all distributions for the input, variables of interest $F_{x_1},\\ldots,F_{x_m}$.\n",
    "2. For each realization draw $p^\\ell_{1},\\ldots,p^\\ell_{m}$, p-values\n",
    "3. Apply the inverse of each distribution to calculate a realization of each variable, $X^\\ell_j = F_{x^\\ell_j}^{-1}(p^\\ell_j),  \\, \\forall \\, j = 1,\\ldots$, $m$ variables.\n",
    "4. Apply each set of variables for a $\\ell$ realization to the transfer function to calculate the ouptput realization, $Y^\\ell = F(X_1^\\ell,\\ldots,X_m^\\ell)$.\n",
    "\n",
    "Monte Carlo Simulation (MCS) is extremely powerful\n",
    "\n",
    "* Possible to easily simulate uncertainty models for complicated systems \n",
    "* Simulations are conducted by drawing values at random from specified uncertainty distributions for each variable\n",
    "* A single realization of each variable, $X_1^\\ell, X_2^\\ell,\\ldots,X_m^\\ell$ is applied to the transfer function to calculate the realization of the variable of interest (output, decision criteria):\n",
    "\n",
    "\\begin{equation}\n",
    "Y^\\ell = f(X_1^\\ell,\\ldots,X_m^\\ell), \\, \\forall \\, \\ell = 1,\\ldots, L\n",
    "\\end{equation}\n",
    "\n",
    "* The MCS method builds empirical uncertainty models by random sampling\n",
    "\n",
    "How many realizations, $L$?\n",
    "\n",
    "The answer is enough! If the MCS computational cost is low then **many** is the right answer. If too few realizations are calculated then the summary statistics and the entire CDF of the output, decision criteria may be incorrect. This is caused by fluctuations due to not enough samples (see the 'Law of Small Numbers').\n",
    "\n",
    "The MCS method is very powerful. You can simulate output distributions that could not be calculated analytically.  \n",
    "\n",
    "#### Limitations\n",
    "\n",
    "The MCS method above assumes:\n",
    "1. **representativity** - the distribution is representative\n",
    "2. **independence** - the variables are independent of eachother\n",
    "3. **stationarity** - all realizations for each variable are from the same distribution \n",
    "  \n",
    "#### Getting Started\n",
    "\n",
    "Here's the steps to get setup in Python with the GeostatsPy package:\n",
    "\n",
    "1. Install Anaconda 3 on your machine (https://www.anaconda.com/download/). \n",
    "2. From Anaconda Navigator (within Anaconda3 group), go to the environment tab, click on base (root) green arrow and open a terminal. \n",
    "3. In the terminal type: pip install geostatspy. \n",
    "4. Open Jupyter and in the top block get started by copy and pasting the code block below from this Jupyter Notebook to start using the geostatspy functionality. \n",
    "\n",
    "There are examples below with these functions. You can go here to see a list of the available functions, https://git.io/fh4eX, other example workflows and source code. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geostatspy.GSLIB as GSLIB          # GSLIB utilies, visualization and wrapper\n",
    "import geostatspy.geostats as geostats    # GSLIB methods convert to Python        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also need some standard packages. These should have been installed with Anaconda 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                        # ndarrys for gridded data\n",
    "import pandas as pd                       # DataFrames for tabular data\n",
    "import os                                 # set working directory, run executables\n",
    "import matplotlib.pyplot as plt           # for plotting\n",
    "from scipy import stats                   # summary statistics\n",
    "import math                               # trig etc.\n",
    "import random\n",
    "from ipywidgets import interactive        # widgets and interactivity\n",
    "from ipywidgets import widgets                            \n",
    "from ipywidgets import Layout\n",
    "from ipywidgets import Label\n",
    "from ipywidgets import VBox, HBox"
   ]
  },
  {
   "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": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# interactive calculation of the random sample set (control of source parametric distribution and number of samples)\n",
    "l = widgets.Text(value='                                      Monte Carlo Simulation Demonstration, Michael Pyrcz, Associate Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "operator = widgets.RadioButtons(options=['add', 'mult'],description='Operator:',disabled=False,layout=Layout(width='230px', height='50px'))\n",
    "\n",
    "L = widgets.IntSlider(min=1, max = 10000, value = 50, description = '$L$:',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "L.style.handle_color = 'gray'\n",
    "\n",
    "uiL = widgets.VBox([L,operator])\n",
    "\n",
    "dist1 = widgets.Dropdown(\n",
    "    options=['Uniform','Triangular','Gaussian'],\n",
    "    value='Gaussian',\n",
    "    description='$X_1$:',\n",
    "    disabled=False,\n",
    "    layout=Layout(width='200px', height='30px')\n",
    ")\n",
    "min1 = widgets.FloatSlider(min=0.0, max = 100.0, value = 10.0, description = 'Min',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "min1.style.handle_color = 'blue'\n",
    "max1 = widgets.FloatSlider(min=0.0, max = 100.0, value = 30.0, description = 'Max',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "max1.style.handle_color = 'blue'\n",
    "ui1 = widgets.VBox([dist1,min1,max1],kwargs = {'justify_content':'center'}) \n",
    "\n",
    "dist2 = widgets.Dropdown(\n",
    "    options=['Triangular', 'Uniform', 'Gaussian'],\n",
    "    value='Gaussian',\n",
    "    description='$X_2$:',\n",
    "    disabled=False,\n",
    "    layout=Layout(width='200px', height='30px')\n",
    ")\n",
    "min2 = widgets.FloatSlider(min=0.0, max = 100.0, value = 10.0, description = 'Min',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "min2.style.handle_color = 'red'\n",
    "max2 = widgets.FloatSlider(min=0.0, max = 100.0, value = 30.0, description = 'Max',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "max2.style.handle_color = 'red'\n",
    "ui2 = widgets.VBox([dist2,min2,max2],kwargs = {'justify_content':'center'})\n",
    "\n",
    "dist3 = widgets.Dropdown(\n",
    "    options=['Triangular', 'Uniform', 'Gaussian'],\n",
    "    value='Gaussian',\n",
    "    description='$X_3$:',\n",
    "    disabled=False,\n",
    "    layout=Layout(width='200px', height='30px')\n",
    ")\n",
    "min3 = widgets.FloatSlider(min=0.0, max = 100.0, value = 10.0, description = 'Min',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "min3.style.handle_color = 'yellow'\n",
    "max3 = widgets.FloatSlider(min=0.0, max = 100.0, value = 30.0, description = 'Max',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "max3.style.handle_color = 'yellow'\n",
    "ui3 = widgets.VBox([dist3,min3,max3],kwargs = {'justify_content':'center'})\n",
    "\n",
    "ui = widgets.HBox([uiL,ui1,ui2,ui3])\n",
    "ui2 = widgets.VBox([l,ui],)\n",
    "\n",
    "def make_dist(dist,zmin,zmax,L):\n",
    "    if dist == 'Triangular':\n",
    "        z = np.random.triangular(left=zmin, mode=(zmax+zmin)*0.5, right=zmax, size=L)\n",
    "        pdf = stats.triang.pdf(np.linspace(0.0,100.0,1000), loc = zmin, c = 0.5, scale = zmax-zmin)* 2 * L \n",
    "    if dist == 'Uniform':\n",
    "        z = np.random.uniform(low=zmin, high=zmax, size=L)\n",
    "        pdf = stats.uniform.pdf(np.linspace(0.0,100.0,1000), loc = zmin, scale = zmax-zmin) * 2 * L\n",
    "    if dist == 'Gaussian':\n",
    "        mean = (zmax + zmin)*0.5; sd = (zmax - zmin)/6.0\n",
    "        z = np.random.normal(loc = mean, scale = sd, size=L)\n",
    "        pdf = stats.norm.pdf(np.linspace(0.0,100.0,1000), loc = mean, scale = sd) * 2 * L\n",
    "    return z, pdf\n",
    "        \n",
    "def f_make(L,operator,dist1,min1,max1,dist2,min2,max2,dist3,min3,max3): \n",
    "    np.random.seed(seed = 73073)\n",
    "    x1, pdf1 = make_dist(dist1,min1,max1,L)\n",
    "    x2, pdf2 = make_dist(dist2,min2,max2,L)\n",
    "    x3, pdf3 = make_dist(dist3,min3,max3,L)\n",
    "\n",
    "    xvals = np.linspace(0.0,100.0,1000)\n",
    "    plt.subplot(241)\n",
    "    plt.hist(x1,density = False,bins=np.linspace(0,100,50),weights=None,color='blue',alpha=0.3,edgecolor='black')\n",
    "    plt.plot(xvals,pdf1,'--',color='black')\n",
    "    plt.xlim(0,100); plt.xlabel(\"$X_1$\"); plt.title(\"First Feature, $X_1$\"); plt.ylabel('Frequency')\n",
    " \n",
    "    plt.subplot(242)\n",
    "    plt.hist(x2,density = False,bins=np.linspace(0,100,50),weights=None,color='red',alpha=0.3,edgecolor='black')\n",
    "    plt.plot(xvals,pdf2,'--',color='black')\n",
    "    plt.xlim(0,100); plt.xlabel(\"$X_1$\"); plt.title(\"Second Feature, $X_2$\"); plt.ylabel('Frequency')\n",
    " \n",
    "    plt.subplot(243)\n",
    "    plt.hist(x3,density = False,bins=np.linspace(0,100,50),weights=None,color='yellow',alpha=0.3,edgecolor='black')\n",
    "    plt.plot(xvals,pdf3,'--',color='black')\n",
    "    plt.xlim(0,100); plt.xlabel(\"$X_1$\"); plt.title(\"Third Feature, $X_3$\"); plt.ylabel('Frequency')\n",
    " \n",
    "    y = np.zeros(L)\n",
    "    ymin = 0.0\n",
    "    if operator == \"add\":\n",
    "        y = x1 + x2 + x3\n",
    "    elif operator == \"mult\":\n",
    "        y = x1 * x2 * x3\n",
    "        \n",
    "    ymax = max(round((np.max(y)+50)/100)*100,100) # round up to nearest hundreds to avoid the chart jumping around\n",
    "    \n",
    "    plt.subplot(244)\n",
    "    plt.hist(y,density = False,bins=np.linspace(ymin,ymax,50),weights=None,color='black',alpha=0.5,edgecolor='black')\n",
    "    plt.xlabel(\"$Y$\"); plt.title(\"Response Feature, $Y$\"); plt.ylabel('Frequency')\n",
    "    plt.xlim(ymin,ymax)\n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.3, hspace=0.2)\n",
    "    plt.show()    \n",
    "   \n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "interactive_plot = widgets.interactive_output(f_make, {'L':L,'operator':operator,'dist1':dist1,'min1':min1,'max1':max1,'dist2':dist2,'min2':min2,'max2':max2,'dist3':dist3,'min3':min3,'max3':max3})\n",
    "interactive_plot.clear_output(wait = True)                # reduce flickering by delaying plot updating    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Monte Carlo Simulation Demonstration\n",
    "\n",
    "* specify the distributions for 3 Random Variables, $X_1$, $X_2$, and $X_3$ and select the operator $y = f(X_1,X_2,X_3)$\n",
    "\n",
    "* observe the distribution of the resulting Monte Carlos Simulation realization histograms of $X_1^{\\ell}$, $X_2^{\\ell}$, $X_3^{\\ell}$, and $y^{\\ell}$\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) | [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy)\n",
    "\n",
    "### The Inputs\n",
    "\n",
    "* **$L$**: number of realizations, **Operator**: addition for $y = X_1 + X_2 + X_3$, multiplication for $y = X_1 \\times X_2 \\times X_3$\n",
    "\n",
    "* **$X_1$, $X_2$, and $X_3$**: distribution type, min and max. Assume mode or mean is centered and 3 st.dev. for Gaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c6da3815eb554b9c8fceaa202d956afa",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                      Monte Carlo Simulation Demonstration, Michael…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6190088f3cae4a24bc78f48562f34d6b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output(outputs=({'output_type': 'display_data', 'data': {'text/plain': '<Figure size 432x288 with 4 Axes>', 'i…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui2, interactive_plot)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This was a basic demonstration of Monte Carlo simulation for uncertainty analysis. A lot more could be done, for example, more complicated transfer functions and a combination of non-parametric and parametric distributions. Also, one could integrate relationships between the variables (we assumed independent here).\n",
    "\n",
    "I have other demonstrations on the basics of working with DataFrames, ndarrays, univariate statistics, plotting data, declustering, data transformations, trend modeling, multivariate analysis and many other workflows available at https://github.com/GeostatsGuy/PythonNumericalDemos and https://github.com/GeostatsGuy/GeostatsPy. \n",
    "  \n",
    "I hope this was helpful,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "#### The Author:\n",
    "\n",
    "### Michael Pyrcz, Associate Professor, University of Texas at Austin \n",
    "*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*\n",
    "\n",
    "With over 17 years of experience in subsurface consulting, research and development, Michael has returned to academia driven by his passion for teaching and enthusiasm for enhancing engineers' and geoscientists' impact in subsurface resource development. \n",
    "\n",
    "For more about Michael check out these links:\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)\n",
    "\n",
    "#### Want to Work Together?\n",
    "\n",
    "I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.\n",
    "\n",
    "* Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I'd be happy to drop by and work with you! \n",
    "\n",
    "* Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!\n",
    "\n",
    "* I can be reached at mpyrcz@austin.utexas.edu.\n",
    "\n",
    "I'm always happy to discuss,\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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}