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",
    "### Monte Carlo Simulation in Python \n",
    "\n",
    "#### Michael Pyrcz, Professor, The 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\n",
    "\n",
    "Here's two dashboards to demonstration of Monte Carlo simulation for subsurface uncertainty modeling workflows. This should help you get started with:\n",
    "\n",
    "* building data science workflows integrate uncertainty sources, i.e, math with distributions\n",
    "* using more advanced data science workflows like bootstrap and Markov chain Monte Carlo (McMC)\n",
    "\n",
    "I have recorded a walk-through of this interactive dashboard in my [Data Science Interactive Python Demonstrations](https://www.youtube.com/playlist?list=PLG19vXLQHvSDy26fM3hDLg3VCU7U5BGZl) series on my [YouTube](https://www.youtube.com/@GeostatsGuyLectures) channel.\n",
    "\n",
    "* Join me for walk-through of this dashboard [09 Data Science Interactive: Monte Carlo Simulation](TBD). I'm stoked to guide you and share observations and things to try out!   \n",
    "\n",
    "* I have a lecture on [Monte Carlo Simulation](https://www.youtube.com/watch?v=Qb8TsSINpnU&list=PLG19vXLQHvSB-D4XKYieEku9GQMQyAzjJ&index=13&t=3s) as part of my [Data Analytics and Geostatistics](https://www.youtube.com/playlist?list=PLG19vXLQHvSB-D4XKYieEku9GQMQyAzjJ) course. Note, for all my recorded lecture the interactive and well-documented workflow demononstrations are available on my GitHub repository [GeostatsGuy's Python Numerical Demos](https://github.com/GeostatsGuy/PythonNumericalDemos).\n",
    "\n",
    "* Also, I have a lecture on [Markov chain Monte Carlo](https://www.youtube.com/watch?v=7QX-yVboLhk&list=PLG19vXLQHvSC2ZKFIkgVpI9fCjkN38kwf&index=34) as part of Bayesian regression in my [Machine Learning](https://www.youtube.com/playlist?list=PLG19vXLQHvSC2ZKFIkgVpI9fCjkN38kwf) course.\n",
    "\n",
    "#### Monte Carlo Simulation\n",
    "\n",
    "Definition: random sampling from a distribution.\n",
    "\n",
    "Procedure: \n",
    "\n",
    "1. Model the representative distribution (CDF), $F_X(x)$\n",
    "2. Draw a random value from a uniform [0,1] distribution (p-value), $p^\\ell$\n",
    "3. Apply the inverse of the CDF to calculate the associated realization, $F^{-1}_X(p^\\ell)$\n",
    "\n",
    "In practice, Monte Carlo Simulation (MCS) refers to the workflow with multiple realizations drawn to buld an uncertainty model. \n",
    "\n",
    "\\begin{equation}\n",
    "x^\\ell = F^{-1}_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(x)$, with cumulative probability, p-value, $p^\\ell$.  \n",
    "\n",
    "#### Monte Carlo Simulation for Uncertainty Modeling\n",
    "\n",
    "It would be trivial to apply MCS to a single random variable, after many realizations one would get back the original distribution. But this general method allows us to do math with distributions, random variables, to calculate uncertainty models! \n",
    "\n",
    "Definitions:\n",
    "\n",
    "* **random variable** - a feature's value at a time or location is unknown and take on a range of possible values, $X_1$, usually defined by a PDF, $f_{X_1}(x_1)$, or CDF, $F_{X_1}(x_1)$ \n",
    "\n",
    "* **data value** - a sample observation at a specific location or time for a feature, $x_1\\left(\\bf{u}_{\\alpha} \\right)$ \n",
    "\n",
    "* **realization** - a random sample drawn from a random variable, $x_{\\alpha}^{\\ell}$ \n",
    "\n",
    "With MCS we can solve problems like this:\n",
    "\n",
    "$Y = X_1 + X_2 + X_3$ given $X_1$, $X_2$, and $X_3$ are random variables with any distribution, parametric or empirical (a set of values known as a reference distribution).\n",
    "\n",
    "#### Monte Carlo Simulation for Uncertainty Modeling General Workflow\n",
    "\n",
    "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 realizations! \n",
    "\n",
    "* If the MCS computational cost is low then **many** is the right answer. \n",
    "\n",
    "* 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 Simulation Perspective\n",
    "\n",
    "There are many uncertainty operations that we can work out analytically, for which we don't need MCS:\n",
    "\n",
    "* $X_1 + c$\n",
    "* $c \\cdot X_1$\n",
    "* $X_1 + X_2$, given $X_1$ and $X_2$ are Gaussian distributed\n",
    "\n",
    "But consider the following:\n",
    "\n",
    "* $X_1 \\times X_2$, given $X_1$ and $X_2$ could have any distribution\n",
    "* what is the probability of winning the game of solitaire?\n",
    "\n",
    "With MCS we simulate the system, sample $L$ outcomes and then summarize!\n",
    "\n",
    "#### Limitations\n",
    "\n",
    "The MCS method above assumes:\n",
    "1. **representativity** - the distributions are representative\n",
    "2. **independence** - the random variables are independent of eachother\n",
    "3. **stationarity** - all realizations for each variable are from the same distribution, e.g., no change over time nor space \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",
    "\n",
    "That's all!\n",
    "\n",
    "#### Load the Required Libraries\n",
    "\n",
    "We will also need some standard Python packages. These should have been installed with Anaconda 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "supress_warnings = True                                   # supress warnings?\n",
    "import numpy as np                                        # ndarrys for gridded data\n",
    "import matplotlib.pyplot as plt                           # for plotting\n",
    "from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks\n",
    "from scipy import stats                                   # summary statistics\n",
    "import math                                               # math operators, e.g., square root and natural log\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\n",
    "cmap = plt.cm.inferno                                     # default color bar, no bias and friendly for color vision defeciency\n",
    "plt.rc('axes', axisbelow=True)                            # grid behind plotting elements\n",
    "if supress_warnings == True:\n",
    "    import warnings                                       # supress any warnings for this demonstration\n",
    "    warnings.filterwarnings('ignore')                  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you get a package import error, you may have to first install some of these packages. This can usually be accomplished by opening up a command window on Windows and then typing 'python -m pip install [package-name]'. More assistance is available with the respective package docs.  \n",
    "\n",
    "#### Declare Functions\n",
    "\n",
    "I just added a convenience function for adding major and minor gridlines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_grid_orig():\n",
    "    plt.gca().grid(True, which='major',linewidth = 1.0); plt.gca().grid(True, which='minor',linewidth = 0.2) # add y grids\n",
    "    plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)\n",
    "    plt.gca().xaxis.set_minor_locator(AutoMinorLocator()); plt.gca().yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks\n",
    "    \n",
    "def add_grid(sub_plot):\n",
    "    sub_plot.grid(True, which='major',linewidth = 1.0); sub_plot.grid(True, which='minor',linewidth = 0.2) # add y grids\n",
    "    sub_plot.tick_params(which='major',length=7); sub_plot.tick_params(which='minor', length=4)\n",
    "    sub_plot.xaxis.set_minor_locator(AutoMinorLocator()); sub_plot.yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks\n",
    "    \n",
    "def sample_dist(mean,stdev,L,dist):\n",
    "    pvalue = np.random.rand(L)\n",
    "    if dist == 'log':\n",
    "        mu = math.log(mean**2/(math.sqrt(stdev**2 + mean**2)))\n",
    "        sigma = math.sqrt(math.log((stdev**2)/(mean**2)+1.0))\n",
    "        z = stats.lognorm.ppf(pvalue,s = sigma, scale = math.exp(mu)) \n",
    "    elif dist == 'norm':\n",
    "        z = stats.norm.ppf(pvalue,loc = mean, scale = stdev) \n",
    "    return pvalue,z\n",
    "\n",
    "def plot_dist(mean,stdev,dist):\n",
    "    p = np.linspace(0.001,0.999,1000)\n",
    "    if dist == 'log':\n",
    "        mu = math.log(mean**2/(math.sqrt(stdev**2 + mean**2)))\n",
    "        sigma = math.sqrt(math.log((stdev**2)/(mean**2)+1.0))\n",
    "        z = stats.lognorm.ppf(p,s = sigma, scale = math.exp(mu))\n",
    "    elif dist == 'norm':\n",
    "        z = stats.norm.ppf(p,loc = mean, scale = stdev) \n",
    "    return p,z\n",
    "\n",
    "def plot_input(mean,stdev,dist,xmin,xmax,L,title,xlabel,seed,ax):\n",
    "    np.random.seed(seed=seed)\n",
    "    pvalue, real = sample_dist(mean,stdev,L,dist)\n",
    "    np.random.seed(seed=seed)\n",
    "    cdf, val = plot_dist(mean,stdev,dist)\n",
    "    ax.plot(val,cdf,color='black',lw=2.0)\n",
    "    \n",
    "    for l in range(0,L):\n",
    "        alpha = max(0.02,(6 - (L-l) )/5)\n",
    "        dhead = 1;lhead = 0.02        \n",
    "        if l == L-1:\n",
    "            ax.plot([0.0,real[l]],[pvalue[l],pvalue[l]],color='red',alpha = alpha,lw=1,zorder=1)\n",
    "            ax.plot([real[l],real[l]],[pvalue[l],lhead+0.02],color='red',alpha = alpha,lw=1,zorder=1)\n",
    "            ax.annotate(format(int(real[l]), ','),xy=[real[l]+(xmax-xmin)*0.01,0.07],rotation=90.0,color='red',alpha = alpha)\n",
    "            ax.scatter(real[l],0.01,s=30,color='red',edgecolor='black',alpha=1.0)\n",
    "        else:\n",
    "            ax.plot([0.0,real[l]],[pvalue[l],pvalue[l]],color='darkorange',alpha = alpha,lw=1,zorder=1)\n",
    "            ax.plot([real[l],real[l]],[pvalue[l],lhead+0.02],color='darkorange',alpha = alpha,lw=1,zorder=1)\n",
    "            ax.annotate(format(int(real[l]), ','),xy=[real[l]+(xmax-xmin)*0.01,0.07],rotation=90.0,color='darkorange',alpha = alpha)\n",
    "            ax.scatter(real[l],0.01,s=10,color='yellow',edgecolor='black',alpha=1.0)\n",
    "\n",
    "    ax.set_ylim([0,1]) \n",
    "    ax.set_xlim([xmin,xmax])\n",
    "    if dist == 'log':\n",
    "        ax.set_xscale('log')\n",
    "        ax.grid(which='both',axis='both',visible=True)\n",
    "        ax.grid(True, which='major',linewidth = 1.0); ax.grid(True, which='minor',linewidth = 0.2)\n",
    "    elif dist == 'norm':\n",
    "        ax.set_xscale('linear');\n",
    "        add_grid(ax)        \n",
    "    ax.set_xlabel(xlabel); ax.set_ylabel('Cumulative Probability'); ax.set_title(title)\n",
    "    return real"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple Monte Carlo Simulation Demonstration\n",
    "\n",
    "Now let's set up our first dashboard."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = widgets.Text(value='                                          Monte Carlo Simulation Demonstration, Michael Pyrcz, Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "L = widgets.IntSlider(min=1, max = 40, value = 1, description = '$L$:',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "L.style.handle_color = 'gray'\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 = 'gray'\n",
    "max1 = widgets.FloatSlider(min=0.0, max = 100.0, value = 90.0, description = 'Max',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "max1.style.handle_color = 'gray'\n",
    "\n",
    "ui = widgets.HBox([L,dist1,min1,max1],kwargs = {'justify_content':'center'})\n",
    "ui1 = widgets.VBox([l,ui],)\n",
    "\n",
    "def make_sample_cdist(dist,zmin,zmax,L):\n",
    "    cdf = np.linspace(0.001,0.999,1000)\n",
    "    pvals = np.random.random(size = L)\n",
    "    if dist == 'Triangular':  \n",
    "        distribution = stats.triang(loc = zmin,c=0.5,scale = zmax-zmin)\n",
    "    if dist == 'Uniform':        \n",
    "        distribution = stats.uniform(loc = zmin,scale = zmax-zmin)\n",
    "    if dist == 'Gaussian':\n",
    "        mean = (zmax + zmin)*0.5; stdev = (zmax - zmin)/6.0\n",
    "        distribution = stats.norm(loc = mean,scale = stdev)\n",
    "    \n",
    "    cdfx = distribution.ppf(cdf)\n",
    "    sample = distribution.ppf(pvals)\n",
    "    \n",
    "    return sample, pvals, cdfx, cdf \n",
    "        \n",
    "def f_make1(L,dist1,min1,max1): \n",
    "    np.random.seed(seed = 73073)\n",
    "    sample, pvals, cdfx, cdf = make_sample_cdist(dist1,min1,max1,L)\n",
    "    \n",
    "    plt.subplot(121)\n",
    "    plt.plot(cdfx,cdf,'--',color='black',linewidth = 3,zorder=10); add_grid_orig()\n",
    "    plt.xlim(0,100); plt.ylim(0,1.0); plt.xlabel(\"$X_1$\"); plt.title(\"Cumulative Distribution Function, $X_1$\"); plt.ylabel('Cumulative Probability')\n",
    "     \n",
    "    for l in range(0,L):\n",
    "        alpha = max(0.02,(6 - (L-l) )/5)\n",
    "        dhead = 1;lhead = 0.02        \n",
    "        if l == L-1:\n",
    "            plt.plot([0.0,sample[l]],[pvals[l],pvals[l]],color='red',alpha = alpha,lw=1,zorder=1)\n",
    "            plt.plot([sample[l],sample[l]],[pvals[l],lhead+0.02],color='red',alpha = alpha,lw=1,zorder=1)\n",
    "            head = plt.Polygon([[sample[l],0.02],[sample[l]-dhead,lhead+0.02],[sample[l]+dhead,lhead+0.02],[sample[l],0.02]], color='red',alpha = alpha,zorder=1)\n",
    "            plt.annotate(int(sample[l]),[sample[l]+1,0.07],rotation=90.0,color='red',alpha = alpha)\n",
    "            plt.scatter(sample[l],0.01,s=30,color='red',edgecolor='black',alpha=1.0)\n",
    "            plt.gca().add_patch(head)\n",
    "        else:\n",
    "            plt.plot([0.0,sample[l]],[pvals[l],pvals[l]],color='darkorange',alpha = alpha,lw=1,zorder=1)\n",
    "            plt.plot([sample[l],sample[l]],[pvals[l],lhead+0.02],color='darkorange',alpha = alpha,lw=1,zorder=1)\n",
    "            head = plt.Polygon([[sample[l],0.02],[sample[l]-dhead,lhead+0.02],[sample[l]+dhead,lhead+0.02],[sample[l],0.02]], color='darkorange',alpha = alpha,zorder=1)\n",
    "            plt.annotate(int(sample[l]),[sample[l]+1,0.07],rotation=90.0,color='darkorange',alpha = alpha)\n",
    "            plt.scatter(sample[l],0.01,s=10,color='yellow',edgecolor='black',alpha=1.0)\n",
    "            plt.gca().add_patch(head)            \n",
    "\n",
    "    plt.subplot(122)\n",
    "    plt.hist(sample,density= False,bins=np.linspace(0,100,20),weights=None,color='red',alpha=1.0,edgecolor='black',zorder=1)\n",
    "    if l > 0:\n",
    "        plt.hist(sample[:-1],density= False,bins=np.linspace(0,100,20),weights=None,color='darkorange',alpha=0.8,edgecolor='black',zorder=5)\n",
    "    \n",
    "    for l in range(0,L):\n",
    "        if l == L-1:\n",
    "            plt.scatter(sample[l],0.1,s=30,color='red',edgecolor='black',alpha=1.0,zorder=l+10)\n",
    "        else:\n",
    "            plt.scatter(sample[l],0.1,s=10,color='yellow',edgecolor='black',alpha=1.0,zorder=l+10)\n",
    "    \n",
    "    plt.xlabel(r\"$x_{\\alpha}, \\alpha = 1,\\ldots,L$\"); plt.title(r\"Monte Carlo Simulation Realizations, $x_{\\alpha}, \\alpha = 1,\\ldots,L$\"); plt.ylabel('Frequency')\n",
    "    plt.xlim([0,100]); plt.ylim([0,10]); add_grid_orig()\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",
    "interactive_plot1 = widgets.interactive_output(f_make1, {'L':L,'dist1':dist1,'min1':min1,'max1':max1})\n",
    "interactive_plot1.clear_output(wait = True)                # reduce flickering by delaying plot updating  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Monte Carlo Simulation Demonstration\n",
    "\n",
    "* $L$ random samples from a parametric distribuiton, $X_1$.\n",
    "\n",
    "#### Michael Pyrcz, Professor, The 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, **$X_1$**: type of distribution, parameterized by the **min**: minimum, and **max**: maximum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2c639cc3c2a54c72a0e2b173262056aa",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                          Monte Carlo Simulation Demonstration, Mic…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b114c3c11b894cdbb05cd604cc9ea82e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui1, interactive_plot1)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Monte Carlo Simulation Demonstration\n",
    "\n",
    "Now let's set up our second dashboard for adding or multiplying distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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, 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='100px'))\n",
    "\n",
    "L = widgets.IntSlider(min=1, max = 10000, value = 1, 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",
    "    \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.7,edgecolor='grey')\n",
    "    plt.plot(xvals,pdf1,'--',color='black',linewidth = 3)\n",
    "    plt.xlim(0,100); plt.xlabel(\"$X_1$\"); plt.title(\"First Predictor 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.7,edgecolor='grey')\n",
    "    plt.plot(xvals,pdf2,'--',color='black',linewidth = 3)\n",
    "    plt.xlim(0,100); plt.xlabel(\"$X_1$\"); plt.title(\"Second Predictor 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.7,edgecolor='grey')\n",
    "    plt.plot(xvals,pdf3,'--',color='black',linewidth = 3)\n",
    "    plt.xlim(0,100); plt.xlabel(\"$X_1$\"); plt.title(\"Third Predictor 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.ylabel('Frequency'); plt.xlim(ymin,ymax)\n",
    "    if operator == \"Add\":\n",
    "        plt.title(\"Response Feature, $y = X_1 + X_2 + X_3$\")\n",
    "    else:\n",
    "        plt.title(r\"Response Feature, $y = X_1 \\times X_2 \\times X_3$\")\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",
    "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    "
   ]
  },
  {
   "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, Professor, The 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": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "355955e223f142bd9c604fda907a69f8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                          Monte Carlo Simulation Demonstration, Mic…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d96c14bb103a47c294964b965fd7ea56",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output(outputs=({'output_type': 'display_data', 'data': {'text/plain': '<Figure size 640x480 with 4 Axes>', 'i…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui2, interactive_plot)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Monte Carlo Simulation Lithium Salar Brine Resource Uncertainty Model Demonstration\n",
    "\n",
    "Now let's buld an uncertainty model for a lithium salar brine resource. How much lithium is available for extraction?\n",
    "\n",
    "* the lithium is disolved in the brine \n",
    "* the brine is in the pore space of a generally flat aquifer bounded by aquicludes below and above, and low permeability rock, pinch outs or fresh water transition around the periphery\n",
    "* we need to calculate the mass of lithium available.\n",
    "\n",
    "The predictor features include:\n",
    "\n",
    "* **area** ($A$) - of the salar brine in units, $m^2$\n",
    "* **average thickness** ($\\overline{th}$) - the average thickness over the area of the salar brine of the bound aquifer\n",
    "* **average effective porosity** ($\\overline{\\phi}_{e}$) - average effective porosity over the bound aquifer \n",
    "* **specific yield** ($s_y$) - the proportion of brine that is receoverable under gravity flow conditions\n",
    "* **average lithium concentation** ($\\overline{[Li]}$} - amount of disssolved lithium, in $\\frac{mg}{L}$. E.g., average lithium concentations, Salton Sea 211 $\\frac{mg}{L}$, Westmorland 48 $\\frac{mg}{L}$, and Brawley 100 $\\frac{mg}{L}$.\n",
    "\n",
    "We will also use the following unit conversions:\n",
    "\n",
    "* $\\frac{1000 L}{m^3}$ - to convert from litres to meter cubed\n",
    "* $\\frac{tonne}{1.0^9 mg}$ - to convert from miligrams to metric tons, tonnes (1,000 kilograms)\n",
    "\n",
    "Now we have the follow transfer function to calculate our lithium resource, with units shown:\n",
    "\n",
    "\\begin{equation}\n",
    "m_{Li} \\ (tonne) = A \\ (m^2) \\cdot \\overline{th} \\ (m) \\cdot \\overline{\\phi}_{e} \\ (\\%) \\cdot \\frac{1}{100\\%} \\cdot 1,000 \\left(\\frac{L}{m^3}\\right)  \\cdot \\overline{[Li]} \\ \\left(\\frac{mg}{L}\\right) \\cdot \\frac{1}{1.0^9} \\left(\\frac{tonne}{mg}\\right) \\cdot s_y \\ (fraction)\n",
    "\\end{equation}\n",
    "\n",
    "We can simplify this lithium resource transfer function to:\n",
    "\n",
    "\\begin{equation}\n",
    "m_{Li} = A \\cdot \\overline{th} \\cdot \\overline{\\phi}_{e} \\cdot \\overline{[Li]} \\cdot s_y \\cdot \\frac{1}{1.0^8} \\ (\\text{tonnes of lithium resource})\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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 Lithium Resources Demonstration, Michael Pyrcz, Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "L = widgets.IntSlider(min=1, max = 100, value = 1, description = '$L$:',orientation='horizontal',\n",
    "                      layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "L.style.handle_color = 'gray'\n",
    "\n",
    "A_mean = widgets.IntSlider(min=1000, max = 10000000, value = 1000000, description = '$\\mu_{A}$',\n",
    "                             orientation='horizontal',layout=Layout(width='230px', height='50px'),\n",
    "                             readout_format=',',continuous_update=False)\n",
    "A_mean.style.handle_color = 'blue'\n",
    "\n",
    "A_stdev = widgets.IntSlider(min=0, max = 10000000, value = 300000, description = '$\\sigma_{A}$',\n",
    "                            orientation='horizontal',layout=Layout(width='230px', height='50px'),\n",
    "                            readout_format=',',continuous_update=False)\n",
    "A_stdev.style.handle_color = 'blue'\n",
    "\n",
    "ui_A = widgets.VBox([A_mean,A_stdev],kwargs = {'justify_content':'center'}) \n",
    "\n",
    "th_mean = widgets.IntSlider(min=1, max = 1000, value = 100, description = '$\\mu_{\\overline{th}}$',\n",
    "                              orientation='horizontal',layout=Layout(width='230px', height='50px'),\n",
    "                              readout_format=',',continuous_update=False)\n",
    "th_mean.style.handle_color = 'darkblue'\n",
    "\n",
    "th_stdev = widgets.IntSlider(min=0, max = 300, value = 60, description = '$\\sigma_{\\overline{th}}$',\n",
    "                               orientation='horizontal',layout=Layout(width='230px', height='50px'),\n",
    "                               readout_format=',',continuous_update=False)\n",
    "th_stdev.style.handle_color = 'darkblue'\n",
    "\n",
    "ui_th = widgets.VBox([th_mean,th_stdev],kwargs = {'justify_content':'center'}) \n",
    "\n",
    "li_mean = widgets.IntSlider(min=0, max = 400, value = 200, description = '$\\mu_{\\overline{Li}}$',\n",
    "                               orientation='horizontal',layout=Layout(width='230px', height='50px'),\n",
    "                               readout_format=',',continuous_update=False)\n",
    "li_mean.style.handle_color = 'red'\n",
    "\n",
    "li_stdev = widgets.IntSlider(min=0, max = 200, value = 50, description = '$\\sigma_{\\overline{Li}}$',\n",
    "                                orientation='horizontal',layout=Layout(width='230px', height='50px'),\n",
    "                                readout_format=',',continuous_update=False)\n",
    "li_stdev.style.handle_color = 'red'\n",
    "\n",
    "ui_li = widgets.VBox([li_mean,li_stdev],kwargs = {'justify_content':'center'}) \n",
    "\n",
    "ui_all = widgets.HBox([L,ui_A,ui_th,ui_li])\n",
    "ui4 = widgets.VBox([l,ui_all],)\n",
    "\n",
    "\n",
    "def f_make4(L,A_mean,A_stdev,th_mean,th_stdev,li_mean,li_stdev): \n",
    "    seed = 73073\n",
    "    np.random.seed(seed = seed)\n",
    "    \n",
    "    phi_const = 0.2\n",
    "    sy_const = 0.4       \n",
    "    \n",
    "    plt_x1 = plt.subplot2grid((3, 3), (0, 0), rowspan=1, colspan=1)\n",
    "    plt_x2 = plt.subplot2grid((3, 3), (1, 0), rowspan=1, colspan=1)\n",
    "    plt_x3 = plt.subplot2grid((3, 3), (2, 0), rowspan=1, colspan=1) \n",
    "    plt_y = plt.subplot2grid((3, 3), (0, 1), rowspan=3, colspan=2) \n",
    "    \n",
    "    np.random.seed(seed = 73073)\n",
    "    \n",
    "    A = plot_input(A_mean,A_stdev,'log',10000,10000000,L,r'Area ($A$)',r'Area ($m^2$)',seed,plt_x1)\n",
    "    th = plot_input(th_mean,th_stdev,'log',1,1000,L,r'Average Thickness ($\\overline{th}$)',r'Average Thickness ($m$)',seed,plt_x2)\n",
    "    phi = np.full(L,phi_const)\n",
    "    li = plot_input(li_mean,li_stdev,'norm',0.0,400.0,L,r'Average Lithium Grade ($\\overline{[Li]}$)',r'Average Lithium Grade (mg/L)',seed,plt_x3)\n",
    "    sy = np.full(L,sy_const)\n",
    "    \n",
    "    resource = (A * th * phi * li * sy) / 1.0e8\n",
    "    last = resource[-1]    \n",
    "    if L > 1:\n",
    "        last_pvalue = resource.argsort().argsort()[-1]/(L-1) \n",
    "        plt_y.scatter(last,last_pvalue,color='red',edgecolor='black',s=30,zorder=1000)\n",
    "        plt_y.plot([0,last],[last_pvalue,last_pvalue],color='red',zorder=200)\n",
    "        plt_y.plot([last,last],[0.0,last_pvalue],color='red',zorder=200)\n",
    "    else:\n",
    "        plt_y.scatter(last,0.0,color='red',edgecolor='black',s=30,zorder=1000)\n",
    "        plt_y.scatter(last,1.0,color='red',edgecolor='black',s=30,zorder=1000)\n",
    "        \n",
    "    #plt_y.plot(np.sort(resource),np.linspace(0.0,1.0,L),color='black',lw=1,ls='--',zorder=200)\n",
    "    plt_y.scatter(np.sort(resource),np.linspace(0.0,1.0,L),color='yellow',edgecolor='black',s=20,zorder=201)\n",
    "\n",
    "    plt_y.set_ylim([0,1]); plt_y.set_xlim(0.1,1000); plt_y.set_ylabel('Cumulative Probability') \n",
    "    plt_y.set_xlabel('Lithium (tonnes)'); plt_y.set_title('Recoverable Lithium Resource')\n",
    "    \n",
    "    method = 'median_unbiased'\n",
    "    \n",
    "    P10 = np.interp(0.10,np.linspace(0.0,1.0,len(resource)),np.sort(resource),left=None, right=None, period=None)\n",
    "    \n",
    "    if L > 1:       \n",
    "        cdfy = np.interp(np.linspace(0.0,1.0,100),np.linspace(0.0,1.0,len(resource)),np.sort(resource),left=None, right=None, period=None)\n",
    "        plt_y.plot(cdfy,np.linspace(0.0,1.0,100),color='black',lw=1,ls='--',zorder=200)\n",
    "    \n",
    "#     cdfy = np.interp(np.linspace(0.0,1.0,100),np.linspace(0.0,1.0,len(resource)),np.sort(resource),left=None, right=None, period=None)\n",
    "#     plt.plot(cdfy,np.linspace(0.0,1.0,100),color='black',lw=1,ls='--',zorder=200)\n",
    "    \n",
    "    plt_y.plot([0,P10],[0.1,0.1],color='black')\n",
    "    plt_y.plot([P10,P10],[0.0,0.1],color='black')\n",
    "    \n",
    "    plt_y.plot([0,np.percentile(resource,50)],[0.5,0.5],color='black')\n",
    "    plt_y.plot([np.percentile(resource,50),np.percentile(resource,50)],[0.0,0.5],color='black')\n",
    "    \n",
    "    plt_y.plot([0,np.percentile(resource,90)],[0.9,0.9],color='black')\n",
    "    plt_y.plot([np.percentile(resource,90),np.percentile(resource,90)],[0.0,0.9],color='black')\n",
    "\n",
    "    plt_y.annotate(r'$Li^{\\ell}$ = ' + str(np.round(last,1)),xy=[last,0.015],color='red',rotation=270,zorder=200)\n",
    "    plt_y.annotate('0.10',xy=[1.1,0.11]); plt_y.annotate('0.50',xy=[1.1,0.51]); plt_y.annotate('0.90',xy=[1.1,0.91])\n",
    "    if L > 1:\n",
    "        plt_y.annotate(r'$F_{\\ell}^{-1}$ = P' + str(int(last_pvalue*100)),xy=[1.1,last_pvalue+0.01],color='red',rotation=0,zorder=200)\n",
    "        plt_y.annotate(r'$P10_{Li}$ = ' + str(np.round(np.percentile(resource,10),1)),xy=[np.percentile(resource,10),0.015],rotation=270,zorder=100)\n",
    "        plt_y.annotate(r'$P50_{Li}$ = ' + str(np.round(np.percentile(resource,50),1)),xy=[np.percentile(resource,50),0.015],rotation=270,zorder=100)\n",
    "        plt_y.annotate(r'$P90_{Li}$ = ' + str(np.round(np.percentile(resource,90),1)),xy=[np.percentile(resource,90),0.015],rotation=270,zorder=100)\n",
    "    else:\n",
    "        plt_y.plot([last,last],[0.0,1.0],color='black',lw=1,ls='--',zorder=200)        \n",
    "    \n",
    "    plt_y.set_xscale('log')\n",
    "    plt_y.grid(which='both',axis='both',visible=True)\n",
    "    plt_y.grid(True, which='major',linewidth = 1.0); plt_y.grid(True, which='minor',linewidth = 0.2)\n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=1.5, top=1.7, wspace=0.3, hspace=0.5); plt.show()\n",
    "    \n",
    "interactive_plot4 = widgets.interactive_output(f_make4, {'L':L,'A_mean':A_mean,'A_stdev':A_stdev,'th_mean':th_mean,'th_stdev':th_stdev,'li_mean':li_mean,'li_stdev':li_stdev})\n",
    "interactive_plot4.clear_output(wait = True)                # reduce flickering by delaying plot updating    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Monte Carlo Simulation for Lithiums Resource Uncertainty Demonstration\n",
    "\n",
    "#### Michael Pyrcz, Professor, The 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, * **area** ($A$) - of the salar brine in units, $m^2$\n",
    "   \n",
    "* **average thickness** ($\\overline{th}$) - the average thickness over the area of the salar brine of the bound aquifer\n",
    "\n",
    "* **average lithium concentation** ($\\overline{[Li]}$) - amount of disssolved lithium, in $\\frac{mg}{L}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d984ca85c1c44208976710e555862138",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                          Monte Carlo Simulation Lithium Resources …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "eed0a2d221db471f9f7a9fadf5e84a0a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui4, interactive_plot4)                            # 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 one could integrate relationships between the variables (we assumed independence here).\n",
    "\n",
    "Note, I have other demonstrations with McMC methods, also!\n",
    "  \n",
    "I hope this was helpful,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "#### The Author:\n",
    "\n",
    "### Michael Pyrcz, Professor, The 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. Professor, Cockrell School of Engineering and 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": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}