{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "

\n", " \n", "\n", "

\n", "\n", "### Monte Carlo Methods 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 three dashboards to demonstrate Monte Carlo methods for data science. This should help you get started with:\n", "\n", "* solving otherwise intractable or computationally expensive problems, numerical integration or optimization by random sampling (drawing realizations).\n", "\n", "whereas Monte Carlo Simulation (MCS) uses random sampling to simulate systems (simulate reality), Monte Carlo methods use random sampling to solve mathematical or statistical problems. \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 [10 Data Science Interactive: Monte Carlo Methods](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 demonstrations are available on my GitHub repository [GeostatsGuy's Python Numerical Demos](https://github.com/GeostatsGuy/PythonNumericalDemos).\n", "\n", "* Join me for walk-through of my Monte Carlo Simulation dashboard [09 Data Science Interactive: Monte Carlo Simulation](https://www.youtube.com/watch?v=5y6eWl4YPSk). \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 Methods\n", "\n", "Definition: the use of random sampling to solve a mathematical or statistical problem.\n", "\n", "Procedure: \n", "\n", "1. Define the interval of possible inputs to the system, $x_{1_{min}} < X_1 < x_{1_{max}}$, $x_{2_{min}} < X_2 < x_{2_{max}}$, $\\ldots$, $x_{m_{min}} < X_m < x_{m_{max}}$ \n", "2. Monte Carlo simulate, draw an independent realization, from this interval, $x_1^\\ell$, $x_2^\\ell$, $\\ldots$, $x_m^\\ell$. \n", "3. Perform the mathematical or statistical calculation, on this realization, $y^{\\ell} = f\\left(x_1^\\ell, x_2^\\ell, \\ldots, x_m^\\ell \\right)$\n", "4. Go to step 2, repeat until enough realizations for an accurate solution\n", "5. Summarize $f(y)$ over the realizations, $y^{\\ell}, \\ell = 1,\\ldots,L$.\n", "\n", "#### Monte Carlo Methods vs. Monte Carlo Simulation (MCS)\n", "\n", "Monte Carlo methods use random sampling to solve a mathematical or statistical problem, while MCS simulates a system to solve a statistical uncertainty problem by summarizing the outcomes over many realizations. Here's some MCS examples to compare with the dashboards that I made below,\n", "\n", "1. MCS from uncertainty distributions, random variables, for average porosity, thickness and area, to calculate and summarize over many realizations of volume of available ground water.\n", "2. Assign a strategy, and have the computer play many hands of solitaire and summaries the probability of winning. \n", "\n", "#### General 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", "#### How many realizations, $L$?\n", "\n", "The answer is enough realizations! \n", "\n", "* If the Monte Carlo method computational cost is low then **many** is the right answer for greater accuracy. \n", "\n", "* If too few realizations are calculated then the accuracy of the output is limited (see the 'Law of Small Numbers'). \n", "\n", "#### Limitations\n", "\n", "The Monte Carlo methods limitations:\n", "\n", "1. **approximative** - an approximation of the solution, but more accurate as $L \\rightarrow \\infty$.\n", "2. **computational complexity** - if the transfer function is expensive then given the large required number of realizations, may be very high computational complexity \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", "import math\n", "from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks\n", "from scipy import stats # summary statistics\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 functions to:\n", "\n", "1. add major and minor gridlines\n", "2. calculate a isotropic variogram between 2 points\n", "3. calculate the azimuth between 2 points to rotate a label" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "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 isotropic_spherical_variogram(x1,y1,x2,y2,r):\n", " dist = np.sqrt((x1-x2)**2 + (y1-y2)**2)\n", " ones = np.ones(len(x1))\n", " gamma = np.where(dist >= r,ones,3/2*(dist/r)-0.5*(dist/r)**3)\n", " return gamma\n", "\n", "def rotation(tail,head):\n", " return np.arctan2(tail,head)*180.0/math.pi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Monte Carlo Method to Calculate the Mean Demonstration\n", "\n", "Now let's set up our first dashboard. This is a simple inferential example, calculate a mean by random sampling from a distribution. \n", "\n", "#### Monte Carlo Method to Calculate a Mean\n", "\n", "1. randomly sample a p-value\n", " * uniform random cumulative probability value\n", "2. sample from the distribution by applying the p-value to the inverse CDF \n", " * $x^\\ell = F^{-1}_X(p^\\ell)$\n", "3. calculate the average\n", "4. go to 1., continue until enough samples for the required accuracy" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "l = widgets.Text(value=' Monte Carlo Method to Estimate Mean 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", " truth_mean = np.average([min1,max1])\n", " counts = np.arange(1,len(sample)+1,1)\n", "\n", " L_all = 40\n", " np.random.seed(seed = 73073)\n", " sample_all, pvals_all, cdfx_all, cdf_all = make_sample_cdist(dist1,min1,max1,L_all)\n", " prop_all = np.cumsum(sample_all, axis=0, dtype=None, out=None)/np.arange(1,len(sample_all)+1,1)\n", " counts_all = np.arange(1,len(sample_all)+1,1)\n", " \n", " ax1 = plt.subplot(121)\n", " plt.plot(cdfx,cdf,'--',color='black',linewidth = 3,zorder=10); add_grid(ax1)\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", " ax2 = 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.vlines(truth_mean,0,10,colors='black',zorder=1000)\n", " plt.annotate('Population Mean',[truth_mean+1,7],rotation=270.0)\n", " \n", " plt.plot(prop_all,counts_all/L_all*10,color='gray',lw=2,ls='--',zorder=999)\n", " plt.plot(prop_all[:l+1],counts_all[:l+1]/L_all*10,color='blue',lw=2,zorder=1000)\n", " plt.scatter(np.average(sample),L/L_all*10,c='blue',edgecolor='black',zorder=1000,label='Sample Mean')\n", " plt.annotate('Sample Mean',[np.average(sample)+1,L/L_all*10+0.2],color='blue',zorder=1000)\n", " \n", " plt.xlabel(r\"$x_{\\alpha}, \\alpha = 1,\\ldots,L$\"); \n", " plt.title(r\"Monte Carlo Arithmetic Average\"); plt.ylabel('Frequency')\n", " plt.xlim([0,100]); plt.ylim([0,10]); add_grid(ax2)\n", " plt.legend(loc='upper left')\n", " \n", " ax2 = plt.gca().twinx()\n", " ax2.set_ylabel('Number of Realizations',color='blue'); ax2.set_ylim([1,L_all])\n", " ax2.spines['right'].set_color('blue'); ax2.tick_params(axis='y', colors='blue')\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": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "ecbc0724abc648bbbcfa979829871983", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(Text(value=' Monte Carlo Method to Estimate Mean Demonstration,…" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "26fdb0a1f38e4051bf2d154f9feb78cd", "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 Method to Calculate a Geometric Ratio of the Area of Circle Inscribed in a Square Demonstration\n", "\n", "Now let's set up our second dashboard. This is an example where we solve a geometric problem, the ratio of circle inscribed in a square by random sampling. \n", "\n", "* All we need to calculate is whether a point is in the circle, easy to do given the circle centroid and radius. \n", "\n", "#### Monte Carlo Method to Calculate Geometric Ratio of the Area of Circle Inscribed in a Square\n", "\n", "1. randomly sample a point within the square $v_1$\n", " * sample independent uniform random values over x and y.\n", "2. check if the point is within the circle with a distance to center of the circle\n", " * inside if less than the circle's radius\n", "3. calculate the proportion of points within the circule\n", "4. go to 1., continue until enough samples for the required accuracy\n", "\n", "From geometry we can calculate the exact solution:\n", "\n", "\\begin{equation}\n", "\\frac{\\pi \\cdot r^2}{\\left( 2 \\cdot r \\right)^2} = \\frac{\\pi \\cdot r^2}{4 \\cdot r^2} = \\frac{\\pi}{4} \n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "title = widgets.Text(value=' Monte Carlo Method to Estimate Ratio Circle Incribed in a Square Demonstration, Michael Pyrcz, Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n", "\n", "l = widgets.FloatSlider(min=0, max = 5, step = 0.1, value = 0.3, description = '$log(l)$:',orientation='horizontal',layout=Layout(width='750px', height='50px'),continuous_update=False)\n", "l.style.handle_color = 'gray'\n", "\n", "ui = widgets.HBox([l],kwargs = {'justify_content':'center'})\n", "ui2 = widgets.VBox([title,ui],)\n", "\n", "def f_make2(l): \n", " L = 100000; radius = 50; cx = 50; cy = 50; seed = 73071\n", " \n", " i = int(10**l)\n", " np.random.seed(seed=seed)\n", " pts = np.random.rand(L,2)*100\n", " pts_l = pts[:i]\n", " dist = (pts[:,0]-50)**2 + (pts[:,1]-50)**2\n", " dist_l = (pts_l[:,0]-50)**2 + (pts_l[:,1]-50)**2\n", " inside = dist < (radius * radius)\n", " inside_l = dist_l < (radius * radius)\n", " \n", " counts = np.arange(1,len(inside)+1,1)\n", " prop = np.cumsum(inside, axis=0, dtype=None, out=None)/np.arange(1,len(inside)+1,1)\n", " truth = math.pi/4\n", " \n", " ax1 = plt.subplot(121)\n", " \n", " plt.scatter(pts_l[inside_l,0],pts_l[inside_l,1],color='red',edgecolor='black',zorder=10)\n", " plt.scatter(pts_l[inside_l==False,0],pts_l[inside_l==False,1],color='gray',edgecolor='black',zorder=10)\n", " plt.scatter(50,50,marker='x',c='black',s=80,zorder=1)\n", " \n", " circle1 = plt.Circle((50, 50),radius,fill=False,edgecolor='black',zorder=1)\n", " plt.plot([52,50+radius-1],[50,50],c='black',lw=1)\n", " plt.gca().add_patch(circle1)\n", " \n", " plt.gca().add_patch(plt.Rectangle((0, 0), 100, 100, fill=False,edgecolor='black',lw=2))\n", " \n", " head = plt.Polygon([[cx+radius-1,cy-1],[cx+radius-1,cy+1],[cx+radius,cy],[cx+radius-1,cy-1]], color='black',alpha = 1.0,zorder=1)\n", " plt.gca().add_patch(head)\n", " plt.xlim([-10,110]); plt.ylim([-10,110]); add_grid(ax1); plt.xlabel(\"X (m)\"); plt.ylabel(\"Y (m)\") \n", " plt.title('Ratio of Circle Inscribed Within Square')\n", " \n", " ax2 = plt.subplot(122)\n", " plt.plot(counts,prop,c='gray',zorder=1)\n", " plt.plot([0,L],[truth,truth],c='red',zorder=3); plt.annotate('Exact Solution',[2,truth+0.01],c='red')\n", " plt.scatter(counts[:i-1],prop[:i-1],c='gray',edgecolor='black',s=30,zorder=9);\n", " plt.scatter(counts[i-1],prop[i-1],c='red',edgecolor='black',s=30,zorder=10); add_grid(ax2)\n", " plt.ylim([0.0,1.3]); plt.xlim([1,100000]); plt.xscale(\"log\"); plt.xlabel(\"Number of Samples\"); plt.ylabel(\"Cumulative Proportion\"); plt.title('Monte Carlo Geometric Ratio')\n", " \n", " plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.3, wspace=0.3, hspace=0.2); plt.show()\n", " \n", "interactive_plot2 = widgets.interactive_output(f_make2, {'l':l})\n", "interactive_plot2.clear_output(wait = True) # reduce flickering by delaying plot updating " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "8a91197dbd02474cb7273fd500eebe88", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(Text(value=' Monte Carlo Method to Estimate Ratio Circle Incribed in a Squa…" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "50789c2a780b41f0afa468fb755d4337", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Output()" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "display(ui2, interactive_plot2) # display the interactive plot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Monte Carlo Method to Calculate Volume Integrated Variograms (Gamma Bar) Demonstration\n", "\n", "Now let's set up our third dashboard. This is an example of an expensive integration that we can approximate with a limited number of samples. \n", "\n", "\\begin{equation}\n", "\\overline{\\gamma}(v_1,v_2) = \\frac{1}{v_1 \\cdot v_2} \\int_{v_1} \\int_{v_2} \\gamma\\left( \\bf{u}_{v_1},\\bf{u}_{v_2} \\right) dv_1 dv_2\n", "\\end{equation}\n", "\n", "For a 2D problem this is a fourfold (quadruple) integral and for a 3D problem this is a sixfold (sextuple) integral, see below: \n", "\n", "In 2D:\n", "\\begin{equation}\n", "\\overline{\\gamma}(v_1,v_2) = \\frac{1}{v_1 \\cdot v_2} \\int_{{v_1}_X} \\int_{{v_1}_Y} \\int_{{v_2}_X} \\int_{{v_2}_Y} \\gamma\\left( \\bf{u}_{v_1},\\bf{u}_{v_2} \\right) dv_{1_X} dv_{1_Y} dv_{2_X} dv_{2_Y}\n", "\\end{equation}\n", "\n", "and in 3D:\n", "\n", "\\begin{equation}\n", "\\overline{\\gamma}(v_1,v_2) = \\frac{1}{v_1 \\cdot v_2} \\int_{{v_1}_X} \\int_{{v_1}_Y} \\int_{{v_1}_Z} \\int_{{v_2}_X} \\int_{{v_2}_Y} \\int_{{v_2}_Z} \\gamma\\left( \\bf{u}_{v_1},\\bf{u}_{v_2} \\right) dv_{1_X} dv_{1_Y} dv_{1_Z} dv_{2_X} dv_{2_Y} dv_{2_Z} \n", "\\end{equation}\n", "\n", "given the semivariogram calculation computational cost, combinatorial of variogram values to average to calculate 1 gamma bar value, and combined with the large number of gamm bar values often required in geostatistical spatial estimation and simulation, this is a too high computational complexity and we need a speed up!\n", "\n", "#### Monte Carlo Method to Calculate Gamma Bar Values\n", "\n", "1. randomly sample a lag vector between volumes $v_1$ and $v_2$\n", " * uniform random value over x and y for each the lag vector tail (in $v_1$) and the lag vector head (in $v_2$).\n", "2. calculate the variogram value over the lag vector\n", " * in this example, I assume a isotropic spherical variogram model\n", "3. calculate the average variogram over all random lag vectors\n", "4. go to 1., continue until enough samples for the required accuracy\n", "\n", "Note, the 'Exact Solution' for the gamma bar was calculated with 100k samples and hard coded in the dashbaord below." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "title = widgets.Text(value=' Monte Carlo Method to Estimate Volume Integrated Variograms, Gamma Bars, Demonstration, Michael Pyrcz, Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n", "\n", "L_dia = widgets.FloatSlider(min=0.25, max = 2.0, step = 0.25, value = 0.3, description = '$log(l)$:',orientation='horizontal',layout=Layout(width='750px', height='50px'),continuous_update=False)\n", "L_dia.style.handle_color = 'gray'\n", "\n", "ui = widgets.HBox([L_dia],kwargs = {'justify_content':'center'})\n", "ui3 = widgets.VBox([title,ui],)\n", "\n", "def f_make3(L_dia): \n", " L = int(10**L_dia)\n", " seed = 73073; r = 200.0; L_all = 100\n", " \n", " gamma_truth = 0.81038 # calculated with 100k samples\n", " np.random.seed(seed=seed)\n", " pts_v1_all = np.random.rand(L_all,2)\n", " pts_v1_all[:,0] = pts_v1_all[:,0]*50.0+25.0; pts_v1_all[:,1] = pts_v1_all[:,1]*50.0+35.0; \n", " pts_v2_all = np.random.rand(L_all,2);\n", " pts_v2_all[:,0] = pts_v2_all[:,0]*100.0+125; pts_v2_all[:,1] = pts_v2_all[:,1]*100.0+10.0; \n", " gamma_all = isotropic_spherical_variogram(pts_v1_all[:,0],pts_v1_all[:,1],pts_v2_all[:,0],pts_v2_all[:,1],r)\n", " \n", " pts_v1 = pts_v1_all[:L,:]; pts_v2 = pts_v2_all[:L,:]; gamma = gamma_all[:L]\n", " \n", " counts_all = np.arange(1,L_all+1,1)\n", " prop_all = np.cumsum(gamma_all, axis=0, dtype=None, out=None)/np.arange(1,L_all+1,1)\n", " \n", " f, (ax1, ax2) = plt.subplots(1, 2, gridspec_kw={'width_ratios': [3, 1]})\n", " \n", " for l in range(0,L):\n", " alpha = max(0.02,(11 - (L-l) )/10) \n", " if l == L-1:\n", " rotate = rotation(pts_v2[l,1]-pts_v1[l,1],pts_v2[l,0]-pts_v1[l,0])\n", " ax1.scatter(pts_v1[l,0],pts_v1[l,1],color='red',edgecolor='black',alpha = alpha,s=20,zorder=150)\n", " ax1.scatter(pts_v2[l,0],pts_v2[l,1],color='red',edgecolor='black',alpha = alpha,s=20,zorder=150)\n", " ax1.plot([pts_v1[l,0],pts_v2[l,0]],[pts_v1[l,1],pts_v2[l,1]],color='red',alpha = alpha,zorder=130)\n", " ax1.annotate(r'$\\gamma(v_1^{\\ell},v_2^{\\ell})$ = ' + str(np.round(gamma[l],2)),\n", " [(pts_v1[l,0]+pts_v2[l,0])*0.5,(pts_v1[l,1]+pts_v2[l,1])*0.5+2],\n", " rotation = rotate,color='red',ha='center',zorder=100)\n", " ax1.annotate(r'$\\overline{\\gamma}(v_1,v_2)$ = ' + str(np.round(prop_all[l],2)),\n", " [5,5],color='red',zorder=100)\n", " else:\n", " ax1.scatter(pts_v1[l,0],pts_v1[l,1],color='gray',edgecolor='black',alpha = alpha,s=20,zorder=50)\n", " ax1.scatter(pts_v2[l,0],pts_v2[l,1],color='gray',edgecolor='black',alpha = alpha,s=20,zorder=50)\n", " ax1.plot([pts_v1[l,0],pts_v2[l,0]],[pts_v1[l,1],pts_v2[l,1]],color='black',alpha = alpha,zorder=30)\n", " \n", " ax1.add_patch(plt.Rectangle((25,35), 50, 50, fill=False,edgecolor='black',lw=2))\n", " ax1.annotate(r'$v_1$',[26,86])\n", " ax1.add_patch(plt.Rectangle((125,10), 100, 100, fill=False,edgecolor='black',lw=2))\n", " ax1.annotate(r'$v_2$',[126,111])\n", " add_grid(ax1)\n", " ax1.set_xlim([0,240]); ax1.set_ylim([0,120]); ax1.set_xlabel('X (m)'); ax1.set_ylabel('Y (m)'); ax1.set_title('Volume Integrated Variograms (Gammabar)')\n", " \n", " counts = np.arange(1,L+1,1)\n", " prop = np.cumsum(gamma, axis=0, dtype=None, out=None)/np.arange(1,L+1,1)\n", " \n", " ax2.plot(prop_all,counts_all,color='gray',zorder=10); ax2.scatter(prop,counts,color='gray',edgecolor='black',s=20,zorder=100)\n", " ax2.scatter(prop[l],counts[l],color='red',edgecolor='black',s=50,zorder=120)\n", " ax2.plot([gamma_truth,gamma_truth],[1,L_all],color='red',zorder=20); \n", " plt.annotate('Exact Solution',[gamma_truth,3],c='red',rotation=270)\n", " ax2.set_xlim([0.5,1.0]); ax2.set_ylim([1,L_all]); ax2.set_xlabel(r'$\\overline{\\gamma}(v_1,v_2)$'); ax2.set_ylabel('Number of Samples'); ax2.set_title('Monte Carlo Gammabar')\n", " add_grid(ax2); ax2.set_yscale('log')\n", " \n", " plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.3, hspace=0.2); plt.show() \n", " \n", "interactive_plot3 = widgets.interactive_output(f_make3, {'L_dia':L_dia})\n", "interactive_plot3.clear_output(wait = True) # reduce flickering by delaying plot updating " ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "scrolled": false }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "a4aa7e0f2b804bf5bcda512607237219", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(Text(value=' Monte Carlo Method to Estimate Volume Integrated Variograms, Gamma …" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "eee0269a75b0439f8979ac4b6a7c91ee", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Output()" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "display(ui3, interactive_plot3) # display the interactive plot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Comments\n", "\n", "This was an interactive demonstration of Monte Carlo methods to solve statistical and mathematical problems. I have other lectures and inteactive dashboards on Monte Carlo Simulation and I wanted to expand our coverage with general Monte Carlo methods here. I also have lectures on Markov chain Monte Carlo (McMC) methods!\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": [] } ], "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 }