04_Stats_and_Montecarlo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Created under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2020 R.C. Cooper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Statistics and Monte-Carlo Models\n",
    "\n",
    "Monte Carlo models use random numbers to either understand statistics or generate a solution [1]. \n",
    "The main element in a Monte Carlo model is the use of random numbers. Monte Carlo methods are very useful if you can easily execute a function lots of time or even in parallel. \n",
    "\n",
    "We can generate random numbers in many ways, but most programming languages have 'pseudo'-random number generators. \n",
    "\n",
    "In Python, we use the numpy library as such\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "x=np.random.rand(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "#Import rcParams to set font styles\n",
    "from matplotlib import rcParams\n",
    "\n",
    "#Set font style and size \n",
    "rcParams['font.family'] = 'sans'\n",
    "rcParams['font.size'] = 16\n",
    "rcParams['lines.linewidth'] = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([3., 1., 0., 0., 1., 0., 0., 2., 1., 2.]),\n",
       " array([0.34837708, 0.41148652, 0.47459596, 0.53770539, 0.60081483,\n",
       "        0.66392427, 0.72703371, 0.79014314, 0.85325258, 0.91636202,\n",
       "        0.97947146]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD9CAYAAABDaefJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAMsUlEQVR4nO3cf4wnd13H8eerPeoKJfGaFhLBstRD6ZESgxdsgtofIKUl4B/U/oGiMdZqExNiUQOiFUGlVdFElNKLmBKLCKWixESltV4B0ytuTVqo2FqEQk2abrmW8qNXKHz8Y+Zw/bLXne/O97vf2/c+H8k3393Z7+x8vjN3z52dmZ201pAk1XDcogcgSZodoy5JhRh1SSrEqEtSIUZdkgrZNe8FnHzyyW15eXnei5GkUm677bYHW2unTDvf3KO+vLzMysrKvBcjSaUkuXcz83n4RZIKMeqSVIhRl6RCjLokFTIo6knOS3JTkvuTPJbkviTvT7J33gOUJA039OqXk4DbgHcAq8CpwOuBg0nOaK1t6iytJGm2BkW9tfZe4L1rpyX5OPCfwIXA22Y/NEnStMYcU/9C//z1WQxEkjTeVFFPcnySE5I8B7gauB/467mMTJI0tWn31G8FHgPuBp4PnNtae2DyRUkuSbKSZGV1dXXUAA8fHjX7tluuJI0xbdRfA5wJvBp4BLghyfLki1pr+1tr+1pr+045ZepbF/w/S0uQbP1jaWnUsCVpIaaKemvtU621W/sTpy8GTqS7CkaSdAzY9InS1trDwD3AntkNR5I0xqajnuTpwHOBT89uOJKkMQZdp57kg8C/A3fQHUv/PuCXgcfxGnVJOmYM/YvSg8BFwOuAE4DPAweAt7bWPjuXkUmSpjb0L0qvBK6c81gkSSN5l0ZJKsSoS1IhRl2SCjHqklSIUZekQoy6JBVi1CWpEKMuSYUYdUkqxKhLUiFGXZIKMeqSVIhRl6RCjLokFWLUJakQoy5JhRh1SSrEqEtSIUZdkgox6pJUiFGXpEKMuiQVYtQlqRCjLkmFGHVJKsSoS1IhRl2SCjHqklSIUZekQoy6JBVi1CWpEKMuSYUYdUkqxKhLUiFGXZIKMeqSVIhRl6RCjLokFWLUJakQoy5JhRh1SSrEqEtSIUZdkgox6pJUiFGXpEKMuiQVYtQlqRCjLkmFGHVJKsSoS1IhRl2SCjHqklSIUZekQoy6JBVi1CWpEKMuSYVsGPUkFya5Psm9SR5NcleStyZ56lYMUJI03JA99V8BvgH8OvAy4CrgUuCGJO7pS9IxZNeA17yitba65vObkxwC3g2cDdw0j4FJkqa34Z72RNCP+Lf++RmzHY4kaYzNHj45q3/+1KwGIkkab+qoJ3kG8GbgxtbaylFec0mSlSQrq6vr7ehLWpTDh3fWche57EUsd8gx9W9JciLwd8DjwM8e7XWttf3AfoB9+/a1MQOUNFtLS5Bs/XLbAkuwk97z4KgnWQI+BJwGnNVau29uo5IkbcqgqCd5EnA98ELgJa21T8x1VJKkTdkw6v216O8BXgy8vLV2cO6jkiRtypA99T8DfgL4XeArSc5c87X7PAwjSceOIVe/nN8/vxG4ZeJx8ZzGJUnahA331Ftry1swDknSDHjvFkkqxKhLUiFGXZIKMeqSVIhRl6RCjLokFWLUJakQoy5JhRh1SSrEqEtSIUZdkgox6pJUiFGXpEKMuiQVYtQlqRCjLkmFGHVJKsSoS1IhRl2SCjHqklSIUZekQoy6JBVi1CWpEKMuSYUYdUkqxKhLUiFGXZIKMeqSVIhRl6RCjLokFWLUJakQoy5JhRh1SSrEqEtSIUZdkgox6pJUiFGXpEKMuiQVYtQlqRCjLkmFGHVJKsSoS1IhRl2SCjHqklSIUZekQoy6JBVi1CWpEKMuSYUYdUkqxKhLUiFGXZIKMeqSVIhRl6RCjLokFWLUJakQoy5JhQyKepJnJnl7kluSfDVJS7I836FJkqY1dE99D3AR8BDw0fkNR5I0xtCof6S19vTW2gXAdfMckCRp8wZFvbX2zXkPRJI0nidKJamQuUQ9ySVJVpKsrK6uzmMRc3f48M5c9iK4rncG1/XW2DWPb9pa2w/sB9i3b1+bxzLmbWkJksUsu23LNbZ5ruudwe28NTz8IkmFGHVJKsSoS1Ihg4+pJ7mw//AH++fzk6wCq621m2c+MknS1KY5UTr5R0fv6J9vBs6eyWgkSaMMjnprbUHnrSVJQ3lMXZIKMeqSVIhRl6RCjLokFWLUJakQoy5JhRh1SSrEqEtSIUZdkgox6pJUiFGXpEKMuiQVYtQlqRCjLkmFGHVJKsSoS1IhRl2SCjHqklSIUZekQoy6JBVi1CWpEKMuSYUYdUkqxKhLUiFGXZIKMeqSVIhRl6RCjLokFWLUJakQoy5JhRh1SSrEqEtSIUZdkgox6pJUiFGXpEKMuiQVYtQlqRCjLkmFGHVJKsSoS1IhRl2SCjHqklSIUZekQoy6JBVi1CWpEKMuSYUYdUkqxKhLUiFGXZIKMeqSVIhRl6RCjLokFWLUJakQoy5JhRh1SSrEqEtSIYOinuR7knwgyReTPJLkb5KcOu/BSZKms2HUkzwZuAl4LvAzwGuA5wD/kuQp8x2eJGkauwa85ueB04Dvb63dA5DkDuC/gF8A/mh+w5MkTWPI4ZdXAgePBB2gtfYZ4F+BH5/XwCRJ0xsS9ecBn1xn+p3A3tkOR5I0xpDDLycBD60z/RCwe70ZklwCXNJ/+uUkd21ueFvmZODBtROSBY1kwctekG9b/1tlB67r9WzJ+t+J/6emWO562+BZm1nmkKgDtHWmHXW4rbX9wP7NDGgRkqy01vYtehw7let/sVz/izfLbTDk8MtDdHvrk3az/h68JGlBhkT9Trrj6pP2Av8x2+FIksYYEvUPAWcmOe3IhCTLwIv6r1WwbQ4VFeX6XyzX/+LNbBuktfUOl695QfcHRrcDjwK/QXd8/S3AU4Hnt9a+PKvBSJLG2XBPvbX2FeBc4G7gL4H3AJ8BzjXoknRs2XBPXZK0fZS9S+OsbkKW5A1JWpKPzWOcVY1d/0lOT3JdkgeTPJrkriSvneeYqxmzDZKcmuTdST6X5KtJ7k7yO97vabgkz0zy9iS39Ouw9ecjh8x7XN+ezyY5nOT2JK8aMm/JqM/qJmT9yeE3Ag/MY5xVjV3/SfYBtwLfAVwMXAC8DTh+XmOuZsw26L9+I/CjwG8CLwf+HHgd8BdzHHY1e4CL6C79/uiU874FeBPwp8D5wEHguiQXbDhna63cA3gt8A1gz5ppzwYeBy6b4vv8E3A1cAD42KLf13Z5jFn/dDsadwIfXPT72M6PkdvgpXQXRLx0YvoV/fxPXvT72w4P4Lg1H1/cr9PlAfM9DXgM+O2J6f8M3LHR/CX31JnBTciSvBp4AfCGuYywtjHr/2y6v4Hw7p/jjNkGJ/TPj0xMf5juh643VxigtfbNTc56Ht02uHZi+rXAGUme/UQzV436qJuQJdkN/DHwa621QzMe204wZv3/cP+8lORgkq8neSDJnyT5zpmOsrYx2+BGultrX5lkb5ITk5xLt/f/ztZdEaf5eR7dnvo9E9Pv7J+fcPtVjfrUNyGb8Ad0l3BeM8Mx7SRj1v9398/vAz4M/Bjw+3S/vv7VrAa4A2x6G7TWDtP9cD1yKOxLdL/6/z3wS7MdptZxEvBw64+5rHFozdePaugNvbajqW5C9q0XJD8C/DTwgnVWqobb1Prn/3Y0rm2tXd5/fCDJ8cAVSfa21rw9xTCb/T+wRPdD9Wl0J1g/B7wQuJzumPqlMxyjvl3Y/P+fslEfcxOyq4F3Afcl+a5+2i7g+P7zR1trj81spDWNWf9f6J9vmJj+YboTdT+A9xwaYsw2+Dm6cxt7Wmuf7qd9JMkXgf1J3tlau31mI9WkQ8DuJJnYsdy95utHVfXwy5ibkJ0O/CLdP/wjjxcBZ/Yfu5eysTHr/8hxw8k9lSN7KZs9+bTTjNkGZwAPrQn6ER/vn08fOTY9sTvpLuf93onpR46lP+H2qxr1MTchO2edx+10J53OAT4w++GWM2b9/wPdSaKXTUw/r39emc0QyxuzDe6n21PcMzH9h/rn/5nRGLW+fwS+BvzkxPSfAj7ZX8V0dIu+lnNO14c+he7M8SfoLt96JV2Y/xs4cc3rnkV3jPDyDb7fAbxOfcvWP/Bb/fTfA14CvJ7uhnLXLPq9bZfHmG0ALNNdzng33R8unQP8aj9thTXXX/vYcDtc2D+uovvt89L+87PWvOZx4F0T810BHAYuozsUdhXdb6mv2HCZi37Tc1yZpwLX9/8QvwT8LRMX/vf/eBvwpg2+l1HfwvVPd6jlsj5KXwPuBd4MPGnR72s7PUZug73A+4HP9z9Q7wb+ENi96Pe1nR79ul3vcWDiNddMzHc83V1x76X7zfUO4MIhy/SGXpJUSNVj6pK0Ixl1SSrEqEtSIUZdkgox6pJUiFGXpEKMuiQVYtQlqZD/BbUtv75EOncoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(x,bins=10, color='b', histtype='bar', edgecolor='w')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `numpy.random.rand(10)` function generates 10 random numbers between 0 and 1. The pyplot function `hist` then displays a histogram of these randomly generated numbers. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Try generating more random numbers and plotting histograms of the results i.e. increase `10` to larger values. \n",
    "\n",
    "What should the histogram of `x` look like if Python is generating truly random numbers?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Examples of Monte Carlo models:\n",
    "\n",
    "Monte Carlo models have a wide array of applications. We are going to use Monte Carlo models in later modules to explore how uncertainty in measurements can be incorporated into computational models. The three main applications for Monte Carlo models are used in three main classes: optimization, numerical integration, and generating population distributions [1]. \n",
    "\n",
    "Here is a brief list of Monte Carlo model use cases in real-world applications:\n",
    "\n",
    "- [Eigenvlaues in supercritical systems](https://link.springer.com/chapter/10.1007%2FBFb0049064)\n",
    "- [average time between failures for reliability](http://www.egr.msu.edu/~mitraj/research/pubs/proc/singh-mitra_em_stdby_ias95.pdf)\n",
    "- disordered materials (physics)\n",
    "- [Calculation of the energy output of a wind farm](http://www.mdpi.com/1996-1073/9/4/286/pdf)\n",
    "- [US Coast Guard rescue missions](https://en.wikipedia.org/wiki/Search_and_Rescue_Optimal_Planning_System)\n",
    "- [Radiation shielding](http://www.sciencedirect.com/science/article/pii/S0920379612000580)\n",
    "- [Predict number of asteroids that hit body of water](https://cneos.jpl.nasa.gov/sentry/intro.html)\n",
    "- [Financial modeling](https://en.wikipedia.org/wiki/Monte_Carlo_methods_in_finance)\n",
    "\n",
    "We will explore Monte Carlo modeling through the use of three examples:\n",
    "\n",
    "1. Calculate the value of $\\pi$\n",
    "\n",
    "2. Calculate the integral of a function\n",
    "\n",
    "3. Propagate uncertainty in manufacturing into uncertainty in failure load"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "\n",
    "# Example 1: Calculate $\\pi$ with random numbers. \n",
    "\n",
    "Assuming we can actually generate random numbers (a topic of philosophical and heated debates) we can populate a unit square with random points and determine the ratio of points inside and outside of a circle.\n",
    "\n",
    "![Unit circle and unit square](../images/MonteCarloPi.gif)\n",
    "\n",
    "![1/4 Unit circle and 1/4 unit square](../images/MonteCarloPi_rand.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The ratio of the area of the circle to the square is:\n",
    "\n",
    "$\\frac{\\pi r^{2}}{4r^{2}}=\\frac{\\pi}{4}$\n",
    "\n",
    "So if we know the fraction of random points that are within the unit circle, then we can calculate $\\pi$\n",
    "\n",
    "(number of points in circle)/(total number of points)=$\\pi/4$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.61262684, 0.8262925 , 0.59955407, 0.43229944, 0.66224665,\n",
       "       0.61154899, 0.34189469, 0.94126304, 0.41459117, 0.6400169 ])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.rand(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def montecarlopi(N):\n",
    "    '''Create random x-y-coordinates to and use ratio of circle-to-square to \n",
    "    calculate the value of pi\n",
    "    i.e. Acircle/Asquare = pi/4\n",
    "    Arguments\n",
    "    ---------\n",
    "    N: number of random points to produce between x=0-1 and y=0-1\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    our_pi: the best prediction of pi using N points\n",
    "    '''\n",
    "    \n",
    "\n",
    "    x=np.random.rand(N,1);\n",
    "    y=np.random.rand(N,1);\n",
    "    R=np.sqrt(x**2+y**2); # compute radius\n",
    "    num_in_circle=sum(R<1);\n",
    "    total_num_pts =len(R);\n",
    "    our_pi = 4*num_in_circle/total_num_pts;\n",
    "    return our_pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean value for pi = 3.188800\n",
      "standard deviation is 0.041097\n",
      "actual pi is 3.141593\n"
     ]
    }
   ],
   "source": [
    "test_pi=np.zeros(10)\n",
    "for i in range(0,10):\n",
    "    test_pi[i]=montecarlopi(1000);\n",
    "\n",
    "print('mean value for pi = %f'%np.mean(test_pi))\n",
    "print('standard deviation is %f'%np.std(test_pi))\n",
    "print('actual pi is %f'%np.pi)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "\n",
    "1. Why is there a standard deviation for the value of $\\pi$ calculated with a Monte Carlo method? Does it depend upon how many times you run the function i.e. the size of `test_pi`? or the number of random points `N`? Alter the script above to discover correlations\n",
    "\n",
    "2. How well does our function `montecarlopi` converge to the true value of $\\pi$ (you can use `np.pi` as a true value)? Plot the convergence as we did in [03-Numerical_error](https://github.uconn.edu/rcc02007/CompMech01-Getting-started/blob/master/notebooks/03-Numerical_error.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 2: Calculate the integral of a function\n",
    "\n",
    "One way to calculate the integral of a function, $f(x)$, is to approximate it as a Riemann sum as such\n",
    "\n",
    "$I=\\int_{a}^{b}f(x) dx \\approx \\Delta x \\left(f(a)+f(a+\\Delta x)+f(a+2\\Delta x)+...+f(b)\\right))$.\n",
    "\n",
    "Another way to integrate this function is with a [Monte Carlo approach](https://en.wikipedia.org/wiki/Monte_Carlo_integration) [2, 3]. We can approximate the integral as such\n",
    "\n",
    "$I=\\int_{a}^{b}f(x) dx \\approx \\frac{1}{n}\\sum_{i=1}^n f(x_i)$. \n",
    "\n",
    "where $x_i$ are uniformly random values of $x$ over the interval $[a,b]$. We were actually doing this in two dimensions in Example #1 when we compared the area of the circle to the area of a square. \n",
    "\n",
    "Visually, this approximation can be represented by the following figure.\n",
    "\n",
    "![Integration approximation](../images/integrals.png)\n",
    "\n",
    "The figure above shows the exact integral on the left as compared to two approximations, the Riemann sum (top-right) and Monte Carlo method (bottom-right). \n",
    "\n",
    "The main benefit in using a Monte Carlo integration over a Riemann sum (or another non-random method e.g. trapezoidal) is that the integration converges as \n",
    "\n",
    "$error \\approx \\frac{1}{\\sqrt{n}}$\n",
    "\n",
    "independent of the number of dimensions in the formula. Non-random methods have convergence rates that depend upon the step size along each axis, so with two dimensions reducing the step size to 1/3 increases n by 4, in three dimensions its 8, in 4D its  16, and so on. We'll do an example with a 1D function to introduce the concept. \n",
    "\n",
    "You can even integrate functions that are not easy to evaluate, such as the [area of light shown by the bat signal](http://leios.github.io/Batman_Montecarlo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise (Discussion)\n",
    "\n",
    "Why does the Monte Carlo method work? Is there a benefit to random numbers as opposed to using an ordered set divided into equally-spaced intervals?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we will calculate the integral  $I=\\int_{1}^{3}f(x)dx$ where\n",
    "\n",
    "$I=\\int_{1}^{3}\\frac{x}{1+x^2}dx \\approx \\frac{1}{n}\\sum_{i=1}^n \\frac{x_i}{1+x_i^2}$\n",
    "\n",
    "Our first objective is to take random numbers between 0 and 1 and turn them into random number between 1 and 3. This process is called _scaling_ or _transforming_ our variable $x_{i}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD9CAYAAAC2l2x5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAfoElEQVR4nO3de1xVVf7/8deSq1wsEVQmA7Rp0LE0k8z8alDmZdT020zeEtRuar+yLGuywKSywBnTccpqrNQGzbLvjPVN+/1yumg2komlppldFLt5B+8OCK3fH1zyiMABzuHA9v18PM4D2Je1PuccfLNde5+1jbUWERFxnia+LkBERLxDAS8i4lAKeBERh1LAi4g4lAJeRMSh/H1dwOkiIyNtXFycr8sQEWk0NmzYcMBaG3W2dQ0q4OPi4sjJyfF1GSIijYYxZldl6zREIyLiUAp4ERGHUsCLiDiUAl5ExKEU8CIiDqWAFxFxqAZ1maQ4S0FBAXl5eRw9epTi4mJflyPSqAQGBhIZGcl5551X6zYU8OIVBQUFfPfddzRv3py4uDgCAgIwxvi6LJFGwVrLyZMn+eGHHwgKCiI4OLhW7WiIRrwiLy+P5s2bExkZSWBgoMJdpAaMMYSEhBAZGcn+/ftr3Y6O4BurdDf/25Z+2Lt1VOLo0aNo2olzx+YfDnm0vU5tzvdoe41VeHg4Bw8erPX+OoIXryguLiYgIMDXZYg0av7+/hQVFdV6fwW8eI2GZUTqpq7/hhTwIiIOpYAXEXEonWSVehc3ZYVP+8/NHOjT/uPi4khKSmLhwoU+raMhy83NpW3btixYsICxY8d6pE1jDNOmTSM9Pd0j7TUGCniRerZs2TKaNWvm6zLkHKCAF6lnXbp08XUJlbLWcurUKQIDA31diiMVFxdjrcXfv36iV2PwIm44fvw47du3p1u3bpw6dap8+cqVK2nSpAlz5851u624uLhaDzucOnWKtLQ04uLiCAwMJC4ujrS0NJeaVq1ahTGGVatWuey7cOFCjDHk5ua61JKcnMz8+fNp3749gYGBrFhR+RDanDlz6NChA02bNqV58+YkJCSwbNkyl23e+7/LGXNDP7rHt6FHhxhuGtSbVSvfLl+/ZOE8Uob0pdclbenZMZbkwX348L133Hr+q1evpnfv3oSHhxMaGkq/fv3YsmWLyzbFxcWkpaURHR1NSEgISUlJbN26tdq2P/nkE4wxvPXWWxXW3XHHHURFRZW/zq+++irXXnstUVFRhIWF0aVLF15++eUK+xljSE1NJTMzk7Zt2xIYGMjnn3/u1nP1BB3Bi7ghNDSUJUuW0L17d6ZOnUpmZib79u1j9OjRDBo0iDvvvLNO7RtjGDNmTLXj8mPGjGHp0qU8/PDD9OzZk+zsbKZPn86OHTt45ZVXatX3Bx98wMaNG5k2bRotW7as9ANqixcvZvLkyTzyyCP06tWLkydPsnnzZvLy8riodJtXFsxjxiMPck2/gUyf/SxNQ0L5cssmfvzh+/J2fvr+e24YmcIFbWIoKipi9bv/j4ljRzD370vpeU2fSutcsWIFQ4YMYeDAgSxatAiAGTNm0KtXLzZv3syFF14IQHp6Ok8++ST33Xcfffv2JScnh8GDB1f7OnTr1o34+HiysrK4/vrry5cXFhaydOlSbrrppvLPduzYsYMbb7yRKVOm0KRJEz788ENuu+02Tp48yYQJE1zaXbhwIe3atWPmzJmEhobyq1/9qtpaPEUBL+KmLl26kJmZyeTJk7nuuuuYOXMmfn5+zJ8/v85t+/n54efnV+U2W7ZsYcmSJS4nCvv27Yufnx9Tp05lypQpdOrUqcZ95+fns2HDBlq3bl3ldtnZ2XTq1IlHHnmkfNmAAQOAkk+yHjt6hKdnPM61/Qcx+4Ws8m3+K6m3SzuTpz5e/v3PP//MlT0T2bXjW17PWlBlwN9zzz0kJiby5ptvli+75ppraNeuHU899RR/+ctfyM/PZ/bs2YwbN46ZM2cCv7xGU6ZMqfa1SElJYfr06Rw+fLh8kq+3336bvLw8UlJSyrd7+OGHXZ5DUlISu3fv5rnnnqsQ8NZaVq5cSdOmTavt39M0RCNSA5MmTaJ///4MGjSIlStX8ve//53IyMg6t1tUVMRLL71U5TYffvghAMnJyS7Ly35evXp1rfru3r17teEOcMUVV7Bx40YmTpzIu+++y4kTJ1zWb8r5hBPHj3HjTWOqbOeLzRu5a+xwrunyGy6Pi6Rr2yg+XvMBuTu+qXSfr7/+mm+//ZZRo0ZRVFRU/ggJCeGqq64qf20+//xzjh8/zrBhw1z2HzFiRLXPD0pey4KCAl5//fXyZVlZWcTHx9OtWzeXekaOHMkFF1xAQEAAAQEBvPjii2zfvr1Cm/379/dJuIMCXqRGjDGkpKRQUFBA586d6d27d/U7eUheXh4A0dHRLsvLwrlsfU2d2V5lRo8ezXPPPce6devo168fERER/P73vy8f0z90qKT/ltGVD0Hs+ekHxo0cwpFD+Ux5bAZ/f+MdXln+Pv+V1JvCgv9Uut++ffsAuPXWW8sDteyxfPny8vladu/eDUCrVq1c9j/z58rExsZy9dVXk5WVVfqcDrFixQqXo/djx47Rp08fNm3aRGZmJmvWrGH9+vXccsstFBQUVGjT3dfXGzREI1IDe/bsYdKkSVx++eV89tlnzJkzh3vuuade+o6IiCiv4aKLLipfvmfPHgBatGgBUD61bGFhocv+lU1a5e7H4Y0xjB8/nvHjx5Ofn8/KlSuZPHkyw4cP54V/vMP5zUv637dnNxe3/+1Z2/j3qvc4euQIf35uAa2iLyhf/p+TJ6vsu+y5ZWRkcN1111VYX3bVT1mY7t27l44dO5av37t3r1vPEUqGaW6//XZ27drFO++8Q2FhIaNGjSpfn52dza5du1izZg09e/YsX17ZnDG+nLJDR/AibrLWMmbMGAIDA/nXv/7FpEmTePDBB9m8eXO99J+YmAiUXMFxusWLFwNw9dVXAyVHoUCFq0vefvttPKV58+YMHz6cYcOGlfdzWUI3QkLD+McrFa8mKVMW5P7+v0xEl7vjGzbmrKuyv/j4eOLi4ti6dSsJCQkVHmXnHjp16kRoaChLly512f/M16wqQ4cOJTg4mMWLF5OVlcXVV1/tcuK5bGjq9Mn08vPzXc4NNBQ6ghdx06xZs3j33Xd5//33iYiIIDMzk1WrVjFy5EhycnLqNM7q7+/PmDFjqhyH79ixIyNHjiQ9PZ2ioiJ69OhBdnY2jz/+OCNHjiwPuejoaBITE8nIyCAyMpKWLVuyaNEivv3221rXBzBu3DjCw8O56qqraNmyJV999RVZWVn07dsXgNCwcO6eMpXMqQ9y37jRDLhhKKGhYWz/4nMCg4K56eZxXNkzEX9/f1In3cHocXdyYN9enpuVQesL2mB//rnSvo0xzJ07lyFDhlBYWMiwYcOIjIxk7969rF27lpiYGO677z7OP/987r33Xp544gnCw8Pp27cv69evr/b8xumaNWvG4MGDmTt3Lrt37+aFF15wWd+jRw+aNWvGnXfeyaOPPsrx48eZPn06kZGRHD7sm+m5K2WtbTCPrl27WnHTtGbuPXzkiy++8Fnf3vDpp5/awMBAm5qa6rL8yy+/tCEhIXbChAlutxUbG2vHjBnjsgyosOxsCgsLbWpqqo2JibH+/v42JibGpqam2sLCQpftvv/+ezto0CB73nnn2VatWtmHHnrIvvDCCxawO3fudKll1KhRbtW9cOFCm5iYaKOiomxgYKCNi4uzkyZNsocPH7abvs8vf8x8fqG95LKuNigo2IaGhdtLLutq/7pgSfn6Pz0738ZddLENDAqy7X7T3mY+86IdfONI+6s2F5Zvs3PnTgvYBQsWuNSwdu1aO3DgQHv++efboKAgGxsba4cPH27Xrl1bvk1RUZFNTU21rVq1ssHBwTYxMdFu3brVAnbatGluPdfly5dbwAYHB9tDhw5VWP/ee+/Zyy67zAYHB9t27drZOXPm2GnTptmSSP0FUOF3pqaq+7cE5NhKMtWUrG8YEhISbE5Ojq/LaBwa+A0/tm3bRocOHXzSt9Q/3fDDe6r7t2SM2WCtTTjbOo3Bi4g4lMbgRTzk559/5udqxpGr+zCTiCfpCF7EQ2655ZYK12if/qjPa+ZFQEfwIh6Tnp7OXXfdVen68PDweqxGRAEv4jFxcXGVTtQl4gsaohERcSgFvIiIQyngRUQcSgEvIuJQCngREYdSwIuIOJRbl0kaY9oADwIJQGegKdDWWpt7xnaVTWzTxVq7sQ51ipO4O4+O1/qv//l5xo4dy6pVq1xueH2uWJ/9EbcNu54Xl77FFVf1rH6H06Snp/Poo4/izTmz4uLiSEpKqvZ+uI2Ru0fwvwaGAfnAmmq2XQhcdcbjq1rWJ+IIU6dOZdmyZb4uwyc6XNKJrDdX0uGSmt8v9rbbbiM7O9sLVZ0b3P2g04fW2lYAxpjbgL5VbPujtfbjOlcm0ggUFBQQFBRU7Xan34GpIXL3edREcXEx1lrCwpvR6fIratVGmzZtaNOmjUfraoistZw6dar8zlSe4tYRvLW28hmURM4R6enpGGPYsmUL/fr1IywsrMLNnSszduzYOn3KddGiRXTu3Jng4GAiIyNJSUkpv/9oGWMM6enpLstyc3MxxrgMP4wdO5Y2bdqQnZ1Njx49aNq0KX/84x8r7dtay+zZs4mPjycwMJDo6Gjuuusujhw54rJd5wub8/SMx3lp7mx+16MzCe1a8vWXX7A++yM6X9ic9dkflW9bXFzMM3+aTu+u7bny4l9x2/DB7PzmKzpf2JznZmWWb1f2mp/5PNPS0vjrX/9K27ZtCQ8PJzExka1bt7pst3LlSgYMGEB0dDQhISFccsklPPXUUxQXF1f5Wp9p9+7d+Pv78/TTT1dYN2PGDAICAti/f3+N+oyLiyM5OZn58+fTvn17AgMDWbFiRY3qcoc3piq4wxjzAFAMfAxMs9ZWN6wj0mgMGTKEW2+9lQcffJAmTWp/nUJSUhK5ubnVjsvPmzeP8ePHM3z4cDIyMvjpp594+OGHWbduHZ9++ilhYWE17vvw4cOMGDGC+++/nyeffLLKu1GlpqaSkZHBnXfeyfXXX88XX3zB1KlT2bRpE6tXr3Z5Dd58fQltYmKZnPY4TZuG0LJVa44dPVKhzeeeyuDFZ2YxZsJEuvdMYtuWTdx9y01u179o0SLi4+OZM2cOhYWFPPDAAwwZMoQvv/wSf/+SWNuxYwe9e/dm4sSJBAcHk5OTQ3p6Ovv37yczM7OaHn4RHR3NddddR1ZWFhMnTqxQR//+/YmKiqpxnx988AEbN25k2rRptGzZ0ivTXHg64BcBy4GfgFjgAeB9Y0wfa+2qs+1gjBkHjAOIiYnxcDkinnf33Xd75Ebbfn5+5WFUmeLiYqZOnUpSUpLLfUXbt29Pr169mD9/PnfffXeN+z527BiLFi1iyJAhVW6Xl5fHrFmzGDNmDM888wwA/fr1IyoqipSUFJYvX87gwYN/2cFanl/8T4JP+4Ox4xvXU3BHDh1i0UvPMzT5Zu59+FEArrr6Gvz9A3jq8TS36g8ICGD58uUu90UdOnQon3zyCT169ABgwoQJp5Vl6dWrF4WFhcycOZMnn3yyRn+cU1JSSE5OZvv27cTHxwOwceNGtmzZwtSpU8u3q0mf+fn5bNiwgdatW7tdR0159DJJa22KtfY1a+0aa+0ioCclYT+9in3mWWsTrLUJZX8FRRqyG264wSPtvPfee3zzzTdVbrN9+3b27dvHqFGjXJb37NmT2NhYVq9eXau+/f39GTRoULXbffzxxxQUFJCcnOyyfMSIEfj7+1fov0dSb5dwP5uvv9zKyRPH6TPov12W9xk4uJI9KurTp49LuF966aUAfPfdd+XLdu/ezfjx44mNjSUwMJCAgADS0tI4dOgQ+/btc7svKHnPw8LCyMrKKl+WlZXFeeed5/IHriZ9du/e3avhDl6+Dt5aexRYAdTuDItIAxQdHV1vfeXl5VXaZ+vWrcvX11TLli3duvlIZf37+/vTokWLCv1HtWxVbZv79+0FIKJFpMvyFpEtq923TEREhMvPZSeI//Of/wAlN18ZPHgwy5cvJy0tjffff5/169eTmprqsp27QkJC+MMf/sDixYux1lJcXMySJUsYOnQowcHBteqzPn6P6mO6YAM0nBu/itTRmSf9vKksyPbs2VNh3Z49e0hI+OVWnEFBQRQWFrpsc/DgwbO26+5zOL3/jh07li8vKiri4MGDtGjRosbtlv0RyDt4wLXWAzU7qq7Kt99+S05ODllZWS7/+3jrrbdq3WZKSgovv/wyH330ESdPnmT37t2kpKTUus/6+D3y6hG8MaYZMBBY581+RJwqPj6eVq1auYy/A6xdu5Zdu3aRmJhYviw2NpYtW7a4bFfXKzO6d+9OUFBQhf5fe+01ioqKXPp318XtO9I0JJR/LX/DZfnK5W/WqdbTnThxAsBlGOfUqVMsXry41m1ec801tGnThqysLLKysoiLi6NXr15e7bOu3D6CN8bcWPpt19KvvzPG7Af2W2tXG2PuB+KBD/jlJOv9QGtg1JntiZzrevfuza5du6och/fz8+Oxxx5j/PjxJCcnk5yczI8//khqaioXX3wxN998c/m2I0aMYPr06TzxxBN0796dNWvWsGTJkjrVGBERwX333UdGRgahoaEMGDCAbdu2kZaWRs+ePRk4cGCN22x2/vkk3zqBF5+ZRUhYWPlVNMteXQSAqcOVSWU6dOhAbGwsqamp+Pn5ERAQwOzZs+vUZpMmTRg1ahR/+9vfOHXqFPfee6/LUbg3+qyrmgzRvH7Gz8+Wfl0NJAHbgRtKH+cBR4B/A7daaz+pW5niKD6YKqAhKi4upqioqNrtxo0bR0hICH/+858ZMmQIYWFhDBgwgD/96U8ul0g+9NBDHDp0iGeeeYbMzEwGDBhAVlYWV155ZZ3qfOKJJ4iKiuL555/n2WefpUWLFowePZqMjIxaXyZ6x+SHsNay7LVFLJk/j0u7dOXxWXMZc0N/wsOb1alegMDAQN544w3uuusuRo8eTUREBLfccgsxMTHcfvvttW43JSWFGTNmAFQ48eytPuvCeHOOh5pKSEiwOTk5vi6jcXB3Phcfhem2bdvo0KGDT/qW+rf5h0N1bmPl8jd44I6bWfA/Kxj7hwEeqMoZqvu3ZIzZYK1NONs63ZNVROrd5s9yWPPeSi7t0pWgoGC++Hwj8+f+hU6XX0GXblf5ujzHUMCL1EHZJXNVqe7DTOeikJBQPl23ltdefpHjx47SvEUUfa//b+55cFq9XqXkdPrNE6mDl19+2eVE59k0pGHQhuLX8R146fXlvi7D8RTwInVw/fXXs379el+XIXJWCniROmjRokWFD/uINBS6ZZ94jYYmROqmrv+GFPDiFX5+fpw6dcrXZYg0akVFRXU6Sa+AF68IDw+vcEMIEamZo0ePlk9mVhsKePGKiIgI8vPzOXDgAIWFhRquEakBay0nTpzgwIED1GUadZ1kFa8ICgoiJiaGvLw8cnNza3ybNGlc9uaf9Gh7245WPaf8uSAoKIhWrVrV6QheAS9eExQURHR0dL3Ony6+8bspnr2faG5mzScxk4o0RCMi4lAKeBERh1LAi4g4lAJeRMShFPAiIg6lgBcRcSgFvIiIQyngRUQcSgEvIuJQCngREYdSwIuIOJQCXkTEoRTwIiIOpYAXEXEoBbyIiEMp4EVEHEoBLyLiUAp4ERGHUsCLiDiUAl5ExKEU8CIiDqWAFxFxKAW8iIhDKeBFRBxKAS8i4lAKeBERh1LAi4g4lL+vC2j00s9zc7vD3q2jrpzyPBqIuCkrPNpebuZAj7Z3rjlX3w8dwYuIOJQCXkTEoRTwIiIOpYAXEXEoBbyIiEMp4EVEHMqtgDfGtDHGPG2MyTbGnDDGWGNM3Fm2CzbG/NkYs9sYc7J0+6s9XbSIiFTP3SP4XwPDgHxgTRXbvQTcDjwCDAJ2A+8YYy6rS5EiIlJz7n7Q6UNrbSsAY8xtQN8zNzDGdAZuAm6x1i4oXbYa2Ao8Bgz2SMUiIuIWt47grbU/u7HZYOAU8Npp+xUBrwL9jDFBtapQRERqxZMnWTsCO621J85YvhUIpGSYR0RE6okn56KJoGSM/kx5p62vwBgzDhgHEBMT48FyGhh353pxSr/iUZ6eS6WhO9eer7d48gjeALaS5ZWy1s6z1iZYaxOioqI8WI6IyLnNkwGfx9mP0puftl5EROqJJwN+K9DWGBNyxvLfAoXANx7sS0REquHJgP9fIAAYWrbAGOMPDAdWWmsLPNiXiIhUw+2TrMaYG0u/7Vr69XfGmP3AfmvtamvtRmPMa8BfjDEBwE7gDqAtMMqTRYuISPVqchXN62f8/Gzp19VAUun3NwNPANOB84FNQH9r7ad1qFFERGrB7YC31lZ5NUzpNieB+0ofIiLiQ5pNUkTEoRTwIiIOpYAXEXEoT05V4FvufiQ//bB365Cz8vRHz3MzB3q0PREn0hG8iIhDKeBFRBxKAS8i4lAKeBERh1LAi4g4lAJeRMShFPAiIg6lgBcRcSgFvIiIQyngRUQcSgEvIuJQzpmLxtPcndtGHMHTc+WINAQ6ghcRcSgFvIiIQyngRUQcSgEvIuJQCngREYdSwIuIOJQCXkTEoRTwIiIOpYAXEXEoBbyIiENpqgKpGXencEg/7N06RHzI01Nb5GYO9Gh7ZXQELyLiUAp4ERGHUsCLiDiUAl5ExKEU8CIiDqWAFxFxKAW8iIhDKeBFRBxKAS8i4lAKeBERh1LAi4g41Lk3F427c6lIg+bpuUBEnEhH8CIiDqWAFxFxKAW8iIhDKeBFRBxKAS8i4lAeDXhjTJIxxp7lcciT/YiISPW8dZnk3cD6034u8lI/IiJSCW8F/DZr7cdealtERNygMXgREYfyVsAvNsYUG2MOGmNeMcbEeKkfERGphKeHaA4DTwGrgSNAF+BhINsY08Vau+/MHYwx44BxADEx+jsgIuIpHg14a+1nwGenLVptjPkQ+ISSE69pZ9lnHjAPICEhwXqyHhGRc5nXx+CttZ8CXwFXeLsvERH5RX2dZDWAjs5FROqR1wPeGJMA/AZY5+2+RETkFx4dgzfGLAZ2Ap8Chyg5yfoQ8CPwtCf7EhGRqnn6KpotwEhgIhAC7AH+CUyz1h7wcF8iIlIFT19FkwFkeLJNERGpHX2SVUTEoRTwIiIOpYAXEXEoBbyIiEMp4EVEHEoBLyLiUAp4ERGHUsCLiDiUAl5ExKEU8CIiDqWAFxFxKAW8iIhDKeBFRBxKAS8i4lAKeBERh1LAi4g4lAJeRMShFPAiIg6lgBcRcSgFvIiIQyngRUQcSgEvIuJQCngREYdSwIuIOJQCXkTEoRTwIiIOpYAXEXEoBbyIiEMp4EVEHEoBLyLiUAp4ERGHUsCLiDiUAl5ExKEU8CIiDqWAFxFxKAW8iIhDKeBFRBxKAS8i4lAKeBERh1LAi4g4lAJeRMShFPAiIg6lgBcRcSgFvIiIQyngRUQcyqMBb4y50BjzP8aYw8aYI8aYfxpjYjzZh4iIuMdjAW+MCQHeB9oDY4AU4GLgA2NMqKf6ERER9/h7sK3bgXZAvLX2GwBjzGbga2A8MMuDfYmISDU8OUQzGPi4LNwBrLU7gX8DQzzYj4iIuMGTAd8R2HKW5VuB33qwHxERcYMnh2gigPyzLM8Dmle2kzFmHDCu9MdjxpjtNew3EjhQw33E2x41oPemodL70sCYGeXf1ua9ia1shScDHsCeZZmpcgdr5wHzatuhMSbHWptQ2/3Fe/TeNEx6XxouT783nhyiyafkKP5MzTn7kb2IiHiRJwN+KyXj8Gf6LfCFB/sRERE3eDLg/xfoboxpV7bAGBMH/FfpOm+p9fCOeJ3em4ZJ70vD5dH3xlh7tmHzWjRU8mGmTcBJII2S8fjHgXCgk7X2mEc6EhERt3jsCN5aexy4FvgKyAIWAzuBaxXuIiL1z2NH8CIi0rA0ytkkNalZw2SMaWOMedoYk22MOWGMsaXnYcSHjDE3GmP+YYzZZYw5aYzZbozJMMaE+7q2c50xpp8x5n1jzB5jTIEx5gdjzFJjjEc+HNrojuBLJzXbBBTwy1j/dCCEkrH+4z4s75xmjEkCXgM2AH5AX6CttTbXh2Wd84wxHwPfAW8CPwBdgHTgS6CHtfZn31V3bjPGjAQuB9YB+4EYYApwIXCptXZXndpvhAF/DyUTl50+qVlbSiY1+6O1VpOa+YgxpklZWBhjbgNeQAHvc8aYKGvt/jOWjQZeBnpba9/3TWVyNsaYeEr++N5vrX2qLm01xiEaTWrWQOlIsGE6M9xLrS/9ekF91iJuOVj69VRdG2qMAa9JzUTqLrH06zafViEAGGP8jDGBxpiLgb8Be4BX69qup+eiqQ+1mtRMREoYYy4AHgPetdbm+LoeAUrG4LuWfv8NJZeX76tro43xCB5qMamZiIAxJoySk61FwM0+Lkd+kQJ0B24CjgD/8sQVaI0x4DWpmUgtGGOCKZk2pB3Qz1r7g49LklLW2m3W2nXW2iVAbyCMkqtp6qQxDtFoUjORGjLGBAD/ALoB11lrP/dxSVIJa+0hY8w3wK/r2lZjPIL31aRmIo2SMaYJJVOH9AaGWGs/9nFJUgVjTCugPfBtndtqhNfBa1KzBswYc2Ppt72BCcD/oeQDHPuttat9Vtg5zBjzHCXvxRPA8jNW/6ChGt8xxiwDPgU2UzL2/hvgXqA10M1a+1Wd2m9sAQ9QOi3BbKAPJSdX3wMm6QM1vmeMqewXarW1Nqk+a5ESxphcKr+t26PW2vT6q0ZOZ4x5EBgGXAQEAt8Dq4AMT+RZowx4ERGpXmMcgxcRETco4EVEHEoBLyLiUAp4ERGHUsCLiDiUAl5ExKEU8CIiDqWAFxFxqP8PyZ8DpT/ZHVoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ri=np.random.rand(100)\n",
    "\n",
    "xi=1+(3-1)*ri # scaling equation for our uniform distribution 0..1 to 1..3\n",
    "\n",
    "plt.hist(xi,label='x_i: our scaled var')\n",
    "plt.hist(ri,label='r_i: our original var')\n",
    "plt.legend();\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Compare the scaled histogram to the original histogram. What is similar? What is different?\n",
    "\n",
    "Make a scaling equation to get uniformly random numbers between 10 and 20. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performing the integration\n",
    "\n",
    "If you have been doing a lot of calculus lately, you may know that \n",
    "\n",
    "$\\int\\frac{x}{1+x^2}dx = \\frac{\\log({x^2+1})}{2}$. \n",
    "\n",
    "If not, Python's [Sympy: Symbolic algebra](https://docs.sympy.org/1.5.1/modules/integrals/integrals.html) can help you jog your memory [4]. It will also give us an exact solution to judge our convergence. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\log{\\left(x^{2} + 1 \\right)}}{2}$"
      ],
      "text/plain": [
       "log(x**2 + 1)/2"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sympy\n",
    "\n",
    "x=sympy.Symbol('x')\n",
    "\n",
    "sympy.integrate(x/(1+x**2),x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0.80471895621705$"
      ],
      "text/plain": [
       "0.804718956217050"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.integrate(x/(1+x**2),(x,1,3)).evalf()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the exact integral that we are approximating is\n",
    "\n",
    "$I=\\frac{\\log{10}}{2}-\\frac{\\log{2}}{2}=0.80471895621705$\n",
    "\n",
    "_Note: We won't be using sympy for general numerical methods. It is helpful to do some algebra for us, but in general we want to use arrays so that our solutions can be scaled up and solved quicker by the computer._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7824846191561594\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'f(x)')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdcAAAE7CAYAAACPAss4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3wUdfrA8c+TIFXqEYqI9CJFiqFJFcSGFOUU0FOKBRTrqVhQj1PPguXn2RBOFPRU9BRFUKwQpaMiIE1EQIrSi3QIPL8/vpNls2ySDdnNbJLn/XrltdnvtGdnZ+eZ73dmviOqijHGGGOiJ8HvAIwxxpj8xpKrMcYYE2WWXI0xxpgos+RqjDHGRJklV2OMMSbKLLkaY4wxUZavk6uIjBARFZFOEY4/wBt/QGwjCyxvrYisjXDcEz6LiHTyykaEjKsikhLNWPOC7KzPeCPOXSKyQkQOet9hL7/jMlnL6HcY7/POj+JpfWUruYpICxEZKyI/i8heETkkIhtEZLKIXCcixWMVqIlPIjLO25ir58KycvXgJ5ddBTwF7AWeBf4JrMhsAm9dqIjsEZFTMxinUdB4C6MedebxpYhIrt1ILyIXi8j/RGS9d4CyR0QWisjTIlInt+KIR3n5wDNUbm9XJ6tQJCOJSCLwHHAzcARIAT4FDgCVgQ7AJcA/gKqxCDSXfAjMBf7wO5AwXgQmAOv8DiSOdfE7gBy4yHu9RFU3ZWO6VOBU4K/AuDDDB3rjRPRbz4tEpATwJnAp7uDkC+AXoDDQCLgFuE1EWqnqAt8CPTnzgTOBbX4HkkfEzfqK9Ac3EpdY5wF9VXVt6Agi0gV4PHqh5T5V3Q3s9juOcFR1G3GwwcQzVf3V7xhyoDJANhMrwDKgDC6JjgseICKFgL8BU4HuOQ8xbo3DJdZPgAHebyVARE7H7cNK5X5oOaOq+8miBcMcF1frS1Uz/QPqA0eBTUC5LMYtFPJ+AKDe66W4WuE+IMUbXhq4F5jhzf8wrmY2GqgUZv7jvPnVAh4AVgMHgeXAjWHGH+GN3wnX7LbIG38D8CiQmFG8YebVGZgMbAUOAb8BbwGNg8api2vaWwjs9Ja11Iv1lDDzXAuszeo7CP0sQWWdvLIRIeNq2jqOZFnhhgFlgcdwG+p+7/MsAV4GTg2aTsP8jfOGV097DzQGpnjz0aDlDAI+9tbnIdwBxCQgOYPvPvRvbVaf0YtjHK5F4rC3rOeB8mHGVVzLTEXgDS+e/V5Z80i+q6B5nQV8ELTNrPS2uxJhtrkMP1cm81dvW3sYOAbUDBne0xvnsrRxw8wjCXjBWyeHvXU0Dqie0XYClPTW3x+4bfw74LwwsYX7C91We3vrdjeuJexH4IZsrOPzvPkuA4pkMW6RoP8vBd7F7UMOeNvlV6GfI/R3hmulmwb8CazJ7HfoDeuF27/twe37vgOuzcbnCzvvSLdTjv8Gw/0NCBovAbgeV4Ha6/3NBi7LIK7mwJfeeDtx23k1b/kaMu4Iju+HrwMWe9tN2n7iNNw2PJ/jv5VVwNNAyexsV9H4Lkifty4A5njrdgtu/1c8ku8ukprrNd6KH62qOzIbUVVTMxjUB9dkNwm38tPGOxN3bmka8D5uI28C3ACcLyLNVXVnmPn9G0gG3sMl/suBl0Wkkqr+I8z4twBdveVPB3oAw3E193sz+0wAIvJ34BncD+pD3E6lKi7hfgf85I16Ga4GMQ33Qy0CdAQe8eLNExeoiIgAn+Ni/gKX/ArjDmoG4JLuXtypggG47+zfwC5vFqHn9moDs4AfgFeB04OGveSN/wWwHfcD7QVcICKdVHWuN95HuBpaT9z3mLaMXWRCROoBM4G/ePP4GbdjuAXoJiKtVXVryGRlvHi343Zc1XDf7dcicqZGULsUkfa4dZiI20434H74w4GuItJRVQ96n+OfuPVYzfs/y88VYhzuAG4A8FBQ+UDcTndKBjEm4Q54a+J2lG/jDhCvwa2btqq6MmSyU3DfVSnvc5UF+gGfikiyqi72xgv3mcD9/tOW/xRwFy6xv4vbgXUFRnvr+Y4IPvtA7/UZVT2U2Yghwx/H7cS/wR3YV8Ztd5+LyOWqOjHMLNoC9+M+/yjcQUaGROROXILYCozHHbxcCrwqIk1U9dYsPltWItlOd+HW/+3eNM8FTb/Qi1OAd4ArcAcp473h3YAPROR2Vf130OdqBnyL27+9h6sQtcf9zjLLEfcA7XD7k6m49Q7ugOUO4GtcQlegFXAn0MHbDo9442a5XYVzkt9FT9zpmkleXOcDNwLlgL6ZLQ+IqOY63fuw50Z6tBXmCCAVaB9meGnC1IZxTVkKPBBSPs4r/x2oGFReDljjLadOmCOmHUDtkPG3445gCoc7Ygkqa4pL4GuA00LiKRQSR5Xg+XllAvzHm2+7kGFricOaK67GpcCzYcYtFbLO0r6T6mHGrc7xI8v7M1h2jTBlZ3rfzVcZHVFG8jlCtt+rQ8of8spfC7PuFHewIEHl//DK74vgu0rAHXkfAzqGDBvrzeehkPIUQo74I1hOoDbqTb82LWagAm4H8u/QcYOmf90rfzCk/BqvfFqY9au4WsopQeX9vfLRkX4mXI1AcQc8RYPKT/HKFGgRwTpY441bM6txI9juKuIOglaFlHcK2i6uDDNd2vARQWW1cNenbCSoFQ53fnyRN36HCOI8Yd4ns52G+20EDRvsTfMSQa15QAlcTfYQQfs+XEJXoFvIfNL2cxpSPsIr3w3UD7P8CgS15gSVP+BN97dsbFc5/i44vp85DLQOKi+KayU9BlTJ6ruL5GrhSt7r76EDROQS7xaR4L/6YebxoarOCC1U1d0avjb8Fu6LOC+DmJ5X1c1B89mBa45NxB1Fhxt/Vcj4H+NWbr0MlpFmMG5neZ+qplsHqpoaEsdGVT0cMo7ijnLJ5PPEqwOhBar6Z+hnjMAfuO/nBKq6JkzZclxSbC8ihbO5rAAROQP3Y/tRVd8MGTwS2Az0C7OMfbiDAQ0qG+e9Jkew6Ha4H/RHqvpNyLD7cTur/hHMJztexx3Nd/be/w2XqF4PN7L3mfviag8jg4ep6hu4Ws25IhLuAsU79XhNAtzvNZXI1k2aobgd2A3qavBpyz6C26mCa/HKStr+aWM2lp3RdrcZmAjUyuDq9+9V9e0IF3EV7uD7SQ1q6VDVvbgECDnfBnK6naYZimvavUNVj6YVquo+XKtbYVyNGG+9nAPMUdVPQuYzguOtkuGMUdUTzoeq6hZvWaFe9l5zut882e/ibT3ecoa3nU7AVZiaZ7XQnF5BeAku+QRbyIknlL/PaAbehVC3Ay1xTXeJQYMrZzDZzEzKmoQZ9mOYsrQfY5mMYvO08F6/yGI8RCQBuBZ35NMQV8uToFEy+jzxZhnu/Op9ItIUd6HIDGBJyA85UotCdsYBIlIbl3DOxZ17CU10f+Hkr95u6r2mhA5Q1YMiMhfX9FOP4037AL+E+bFHur1ktdzNIrICaCIiJVV1TwTzi8T7uCvKB+Ca1wbiaqoZ3X5TH3ckPkfDN6d+g/scTYH1QeW7NOSCRlVNFZHNRLZu0rTCnWa5ybVKpnNKUIwxISKVgPtwzX5VcesiWGVcbS9YhvuxMDLcBnAHjsHjnKycbqd4t082wjXt3h/mu0jyXtO+i7T96+zQEVV1o4isw51mCCezPHA5Lpc0xZ1qCK745XS/ebLfRU7yRkTJdTNuxZ6GO18VoKpDgCHgOjng+FFAqC3hCkXkCtyRwB7gM9zGnFZbuh3Xph9O6DmytDgh/BWB4a4ATjvCSgwzLFhp4FAGNexQLwA34c4hTeT4RVplgNvI+PPEFW9n2Rl3kcFlwMXeoPUi8qiqjsnmLDP6/uvgLmIoiTtH/SHuXO4x3PmvJuRsnaVtC5szGL4pZLw0J2wv3jqBrLeXSJfbxBsvKslVVfeJyHtAXxF5HbfDvD2TSaK2bjypRLZu0pTD7X8y2meAa5bMyibc6YcquIuTsiQi5XDb3em4g/KpuM91DNfS0ZHw213Y7TgDGa5fVd0tIofI+dXLOd1OwSUywbV6RPJdpJ1nDrcPBreOMkquGe0H7sa1nmzBfRcbcRc84cWU0/3myX4XOckbESXXObiNrRPHs3x2ZVTb+QcumTbXoNsovBPswzKZXxIhiR53vgTc0XA07QJqi0i5zBKsiFTEnexeBLRR1QNBw1rhkqufjnG8RhCqFCHrTd1FPjeKyFDcjrorbmc9WkS2quqH2Vh2Rt//7bgDjytV9Z3gAd46C9cKkR1pn6liBsNjtc34tdzXcVdfv4E7x/RWJuP6FWPw8g+o6ulZjpm52bjkei4RJldc61JVXJNqutsHRWQUbn8XTnZabYLXb7qduoiUxiWMWK3b7EiLYZaqtotg/LSDwaQMhlfIZNoT1p93u9gDuNOOTTToNipvn5pZwo+UL99FJOdcx+N2zIO9I75oqgUs1xPvT2wGFMtkunAbQVrZomgEFuQ77/X8LMargTsC/Co4sXraRjmmk7ELqOh1CBIgItVwR69hqeoxVV2sqs9w/Hx2j6BR0s7RZKfWkqaW9zo5JKZihD+nkd1lpTWJdggdICJFcE2TBznxQC2nMltuBVxL0OooNgkDoKozcZ0nVAGmaMj9niF+xn321hmc106LPSe/p6MQ6IQm1HygSgbndLMj7Zzy373vNENBnzOj7U6ANjmMJ02G2wDHk3du9ph1lDC/G28bXAE0yqiXrxBp28M5oQNEpApwRjbjKo87uJ8TZnvNaL+Z2XYVji/fRZbJ1TsB/Rwu63+awYl+cM2n2bUOqOPtcAAQkVK4K+Ayc6t3VJM2TTngbtxKfyfDqU7OGNzBxWMikq7tX0QKBcWe1nNSGwk6cSEidXHndvz2A67mGrjgS0ROwV2eno6I1MjgwrS0dR588JBWm69yEjGlrbPAj8hbd48R/gg4W8tS1XW4c4dni0joxTF34c7lTDiJC7SyMhNXi7pUREJ3Qo/izu+9EeVlprkCd4vBXZmN5J1nfRe3DtLd8iIiV+EOblK8dXiyMvu+XvBex3q1h3S8bbB6VgtQ1a9w55sbAO+FqwCISGUR+S/HE8IJ253ndnLeWpLmbdz+6G7vlqe0WEpw/BaSWG0D4ewAymdwAPICbv/9soiEnntGRBqm7ee88+1zcfu5biGjjiD71/FswWu99A6q05ZZGbcfCCe7+xxfvotIV8QwXNV5KLBSXKfwS3BHvhWB1riN+w8gO73kvIhLpAtE5ANvGRfhmghOuDo5yAJgkXeOKe0+1yrAw6r6SzaWnyVVXSgiw3BJaIWITMR9ztNwV7E9DTynqr+LyIe4Hdt3IjLdG6cH7n7H3tGM6yS8hLvY5TUR6YqryXbBNfOEXjDUBPjQu+BnKe4HUAN3HnQ/8ErQuNNxO/JXROR9b/hPYa4kDOcV3IU3E0XkXdw5jvbeslJwpyKCzcVtc7d5B2HbgN2qOoqM3YhLdm97F02sxCWOC3C3cdwTQZzZoqrHRGQQ7vzRNG873Yg7Sm6DO9AZmckscrLshUR+FD7Mi+kJETkX97uqi9uGt+PWXU5Mx3XL+K6IfIa7Snqmqs5U1U9F5HHcgecqEfkcdxtMEu5WrNbAlZx4UVE4A3D7sl7AWm9ev3plDXFNxgm4e9XBdZV4D/Ci97k34K6ubY27eC80aWSbqq4SkfuBJ4GfROR/HL+3sgbwcpgryWNpOu4zThKRWbjTBp+o6k+4uxnOAa4GOonINI7f+3sW7mKfNhw/Z3oz7j7XD71tO+0+1+q4DiIaRxqU91t5BXeA96OIfII7H3+Jt4xwd3NkuF1lsAx/vous7tXR9PcQtQJewzU97fM+1EbcBnkDXs89euL9QgMymJ/gEvYy3NHLelyyLUn4exbHefOrBTyI2zkewjVr3BRm/iMIuTc0s2GZxYtrFv6M4z0vrcX9SBsGjVMSV8v/jeO9M92K+wIVr0eSoPFP+IyZrPtw8XYiwvtcgz7D994624w7uDlhXeMu9HgCd4/bFu+zrMadIjgzzHzvw+3MjgR/ToJ6aMrkc3XBnTfbg9uhfwDUIYP7Z3EHKz9424uGxB12feIusHiD9L2AvQhUCDNu2HWX1bAMxm+Ku7Btu7fcVbij8VPDjJtCDu5zPdlxcS0EL3rr5LC3jt4g/H2gGW6v4YbhWkqewf2uUzPYVi/C7T+2ecvfiGttuJMwPWhl8RkvxtVi13vb+D7czv5Zgu5z98Zthus4YyfuoO5z3J0BI4jwdxbJcNwFgTNxF+rt97bd67PxmcLOO7vbKa7pdaz3/R4lzH4Od8vKdG+dHPK2ic9xB1klQsZNxl2EuM8bfyLu9/4T7oA3eNwT1mnI8MK4+85X4fY1v3rTFM7gs2S4XUXjuyDzPJDhsNC/tBvO8wQRGYe7H6mGhunf2BhjjD+8c7ZbcLfstfQ7Hr/l6+e5GmOMiS4ROUVE/hJSloBrdi2G6y6wwMu3j6EyxhgTE6WBDSLyBe4ahhK4uzUa4U7RZXVBaoFgydUYY0x27MXdAtUFd7FYEdz5z+eAR9R1K1jg5alzrsYYY0xeYOdcjTHGmCiz5GqMMcZEmSVXY4wxJsosuRoTRERURGLS56uI9BeRH0Rkv4hsF5EPMuhmMqPpB3jxZfY3Nsx054nIZyKywVv2LyLyqrjH/QWPd4qI9BaRN0RkhTfubhGZKSLXZBDTiExisQtbTIFlVwsbkwvk+CMZ1+AeAl0W97DyriLSVl03dFlZyPG+UEMNwD02LN1zh0XkNtxVnNtxj/TbgbtlYhDQx1v2Ym/0Wrgejv7EPRP2I9zzdC8Dxnvjhj6/Oc14TuyqMNp9NhuTZ9jVwsYEERHFPdw9pw+yDp5nfVxf3CuBVuo9Dcfr1P9b3BNB2udg/mVx/UPvA05T7+Hn3oMZtuG6iGukqn8ETTME16fsOFUd6JVVAboD4zX9IxOTcE+xqe7FPz9o2AjcQcO5qppysp/BmPzGmoWNib0BuMd9/UuDHjOnqrNxjz1rl53m4TCuwt1r+HZaYvX8Bden7I/BidWT9mCFwFNCVHWjqr6iIY9MVPds3zHe23CP7TLGhLBmYWNiL+2ZkV+GGfY57mkuHXC925yMQd7rayHlm3HNwc1EpKKqBj8o+mLvdVqEyzjivaZmMLy9iLTEdQq/HPg6JNEbU6BYcjUmC95zRQdkY5K1qjou6H0dYI+qbgkzbtojEmuHGRZJbE1wT3j5UVV/DB6mquqdcx2He9TWh7gnmDTEPXLvVY4/VzWzZSTiasfgnoQSzsMh738XkatVNdLkbUy+YsnVmKxVx51XjNQ3uISWphTHn4UZ6k/v9YQHhkfoWu81tNYKgKq+JSLbgLdwj4VMMw94Q1WPhJsuxEO4x+eNV9UlIcMW4p5U9Q2upnw67kKt+4EpItJCVZdG+mGMyS/snKsxWVDVFFWVbPx1CjebaMclIoVxDxQ/hEue4ca5Dnde93XcQcKpuAdjC+5B7pk+GFxErsY9O3kxcEvocFX9SFXfUNXfVPWgqq5S1UeB23BPSLn3JD+eMXmaJVdjYu9PMq6ZlvJed5/EfHvhLlr6SFV3hg70LpIaBUxW1bu9BLhPVefgHjp/BHgio5mLyF9xSXkF0DX4YqwIjMedn22bjWmMyTesWdiYLEThnOsvQGsRqRDmvGsd73XVSYSWdiHTCR1HeLrifuMpoQNUdbOILAeaikhhVU13T6qIXAq8A6wGumRwvjhDqnpYRPYAxbMznTH5hSVXY7JWnZydc/0GaI1LdqHNtxd4r99mJyAROd2b3zpchw/hFPFey2cwPAk4RsgVwCLSA3jXm/e5YW7jiSS+2riOMuZld1pj8gNrFjYmC1E45zoOd4vKcBEpmVbodSLRHZilqiuCyk8RkfoiUiuTsAbgfr/jVPVYBuPM9l5vEJHKwQO8c7FVgRnB03vnYP8H/I5LrBszCkBECotIcpjyMsB/vLcTMvkMxuRb1kOTMUFi0UOTN98RHO/+cCJQBuiHS7rpuj/0mqHXAL+pavUw8xJcU3NNoKaqrs1kuROAPrhzuh8CW3FX/nYF9gMdVfV7b9z6uKt/iwBveDGEWqiqH3njl8Hd2rPQ+9sMnAZchKstfwb0iPCKZGPyFWsWNiYXqOoIEVkL3AoMBQ7gOpAYrqrLszm7Drh+gKdlllg9VwEzgGuAvwJFcbcFvQU8pqrLgsatxPGm5LAd9eMuVPrI+/8A8Dyuybsbrhl4H/AT8F/gP5nUqo3J16zmaowxxkSZnXM1xhhjosySqzHGGBNlllyNMcaYKLPkaowxxkSZXS0conz58lq9enW/wzDGmDzlhx9+2KaqSVmPWTBYcg1RvXp1vv/+e7/DMMaYPEVEfvM7hnhizcLGGGNMlFlyNcYYY6LMkqsxxhgTZZZcjTHGmCiL2+QqIlVF5H0R2S0if4rIRBE5I8JpNYO/qHbGbowxxoQTl1cLi0hxYBpwCOgPKPAoMF1EzlLVfRHMZhwwOqRsZTTjNMYYY8KJy+QKXI97nFY9VV0FICKLcY/ZGgw8G8E8Nqrq3NiFaIwxxoQXr8m1BzA3LbECqOoaEZkF9CSy5GrixNGjR/nll19YuHAhixYtYvHixezatYuSJUtSqlSpdK/ly5enZcuWNG/enMKFC/sdujHGnJR4Ta4NgUlhypcCl0c4jxtF5G7cw6jnAv9Q1RlRis9kYdmyZYwaNYr58+fz008/ceDAgWxNX6xYMVq1akX79u1p164dbdq0oWTJkjGK1hhjoitek2s5YGeY8h24BzJn5b/AFOB3oBpwNzBNRLqqakroyCJyA3ADwBlnRHTNlMnAggUL+Ne//sXEiRNzNJ8DBw6QkpJCSkoKAImJiXTt2pUBAwbQs2dPihYtGoVojTEmNuLyYekichh4RlXvCyn/F3CPqmbroEBESgJLgPWq2i6zcZOTk9W6P8y+2bNn869//YtPP/007PDKlSvTtGlTmjRpQpMmTahcuTJ79+5lz549gb8///yT1atXM3PmTNasWZPhssqUKUO/fv0YOHAgycnJiEisPpYxJkIi8oOqJvsdR7yI15rrTlztNVRZwtdoM6Wqe0TkE+DanAZm0lu8eDG3334706dPP2FY9+7duemmmzj77LNJSspef94bN25k5syZzJgxg5kzZ7Jo0aLAsF27djFq1ChGjRpFgwYNGDJkCNdeey3FixfP8ecxxphoiNea6zSgcGgtU0RScDF3PIl5jgIGqmqm7YlWc43cxx9/TL9+/di/f3+gTES4/PLLuf/++2nSpEnUlrV27VreeOMNxo0bF7ZWm5SUxJ133smNN95IqVKlorZcY0xkrOaaXrx2IvEx0FpEaqYViEh1oK03LFtEpBTQDZgXpfgKNFXlmWeeoVevXoHEmpiYyDXXXMOyZct49913o5pYwT2t6KGHHmLVqlWkpKTQv3//dDXVrVu3cu+991K9enVGjBjBjh07orp8Y4zJjnituZYAFgEHgAdwnUg8ApQEzlLVvd541YBfgYdV9WGv7C6gHjCd4xc0pZV1yeqKYau5Zu7IkSPcfPPNjBkzJlBWo0YNJk+eTMOGDXM1lj179jB+/HhGjhzJ+vXr0w079dRTGTp0KPfeey9lypTJ1biMKYis5ppeXNZcvR6YOuN6VHoTeAtYA3ROS6weARJJ/zl+BhoAzwNf4u6JXQO0s1txcmbXrl1cdNFF6RJr27ZtmTdvXq4nVoCSJUty8803s2rVKl599VVq164dGLZ3716efPJJateuzQsvvMCRI0dyPT5jTMEVlzVXP1nNNbzVq1fTrVs3VqxYESi76qqrGDt2LEWKFPExsuNSU1N57733eOyxx1i6dGm6YXXr1mXkyJH06NHDri42Jgas5ppeXNZcTXzZvHkzHTp0SJdYH374Yd588824SawAhQoV4sorr2Tx4sVMmDCB6tWrB4atXLmSXr160blzZ3744Qf/gjTGFAiWXE2mUlNT6du3Lxs3bgSgSJEivPPOOzz44INxWwNMSEigT58+LF++nJEjR1K6dOnAsJSUFFq0aMHgwYPtoidjTMxYcjWZGj58eKCXJBHho48+om/fvv4GFaGiRYty9913s2rVKm6++WYSExMBd7XzmDFjqFevHuPGjcNOjRhjos2Sq8nQhx9+yMiRIwPv//nPf3LhhRf6GNHJKV++PC+88AJLly6lW7dugfJt27YxcOBAOnbseMI5WmOMyQlLriaslStX0r9//8D7iy++mOHDh/sYUc7Vq1ePyZMn8+GHH1K1atVA+YwZM2jatCnDhg1j375IHhVsjDGZs+RqTrBv3z569+7Nnj17ANeBw5tvvklCQt7fXESEXr16sXz5coYNG0ahQq4H0NTUVJ566ikaNWrEV1995XOUxpi8Lu/vLU1UqSqDBw9myZIlgLuA6YMPPqBcuXBdPeddJUqU4Mknn+THH3+kffv2gfK1a9fStWtXBg0axM6d2e7G2hhjALvP9QQF/T7Xl19+maFDhwbejx07lkGDBvkYUXQ9M+LEMlXl+0Xj+fjzv3Pg4PGEWvLUSlx28Us0PvOyQNmdYaY3xth9rqGs5moCli1bxu233x54f+211+arxJoREaFF0wEMG7qcsxpcHijfs3cT49/rzfj3evPnnj98jNAYk9dYcjUBw4YNC3QT2KxZM1588UWfI8pdJU+tyDWXv8eAPh9S6tTKgfKflk/kqZcbsmDxW3bbjjEmIvH6PFeTy77++ms++eQTwNXkXn/9dYoWzfTpfPnG7JTQkl50aNKJJb/ezdpNrwJw4OBO3v7wbyw68z16XzKaUqdWAqyZ2BgTntVcDceOHeOuu+4KvB8wYEDUHxmX1xQuVIbm9f5Du7O+onjR6oHypT9/zFMvNbBarDEmU5ZcDf/9739ZuHAhAMWKFeORRx7xOaL4UaFsF7qcvZgap90YKEurxY5/7zI2bdrkY3TGmHhlybWAO3DgQLrOIe666y6qVKniY0Tx55RCJWlW52XanfUVZUtXC5QvWfERDRs25N133/UxOmNMPLLkWsA999xzbNiwAYAKFSpw9913+xxR/KpQtgt33vgTbc4eEijbsWMHffv2pU+fPmzbtl3DjAcAACAASURBVM3H6Iwx8cSSawG2ZcsWHn/88cD7hx9+mJIlS/oYUfwrWqQkvS8ZxQ1Xf0mZ0mcEyt977z0aNWrExx9/7GN0xph4Ycm1APvnP/8Z6OLwzDPP5Nprr/U5oryjbs3zuOvGn9Kts82bN9OzZ08GDBjArl27fIzOGOM3S64F1IoVKxg9enTg/VNPPRXoZ9dEpmiRUrz66qtMmTKFypWP3xc7fvx4GjduzBdffOFjdMYYP1lyLaDuuecejh49CkDnzp25+OKLfY4o7+rWrRtLlizhyiuvDJRt2LCBCy64gJtuuom9e/f6GJ0xxg+WXAugb7/9NnBuUER4+umnERGfo8rbypUrx1tvvcX7779P+fLlA+WjRo2iSZMmzJgxw8fojDG5zdoBC6Bnnnkm8P/VV19Ns2bNfIwm70v/MIDeDB3Qng+mDGbJio8AWL16NR07dKRDm79zYedHOaXQ8Z6vrIcnY/Inq7kWMGvXrmXy5MmB93n9AejxqGSJCvS/YiL9er1B0SKlAVCUb+Y8w/+Nbs763wvuU5eMKSgsuRYwr7zySqDbvvPPP5+6dev6HFH+JCKc3eRq7rppCXVrnR8o37JtOS+82prPpj1I6tHDPkZojIklS64FyMGDB3n11VcD74Of22pio0yp07n+qs/o3e0VCp9SAoBjepSvZjzKv//TItDtpDEmf7HkWoC8++67bN++HYBq1arRrVs3nyMqGESENsmDufPGxdSs1jFQ/sfmxbRo0YJHHnkk8Kg/Y0z+YMm1AHnppZcC/994440kJib6GE3B85eyNRnSfxo9L3iOUwoVAyA1NZWHHnqINm3asGTJEp8jNMZEiyXXAmL+/Pl89913ABQpUsR6Y/JJgiTQvvVt/H3IQqqd3iZQ/sMPP3D22Wfz2GOPkZqa6mOExphosORaQATXWvv27ZvuXkyT+5L+UpehA2cwcuRIChcuDMDhw4cZPnw4rVu3tlqsMXlc3CZXEakqIu+LyG4R+VNEJorIGVlPecJ87hMRFZGZsYgzL9i6dSsTJkwIvLcLmeJDQkIid999Nz/++CMtWrQIlFst1pi8Ly6Tq4gUB6YB9YH+wNVAHWC6iJTIxnxqAsOBLbGIM68YO3Yshw+72z5atmyZbkdu/NegQQNmz57NE088YbVYY/KJeO2h6XqgJlBPVVcBiMhi4BdgMPBshPMZBbwF1CN+P2tMHT16lFGjRgXeW601fsxOAUakvStEIe7h1uu68+6kgazfOB9wtdimTZrTpcMDdG53L4USXfK1np2MiW9xWXMFegBz0xIrgKquAWYBPSOZgYhcCTQH7otJhHnElClTWLduHQDly5fniiuu8Dkik5lKSQ24edAsLu7yBIleIj167AhfpPyDf49Jtt6djMkj4rU21xCYFKZ8KXB5VhOLSFng/4BhqrqjIHdKH3wh03XXXUfRokUzGdv4ZXZK8LtCFOUezm3WgwU/D2LHnrkA/LHlJ/79n1bUrXo3R478g1NOKRaYwmqyxsSXeK25lgN2hinfAZSNYPqngJXAuEgWJiI3iMj3IvL91q1bIw4y3v388898+eWXACQkJDBkyBCfIzLZUarEmXRsNpPGtZ4lMSEtkR5j5foneXZ0U9asK7DX6BkT9+I1uQJomLIsq6Ai0h64BrhR0zrRzWpBqmNUNVlVk5OSkrIZZvx65ZVXAv93796datWq+RiNORkiidQ5/Q66JP9EUplzA+Vbt6/kpdfb88EnN3Hg4G4fIzTGhBOvyXUnrvYaqizha7TBRgNjgQ0iUkZEyuCavxO990WiG2p8Onr0aLrbb2666SYfozE5dWqxWrQ76yua1RlNocSSgfI534/iqZcbMGlSuLMoxhi/xGtyXYo77xqqAbAsi2nPBIbgknDaX1ugtff/jdELM37NmjWLTZs2AVChQgW6dOnic0Qmp0QSqHHaDZyXvJQz614SKP9zz+/06tWLv/71r/zxxx8+RmiMSROvyfVjoLV3nyoAIlIdlyQ/zmLac8P8LQKWeP+/H/1w48///ve/wP+XXXaZ9SOcjxQvWpVBfT/mb70ncGqJCoHyDz74gDPPPJMxY8Zw7NgxHyM0xsRrcv0PsBaYJCI9RaQH7urh9bhmXwBEpJqIpIrIQ2llqpoS+gfsAnZ77zfk6ifxwbFjx/jggw8C7y+/PMsLrE0eIyI0bdSHYUOX07LpoED57t27GTx4MB07dmTp0qU+RmhMwRaXyVVV9wGdcVf8vonrCGIN0FlV9waNKkAicfo5/DJr1qxA82BSUhIdOnTwOSITK8WLleOKnmP5+uuvqVWrVqB85syZNG3alPvuu4/9+/f7GKExBVPcJiVVXaeqvVW1lKqWVNVeqro2ZJy1qiqqOiKLeXVS1XaxjDeehDYJFyoUr7czm2jp3LkzP/30E/fdd1/g+05NTeWJJ56gYcOGfPrppz5HaEzBYnvdfMaahAuW2SnePyMAivGXwo9x+/VX8cEnQwL3wa5du5Zu3bpxVoO/0vOC5yhdqgpgHU8YE0txW3M1J2f27Nn8/vvvgOvusGPHjj5HZHJbpQoNuXHAN1ze/VWKFzt+R9viZe8z8qX6pMx+mqNHj/gYoTH5nyXXfMaahA24h7K3an4tw4au4Owm1wTKDx3ey5Qv7+bZ0U1JSUnxL0Bj8jnb8+Yjx44d4/33j99pZE3CBdPslOB3SVQrM57iTQay8Jeh7NnvbhPfvHUZ5557Ls0aX0n3rk9TqmRlwJqKjYkWq7nmI3PmzEnXJNypUyd/AzJxI6lMJ7qcvZBGNZ8iMeH4I5F//OltnnyxHt/Medaaio2JIkuu+Uhwk/Cll15qTcImnYSEU6hb9S66tlzB6Ul9AuWHDu9h8hd38swrZ/HFF1/4GKEx+Ycl13zCmoRNpIoXOZ2WDSYw+OqvqFC+fqB8y7YVXHDBBfTq1YvVq1f7GKExeZ8l13xi7ty5bNy4EYC//OUvnHvuuVlMYQq6OjW78Pchi7ik69MUKXz8YQCTJk2iQYMGDB8+nL1792YyB2NMRiy55hPWJGxORqHEwnQ6507uuWUlyU0HBMoPHTrEY489Rr169Rg/frz1VWxMNllyzQesSdjkVKlTK9G35+vMnTuXli1bBsp///13BgwYQIsWLfj22299jNCYvMWSaz4wb948NmxwzyMoV66cNQmbk9aqVSvmzJnDuHHjqFSpUqB8wYIFdOzYkd69e/Prr7/6GKExeYO1HeYDoU3Cp5xyio/RmLxodor3zwhwx9z9uWVgb6bPfJKUOU+TmnoQgIkTJzLpo8m0bXkL53UYHugByu6PNSY9q7nmcarKxIkTA++tSdhES5HCp3Jh50e45+afad74qkD50WNH+Hbuszz+fC1SZj/NES/xGmOOs+Sax/3888/89ttvAJQqVYrOnTv7HJHJb8qWPoMrL/svt1w7l+pVzwmUHzi4iylf3s2TL9bjzTfftIuejAliyTWPC77pv0uXLtYkbGKm2umtGDpwJlf/9T3Kl6sdKN+1ex3XXHMNZ599tnVCYYzHkmseF7wzO//8832MxBQEIkKThpdz101L6XXRC5QonhQYtnDhQi644AK6dOnCvHnzfIzSGP9Zcs3DDh06xPTp0wPvLbma3FIosTDtWt7Mfbeu4rz2D1CsWLHAsGnTptG6dWsuvfRSli5d6mOUxvjHkmseNmfOHPbv3w9ArVq1qFmzps8RmYKmaJFSXNj5EVatWsX1119PYmJiYNhHH31E48aN6d+/P2vXrvUvSGN8YMk1D7MmYRMvTjvtNMaMGcOyZcvo0+f4QwFUlTfeeIO6desydOjQQBedxuR3llzzMEuuJt7UrVuXCRMmsGDBAi666KJA+ZEjR3j55ZepVasWt912G5s2bfIxSmNiz5JrHrV161YWLFgAQGJiovXKZOJKs2bN+PTTT/n2229p27ZtoPzQoUM8//zz1KxZk7vuuoutW7f6GKUxsWPJNY/6+uuvUVUAWrduTenSpX2OyJgTtW/fnhkzZvD555+n67P4wIEDPPPMM9SoUYN77rmHLVu2+BilMdFn3R/mUdYkbPKCZ0YACHA+l1/YlUa1PuXzlIfY+Idrddm3bx8jR47k/559gTbJN9LpnLsoVbIyYF0qmrzNaq55kKpacjV5jojQoG43br/+ewb0+ZDKFc8KDDuSeoBv5z7LY/+uwUdTb2XXnxt8jNSYnLPkmgctX748cNVlmTJlSE5O9jkiYyInIjSq34s7Bv9I/ysmclqlpoFhqUcPMXP+Czz+fC2GDBnC6tWrfYzUmJNnyTUPCu3y0B6MbvKiBEmg8ZmXcscNCxjY92OqntYiMOzo0cOMHj2aOnXqcOWVV7J48WIfIzUm+yy55kHByfWCCy7wMRJjck5EaFivO7deN4/rrppKtdPbBIYdO3aMd955hyZNmtC9e3dmz57tY6TGRM6Sax5z6NAhUlJSAu+7du3qXzDGRJGIUL/2hdw8aBZDrpl2wrY9ZcoU2rZtS4cOHZgyZYo9hcfEtbhNriJSVUTeF5HdIvKniEwUkTMimK6aiEwSkd9E5ICIbBORFBG5KKtp84KZM2dy4MABwN2wX716dX8DMibKRITaNc7liy++4LvvvuOyyy5DRALDZ8yYQffu3WnUqBFjx47l0KFDPkZrTHhxebJORIoD04BDQH9AgUeB6SJylqruy2TyU4FtwAPABqAUcD3wqYj0VtWJmUwb9+wqYVMQzE4BRgAkc07jD6hdeQXTZj3JgsX/5dixVMBd2Hfddddxx+0P0K7lrbRJHkLxYmUBu43H+C9ea67XAzWBXqr6kapOAnoA1YDBmU2oqktV9VpVfVNVp3vT9sIl2oGxDjzWLLmagqhC+fr07fk699+6mo5t7qJI4ZKBYXv2bmLqtPt59P+q8uGnt7B1+y8+RmqME6/JtQcwV1VXpRWo6hpgFtAzuzNT1VRgN3AkahH6YPPmzSxcuBCAQoUK0alTJ38DMiaXlSldle7nP8UDd6yn23kjKV2ySmDY4SP7mPXdi4x8sR49evRg2rRpgV7MjMlt8ZpcGwJLwpQvBRpEMgMRSRCRQiJSSUQeBOoCL0Uxxlz31VdfBf4/55xzKFmyZCZjG5N/FStamnPb3s19t62mb6/xVK7QODBMUSZPnkyXLl1o2rQpr7/+OgcPHvQxWlMQxWtyLQfsDFO+Aygb4TxG4mqqfwDDgL6q+nV0wvOHNQkbk16hxMIkN7mGvw9ZxA1Xf8mZdbqlG7548WIGDRpE1apVuf/++1m3bp1PkZqCJl6TK7iLmEJJmLKMPAe0ALoDU4G3ReSScCOKyA0i8r2IfB+vT+mwLg+NyZiIULfmeVx75RSGDV3BTTfdRPHixQPDt23bxuOPP06NGjXo3bs306dPtyZjE1Pxmlx34mqvocoSvkZ7AlXdoKrfq+oUVb0CmAs8ncG4Y1Q1WVWTk5KSTjroWFq+fHngGZjlypWjefPmPkdkTHyqUL4eL730EuvXr+eJJ56gatWqgWHHjh1j4sSJdO7cmcaNGzNq1Cj+/PNPH6M1+VW8JteluPOuoRoAy05ynt8DtU86Ip/NnDkz8H+HDh1ITEz0MRpj4l+5cuW45557WL16NR9++CGdO3dON3zp0qXcdNNNVKlShSFDhrBo0SKfIjX5Ubwm14+B1iJSM61ARKoDbb1h2SIiCUA74NcoxZfrgpNru3btfIzEmLylUKFC9OrVi6+//jqQUEuUKBEYvnfvXkaPHk3Tpk0555xzeOONNwIdtRhzsuI1uf4HWAtMEpGeItIDmASsB0anjeT1xpQqIg8FlY0QkedFpI+IdBSRPsBnQEvgH7n6KaLIkqsxOdegQQNeeuklNm7cyAsvvECDBulvPpgzZw79+/enSpUq3HbbbSxZEu6mBWOyFrPk6t0C01xE2opIPREpHOm0Xg9MnYGVwJvAW8AaoLOq7g1eDJBI+s+xAGgEvAB8gbtq+CDQXlUn5OQz+WXjxo2sWbMGgGLFitGsWTOfIzImbytdujQ333wzS5Ys4ZtvvqFv376ccsopgeE7d+7k+eefp3HjxrRp04axY8eyd+/eTOZoTHpR7f5QRJKB64ALgaohgw+LyHfAO8BbqprpVQSqug7oncU4awm5glhVP+Ykmo7j2axZswL/t2rVisKFIz5OMcZkQkTo0KEDHTp0YMuWLbz22muMGTMmcDALMHfuXObOncsdd9xBv379GDhwIK1atUrX37ExoaJScxWRZBFJAeYDbYDJuC4MewEXAP1wPYVuA54A1ovIcBEpGo3l53fWJGxM7FWoUIF7772XVatW8eWXX3LFFVekq83u2bOHMWPG0KZNGxo2bMhTTz0VuILfmFDRahb+BlgINFTVJqp6i6q+pqqTVfUrVX1PVZ9U1cuAisANwGXA3VFafr4WXHO15GpMbCUkJHDeeefx7rvvsnHjRp5++mnq1auXbpzly5czbNgwTj/9dLp3787EiRM5fPiwTxGbeBSt5FpLVW9X1eVZjaiqB1X1XVU9GxgTpeXnW3v27An0J5yQkECbNm2ymMIYEy1JSUnceeedLF++nBkzZjBo0KB0VxofPXqUKVOm0Lt3bypXrszQoUOZN2+edVBhopNcVfWk2kZUdXM0lp+fzZ07N/BQ6LPOOotSpUr5HJExBY+I0K5dO8aOHcumTZt4/fXX6dChQ7pxduzYwcsvv0zr1q2pX78+jz76KGvXrvUnYOO7qF8tnNVDyUXEmoKzwc63GhNfTj31VAYMGMA333zDL7/8wvDhwznjjDPSjbNy5UoefPBBatSoQYcOHXjllVfYvn27TxEbP8TiVpxPROT/Qm+98W7N+Qp4LAbLzLcsuRoTv2rXrs2jjz7KmjVrmD59OgMHDjzhaVUzZszgxhtvpFKlSnTv3p0JEyawf/9+nyI2uSUWyXUI7krhuSJSD0BEugOLcY996xKDZeZLR44cYe7cuYH3bdu29TEaY0xGEhIS6NSpE6+99hqbNm3i7bff5sILLyQh4fguNjU1lSlTptCvXz8qVKjA3/72NyZPnsyhQ4d8jNzEStSTq6qOwT2NJhH4QUQmAh8B3wJNVPXbaC8zv1q4cGHgCLd69eqcfvrpPkdkjMlK8eLF6devH1OnTuX333/n+eefp1WrVunG2bdvH2+99RY9evSgYsWKDBo0iM8//5wjR474FLWJtpj00ORdNTwQl2B74TrN76uqET3RxjjWJGxM3laxYkVuueUW5s6dy6pVq3jkkUeoX79+unF2797N66+/zoUXXkjlypUZPHgwX331FampqT5FbaIhJslVRG7A1VRXAPfjuiOcIyJ59qk0fghOrtYkbEzeVqtWLR544AGWLVvGjz/+yL333kv16tXTjbN9+3bGjBlD165dqVSpEtdff73VaPOoWFwt/AHwCvAq0EpVn8B1ml8E+FFEBkZ7mfmRqlrN1Zh8SERo2rQpjz/+OKtXr2bevHn8/e9/p0qVKunG2759O6+++ioXXnghFStWZODAgUyePJmDBw/6FLnJjljUXNsBl3idShwGUNWluPOw43BJ12Rh1apVbNmyBYAyZcqc8PQOY0zeJyK0bNmSZ555hnXr1jFjxgxuu+22ExLtzp07GTduHD169CApKYkrrriCd955xx70HsdikVybqOqnoYWqekhVbwF6xGCZ+U5ok3DwVYfGmPwnISGBdu3a8dxzz7Fu3Tpmz57NHXfcQdWq6Z+BsnfvXv73v/9x5ZVXkpSUxMUXX8yYMWP4448/fIrchBOLq4Uz7a1JVT+J9jLzI2sSNqbgSuvq9Nlnn+W3335j7ty53HPPPdSpUyfdeIcPH2bq1KkMHjyY0047jVatWvHYY4+xZMkS64LRZ9F6Ks6lJzFNZRFpHY3l50eWXI0x4JqOW7VqxRNPPMHPP//MkiVLePjhh8M+13n+/PkMHz6cxo0bU6dOHe644w6+/vpre6iAD6JVc31JRBaJyBARKZfZiCLSXkTGAKuAs6K0/Hxly5YtrFy5EoDChQuTnJzsc0TGmHggIjRs2JAHH3yQBQsWsHr1ap577jnOPfdcEhMT043766+/8txzz3HeeedRvnx5Lr/8csaPHx+4lsPEVrQell4HuBN4GHhBRJYDi4CtwCGgLFATSAZK427T6aqqs6O0/Hwl+BFzLVq0oGhRe+ytMeZENWrU4LbbbuO2225jx44dTJ06lUmTJjF16lT27t0bGG/Pnj28//77vP/++4GLqC666CIuvvhizj77bLumIwai9VScfar6MFAFuArXacTZwCDgDqA7rkOJf+Oe+XquJdaMWZOwMSa7ypUrx1VXXcV7773Htm3b+Oyzz7jllluoUaNGuvFUlXnz5jFixAhatmxJpUqVuOaaa5gwYQI7duzwKfr8J1rnXHeISHNVPQJcCDyiqg1UtYyqFlXVKqraRVX/qaororHM/MySqzEmJ4oUKcIFF1zA888/z6+//srSpUt58sknad++/Qm11K1bt/Lmm2/Sr18/kpKSuO+++3yKOn+JVltACVwnEQADgKQozbfA2bdvHwsWLAi8P+ecc3yMxhiT14kIDRo0YNiwYXz77bds3bqVCRMm0L9/fypUqJBu3GPHjlGzZk2fIs1fonXO9TfgehFJS7DNRCTDE4XWeX/G5s+fH+hTtGHDhpQrl+n1YcYYky3lypWjT58+9OnTh2PHjrFgwQKmTp3Kp59+yvz587nookwfyW0iFK3k+gQwGugPKPByBuOJNzwxg+EF3rx58wL/W63VGBNLCQkJJCcnk5yczIMPPsjOnTspW7as32HlC1FJrqr6mohMxT2vdTpwK7A8GvMuaL777rvA/y1btvQxEmNMQWOJNXqiVXNFVf8A/hCR8cAnqromWvMuSObPnx/4v0WLFj5GYowx5mRFLbmmUVV76s1J2rRpExs2bACgWLFiNGzY0OeIjDHGnAy7cziOBDcJN2/enEKFon7sY4wxJhdYco0j1iRsjDH5gyXXOGIXMxljTP5gyTVOqGq65Go1V2OMybviNrmKSFUReV9EdovInyIyUUTOiGC6ZBEZIyIrRGS/iKwTkbdEpEZW0/ppzZo1gX49y5YtS61atXyOyBhjzMmKy+QqIsWBaUB9XMcUV+OevDNdREpkMXlfoCHwPHARcC/QHPheRKrGLOgcCj3fKiI+RmOMMSYn4vVy1Otxj6irp6qrAERkMfALMBh4NpNpn1TVrcEFIjILWOPN96GYRJxD1iRsjDH5R1zWXIEewNy0xArgdUoxC+iZ2YShidUr+w33bNkqUY4zaiy5GmNM/hGvybUhsCRM+VKgQXZnJiJnAhWI0y4ZU1NT+eGHHwLv7UphY4zJ2+I1uZYDdoYp3wFkq/NLESkEvIKruY7NYJwbROR7Efl+69YTKr4xt3z5cvbv3w9AlSpVqFy5cq7HYIwxJnriNbmCe3pOqJO5yudF4Bzgb6oaLmGjqmNUNVlVk5OScv9RtNYkbIwx+Uu8XtC0E1d7DVWW8DXasETkceAGoL+qfhGl2KIu+EphaxI2xpi8L16T61LceddQDYBlkcxARIbjbsO5VVXfjGJsUWc1V2OMyV/itVn4Y6C1iNRMKxCR6kBbb1imRORW4FFguKq+EKMYo+LgwYMsXrw48D45OdnHaIwxxkRDvCbX/wBrgUki0lNEegCTgPXA6LSRRKSaiKSKyENBZX2B54DPgGki0jroL9tXGsfawoULSU1NBaBu3bqUKVPG54iMMcbkVFw2C6vqPhHpDPwf8CbuQqavgdtVdW/QqAIkkv4g4UKv/ELvL9g3QKcYhX1SrEnYGGPyn7hMrgCqug7oncU4awm5glhVBwADYhVXtFlyNcaY/Cdem4ULDEuuxhiT/1hy9dHu3btZsWIFAImJiTRr1szniIwxxkSDJVcfBXd52LhxY4oVK+ZjNMYYY6LFkquPrEnYGGPyJ0uuPrKemYwxJn+y5Oojq7kaY0z+ZMnVJ5s2bWL9+vUAFCtWjIYNw/X2aIwxJi+y5OqT4Fpr8+bNKVQobm85NsYYk02WXH1iTcLGGJN/WXL1yYIFCwL/W2f9xhiTv1hy9cmPP/4Y+N86jzDGmPzFkqsPNm/ezO+//w64i5nq1avnc0TGGGOiyZKrD4JrrU2aNCExMdHHaIwxxkSbJVcfWJOwMcbkb5ZcfWDJ1Rhj8jdLrj4ITq7Nmzf3MRJjjDGxYMk1l+3evZtVq1YBUKhQIRo1auRzRMYYY6LNkmsuW7RoUeD/Bg0aUKRIER+jMcYYEwuWXHOZnW81xpj8z5JrLrPzrcYYk/9Zcs1lwd0eWs3VGGPyJ0uuuejgwYMsW7Ys8L5JkyY+RmOMMSZWLLnmoiVLlnD06FEAateuTalSpXyOyBhjTCxYcs1Fdr7VGGMKBkuuucjOtxpjTMFgyTUX2W04xhhTMFhyzSVHjx5l8eLFgfeWXI0xJv+y5JpLfv75Zw4cOABAlSpVqFChgs8RGWOMiZW4Ta4iUlVE3heR3SLyp4hMFJEzIpz2MRH5QkS2i4iKyIAYh5slO99qjDEFR1wmVxEpDkwD6gP9gauBOsB0ESkRwSxuAYoBU2IWZDbZ+VZjjCk4CvkdQAauB2oC9VR1FYCILAZ+AQYDz2YxfWlVPSYitYFrYhpphCy5GmNMwRGXNVegBzA3LbECqOoaYBbQM6uJVfVYDGPLNlW15GqMMQVIvCbXhsCSMOVLgQa5HEuOrV27ll27dgFQtmxZqlWr5nNExhhjYilek2s5YGeY8h1A2WgvTERuEJHvReT7rVu3Rnv2J9RaRSTqyzDGGBM/4jW5AmiYsphkJVUdo6rJqpqclJQU9flbk7AxxhQs8Zpcd+Jqr6HKEr5GG9csuRpjTMESr8l1Ke68a6gGwLIw5XEt+B5X67DfGGPyv3hNrh8DrUWkZlqBiFQH2nrD+q4ArwAACUdJREFU8ozNmzfzxx9/AFC8eHHq1q3rc0TGGGNiLV6T63+AtcAkEekpIj2AScB6YHTaSCJSTURSReSh4IlFpKOI/BW40CtKFpG/emW5KrhJ+KyzziIxMTG3QzDGGJPL4rITCVXdJyKdgf8D3sRdyPQ1cLuq7g0aVYBETjxI+CfQMej9UO8vbZpcY+dbjTGm4InL5AqgquuA3lmMs5YwyVJVO8Umquyz863GGFPwxGuzcL6xcOHCwP9Nmzb1MRJjjDG5xZJrDO3du5dff/0VgMTERBo1auRzRMYYY3KDJdcY+umnn1B1fWHUr1+fokWL+hyRMcaY3GDJNYaCm4SbNGniYyTGGGNykyXXGFq0aFHgf0uuxhhTcFhyjaHg5GoXMxljTMFhyTVGjh49yuLFiwPvreZqjDEFhyXXGPn111/Zv38/ABUrVqRixYo+R2SMMSa3WHKNETvfaowxBZcl1xixziOMMabgsuQaI1ZzNcaYgsuSa4xYcjXGmILLkmsMbN++nQ0bNgBQpEgR6tWr53NExhhjcpMl1xgIrrU2atSIQoXi9uFDxhhjYsCSawxYk7AxxhRsllxjwK4UNsaYgs2SawxYzdUYYwo2S65RdvjwYZYtWxZ4f9ZZZ/kYjTHGGD9Yco2yFStWcOTIEQCqV69OmTJlfI7IGGNMbrPkGmX2DFdjjDGWXKPMzrcaY4yx5Bpl9gxXY4wxllyjSFWtWdgYY4wl12j6/fff2b59OwAlS5akevXq/gZkjDHGF5Zcoyj0fGtCgq1eY4wpiGzvH0XWJGyMMQYsuUaVXSlsjDEGLLlGlV0pbIwxBiy5Rs2+fftYuXIlAAkJCTRq1MjniIwxxvglbpOriFQVkfdFZLeI/CkiE0XkjAinLSoiT4nIHyJyQETmiEiHWMa7ZMkSVBWAunXrUqxYsVguzhhjTByLy+QqIsWBaUB9oD9wNVAHmC4iJSKYxVjgeuAh4BLgD+BzEYlZW601CRtjjElTyO8AMnA9UBOop6qrAERkMfALMBh4NqMJRaQJcCUwSFVf98q+AZYCDwM9YhGwXSlsjDEmTVzWXHEJcG5aYgVQ1TXALKBnBNMeAd4NmjYVmABcICJFoh+uXSlsjDHmuHhNrg2BJWHKlwINIph2jaruDzNtYaB2zsNL79ixYyxevDjw3pKrMcYUbJJ2EU48EZHDwLOqem9I+aPAvaqaYXO2iHwBlFLV1iHl5wFfAh1UdUbIsBuAG7y39YCfcxB+eWBbDqYvaGx9ZY+tr+yx9ZU9OVlf1VQ1KZrB5GXxes4VIFzWlwimk+xOq6pjgDERxpX5wkW+V9XkaMyrILD1lT22vrLH1lf22PqKnnhtFt4JlAtTXtYblpkdmUybNtwYY4yJmXhNrktx505DNQCWRTBtDe92ntBpDwOrTpzEGGOMiZ54Ta4fA61FpGZagYhUB9p6w7Ka9hTg8qBpCwF9gC9U9VC0gw0RleblAsTWV/bY+soeW1/ZY+srSuL1gqYSwCLgAPAA7hzqI0BJ4CxV3euNVw34FXhYVR/+//buLVSzMY7j+PeXzCaUXKBQGoc7JJKkmNEMKdxMiqKUC6dCyPGCMibCKBdyJTmmJEWUQzY3YxxySJjIyFzIOEXIMP4u1trZ2HtmzFrttV77+6m3993v3rt+7661/+t51nqe/6zffxw4FbgG+Ay4mGYziROq6u0F/CiSpEVolCPXqvoJWA5sAB4CHqEpkstnCmsrwC78+3NcADwA3Ao8CxwEnGZhlSQthFGOXCVJmmSjHLlOki4NBhajJAcmubdtpvBzkmqvp2sOSVYleTLJ520Tio+TrEmy19DZxijJqUleTvJlkl+TbEryRJLtbT6jVpLn2+Py1qGzTDKLawc9NBhYjA4FzqZZUvXadn5WcDWwFbgBOA24j+YegheSePz+2z7AW8BlwErgepqVB+vaezS0DUnOAdxirgdj3kRiEux0g4FF7NWq2g8gyYU0/wA1vzOqavOsr6eTfAs8CJxMc3KnVlU9Bjw2+70k64GPgFXAXUPkmgRJ9gbWAlcCjw4cZ+J55ttNlwYDi1JV/TF0hknyj8I64432+YCFzDLBvmmffxs0xfjdAXzQnqCoI4trN10aDEg766T2+cNBU4xYkl2SLElyGHA/8CVNZyzNIcmJwPnAJUNn+b9wWribfZh7O8Zv+Wu7Rak3SQ6g6Uv8YlW9OXSeEXsdOKZ9/QnNMr6vBswzWkl2pTkBubOqujQt0SyOXLvb2QYD0n+SZE/gaeB3mrXcmt95wPHAucAPNDeAHTxkoBG7FtgdWD10kP8TR67ddGkwIO2wJLvRbO25FDipqjYNHGnUqmpmyvz1JM8BG4HrgIsGCzVC7bLBG4ELgakkU7O+PdXe5PRjVW0dJOAEc+TaTZcGA9IOaaftngSOA06vqvcHjjRRqup7mqnhQ4fOMkJLgd2Ah2kGBDMPaJaBfQccMUy0yWZx7aZLgwFpu9q1rI8ApwBnVdW6gSNNnCT70axF/3ToLCP0DrBsjgc0BXcZdhLbKW5/2MGONhjQ3yVZ1b48hWaa7hJgM7C5qqYHCzZCSe6j+RutBp75x7c3OT38d0meAt4G3qO51no4zbrN/YHjqmrDgPEmRpICVlfVTUNnmVQW147aaxZrgRU0NzK9BFxRVRuHzDVm7YE7l+mqOnkhs4xdko3AfDsL3VJVNy9cmvFLci3NDmCHAEuAL4BXgDUekzvO4tqdxVWSpJ55zVWSpJ5ZXCVJ6pnFVZKknllcJUnqmcVVkqSeWVwlSeqZxVWSpJ5ZXCVJ6pnFVZKknllcpZFKskeSj5KsbzvjzLy/MskfSS4dMp+k+bn9oTRiSY4G1gFrq+q6JPvSbEq/vqrOHDadpPlYXKWRS3IlcBewkqbH5hHAUVX19aDBJM3L4iqNXJIAzwLLaTq9rKiql4ZNJWlbvOYqjVw1Z8APAVPAuxZWafwsrtLIJdkfuIemCfhRSS4fOJKk7bC4SiPWTgk/CGwBVtAU2duTHDloMEnb5DVXacSSXAXcASyvqukkS2juHp4Cjq2qXwYNKGlOjlylkWqX4dwGrKmqaYCq2gKcAxwM3D1cOknb4shVkqSeOXKVJKlnFldJknpmcZUkqWcWV0mSemZxlSSpZxZXSZJ6ZnGVJKlnFldJknr2JyX0Vfrsz2CuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f= lambda x:x/(1+x**2)\n",
    "\n",
    "n=10\n",
    "ri=np.random.rand(n)\n",
    "\n",
    "xi=1+(3-1)*ri # scaling equation for our uniform distribution 0..1 to 1..3\n",
    "# we scaled our variable by a factor of 2 so we scale the bins accorcdingly\n",
    "bin_scale = 3-1 \n",
    "\n",
    "my_I=sum(f(xi))/n*bin_scale\n",
    "print(my_I)\n",
    "\n",
    "plt.bar(xi,f(xi),width=bin_scale/n,color=[0.1,0,1,0.5])\n",
    "plt.plot(np.linspace(0,4),f(np.linspace(0,4)),color='k')\n",
    "plt.title('Graphical illustration of Monte Carlo'+\\\n",
    "          ' integration\\n I={:.4f}'.format(my_I))\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('f(x)')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def montecarlo_integrate(f,a,b,N):\n",
    "    '''Estimate the integral of a function f(x) from a->b with N \n",
    "    randomly chosen points\n",
    "    Arguments\n",
    "    ---------\n",
    "    f: a function of the form f=lambda x: ...\n",
    "    a: the lower bound of the integration\n",
    "    b: the upper bound of the integration\n",
    "    N: number of random points to produce between xi=a-b\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    our_I: the best prediction of the integral using N points \n",
    "    '''\n",
    "    ri=np.random.rand(n)\n",
    "\n",
    "    xi=a+(b-a)*ri # scaling equation for our uniform distribution 0..1 to 1..3\n",
    "    # we scaled our variable by a factor of 2 so we scale the bins accorcdingly\n",
    "    bin_scale = b-a\n",
    "\n",
    "    our_I=sum(f(xi))/n*bin_scale\n",
    "    \n",
    "    return our_I\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above function, we have defined a general Monte Carlo integration function. Now, we can create a for loop to try the function multiple times and get an average and standard deviation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean value is     I = 0.81410075695566\n",
      "standard deviation is 0.028138512733038436\n",
      "actual            I = 0.80471895621705\n",
      "error in mean       = -0.009381800738609947\n"
     ]
    }
   ],
   "source": [
    "test_I=np.zeros(10)\n",
    "for i in range(0,10):\n",
    "    test_I[i]=montecarlo_integrate(f,1,3,500);\n",
    "\n",
    "print('mean value is     I = {}'.format(np.mean(test_I)))\n",
    "print('standard deviation is {}'.format(np.std(test_I)))\n",
    "print('actual            I = 0.80471895621705')\n",
    "print('error in mean       = {}'.format(0.80471895621705-np.mean(test_I)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example 3: Determine uncertainty in failure load based on geometry uncertainty\n",
    "\n",
    "In this example, we know that a steel bar will break under 940 MPa tensile stress. The bar is 1 mm by 2 mm with a tolerance of 10 %. What is the range of tensile loads that can be safely applied to the beam?\n",
    "\n",
    "$\\sigma_{UTS}=\\frac{F_{fail}}{wh}$\n",
    "\n",
    "$F_{fail}=\\sigma_{UTS}wh$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAElCAYAAAAskX9OAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZgcVb3G8e8LSVgFEllkD/uu6I0simwiOwFlEUVBFvEiCrJdZRFDAFFRZBeQ3aCggBBQISAEBMMSNChhEwibLCYEgUAgBH73j3OadDrdPTUz1TPTmffzPP3UdNWpqtM1M/Xrs9Q5igjMzMzKME9vZ8DMzOYeDipmZlYaBxUzMyuNg4qZmZXGQcXMzErjoGJmZqVxUGlDkk6SFJI2qVq3VV53XG/mrQhJA3Jeb+3tvED969kk7ao57YU9kTebu/W1/4UyOKi0iKSh+Y+l0eu/vZ1H6ztysDpB0o2SXsx/I+O7eKwhkk6R9JCkaZJelnS3pL0lzVsn/YKSTpT0mKR3JE2RdK2kdbv/yZrmc938Ob/VzePMK+kISQ9LelvSS5IulrR0J46xhqTjJN2V939H0tOSzpe0fJ30y0s6XNKfJT0vaYakf0saJWnt7nyeDvJ5V75mr0tavM72SpDq0t9OGQb01on7kceAK+usf7sbxzwdGAU8041jWNc8A6wFlP2lYHPgeGAm8Ajwka4cRNKHgQeAFYHbgD8CCwM7A5cBmwH7V6VfEBgLfBL4G3AOsCSwG7CNpC0j4t6u5KWA4Xk5upvHuRT4CvAgcAawArA3sJWkDSLipQLHOBnYFZgAXA28BWwIHAjsJunTEfFoVfpDgSNI/99/BKYC6wF7AbtK2joi/tLNz9XMh4BjgcNaeI6uiQi/WvAChgIBXNdD59sqn++43v7sBfI6IOf11t7OS87PSTk/m/RiHlYh3cTmr7o+47twnGPyviNr1i8CPJW3LVe1/ui8bhQwT9X6DUgB7qHq9U3OOwqY2cm83gNM6OZ12zbn/xZgQNX6vfL6SwoeZ19g3Trrv5ePM7pm/a7Ap+qk/1JO/4+C5+3U/wJwV07/JOmL6QoNjtfpv52yXq7+6gMkzSfpUEm35iL0DEkv5KL0KnXSF2oDaFb/32hbLso/IenDki7IVTHv17TfrCbpsqq8Pifp7HrF8c6StH6uepmSqyAelzQyf6OuTXtAri56JqedLOk6SZ9ocOwlJV2Sjz1N0p2SPtPJ/DW6bmtKuqIqLy9LGifpiCLHjYgnI+LeiOhOCRZg5bz8Y83xXwcq35yrf0875+UJEfF+Vfr7gD8A6wCdukZFSFqKFLi6W0o5IC+Pi4iZlZURcQXwT2BPSQt3dJCIuCQiHqqz6afAO8CmNemviYi/1jnOb0jBez1JixX/GLOTNDD/PYWkH9dJcjwwH3BCV8/RKg4qfcMSwM8Akf7Jfg7cD+wJ3CdpxR7Oz/ykKpFPA78DzgfeAMjB5YGct7tJVXETgG8Cf+3mP9JmwF+BHUg3xZ8DrwLfB/4sab6aXc4FhgBjctrbSN9c75Y0rObYHwLuBL4G/AM4C5hC+oa7QVfznI+9PHAv6QZ9F3Aa8HvgPaqqmnrIw3m5bfXK/Pk3AV4EqqtxlsrLp+scq7Ju89JyN8tOzPp7747NgNeB++psG0P6W+7O7zdIJbaZHSWs8m5evteVE0paALgO+DLwfxHx3TrJrgfGAV9tZRtOV7hNpfXWlDSizvorY1Yd7RRg+Yh4sTqBpC1JN72jgf9taS5ntywwHtgjImZU5Wc+4Dekf5pPRsRjVdu+BPwaGAF8p7MnVGpAvpT07WuziLgrrxdwCbAPqQ77h1W7rRERk2qOsw7pBn8Ss99YjwbWAH4eEYdXpT+QFDS7Y3dS9dIOETFbCSG3cfSkC4CvAj+Q9GlSO8nCpPaL94Bda0pDU0hVtSsCT9Qca2hert6CfA4HXiB9QekSSYuSSl1/j1z3U+Nfebkq6QtHV+wGLATcUDBPG5L+zu6JiDc6e7L8pexGYCNg/4i4uEnyo0lf/k4GPt/Zc7VMb9W7ze0vZrWpNHrtUvA4DwNP1Kybow2AOm0qpH+mAC6sc9y624Dn8/o16uyzR972nQZ5nQC8VOAzzVGPDGyZ1/22TvplSFUQjxW8Zn8EpjN7HfuzwDRgcE3aeUg3n0JtKvWuG3B4XrdFSX873aoXJwWR39X8vb1NCvgL1KQ9MW+/jNnbVIaRvjwE8IcC5yzcpgIsALwJnNfN67RCzt/YBtsr7RtHdfH4ywEvkRrtVy+QfpH8//oe6YtRZ37Xt5I6ZzyYf1efb5C+0qaycH7/p/x+wzL+dsp4uaTSetdHxC4dJZL0SeAoUpXTksxeinyrRXlrZFpUlUKqbJiXH21Q+poPWErSYhHR2d5R6+fl2NoNEfGCpH8B60haICKmA0handQwvTmwNDCoZtfBwGRJg4Hlgfsi4tWaY78vaRwpWHTVaNK3xRsk/ZZUurwrIp7rxjG7JLdV3EDqHfQ5UqltEVLvqJOBzSVtEfkORKp23YPUW2odSXcwq/fXI6QeTe/VnON5Umm23vnrlRiWj4jnq95/DliQqqqvXCrftGa/2yLizmYfNy9Ln78jl4JuJF2LvSPi8Q7SDyIF8rWAYyPijk6ecmlSwFgS2D4iipasjga2AU4hfTHrdQ4qfUBuS7iFVG97M6ka4i3SP8t+pG/qPek/DdYPyct9O9h/ITrf5XaRvHy5wfaXSI3GiwDTJa1JumEuSPqWdy2p3SeAL5BuhpU2mEXzstHnanTOQiLiCUmfIpUE9iRfH0n3AkdExN3dOX4nnUHqHrx2RDyS170B/FjSMsAhpKqn63Pe/5vzfkJe/21SaXUE8BxwBTC55hynMev3VfEFYF1gZJ08vV7zfjippFJ949yS1EW22kxSO1gjr+Xlog22L1KTrpDc/nQT8DHgmxExqoP0A4CrgK2Bn0TED5ulb2AZ0ud4CPh70Z0iYkL+IvNFSVvT9Wq+0jio9A3HAANJ1S+zNThK2qsbx6305pnjgTfmvClUa/TNr3Jz2CIixnY1Ux0ce6kG25eqSXcY6TPsERG/q06Y2xLWq1pVuaks2cGxuywi/g7sLKnSMDwcOBj4k6S1IuLf3T1HQdsB/6kKKNXGkoLKx8lBBSAiXgG+lV8fkFQJELO1e0TEabUHlrQqKZCNaJa53Ea2AzAmqtp2IuI4oFOjQeSAOAVYRZKqSl8Vq+VlbVtRs/wtTKo+3YhUzfuLDtIPILUz7gKcHvUb1Yt4gFTS+QUwRtJWEVE0GB5H6uL8Q2DjLp6/NO791TesArxcJ6Asw6zG0q6olBbqVVV8vAvHq+Rvo65lp6kJeVlbBYLSk9GrA49Xqr5I1wxqGlBz1+P1q9flKq/nSNU7g2vSz0OJ/4gR8XZE3BkRRwI/JlVDbV7W8QsYBCyab3a1lsjLdzo6SO44sQepXeXa8rLHhqS2g+72+qq4g/Tlol4Pr61J7RP1eobNQdJCpG7Um5B6XZ3RQfp5SSW53YBzIqJbDyJGxPmkwD6MFFgalcBq93sCuAj4H1KnkV7loNI3PAsskdsIgA96Wp1NN0qTETGV1C10s+puybne/ZguHPJaUo+dY2q77ObjLph7v3TFHaSn1Xerc4yTSTfLy6vWPZuXn646/zykG3m952VGkarlvl+z/gC6156CpA0a9PKqlICm19nWLZIWzc/G1D55P45U7Xd0TfrFSB0KIF3r6m2L1LyfBziV1IvpzCj2RHpRw0kl6D+UdLyL8vKk6kCaS/jrAVdFxLSq9Qvm6zbb0Cu5G++NpC81x0bEqc1Omq/Rr0iB9wJStWG3RcS5pMCyAXBz7e+miRNIf2f1qh97lKu/+oazgS2AcZKuIv3TbUUK+v8EutMP/fT8ukfS1aSeN8NJD8Kt3GzHWhExXdLupOqBeyWNIfV2GQisRHpm4E5gx85mMiLek7RfPvbYfB1eJF2XDUnfNn9atct5pMbl0ZKuJFWLbUZqkL+TOUs8p5Dq/Q+TtD7pae41ge1J7Vmf62yeq+wNHCDpdma1h32C9Dt8hNRDpylJSwI/qbzNy6GSLs0/v1xTtbI78EvSTfWAqvXfI1Vzjcx17JWG+p1J1X+/iTkf2ntJUqUtbwDpG/7apN9FbTtHdw0HxkVEbTtNl0TEnyRdQXqCfnz+HMuTbvbPUxNcgU+Rft9/Jv1+Kn5JKlE+DQxq0BHltEgPkUK6eX+JNDzLS6Qu3LXpL46IZ2tXFvhM5+SgdSZwk6RtooPuyRHxoqQzga5Wv5Wnt7qdze0vOjlMC+kPdALphvQScDHpJnAXNV01KdiluGrbEaSnfGcAj5Pq1VejcZfiJzrI64qkBw+fIlWlTCU9UHgGMKzAZ204NAXpZnwd8Eo+9r/y512oTtrPkb6ZTyM9b/E7UrXYKGqGI8nplyI9C/MKqaH4L6SnxQsP00L9LsUbk76tTiS137yRfz4RGFLw9185bqNXbbfyA+r9/vK2dfI1eJ5UfTWNFJQPps6QK6Qb6uP5b+8NUsD9Rr20TfLfYZdi0hePIFUtlfm/Ni9wJCmAv03qeHEJsEydtJX/k1tr1le66jZ7VQ9vM6pA+iJ/T83+Fw7N2+5iVhfi2boU16QfTHpYuFe7FCtnxsyspSR9hzTywVox++CMNhdxUDGzHiHpz6QBEFfrMLG1LQcVMzMrjXt/mZlZaRxUzMysNL3SpVjS9qSuj58gdZ99nNQj5La8fTCpn/wupC6w44DDIuKfNceZn9TD5ivAYqTeU9+N5uMFzWbxxRePoUOHdvcjmZn1Gw888MCUiFii3rYeDyqSvkF6LuNsUkCYh/QE9IJ5e2WOhZVIDxS9Suprfruk9WP2gekuIg35cBSpe+vBpAeGNo6ICRQwdOhQxo/vtemczczajqSGU5n3aFCRNJT0IN5REXF61aabq34eThomYcuIuD3vNw6YBPwf6RkLJH2MNInNfhFxSV53B+n5gJHMmv/azMx6SE+3qexHqu46r0ma4cALlYACEGlgtRuYNfVpJd27pNFBK+lmAlcC22jOWQLNzKzFejqobEKaynRPSU9Kmqk0H/rBVWnWIQ3/XGsisIJmzTe9DjApImrnGplIGieqW+M5mZlZ5/V0UFmGNDzIqcCPSGMM3QKcLenQnGYIqR2l1tS8HFww3ZA624A0hayk8ZLGT55cyhBEZmZGzweVeUhDgX8jIn4ZEbdFxEGkCXGOzo30ov58HrWjtRVNN4eIuCAihkXEsCWWqNuBwczMuqCng8oreXlLzfoxpMH+liaVNOqVMiollErppKN0U+tsMzOzFurpoDKxwfpK6eL9nGadOmnWBp6NWXMjTARWypMy1aabQSdmezMzs3L0dFD5fV5uU7N+G+D5SJMBjQaWzfO2Ax9MIrQTs88WN5o0j8fuVekGAF8kTVXa4ex2ZmZWrp5++PGPwO3A+ZIWJz2wuBupwX7fnGY06Qn6UZKOYtbDj2LWJEZExIQ8kdPpkgaSnmM5iPTQZHfmdTczsy7q0aASESFpF9IsfCeQ2j8eBfaKiF/nNO9L2pE0y9+5wPykILNFRDxXc8h9SVPNnkQapuVBYNuI+FtPfJ6h3+v6jKhP/2iHEnNiZtY39PgwLZGm4zw4vxqlmUp6UHK/Do41nTTv9uHN0pmZWc/wHPVtyCUkM+urPPS9mZmVxkHFzMxK46BiZmalcVAxM7PSOKiYmVlpHFTMzKw0DipmZlYaBxUzMyuNg4qZmZXGQcXMzErjoGJmZqVxUDEzs9I4qJiZWWkcVMzMrDQe+t46xcPum1kzLqmYmVlpXFLpJd35xt8fuYRk1h5cUjEzs9I4qJiZWWkcVMzMrDQOKmZmVhoHFTMzK417f/Uz7nVmZq3koGI9xgHNbO7n6i8zMyuNSypmfZQf+LR25KBic73uVrv5Bm1WnIOKmc3GQdi6w20qZmZWGpdUzKxUbgvq3xxUzKzP6K1u5w5m5XFQMWshP5tj/Y3bVMzMrDQOKmZmVhoHFTMzK43bVMys3/OzOeVxUDGbC7mDgPUWBxWzDvgGbVZcoaAiaWdgSERckt+vCFwJrAvcDHwtIqa1LJdmZjaHvvigadGG+uOAJarenwYsB1wAbAqMKDdbZmbWjooGlVWAfwBIWgDYHjg8Io4AjgE+35rsmZlZOykaVOYHpuefP0WqNhuT3z8GLFNyvszMrA0VDSpPA5vkn3cGHoiI1/L7JYHX6u1kZmb9S9HeX+cDP5X0eWB94KCqbRsDD5edMTMzaz+FgkpEnCFpCrARcGZEXF61+UPApS3Im5lZW+iLvbB6S9EuxSsAv42IK+ps/hawdKm5MjOztlS0TWUS8PEG2z6at5uZWT9XNKioybaBwPsl5MXMzNpcw+ovSYsBQ6pWLStp5ZpkCwD7AC+1IG9mZnO9uW0YoGZtKocCPwAiv65ukE45nZmZ9XPNgsp1pOdTBFwMnAQ8WZPmHeDhiPhHS3JnZmZtpWFQiYgHgQcBJAVwY0S80lMZMzOz9lP0OZXLWp0RMzNrf4XnU5G0D/AlYAXSWGDVIiJWKTNjZmbWfoo+/Ph94ATgIWACqS3FzMxsNkVLKvsDZ0TEYa3MjJmZtbeiDz9+GLihlRkxM7P2VzSo3AF8rBUZkHSTpJB0Us36wZIulDRF0puSbpW0Xp3955d0qqQXJU2XNE7Spq3Iq5mZNVc0qHwH2FfS3pIWlzRP7asrJ5f0JeoEK0kCRgPbAt8GdiUNB3O7pOVqkl8EfB04HtgReBG4WdL6XcmTmZl1XdE2lcfz8pIG26MTxwI+GAbm58BhwK9rNg8nTQq2ZUTcntOPIw1c+X/AIXndx4AvA/tFxCV53R3ARGBkPo6ZmfWQooFgJClwlOknwMSI+I2kekHlhUpAAYiI1yTdQJp58pCqdO8CV1WlmynpSuB7kuaLCPdUMzPrIUUffhxR5kklbQLsTeN2mnVI3ZdrTQT2lrRwREzL6SZFxFt10g0CVs0/m5lZD+hSW0h3SBpInp44Ih5rkGwI8Gqd9VPzcnDBdEPqbEPSgZLGSxo/efLkYhk3M7MOFX348fgOkkREnFjwnN8lDZl/crNTUr+6rXZel6LpZhMRFwAXAAwbNqzsaj0zs36raJvKiCbbKjflDoNKnpb4WOAAYD5J81Vtni833r9BKmnUK2VUSiiV0slU0rAxjdJNrbPNzMxapFD1V0TMU/siPRD5NVLbx6oFz7cyadywUaTAUHkBHJl/Xo/UDrJOnf3XBp7N7SnkdCtJWrBOuhnAEwXzZWZmJehym0pEvBoRlwOXAucU3G0CsEWdF6RAswUpEIwmzTS5WWVHSYsAO+VtFaNJz6/sXpVuAPBFYIx7fpmZ9axOPVvSwIMUqPoCiIj/AmNr16dnHXkmIsbm96OBccAoSUeRSjBHk9pKflJ1vAmSrgJOzx0AJgEHASsBe3X5E5mZWZeU0ftrR6DULlQR8X4+7i3AucDvgfeALSLiuZrk+5IeyjwJ+AOwPLBtRPytzDyZmVnHivb+urjO6kHAuqQ2kG7NUR8Rc/TWioipwH751Wzf6cDh+WVmZr2oaPXXlszZdfdt4BngdMAzQ5qZWeEn6oe2OB9mZjYX6PEn6s3MbO5VOKhIWlrSTyXdL+lJSfdJ+omkj7Qyg2Zm1j4KBRVJq5OeMTkEmAbcB7wJHApMkLRay3JoZmZto2hD/Y+B14ENI+LpykpJKwJj8vYvlJ47MzNrK0Wrv7YAvl8dUAAi4hnSuGBb1NnHzMz6maJBZRBpoMd63sjbzcysnysaVCYA366diz7PJf/NvN3MzPq5zkwnfCPwSB5r60XgI6SBHFcDdmhN9szMrJ0UffjxJkk7ksbXOpZZk2M9AOwYEWNal0UzM2sXhUcpjoibgJvy3CWDgVfrzA1vZmb9WNEBJQcCgyLizRxI3qrathAwIyLebVEezcysTRQtqVxImgzry3W2nU+aZbHpaMJmZjb368xzKtc32DYa+Gw52TEzs3ZWNKgsCfynwbbJwFLlZMfMzNpZ0aDyH9JkXPWsB7xSTnbMzKydFQ0qNwLfl/TR6pWS1iN1Mb6h7IyZmVn7KdpQfzzwOeABSfcDzwPLAhsAk4DjWpM9MzNrJ4VKKhExBfgkcArpwcf18/Jk4JN5u5mZ9XOdefjxv6QSy/Gty46ZmbUzTydsZmalcVAxM7PSOKiYmVlpHFTMzKw0DipmZlaaTgUVSfNIWlfSZnl0YjMzsw8UDiqSDgZeAh4EbgPWyOuvk3RIa7JnZmbtpFBQkfR14AzgOuCLpAcfK/4C7Fp+1szMrN0ULakcDvwsIg4Efl+z7VFyqcXMzPq3okFlJeDmBtveBBYrJztmZtbOigaVKcDQBtvWAP5dSm7MzKytFQ0qNwDHS1q5al1IWhw4jNTWYmZm/VzRoHIc8A7wEHArEMCZwCPAe8DIluTOzMzaStGh718BhpGGvh8IPEka4fhsYOOIeK1lOTQzs7bRmaHv3wBOzC8zM7M5FH1O5TRJ67c6M2Zm1t6KtqnsS5pK+CFJR0latpWZMjOz9lQ0qCwF7AE8Qar+ekbSrZK+6jHAzMysomhD/YyIuCYidgGWBg4BFgAuA16W9KsW5tHMzNpEp4e+j4hXI+LciPg0sAXwKvDl0nNmZmZtp3Dvr4pc3bUb8BVgc2AmcE252TIzs3ZUtPfXPJK2lXQF8DJwMTAf8E3gIxGxRwvzaGZmbaJoSeUFYAlSQ/2PgV9FxNOtypSZmbWnokHlGuDyiLi3lZkxM7P2ViioRMTBrc6ImZm1v4ZBRdKmwN8iYlr+uamIuLPUnJmZWdtpVlIZC2wE3Jd/jgbplLfNW2bGzMys/TQLKlsAD+eft6RxUDEzMwOaBJWIuKPq57E9khszM2trRZ9TeUrSxxpsW1fSU+Vmy8zM2lHRYVqGkh52rGd+YMVScmNmZm2tM2N/NWpTGQb8t4S8mJlZm2vWpfgw4LD8NoAbJM2oSbYAMAS4sjXZMzOzdtKs99dTwJ/zz/sA44HJNWneIfUQu7D8rJmZWbtp1vvreuB6AEkAIyNiUg/ly8zM2lDRYVr2bXVGzMys/RWeT0XSIGA7YA1Sj69qEREnlpkxMzNrP4WCiqRlgLtIXYuDNDQLzN4jrMOgImk34EukHmNLAs8C1wI/jIg3qtINBk4FdiF1BhgHHBYR/6w53vz5vF8BFgMmAN/1OGRmZr2jaJfiU0mN9CuQAsqGwMrAyaQ5VlYueJwjgfeAY4BtgV8ABwG3SJoHQKkBZ3Te/m1gV2AgcLuk5WqOdxHwdeB4YEfgReBmSesXzI+ZmZWoaPXXZ0gB4YX8/v08SdfxkuYFzgR2LnCcnSKiugfZHZKmApeRpia+DRgObAJsGRG3A0gaB0wC/g84JK/7GPBlYL+IuCSvuwOYCIzMxzEzsx5UtKTyYeCFiHgfeBMYXLXtNlJA6FBNQKm4Py+Xzcvh+Vy3V+33GnADsweu4cC7wFVV6WaSnpnZRlKjEQDMzKxFigaV54HF889PAltXbdsAeLsbedgsLx/Jy3WAh+qkmwisIGnhqnSTIuKtOukGAat2I09mZtYFRau/bifd/K8DzgfOye0W7wLb5HWdJmlZUlXVrRExPq8eAjxdJ/nUvBwMTMvpXm2SbkiT8x4IHAiwwgordDrfZmZWX9Ggchz5Jh0Rv5A0APgisCDwE1Jg6JRc4rgemAlUPwdTmfRrjl3qvC+Sbg4RcQFwAcCwYcM8T4yZWUmKPvw4BZhS9f4s4KyunjR3BR5N6jW2WUQ8X7V5KvVLGZV2nFer0tUrZgyu2m5mZj2oM6MUl0LSQOAaUlvM9rXPnpDaRNaps+vawLMRMa0q3UqSFqyTbgapq7OZmfWgZqMUX9yJ40RE7N9RovwsyhXAZ4EdIuKeOslGA/tK2qwy+6SkRYCdgF/XpDsB2J3UJZmqarkxEfFOJ/JvZmYlaFb91Zl56YumO4cUBE4G3pS0UdW253M12GjSE/SjJB1Fqu46mtRW8pMPThgxQdJVwOm59DOJ9CDlSsBeBfNjZmYlajZK8dAWnG+7vDw2v6qdAIyIiPcl7Qj8FDiXNM7YOGCLiHiuZp99SQHqJNIwLQ8C20bE31qQdzMz60DhASXLUDRQRcRUYL/8apZuOnB4fpmZWS8r3FAvaSFJh0i6WtLtklbL6/eUtGbrsmhmZu2i6CjFywNjgeWAR4F1gQ/lzVsAWwEHtCB/ZmbWRoqWVH5Gmjp4NeB/mP0BwzuATUvOl5mZtaGibSqfAw6MiGfzqMTV/s2swSDNzKwfK1pSGQS80WDboqQxwMzMrJ8rGlT+QZosq57tgAfKyY6ZmbWzotVfpwJXp0kZP3iqfW1JOwP74wmxzMyM4gNKXivpm8CPmPXsyOWkKrFvRcRNLcqfmZm1kcIPP0bEeZJ+BWwMLAm8Avw1Ihq1tZiZWT/TYVCRNIg0Ze/PI+JO4NaW58rMzNpShw31ETGD9HBjjw+Tb2Zm7aVooLgb2KjDVGZm1q8VbVM5ArhO0jTSPPUvUjPcfUS8X3LezMyszRQtqfwTWAU4A3iGNLPiu1WvGS3JnZmZtZWiJZWRFJ+Iy8zM+qmiz6mMaHE+zMxsLuAeXWZmVhoHFTMzK42DipmZlcZBxczMSuOgYmZmpXFQMTOz0jiomJlZaRxUzMysNA4qZmZWGgcVMzMrjYOKmZmVxkHFzMxK46BiZmalcVAxM7PSOKiYmVlpHFTMzKw0DipmZlYaBxUzMyuNg4qZmZXGQcXMzErjoGJmZqVxUDEzs9I4qJiZWWkcVMzMrDQOKmZmVhoHFTMzK42DipmZlcZBxczMSuOgYmZmpXFQMTOz0jiomJlZaRxUzMysNA4qZmZWGgcVMzMrjYOKmZmVxkHFzMxK46BiZmalcVAxM7PSOKiYmVlpHFTMzKw0DipmZlYaBxUzMytN2wcVSctLulrSa5Jel3StpBV6O19mZv1RWwcVSQsCtwFrAvsAXwVWA26XtFBv5s3MrD8a0NsZ6KavAysDa0TEEwCS/gH8C/gGcFov5s3MrN9p65IKMBy4pxJQACJiEnA3sHOv5crMrJ9q96CyDvBQnSnH/L4AAAtuSURBVPUTgbV7OC9mZv1eu1d/DQFerbN+KjC40U6SDgQOzG+nSXqsBXlrtcWBKb2diTbna9h9vobd1yvXUD/u1u4rNtrQ7kEFIOqsU9MdIi4ALmhNdnqGpPERMay389HOfA27z9ew++a2a9ju1V+vkkortQZTvwRjZmYt1O5BZSKpXaXW2sDDPZwXM7N+r92DymhgI0krV1ZIGgp8Om+bm7V19V0f4WvYfb6G3TdXXUNF1GuSaA/5AccHgenAcaT2lROBDwEfjYhpvZg9M7N+p61LKhHxJrAl8DjwK+AKYBKwpQOKmVnPa+uSipmZ9S1tXVKZ20haTtJZksZJektS5DaiovsvK+liSS9JekfSJEmntC7HfU93rqGkD0s6Q9JTkqbn63e2pCVam+u+Q9Jukq6R9Ey+Bo9JOkXShwrsO7+kUyW9mPcdJ2nTnsh3X9LVayhpmKQLJD2a/3aflXSFpJV6Ku9lcFDpW1YF9iB1h/5LZ3bMN877gNWBQ4CtgRHAzDIz2Aa6dA0lidS548vAqcB2efklYHTe3h8cCbwHHANsC/wCOAi4RVJH94uLSOPxHQ/sCLwI3Cxp/dZlt0/q6jXck9Sb9UzS39/3gE8A4yUt39Iclyki/OojL2Ceqp8PIHU8GFpw35tIQWVgb3+OdryGpGAcwIE16/83r1+jtz9bD12/Jeqs2ztfgy2b7PexnGbfqnUDgMeA0b39udrkGtbbb0XgfWBkb3+uoi+XVPqQiHi/K/tJWgXYBjgrIt4tN1ftpavXEBiUl6/XrP9vXvaL/5WImFxn9f15uWyTXYcD7wJXVR1rJnAlsI2k+UrLZB/X1WtYb7+IeAaY3Gy/vqZf/KP0A5/Oy+mSbsntKa9KulzSh3s1Z+1jInAn8P1ct72wpA1IVTl/iohHejd7vWqzvGx2DdYBJkXEWzXrJ5IC9qqtyFgbKXIN5yBpLWDJzu7XmxxU5g7L5OXFpO7V2wHfBXYg1Wn799yBSHUN25Oqa+4H3gDuBZ4Cdu3FrPUqScsCI4FbI2J8k6TNBnetbO+XOnENa/cbAJxHKqlc1KLslW5uGFDSZn05GBsRB+efb5P0Grn6AfhTr+SsvfwS2IjUjvIIsBZwAnC1pJ26UbXWliQtDFxP6uyxb0fJ6cLgrnO7Tl7DWmcDnwJ2iIi2GcvQQWXu8Epe3lKzfkxefhwHlaYk7UDq6bVVRPw5r75T0lOk67gT6ebQL0ian9QbbmVgs4h4voNdpgIr1Fk/uGp7v9KFa1i97ymk6Tn2iYgxHaXvS1wtMneYmJeNnmTtV9+wu2i9vLy/Zv19eblWD+alV0kaCFwDbABsHxH/LLDbRGAlSQvWrF8bmAE8Mecuc68uXsPKvseSuhMfGhG/alEWW8ZBZe5wD/ASqU98tcr72hulzemlvNygZv2GefnvHsxLr8ntb1cAnwV2joh7Cu46GhgI7F51rAHAF4ExEfFO2Xntq7pxDZF0CHAScGxEnNWiLLaUq7/6GEm75R//Jy+3kzQZmBwRd+Q0M4HLImJ/SF03JX0PuFTSecC1pN42JwNjgdt68CP0uq5cQ9I1Oxm4XNKJwKPAmsAPgOeA3/dU/nvZOaTAcDLwpqSNqrY9HxHPS1oReJL07MRIgIiYIOkq4PT8LX0S6YG/lYC9evQT9L4uXUNJewKnk545u61mv9cjoj2m8+jtB2X8mv1FqsKq9xpbk+bSOvt+FXgIeIf0NPNZwMK9/Zna5RoCy5N62UwC3s7LXwLL9vZn6sFr93ST6zcipxla/b5q3wWA00ilvrdJvec27+3P1C7XELi0yN9uX395QEkzMyuN21TMzKw0DipmZlYaBxUzMyuNg4qZmZXGQcXMzErjoGJmZqVxULE+TdJ+kv4laYak/3a8x2z7hqQRVe9HSOpTfehr89hD5xyaz/u1gunPknRD1fvN8/5bdbDf2Jzu8jrbDqid6lnS9ZLOKfxBrE9yULE+S9IywAXAX4EtgaY3sTo2Bi4sO1/9SZ4A7huk0Zq7ai9JaxdINwL4uqTVu3Eu62UOKtaXrQbMSxpO5a7oxFwUABFxT3RiZNjO6iezGX4HeLCz177KBGAKcGJHCSPi7zn9d7p4LusDHFSsT5J0KWncMoA/56qSS/O2PSXdJmmypGmS/i5pnzrHaFq11KgaqKp6Z/OqdWMl3SVpp3y+d4Bv5m0DJB0t6dE86+YLkn6Whz7vymffVtI4SdMlvSbpOklr1KTZWtIfJb0o6S1JD0k6QtK8NekWlHSupFfytRoNLFcwH/MBXwF+XSDtyrma8m5Jg6s2vQn8EPiCpGEFTnslqWSzQJE8Wt/joGJ91YnAIfnng0lVWZVvuysDV5MGKtwFuAG4UNL/tjhPqwNnksZU2waozLsyCjiOdPPdATgF2J80Um2nSNoW+AMwjTTC70HAusBdeQbBipXz+ffL57yMVH10cs0hzwcOII3J9QXSzJYdBolsI2Ax4C8d5PnjpCrKR0jz0dROKHUe8Cxp9N2O3AksQvp9WxvyKMXWJ0XEk5Iq83I/HFXDh0fEDys/52HGxwJLk27A57UwW4sDW0fEhKrzf4Z0898nIioN0rdKmgqMkrR+dfoCTiJNYbxdRMzM5xhHmib6COBwgIj44HNKEunGPwg4UtIxEfF+Lt18mTSM+o9y8jFKsxEWCcAbkQYz/EejBJI+SxrB+XfAgRHxXm2aiHhH0khS4N80Iu5scs4HSfP/bEQ/G117buGSirUdSatJ+o2kfwPv5tcBwBrN9+y2p+sEiG1Jk1Bdk6vBBuR5RCqz9W1a9OCSFgI+AVxVCSgAETEJuBvYrCrt0pLOl/RMPv+7pIC0GLBkTrYh6X/8tzWnurJglpYhDbk+o8H23YE/AudExP71AkqVS0mBsbYkNZuIeBd4LZ/b2pBLKtZW8rfsW4C3SLPjPUm6qR5EqgpqpRfrrFuSVEKY1mCfD3fi+INJc7rXO89LwIrwQelsNOnGO4I098t0UlXgsUClLWfpvHy55li17xuZnzSNQiO75vNe0tGBIuI9SccDV0raroPk00nD6FsbclCxdrMx6eb6mYi4q7Iylw466+28HFSzvlEgqPeMyyv5OJ9psM8LncjPq/kcH6mz7SP5XACrAMOAr0bEqEoCSTvV7FMJTkuRqtSoel/EK8yaY76eA4EjgbGStoyIRzs43m9JXwROonk15RBSjzFrQ67+snZTmQP93cqK3Nto5y4c62XSN/F1a9bv0Ilj3ET6Rr9oRIyv8yocVCLiTeABYPfqXlx5lsBPAXfkVfWuwUDmnGHxXlL7xB416/csmKVHgYGSGvUWe53UYeEp4HZJazU7WKTJm44jVfHtWi+NpI+QrudjBfNofYxLKtZu/kq6mZ0j6QfAQqQb1RRg0c4cKCJCaQrc/SU9TrqR7QBs3oljjJX0G+BqSacB95Fu5EOB7YHvRsTjncjW90m9v26UdC6wMOnBw9eAn+U0jwDPACdLeo8UXA6rk7fHJP0aGJmrzO4HPpfzVUSlQX0DoO7zPhHxRlWPtdslfTYiJjY6YET8QdLdpGBUz4Y157Y245KKtZWImAx8nvRQ5NWk7rsXkrr1dsWhpPnpRwBXkb4lf7uTx/hK3n834Pqcr28B/6J4+wUAEXETKbAtRqouOo8URDaplHpyw/kupHaWy0lzot8J/KjOIb9BmiL5SFIvrTVJPcKK5OVpUpCsrVarTTeNFKgeJs2tXlvyq3Vsk207Ag9ExBNF8mh9j6cTNrOG8oOhZwBLR8RbLT7X/KR2oCMj4qJWnstaxyUVM2vmV8C/yaMHtNg3gP+QHuS0NuWgYmYN5WdP9iN14W61d4CvVT+jY+3H1V9mZlYal1TMzKw0DipmZlYaBxUzMyuNg4qZmZXGQcXMzErz/z00/+E2CiW6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N=10000;\n",
    "r=np.random.rand(N,1);\n",
    "wmean=1; # in mm\n",
    "wmin=wmean-wmean*0.1;\n",
    "wmax=wmean+wmean*0.1;\n",
    "hmean=2; # in mm\n",
    "hmin=hmean-hmean*0.1;\n",
    "hmax=hmean+hmean*0.1;\n",
    "\n",
    "wrand=wmin+(wmax-wmin)*r;\n",
    "hrand=hmin+(hmax-hmin)*r;\n",
    "\n",
    "uts=940; # in N/mm^2=MPa\n",
    "\n",
    "Ffail=uts*wrand*hrand*1e-3; # force in kN\n",
    "plt.hist(Ffail,bins=20,)\n",
    "plt.xlabel('failure load (kN)')\n",
    "plt.ylabel('relative counts')\n",
    "plt.title('Failure load is {:.2f}+/- {:.2f} kN'.format(np.mean(Ffail),np.std(Ffail)));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Normally, the tolerance is not a maximum/minimum specification, but instead a normal distribution that describes the standard deviation, or the 68% confidence interval.\n",
    "\n",
    "So instead, we should generate normally distributed dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAESCAYAAADuVeJ5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAf60lEQVR4nO3debgcVbnv8e8PwygCiUwCYpAgnORBUQPCFSVwkVGDyiQCikxXjqJHgSMBB2QQHEHhoCAIajyCUYQgKIiEMBjAgIAEBIOJyBwIUwyE6b1/rLVN0enee/W0eyf5fZ6nn9pd9Vb16trd/dYaqkoRgZmZ2UCW6XUBzMxs8eCEYWZmRZwwzMysiBOGmZkVccIwM7MiThhmZlbECWOIkXSipJC0dWXe9nneF3tZthKShuWyXtXrskD9/dlP7Kgce85glM2WbEPtu9AJThgtkDQyfxAaPZ7qdRlt6MiJ6KuSfiPp4fwZmd7itkZIOlnSnZLmSXpU0g2SPibpNXXiXyvpC5Juk/S0pCckTZf0GUkrtP/uGpbz/fl9vr/N7ayYk/7fJS2QdL+k0ySt2sQ23iHpa5JulvS4pOcl3SvpW5JG1InvO8jo7zGhnffVoJzX520/I2n1Osv7ElBLn51OGNarF15C3ANcUGf+821s8zRgIvCPNrZhrfkH8B9ApxP+OODLwEvA3cDarWxE0uuBW4A3AVcDlwMrA7sBPwa2AQ6qxC8HXAu8A7gZOAdYFtgJ+C7wAUk7RHfO3h0PzAf+0OoGJA0jvcdxwPXAL4DRwGeBbSRtHRH/KtjU2aR9MI303Xolb/MI4MOStoyIxyrxV5P+V/V8HngdcGWz76cJrwOOBT7XxddoTUT40eQDGAkEcPEgvd72+fW+2Ov3XlDWYbmsV/W6LLk8J+bybN3DMmwIvAtYobJ/prewnWPyusfXzF8F+Htetl5l/kfzvJ/UxC9L+vEM4N0Fr3s9MLOJcgp4qN3vB/DJXMbzauYfm+d/pXA7hwMj65Tx+3k73yvczltz/O2F8U19F/J+DuA+0kHn+g221/Rnp1MPN0l1maTlJX1W0lWSHpT0gqSHJE2UtGGd+KI29/7a2xstk/SApJmSXi/p7Nw88kpNf8lGkn5cKes/JZ1Rr4rcLEmbSbooNwssyM0Cx0taqU7swbkJ5x85do6kiyW9o8G215R0Xt72PEnXSnpPk+VrtN82kfSzSlkelTRN0hEl242I+yLipohop+YJ8OY8vbxm+88A1+WnqxfEv8jCI+Q12ixTPZsDbwAmt7mdg0m1gWNq5n8TeAI4SJIG2khEnB4Rs2vmBfC1/PS9heU5ME9/VBhfl6Rl8+cpJH29TsiXgeWBr7bzOt3ghNF9awDfJh3RTAZOBf4EfAS4WdKbBrk8KwDXAO8GJgFnAc8C5MRxSy7bDaTmsduA/wT+KGm1Vl9U0jbAH4FdST9gpwJPAl8C/iBp+ZpVzgRGkH7YTiU1E+wE3CBpbM22X0dqejkAuAM4HXgc+D2wRatlztt+I3ATqdnneuA7wK+Bl6k0/wySu/J0p+rM/P63Bh4G/loQPwx4H6nJ6KYulHM86Yf+N61uQNJrSc1IMyLi4eqyiHiB9Bl+I7BB68XkxTxt1PxULc9ywL7AC6RmrZZIWhG4mFT7+++I+EKdsEtINcD9JY1u9bW6wX0Y7dlE0nF15l8QEX1f3MeBN9Z+6CVtR/pBm0Cqeg+WdYHpwF75i9dXnuWBn5O+RJtHxD2VZfsA/wscB/xXsy+o1Bl7PumoaZuIuD7PF3Ae8HFSe/LXKqttHBGzarYzhvQDdyKv/hGcAGwMnBoRn6/EH0pKiO3Yk9Tks2tEvOpIPfcpDKazgf2Br0h6N3ArqQ9jPCmB7V5Ti7kYuAL4uKSNSQcBfX0YI4B9aj+XHTIeuCle3S/QrI1IB1kzGyz/W56OIjXHteLjefr7gtjxpNrbpIh4opUXywdcvwG2BA6KiP5qKhNISfEk4EOtvF5X9KotbHF+sLAPo9Hjg4XbuYuatmHqtLlTpw+D9EUJ4Jw62627DHggz9+4zjp75WX/1aCstwGPFLynRdptge3yvF/UiV8HWADcU7jPLgeeA4ZV5t0PzAOG18QuQ/phKerDqLffSJ2cAWzboc9OW+3QpAQxqebz9jwpma/Y4PXOqIl/mdTpPaLwNYv7MEgd8gEc3eZ+em/ezvkNlk/Iy/dscftvA/4FPAqsXvi5C2CnFv7XV5EGOtye/1cf6mc/B7Byfv7b/PxdnfjsdOLhGkZ7LomIDw4UJGlz4ChSM9CavLpmN79LZWtkXlRqDxXvytO3Nqg1LQ+sJWm1iGh2FNFmeXpN7YKIeEjS34AxklaMiOcAJL2F1HY9jtQevlzNqsOBOZKGk5ombo6IJ2u2/YqkaaRE0KrJpKO8SyX9gnQ0en1E/LONbbZE0lrApaRRNO8j1bZWAfbLZRwnadvIvy65WWcS6cdxd1Kz3nKk5rVTgR0kjY080ig3Vb1IA5JqR1O9HBG1vyG75em/+y8kfZjUYVx1UUTc0d/bzdOOj+DKzYyXAq8BPhIRjw8Qvw6wA+mAq5XRUW8gJYM1gV0i4urC9SYAOwInkw66es4Jo8ty2/3vSe2kV5Cq2PNJX4QDSUfYg6lRM0HfePRPDLD+a2l+2Okqefpog+WPAGNy3HOSNiH9GK5EOjq7iNTPEsCHgU1JCQygbzx+o/fV6DWLRMRMSf+HdAT/EfL+kXQTcERE3NDO9pv0XVKH8uiIuDvPexb4ev5R+wyp6eSSvOxLwM6kH6nfVrbzQ0krk/pjDs7bhdTvUK+j9UBSzeZ7NfNfqRM7HrgvIu6qzPswqf2/aiapv6mRp/O00fkWq9TEFZH0BtJQ37WBPSJiSsFqB5CSy48jot57Hsg6pPdxJ/Dn0pUi4rZ8kLK3pB1ICb+nnDC67xhSu/HWEXFzdYGk2i9RM/o+uIucrMXCL1M9jY7YnsnTbSPimlYLNcC212qwfK2auM+R3sNeETGpGpjb7jetzOr7wVhzgG23LCL+DOymdKLbFqQfxU8Bv5X0HxHxYLuvUWhn4LFKsqi6hpQw3s7ChLFzZVm9eHI8kGpkpMT4KpK2B9aOiEWW1cStSmpKOqM6PyL2I9WCmjGT9FltVDvcqBJXRNLawBRSR/lHIqJ0FNcBuSznlb5WjVtINb3vA1dK2j4iShPdF0m1w68BW7X4+h3jUVLdtyHwaJ1ksQ6pL6RVfUf569ZZ9vY68wbSV74tWytOv27L00WGL+YjvrcA9/Y1R5H2GaRmg2rsSixs3gIgN0P9k9SkNbwmfhk6+CWLiOcj4tqIOBL4OqlpaFyntl9gOWDV3HRUq2947ILKvL5aWL0h0fXi27Uz6eCo3eG0RMQ8Uqf+mPwZ+bc8YmkcqYlo1qJrL0rSmqQj9FHA/hHxq8L13ktKTlMj4r7iN1AjIs4CPg2MJSWNojPVI2ImcC7wTtIAjJ5ywui++4E1cps88O8RSWfQRg0vIuYCs0lnvP57aG5u564dt17iItLJVsfUDlvN211J0rsWXa3IVNJZ1HvU2cZJpB/Cn1Tm3Z+n7668/jKkH+l6P34TSU1lX6qZfzDt9V8gaYsGo6H6ai7P1VnWFkmr5nM/as8In0ZKAhNq4lcjdc5D2td9/pinX66er5BrSsfUiW/XeNJQ6es7tL1zSb9RX6uZfxTweuDcvv4aaLzflM4hupo0ku6AiKh3dYZGOnLuBUBEnElKGlsAV0jqryWg6qukz9nx7ZahXW6S6r4zgG2BaZIuJDUlbU/6IvyFdKmDVp2WHzdK+iWwIulLex0LT9oqEhHPSdqTNBrkJklXkkZxLUuqwm9DOteh6WsDRcTLkg7M274m74eHSfvlXaTazbcqq/wA+BgwWdIFpKaqbUid29eyaE3lZFI7+eckbQbcCGwC7ELqP3pfs2Wu+BhwsKQpLOx/egfpf3g3aSRLv/LR7Tf6nubpSEnn578fjVePx98T+CHpB/PgyvyjSU1Jx+c27b5O791ITXI/j4g/VuJPAD6QtzFW0lRSct6ZVLu9jnS5jbblWs9OwOURMeB5DYV+COwNHCBpFCkRjSG9p9tJJ/BVNdpvF+f1bgNG1RnU8UpELPJjnPt59iB9/n7Z7psBiIj/yQc/3wN+J2nHiHh2gHUelvQ9oN45G4OrV8OzFucHTV4aBNiH9GGdT+rg/RHpC3498FJNbNGw2sqyI0jj0F8A7iW1Y29E42G1/Q6PJA2LPDNvcwEwl9Q5+V1gbMF7bXg5BNIP7cWks3QXkIa8ngi8tk7s+0hH1PNI57JMIjVVTaTmEhg5fi3SuR5PkIZLXge8p97+7Kfs9YbVbkU6/2EGqb/k2fz3CZQPS+3bbqNH7dDqg+v9//KyMXkfPEAa0TSPlHA/BSxTJ/6NpATc9xmZn/+fxwIrFJZ/wGG1LBw6vVeHv2srkWqhs/Jn5oH8WVy1Tmzd/cbC4eSNHi81eO2D8vKzWix7f9+Fz+Zl17NwGO2rhtXWxA8n1d56OqxWuTBmZi2TdBrpigCrR7pUiS2BnDDMrG2S7iMNp92h12Wx7nHCMDOzIh4lZWZmRZwwzMysyBI9rHb11VePkSNH9roYZmaLlVtuueXxiFjkXilLdMIYOXIk06f37Pa3ZmaLJUl1bxHtJikzMyvihGFmZkWcMMzMrIgThpmZFXHCMDOzIkUJQ9J6kk6XNE3SfEkhaWSduGjw2KwmbhlJEyTNlvS8pNsl7d7gtQ+R9FdJCyTdI+mTrbxRMzNrT2kNYxSwF+lqidcNEHs+6Qqf1ce9NTEnkO7sdQbpUss3ApMk7VINknQIcBbwK9KlkycBZ0o6rLDcZmbWIaXnYVwbEWsBSDqYdEP0Rh6MiBsbLcz3BjgSOCUi+u6BMCVf7/4U0j0T+q6vfxLw04g4thK3DnCCpHMiouEN683MrLOKahjR2o3PG9mRdBOXiTXzJwKbStogP9+KdBvJ2rifku62tXUHy2RmZgPoxpneh0k6CniZ1NT0lYioNmONId0Ipfbm7TPydDTpZilj8vM7+4mb0qlC29Jp5NGXtbzu7FN27WBJzIa+To+Smki6icr2wKGkmsDVksZVYkYAT8Wi11WfW1lenT45QNyrSDpU0nRJ0+fMmdP8OzAzs7o6WsOIiP0rT6+TdAmphnAiC5uQRLrNYC01eN7UDTsi4mzSLTUZO3asb/ZhZtYhXT0PI9LNzS8DNq/MngsMl1SbIIZXllentTWJETXLzcxsEAzGiXu1NYoZwPLAhjVxo/P0rkocLOzLaBRnZmaDoKsJQ9IqwK7ATZXZvwNeAPatCd8PuDMiZuXn04DHG8TNBW7oeIHNzKyh4j4MSXvkP9+ZpztLmgPMiYipko4ENiaNXHoIeBPpfIu1qfzoR8Rjkk4FJkh6FrgV2BvYDtitEveipC+RTtR7ELgqxxwIHB4RL7Tyhs3MrDXNdHpPqnl+Zp5OBcYB9wAfyo9VgWdItYCDIuLmmnWPBeYBnyUllHuAvSLi0mpQRPxAUgBHAEcB9wOfjogzMTOzQVWcMCKitpO6dvmlwKX9xVRiXyaNnDqxIPYs0uVBzMysh5boW7SadZNP+rOljS9vbmZmRZwwzMysiBOGmZkVccIwM7MiThhmZlbECcPMzIo4YZiZWREnDDMzK+KEYWZmRZwwzMysiBOGmZkVccIwM7MivvigLfbauQigmZVzDcPMzIo4YZiZWREnDDMzK+KEYWZmRZwwzMysiBOGmZkVccIwM7MiThhmZlbECcPMzIo4YZiZWREnDDMzK+KEYWZmRZwwzMysiBOGmZkVccIwM7MiThhmZlbECcPMzIo4YZiZWREnDDMzK+KEYWZmRZwwzMysiBOGmZkVccIwM7MiThhmZlbECcPMzIo4YZiZWREnDDMzKzKs1wUwWxqNPPqyltedfcquHSyJWTnXMMzMrIgThpmZFXHCMDOzIk4YZmZWxJ3e1nPtdACb2eBxDcPMzIo4YZiZWREnDDMzK+KEYWZmRZwwzMysiBOGmZkVccIwM7MiThhmZlakKGFIWk/S6ZKmSZovKSSNrBO3gqRvSnpY0nM5/r114paRNEHSbEnPS7pd0u4NXvsQSX+VtEDSPZI+2eybNDOz9pXWMEYBewFPAtf1E3cucAjwZeD9wMPAFZI2q4k7ATgOOAPYGbgRmCRpl2qQpEOAs4BfATsBk4AzJR1WWG4zM+uQ0kuDXBsRawFIOhjYoTZA0tuAjwIHRsR5ed5UYAZwPDA+z1sTOBI4JSK+lVefImkUcApweY4bBpwE/DQijq3ErQOcIOmciHix2TdsZmatKaphRMQrBWHjgReBCyvrvQRcAOwoafk8e0dgOWBizfoTgU0lbZCfbwWsUSfup8Drga1Lym5mZp3RyU7vMcCsiJhfM38GKUGMqsQtAGbWiQMYXYkDuHOAODMzGwSdTBgjSH0cteZWlvdNn4qIKIijzjZr415F0qGSpkuaPmfOnKKCm5nZwDqZMATUJoG++a3G0SC2oYg4OyLGRsTYNdZYo5lVzcysH51MGHOpf9Q/vLK8bzpcUm2CqBdHnW2OqFluZmaDoJMJYwawgaSVauaPBl5gYZ/FDGB5YMM6cQB3VeJgYV9GozgzMxsEnUwYk4FlgT37ZuShsXsDV0bEgjz7d6QEsm/N+vsBd0bErPx8GvB4g7i5wA0dLLuZmQ2g+BatkvbIf74zT3eWNAeYExFTI+I2SRcCp0laFpgFHAZsQOVHPyIek3QqMEHSs8CtpKSyHbBbJe5FSV8inaj3IHBVjjkQODwiXmjtLZuZWSuauaf3pJrnZ+bpVGBc/vsTpJPtTgRWA24HdoqIW2vWPRaYB3wWWBu4B9grIi6tBkXEDyQFcARwFHA/8OmIOBMzMxtUxQkjImo7qevFPAd8Pj/6i3uZlFROLNjmWaTLg5iZWQ/5arVmZlbECcPMzIo4YZiZWREnDDMzK+KEYWZmRZwwzMysiBOGmZkVccIwM7MiThhmZlbECcPMzIo4YZiZWREnDDMzK+KEYWZmRZwwzMysSDP3wzCzIWDk0Ze1tf7sU3btUElsaeMahpmZFXHCMDOzIk4YZmZWxAnDzMyKOGGYmVkRJwwzMyvihGFmZkWcMMzMrIgThpmZFfGZ3tYR7Z59bGZDn2sYZmZWxAnDzMyKOGGYmVkRJwwzMyvihGFmZkWcMMzMrIgThpmZFXHCMDOzIk4YZmZWxAnDzMyKOGGYmVkRJwwzMyvihGFmZkWcMMzMrIgThpmZFXHCMDOzIk4YZmZWxAnDzMyKOGGYmVkRJwwzMyvihGFmZkWcMMzMrIgThpmZFXHCMDOzIk4YZmZWxAnDzMyKOGGYmVkRJwwzMyvihGFmZkWcMMzMrIgThpmZFRnWyY1JGgdMqbPo6YhYrRI3HPgm8EFgRWAa8LmI+EvN9lYATgD2A1YDbgO+EBHXdrLcZkuTkUdf1vK6s0/ZtYMlscVNRxNGxWeAP1Wev9T3hyQBk4ENgMOBJ4EJwBRJm0XEA5X1zgV2BY4C/g58CrhC0lYRcVuXym5mZnV0K2HcHRE3Nlg2Htga2C4ipgBImgbMAv6blGyQ9Dbgo8CBEXFenjcVmAEcn7djZmaDpBd9GOOBh/qSBUBEPA1cCuxWE/cicGEl7iXgAmBHScsPTnHNzAy6lzB+JullSU9I+l9J61eWjQHurLPODGB9SStX4mZFxPw6ccsBozpeajMza6jTTVJPA98GpgLPAG8HjgGmSXp7RDwGjABm11l3bp4OB+bluCf7iRtRrwCSDgUOBVh//fXrhZiZWQs6mjAi4s/Anyuzpkq6FriZ1DfxRUBA1FlddZ6XxNWW4WzgbICxY8fWW9/MzFrQ9T6MiLgVuBfYPM+aS/3awfA8fbIwbm6dZWZm1iWD1eldrS3MIPVP1BoN3B8R8ypxG0haqU7cC8DMbhTUzMzq63rCkDQWeAtwU541GVhX0jaVmFWAD+RlVOKWBfasxA0D9gaujIgFXS66mZlVdPpM75+Rzqe4FXiK1Ok9AXgQOD2HTSad2T1R0lEsPHFPwDf6thURt0m6EDhN0rJ5u4eRTvjbt5PlNjOzgXV6lNSdwD6kM7hXAh4BLgK+EhGPA0TEK5LeD3wLOBNYgZRAto2If9Zs7xPAScCJpEuD3A7slPtFzMxsEHV6lNTJwMkFcXOBA/Ojv7jngM/nh5mZ9ZCvVmtmZkWcMMzMrIgThpmZFXHCMDOzIk4YZmZWxAnDzMyKOGGYmVkRJwwzMyvSrVu02mJo5NGX9boIZjaEuYZhZmZFnDDMzKyIE4aZmRVxwjAzsyJOGGZmVsQJw8zMijhhmJlZEScMMzMr4oRhZmZFnDDMzKyILw1iZsXauXzM7FN27WBJrBdcwzAzsyJOGGZmVsQJw8zMijhhmJlZEScMMzMr4oRhZmZFnDDMzKyIE4aZmRVxwjAzsyJOGGZmVsQJw8zMijhhmJlZEScMMzMr4oRhZmZFfHnzJUw7l582M+uPaxhmZlbENQwzGxS++dLizzUMMzMr4oRhZmZFnDDMzKyIE4aZmRVxwjAzsyJOGGZmVsQJw8zMijhhmJlZEScMMzMr4oRhZmZFnDDMzKyIryVlZkNeu1dh9rWoOsM1DDMzK+KEYWZmRZwwzMysiPswhhjfMc/MhirXMMzMrIgThpmZFXGTlJkt8Xx72M4Y8jUMSW+U9EtJT0t6RtJFktbvdbnMzJY2QzphSFoJuBrYBPg4sD+wETBF0mt7WTYzs6XNUG+SOgR4M7BxRMwEkHQH8Dfg/wHf6WHZGvJIJzNbEg31hDEeuLEvWQBExCxJNwC7MUQThpktOdz/sdBQTxhjgEvqzJ8B7NnNF3Ytwcza1avfkW4lqqGeMEYAT9aZPxcYXm8FSYcCh+an8yTd08brrw483sb6SxPvq+Z4fzXH+6sJ+nrb++tN9WYO9YQBEHXmqWFwxNnA2Z14YUnTI2JsJ7a1pPO+ao73V3O8v5rTrf01pEdJkWoXI+rMH079moeZmXXJUE8YM0j9GLVGA3cNclnMzJZqQz1hTAa2lPTmvhmSRgLvzsu6rSNNW0sJ76vmeH81x/urOV3ZX4qo10UwNOST824HngO+SOrPOAF4HfDWiJjXw+KZmS1VhnQNIyL+BWwH3Av8FPgZMAvYzsnCzGxwDekahpmZDR1DuobRDZLWk3S6pGmS5kuK3C9Ssu4ykiZImi3peUm3S9q9uyXunTb31ewcX/v4YHdL3TuS9pD0K0n/kPScpHsknSzpdQXrriDpm5IezutOk/TewSh3L7S5r+p9rkLSZoNR9l6QtKOkqyU9ImmBpAck/ULS6IJ1h0s6R9Ljkv4l6SpJm7ZSjqUuYQCjgL1Iw3Kva3LdE4DjgDOAnYEbgUmSdulkAYeQdvYVwBXAVjWPqR0r3dBzJPAycAywE/B94DDg95IG+q6dS7p22peB9wMPA1cswT+C7ewrgPNZ9LN1b1dKOjSMAG4BPg3sAEwgjSC9UVLdk+wAJIk0QGgn4HBgd2BZ0gVc12u6FBGxVD2AZSp/H0zqSB9ZsN6awALgqzXz/wDc0ev3NZT2VY6fDUzs9XsY5P21Rp15H8v7bbt+1ntbjvlEZd4w4B5gcq/f11DaVzkugBN7/R56/QA2zvviiH5idssx21bmrUq6Wsb3mn3Npa6GERGvtLjqjsBywMSa+ROBTSVt0FbBhqA29tVSKSLm1Jn9pzxdt59VxwMvAhdWtvUScAGwo6TlO1bIIaKNfWULPZGnL/YTMx54KCKm9M2IiKeBS0nJpClLXcJowxhSDWNmzfwZeTpgW+JS6AO572OBpBuX5P6LfmyTp3f3EzMGmBUR82vmzyAdpIzqRsGGoJJ91eew/Lman9v239PNgg0Vkl4jaTlJGwFnAY+QDiwaGQPcWWf+DGB9SSs38/pOGOVGAE9FrtNVzK0st4UuJbWZ7gjsCzwP/FrSfj0t1SCStC5wPHBVREzvJ7S/i2z2LV+iNbGvINXq/xPYnnSh0dcDV0sa19VCDg03kQ5c7wXeSmq+e6yf+IE+W3Uv4trI4nDxwaFCNHkhxKVZRBxefS7p16RBAiezaLPeEicfuV0CvAR8YqBwluLPVpP7iojYv/L0OkmXkI6iTwS27kohh479gVVIN5Y7kjRIYOuImN0gvqOfLdcwys0FhudRB1XDK8utgYh4GZgErCfpDb0uTzdJWoE0MuXNwI4R8cAAq8yl8UU2+5YvkVrYV4uIiGeBy4DNO1y8ISci7o6ImyLi58D/BVYGju5nlYE+W01dxNUJo9wMYHlgw5r5fX0XvhjiwPqS7RJ7tqikZYFfAVsAu0TEXwpWmwFskO9hXzUaeIFF+82WCC3uq4abYwn+XNUTEU+RPhv99XH1dwHX+6PJK2Y4YZT7HenLu2/N/P2AOyNi1uAXafEhaRjpLon3R8QjvS5PN+TzB35GOvLbLSJuLFx1Mmls/L/vIpn3197AlRGxoNNl7bU29lW9ba0C7Epq319qSFoL2AS4r5+wycC6kvoGFPTtrw/QwgVcl8o+DEl75D/fmac7S5oDzImIqTnmJeDHEXEQQEQ8JulUYIKkZ4FbSV/o7WhheNriopV9JWkf0j65HPgnsBbwqbyNfQax+IPtf0g/+icB/5K0ZWXZAxHxQD7J6j7g+Ig4HiAibpN0IXBaPuqeRTqJbQMWPUBZUrS0ryQdSTr/YArwEOnOcEcCa7Pk7qu+PsBbgTuAZ4C3AJ8j9ft8O8dsQzov7MCI+EledTIwDZgo6ShSE9QEUo3sG00XpNcnn/TohJdo8LimJub8mvVeQ7pq7j9IIxXuAPbo9fsZavsK2BK4GniUNEb8aeAqUht1z99TF/fV7H7213E5ZmT1eWXdFYHvkIZJPk86Wh7X6/c01PYV6cj4BtLtR18knYswGdii1++py/vrC6QzvZ8C5pNO6jyLyom0wLi8vw6oWXcE8CNSf8Z8UlJ5Wyvl8MUHzcysiPswzMysiBOGmZkVccIwM7MiThhmZlbECcPMzIo4YZiZWREnDDMzK+KEYWZmRf4/nJ+iY8c11WgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N=10000;\n",
    "wmean=1; # in mm\n",
    "wstd=wmean*0.1; # standard deviation in mm\n",
    "hmean=2; # in mm\n",
    "hstd=hmean*0.1; # standard deviation in mm\n",
    "\n",
    "\n",
    "wrand=np.random.normal(wmean,wstd,size=N);\n",
    "hrand=np.random.normal(hmean,hstd,size=N);\n",
    "uts=940; # in N/mm^2=MPa\n",
    "\n",
    "Ffail=uts*wrand*hrand*1e-3; # force in kN\n",
    "plt.hist(Ffail,bins=20)\n",
    "#plt.xlabel('failure load (kN)')\n",
    "#plt.ylabel('relative counts')\n",
    "plt.title('Failure load is {:.2f}+/- {:.2f} kN'.format(np.mean(Ffail),np.std(Ffail)));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this propagation of uncertainty, the final value of failure load seems to be independent of wheher the distribution is uniformly random or normally distributed. In both cases, the failure load is $\\approx 1.9 \\pm 0.25$ kN.\n",
    "\n",
    "The difference is much more apparent if you look at the number of occurrences that failure will occur whether the dimensions are uniformly random or normally distributed. \n",
    "\n",
    "For the uniformly random case, there are approximately 500 parts out of 10,000 that will fail at 1.9 kN. \n",
    "\n",
    "For the normally distributed case, there are approximately 1500 parts out of 10,000 that will fail at 1.9 kN. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Where does a normal distribution come from?\n",
    "\n",
    "\"Everybody believes in the exponential law of errors: the experimenters, because they think it can be proved by mathematics; and the mathematicians, because they believe it has been established by observation\" [5].\n",
    "\n",
    "In the previous example, we drew dimensions from uniformly random distributions and normally distributed random distributions. Why do we use the \"normal\" distribution to describe data with a mean and standard deviation? There are exact statistical methods to derive the normal distribution, but let's take a look at a Monte Carlo approach. \n",
    "\n",
    "Let's say there are 10 different independent factors that could change the dimensions of the steel bar in question e.g. which tool was used, how old the blade is, the humidity, the temperature, and the list goes on. \n",
    "\n",
    "Let's consider one dimension. \n",
    "Each of these factors could change the dimensions of the part, let's use a uniform scale of -1/2-1/2.\n",
    "If the effect is 0, the dimension is exactly as specified. If the effect is -1/2, the dimension is much smaller. Conversely, if the effect is 1/2 the dimension is much larger. Now, we use a Monte Carlo model to generate 10 effects on 10,000 parts as shown in the next block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "factors = np.random.rand(10000,10)-1/2 # each row represents a part and each column is an effect (-1/2-1/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we have created 10,000 parts with 10 uniformly random effects between -1/2-1/2. \n",
    "\n",
    "We sum the effects and look at the final part distribution. The x-axis is labeled \"A.U.\" for arbitrary units, we are just assuming an effect of -1/2-1/2 for each of the 10 factors.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'number of parts')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEQCAYAAAB4JulQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAdSUlEQVR4nO3debRcVZn+8e/DlDA2CWDTRGKCaC8ZDNqRBkVGf0ATCCIiMjRDWmgRVERmFCIkTEEZmkZlAQk/oIVGpA0NmpAmhNYOCoQgBEGBJJAmjEEgEAMhb/+xd0mlUvfm3JtTValbz2etWqfq7F2n3rrr5r7ZZ0+KCMzMzMqwWqsDMDOzvsNJxczMSuOkYmZmpXFSMTOz0jipmJlZadZodQCttvHGG8eQIUNaHYaZWVt56KGHXomITWrPd3xSGTJkCA8++GCrwzAzayuS5tY779tfZmZWGicVMzMrjZOKmZmVxknFzMxK46RiZmalcVIxM7PSOKmYmVlpnFTMzKw0TipmZlaajp9Rb1aWIaffWajenAtHNDgSs9ZxUjFrMicf68t8+8vMzErjpGJmZqVxUjEzs9I4qZiZWWmcVMzMrDROKmZmVhonFTMzK42TipmZlcZJxczMSuOkYmZmpXFSMTOz0jipmJlZabygpNkqygtPWjtyUjFbgaJ/3M3Mt7/MzKxETipmZlYaJxUzMyuNk4qZmZXGScXMzErjpGJmZqVxUjEzs9I0PalI+oykyZJekvSGpBmSRtXU6S9pnKT5khZJmi5p5zrXWk3SGZLmSPqzpEckHdi8b2NmZtWamlQkfRyYAqwJHAMcCDwAXCvpuKqq1+bys4F9gfnAJEnb1VzyPGA0cCXwD8D9wK2S9mng1zAzsy40e0b9l4HVgf0iYmE+d7ekYcARwA/z80OBURExHkDSNGAWcC4wMp/7AHAycGFEXJKvNVXSlsCFwF1N+k5mZpY1O6msBbwLLKo5/ydgQH4+Mte5pVIYEUsk3QycLqlfRCwG9srXu7HmWjcC10kaGhGzG/AdrI/w8itm5Wt2n8qEfLxC0maSNpR0DLAHcGku2xqYHRFv17x3FimJbFlVbzHwVJ16AFuVGbiZma1YU1sqEfGYpF2B24Gv5dPvAl+NiJvz64HAa3XevqCqvHL8U0TECuotR9KxwLEAgwcP7slXMDOzbjS7o/4jwG2k1sR+wOeAHwE/knRYpRpQmygq52tfF6m3nIi4OiKGR8TwTTbZpGj4Zma2As3uUzmf1DLZNyLezef+S9JGwOWSfkJqadRrPlT6XBZUHQdIUk1rpbaemZk1SbOTyrbAI1UJpeK3pBFfHyC1Yg6QtE5Nv8pWwDu834cyC+gHfJhl+1UqfSmPlxy72SrJm3nZqqTZHfUvANtJWqvm/N8Dfya1LiaS5rEcVCmUtAZwMDA5j/wC+CUpyRy27KU4HHjMI7/MzJqv2S2VK4FbgTskXUUaWjwSOAS4NCLeAWZKugW4TNKawGzgOGAoVQkkIl6SdClwhqQ3gRmkxLM7sH8Tv5OZmWXNHv310zzb/TTgGqA/8DRwPPDjqqpHA2OBMcCGwCPA3hExo+aSZwELgW8CmwJPAl+KiDsa+T3MzKy+XicVSQNJrYfHqm5JrVBE/AL4xQrqLAJOyo/u6r1HSjxjin6+mZk1TqE+FUnfkXRB1eudgTmkDvY/5qHCZmbW4Yp21B8OPFP1+mLSLanPAy+SFnY0M7MOV/T21yDgjwCSNgE+BewREffmkVxXNCg+MzNrI0VbKu+R1t0C2Jk0/PfX+fXLdLMkipmZdY6iSWUWcLik9YBRwLSqCYybAy81IjgzM2svRW9/nQv8nDRP5F3SsvMV+5DmiJiZWYcrlFQiYpKkjwGfBGZGxNNVxfcBMxsRnJmZtZeiQ4qPAN6IiNtqEgqkGfJ/W3pkZmbWdor2qYwnLdxYz9BcbmZmHa5oUuluj5J1gSUlxGJmZm2uyz4VSduR+lAq9pO0TU21tYEvk+ewmJlZZ+uuo35/4Jz8PEiLN9bzKvBPZQZl1ltF9xYxs8boLqlcBkwg3fp6BjiA5Ud5LQZerLNPvJmZdaAuk0pEvA68npdhOQl4JiLmNi0yMzNrOyvsqM8bZ13A+3u/m5mZ1VV09NfvgS0aGYiZmbW/oknlbOC7krZtZDBmZtbeiq79dRqwHvCwpDnAfNKIsIqIiF1Kjs3MzNpM0aTyHvB4IwMxM7P2V3RByV0bHIeZmfUBRftUzMzMVqjo7S8AJA0APgL0ry2LiPvKCsrMzNpToaQiqT9wHfAlul5ccvWygjIzs/ZUtKXyXWBX4EjgBuB40j71RwF/A3yzAbGZWYl6si7anAtHNDAS68uK9qkcSNpS+Ob8+jcRMT4PI34E2LsRwZmZWXspmlQGA7Mi4j3SHvXrVpVdBxxcdmBmZtZ+iiaVV0mTHwGeA4ZVlW1M2lfFzMw6XNE+lfuBTwC/AG4DzpO0PmnHx28Dv2pMeGZm1k6KJpWLSLfAAMYAW5L6WFYnJZzjyg/NzMzaTdEZ9Q8CD+bnbwIHSuoH9IuINxoYn5mZtZEeTX6sFhGLSTs/mpmZAT1YpkXSRyRdL+kPkt7KxwmStmxkgGZm1j4KJRVJu5Lmo+xL6kO5Kh/3Ax6V1KNl7yXtI+k+SQslvSHpQUm7V5UPkHSNpFdyAptSby8XSf0ljZM0X9IiSdMl7dyTWMzMrDxFb399H3gY2CsiFlZO5hFgk3P58CIXkvTPwJX5cR4psW0HrJPLBUwEhgJfB14DzgCmStouIuZVXe5aYARwCvAMaab/JEk7RsTMgt/NzMxKUjSpbAUcXJ1QIHXaS7oI+EmRi0gaAlwGnBIRl1UVTap6PhLYCdg9Iqbm900HZgOnAt/I54YBhwKjImJ8PjcNmEUamTay4HczM7OSFO1TmQes1UXZWsD/FrzOKGAp8KNu6owEnq8kFICIeB24A9i/pt67wC1V9ZaQlpLZK49OMzOzJiqaVC4CvidpUPXJ/Poc4PyC19kJeAL4sqSnJS2R9JSk46vqbA08Vue9s4DBktarqjc7It6uU28t0lwaMzNroqK3v3YB1geelnQ/8CLw18AO+fmuuTMf0n71R3Zxnc3yYxxwJvA0cBBwpaQ1IuJyYCAwp857F+TjAGBhrvdaN/UGdvVlJB0LHAswePDgrqqZmVkPFU0qO5H2qZ8PfCg/yK8BPltVN7q5zmqk5HRURPwsn7sn97WcIekK0n4t9a5Ru49L0XrLiYirgasBhg8f3l28ZmbWA0Vn1A8t6fNeJe0ceXfN+cmk5fP/htTSqNfKGJCPldbJAt5fOqZevQV1yszMrIGavUf9rC7OV1oXS3OdrevU2Qp4tmoE2ixgqKR16tR7B3hqJWM1M7MeanZSuT0f96o5vxcwLyJeIM1RGVQ9oVLSBqSJlhOr3jMRWJPUJ1OptwZpb5fJeRkZMzNrol6v/dVLdwFTgR9L2pg0YfGLwJ7A0bnORGA6cKOkU3h/8qOAiysXioiZkm4BLpO0Jmkey3GkSZOHNefrmJlZtaYmlYgISZ8HLgC+R+r/eAI4LCL+LddZKmlf4BLScjD9SUlmt4h4ruaSRwNjScvxb0je2jgiZjTj+5iZ2bKa3VIhL5V/fH50VWcBaaLkqBVcaxFwUn6YmVmLddmnImmGpK3z87Mlbda8sMzMrB1111G/Le/vS38O8MHGh2NmZu2su9tfzwOfl/QiqZN8U0ldTj+PiGfLDs7MzNpLd0nlx6QO8FNJM9dv76YupP3qzcysg3WZVCLifEl3kyYTjieN2HqmWYGZmVn76Xb0V0Q8ADwg6Sjghoh4oilRmZlZWyq69tdujQ7EzMzaX+FlWiRtK+mnkl7O+6C8JOnf6+0db2ZmnalQS0XSp4BpwCLSMiovAJuS1uMaIWnniHioYVFaxxty+p2tDqGjFP15z7lwRIMjsXZTdEb9BaTdGPeIiDcrJyWtD0zJ5XuWH56ZmbWTore/dgAuqE4oAPn1RcCOZQdmZmbtp2hSWdHuiN490czMCieV3wBn5ttdfyFpXeA04P6yAzMzs/ZTtE/lTOBeYK6k/yTtTb8pMAJYG9i1EcGZmVl7KTpP5beSdgDOJu3SOJC0B/w9wHkR8WjjQjQzs3ZReD+ViPgdaZdGMzOzupq9R72ZmfVhTipmZlYaJxUzMyuNk4qZmZXGScXMzEqzwqQiaS1JCySNbEZAZmbWvlaYVCLiHWAJ8OfGh2NmZu2s6O2v/8BzVMzMbAWKTn78BXCFpJ+SEsx8ahaRjIh7So7NzMzaTNGkcls+fiE/KgJQPq5eYlxmZtaGiiYV71FvZmYrVHRByWmNDsTMzNpf4QUlASRtTNoFciPgjohYIKk/8E5ELG1EgGZm1j4Kjf5SMg6YB0wErgOG5OKfA2c1JDozM2srRYcUnwGcAJwL/D2pc77iDmDfkuMyM7M2VDSpfAU4NyLOB2bUlD0FfLi3AUj6paSQNKbm/ABJ10h6RdJbkqZI2rbO+/tLGidpvqRFkqZL2rm38ZiZWe8VTSqD6Hof+neAdXvz4ZIOAYbVOS/Sbba9ga8DBwJrAlMlfbCm+rXAMaRdKfclzaGZJGm73sRkZma9V7Sj/n+BbYCpdcqGAbN7+sGSNgQuBb4F/FtN8UhgJ2D3iJia60/Pn3Mq8I18bhhwKDAqIsbnc9OAWaRbdV6vbBU35PQ7Wx2CmZWoaEvlVuBsSZ+pOheSPgp8G7i5F599MTArIn5Sp2wk8HwloQBExOuk/pv9a+q9C9xSVW9JjmcvSf16EZeZmfVS0aQyGngCuA/4Yz53K/Bofn1hTz5U0k7AEcDXuqiyNfBYnfOzgMGS1quqNzsi3q5Tby1gy57EZWZmK6fo5MdFknYl3Wrai9Q5/ypwHnBTbh0UImlN4MfAJRHxZBfVBgJz6pxfkI8DgIW53mvd1BtYNC4z67mity/nXDiiwZHYqqLw5MeIeA+4IT9WxmnA2sDYbupU1hOrd7439ZYtlI4FjgUYPHhwd1XNzKwHejqjfktge9JosHnAAxHxVA/eP5g0UfIrQL+aPo9+ufP+TVJLo14rY0A+VlonC4B6WWFAVflyIuJq4GqA4cOH10tKZmbWC4WSSl6K5SrgH1l2NeL3JF0PHB8RiwtcagugP3BjnbKT8+MTpD6RPevU2Qp4NiIW5tezgAMkrVPTr7IVaahz4YRnZmYrr2hH/SXAYcA5pM7v9fNxNCnRjCt4nZmkFY9rH5ASzW6kRDARGCRpl8obJW0A7JfLKiaS5q8cVFVvDeBgYHLBRGdmZiUpevvry8D38oz6imeAsWmeIt8izx3pTkT8Cbi39ny+xtyIuDe/nghMB26UdArpdtcZpL6Si6uuN1PSLcBleQDAbOA4YCgpCZqZWRMVban0A37bRdlvSMN3S5NXPN4XuJt02+124D1gt4h4rqb60cB4YAxwJ7A5sHdE1C4nY2ZmDVa0pTKF1McxpU7ZnsBKbSUcEcuN1oqIBcCo/OjuvYuAk/LDzMxaqMukImmLqpc/AG6QtC5p0uOLwF8DXwL2AQ5vZJBmZtYeumupPMWyc0BE6q/4as05gGl4j3ozs47XXVI5umlRmJlZn9BlUomI65sZiJmZtb+io7/MzMxWqPAyLZL2Jk0y3Jw0K75aRMQuy7/LzMw6SaGWiqRTgbtIc0fWJc0ZqX4sbVSAZmbWPoq2VE4gLVd/Ql6t2Kxb3tHRrDMV7VPZALjVCcXMzLpTNKlMAnZoZCBmZtb+enL763ZJAUymzm6LEfFMmYGZmVn7KZpUgrR51ljSwo31eEa9mVmHK5pUJgCfBi4FniBtgGVmZraMokllV9LIrwmNC8XMzNpd0Y76V0grE5uZmXWpaFK5AviaJC/rYmZmXSp6+2sAsA3wuKS7WX70V0TEOaVGZmZmbadoUjmr6vlH65QH4KRiZtbhCiWViPBtLzMzWyEnCzMzK03hpe/NzHqr6AKjcy4c0eBIrNEKJRVJS1l2v/rlRIRn1JuZdbiiLZVzWT6pbATsCfQjzbg3M7MOV7SjfnS985JWB+4AXi8xJjMza1Mr1VGf91e5CjixnHDMzKydlTH6qx8wsITrmJlZmyvaUT+4zum1SLPsLwQeLDMoMzNrT0U76udQf/SXgKeB48sKyMzM2lfRpDKK5ZPKn4G5wAPeu97MzKD46K8JDY7DzMz6AC/TYmZmpSm8TIukI4FDgMFA/5riiIgPlxmYmZm1n0ItFUnfBcYDmwEzgWk1j/sKXueLkm6TNFfSIklPSrpA0vo19QZIukbSK5LekjRF0rZ1rtdf0jhJ8/P1pkvauUgsZmZWvqItlX8CLo+Ib63k550MPAucCcwDPgGMBnaT9OmIWCpJwERgKPB10oZgZwBTJW0XEfOqrnctMAI4BXiGNAptkqQdI2LmSsZqZmY9VDSpbERajmVl7RcRL1e9niZpAXA9sCtwDzAS2AnYPSKmAkiaDswGTgW+kc8NAw4FRkXE+HxuGjCLtFbZyBLiNTOzHijaUT8NGLayH1aTUCoeyMdB+TgSeL6SUPL7Xicltf2r3jcSeBe4pareEuBmYC9J/VY2XjMz65miLZUTgZ9JehW4C1hQWyEilvYyhl3y8ff5uDXwWJ16s4AjJK0XEQtzvdkR8XademsBW+bnVqKi+2KY9Yb3XWl/RZPKH/JxfBfl0YNr/YWkQaRbVVMiorLUy0DSDP5alUQ2AFiY673WTb0u1yOTdCxwLMDgwfVWoDEzs95Ymf1UVoqk9YCfA0uAo6uLuvgs1XldpN5yIuJq4GqA4cOHl/q9zMw62Urtp9JbkvqTRnhtAexSM6JrAfVbGQPy8bWqevWaGQOqys3MrImaPqNe0prAbcD2wD4R8WhNlVmk/pJaWwHP5v6USr2hktapU+8d4KnyojYzsyKamlQkrQbcBOwB7B8R99epNhEYJGmXqvdtAOyXy6rrrQkcVFVvDeBgYHJELC7/G5iZWXd63Lm+kv6VlATGAm9J2qGqbF6+DTYRmA7cKOkU3p/8KODiSuWImCnpFuCy3PqZDRxHmjR5WDO+jJmZLavZt7/+IR/PIiWO6sdX4C9Dk/cF7iZtVXw78B6wW0Q8V3O9o0kj0sYAdwKbA3tHxIzGfg0zM6unqS2ViBhSsN4C0h4uo1ZQbxFwUn6YmVmLeel7MzMrjZOKmZmVxknFzMxK46RiZmalcVIxM7PSOKmYmVlpnFTMzKw0zZ5Rb6so75NiZmVwS8XMzErjpGJmZqVxUjEzs9K4T8XM2o73sl91uaViZmalcVIxM7PSOKmYmVlpnFTMzKw0TipmZlYaJxUzMyuNk4qZmZXG81T6OK/pZZ3M81mazy0VMzMrjZOKmZmVxknFzMxK46RiZmalcVIxM7PSOKmYmVlpnFTMzKw0TipmZlYaT35sQ57QaGarKicVM+t4PfmPmmffd89JZRXiFoiZtTv3qZiZWWn6RFKRtLmkn0p6XdIbkn4maXCr4zIz6zRtf/tL0jrAPcBi4EgggDHAVEkfj4i3WhmfmfUtXvm4e22fVIBjgC2Av42IpwAk/Q74I/DPwA9aGJuZWUfpC0llJHB/JaEARMRsSb8G9mcVSCrugDfrPJ3aoukLSWVr4Od1zs8CDmrkBztZmNnK6mvJpy8klYHAa3XOLwAG1HuDpGOBY/PLhZKebFBsq7qNgVdaHUSH8c+8+frEz1wXtTqC5Xyo3sm+kFQgdc7XUpeVI64Grm5cOO1B0oMRMbzVcXQS/8ybzz/z5uoLQ4pfI7VWag2gfgvGzMwapC8klVmkfpVaWwGPNzkWM7OO1heSykRgB0lbVE5IGgJ8JpdZ1zr+FmAL+GfefP6ZN5Ei6nVHtA9J6wKPAIuA75D6V84D1gc+HhELWxiemVlHafuWSp4xvzvwB+AG4CZgNrC7E4qZWXO1fUvFzMxWHW3fUrGVJ+mjki6X9DtJCyXNlzRR0rBWx9YXeMHT5pL0RUm3SZoraZGkJyVdIGn9VsfWCdxSMSSdQJoMej0wA9gQOBX4BPCZiHioheG1tbzg6SOkBU8rfX5jgHVIfX5e8LRkku4HniWttDGP9Hs8GngC+HRELG1ddH2fk4ohaWPg1aj6ZZD0V8Ac4I6IOKJVsbU7Sd8krT9XveDpUNKCp6dGRMvXputrJG0SES/XnDuC9J+mPSLintZE1hl8+8uIiFei5n8XEfE6afDDoNZE1WfUXfAUqCx4aiWrTSjZA/no3+cGc1KxuiQNBLYBft/qWNrc1sBjdc7PIk3QtebYJR/9+9xgTirWlX8hrZ92WasDaXM9XvDUyiVpEHAuMCUiHmx1PH2dk0ofJOlzkqLA494u3n8GcChwQvVtG+u1Hi14auWRtB6pw34JcHSLw+kIfWWVYlvW/wAfK1Dv7doTkr4KnA98JyKuKzuwDuQFT1tEUn/SUk1bALtExLwWh9QRnFT6oIh4mzR8skck/SNwFfD9iBhbemCdyQuetoCkNYHbgO2Bz0XEoy0OqWP49pcBIOkAYDxwTUSc3Op4+hAveNpkklYjLde0B7B/RNzf4pA6iuepGJJ2BiaT/ud8AlA9OWxxRDzcksD6AC942nySfgh8FRgL/GdN8TzfBmssJxVD0mjgnC6K50bEkOZF0/fkJVkuBf4fqYP+v4ATI2JOK+PqqyTNoYutboHvRcTo5kXTeZxUzMysNO5TMTOz0jipmJlZaZxUzMysNE4qZmZWGicVMzMrjZOKmZmVxknFbAUknSnpWUlLJM3M5zbNWy4vyItznljyZ54o6Qu9eN93cjw/68F7jsrv2bKL8pA0pqexWGfy2l9m3ZC0PWlm9jjgP4A3c9HZpD06jgLmk3bJLNOJwK+Awskhq+zSOULSRhHxarlhmXXPLRWz7lVWe/5RRPxP1cKEHwMeiYjbI+L+iHihRfH9haRPAx8B7gLWAg5pbUTWiZxUrGNJGpZvYb0maZGkX0v6bFX5vcCE/PLpfBtogqQAdgU+W7U3zZD8nqGSbpL0sqTFkmbmxTrrffbtkl7Nn/1k3semepmRw6quP6H2GnUcCbwHHAM8x/utFrOmcVKxjiTpk6R9ZwaS/ggfCLwKTJH0d7na14AL8vMvADuS1kjbEfgd8HB+viMwX9LmwG+AYcC3SPvTzwBukzSy6rO3B6YDH871RgA/AD6YqxwAvABMqrr+eSv4Pv2BLwF3R8TzwI3ApyQV2VfHrDTuU7FONQ54Ftg9It4BkDSJtJ/8d4HPR8Tjkp7J9R+uWgByrqQ3gSXVy6rnhTlF2hCq0pcxKSebc3l/qftLSAlsh7z3DcA9letExMOSFgOv9GDZ9v2BDYH/n19fD5xBar2cXvAaZivNLRXrOJLWJnWy3woslbSGpDVICWEKsHMvL703qT/j9co183UnAcMkbSBpHdJeKjdVJZQyHAm8QRpMQEQ8SWo1HZ73FzFrCv+yWScaCKxOapG8W/M4ARjQyz/EHyD1Y9Rec1wu34i0jfBqQGl7ekjaFNgTuBPoJ2lDSRuSdj4cRNqsqjtL8nH1OtdevaaOWbd8+8s60Z9IG5H9K+/fLlpGRCytd34FXgX+G7ioi/LnSX+4l5L+2Jfl8HzdQ6g/4utI4O5u3v9SPm4GPFlTtlk+vrgyAVrncFKxjhMRb0n6b1KH+oxeJpB6fknqVJ8VEYu6qiTpV6TbUud2U28xsHbBzz0CmEuaM1PrNOAASetHxJt1yiHdJltEGqwwtabswHy8t2As1uGcVKxTnQTcR+pIv5Y0gXFj4JPA6hHRm87ts4HfAvdJupI0IXIAsA2wRUSMyvVOBqYB0yV9n3QrbAtgu4j4eq7zOGnI8r6kkWCv1NspMo9i2xYYHRH31invT+rr+SIwXtKHgKeBcyPiXICIeF3SWGCMpH7AHcA7pH6nbwMTIuL3Vdc8ArgO2CMipvXi52R9mPtUrCNFxAzgU6RbVlcAk4HLSX+g7+vlNZ8FhpP2pD+fdMvph6Q/ztWjux4gddY/B/wLqXP/FJbtZzmDdCvq34EHgNFdfOyRpH3vJ3RRPjl/zpH5tUi3ypb5tx8RY4Gjge2Am0kd/iOBs4Cv1FxztXwNdfGZ1sG8nbCZmZXGLRUzMyuNk4qZmZXGScXMzErjpGJmZqVxUjEzs9I4qZiZWWmcVMzMrDROKmZmVpr/A/wt4URURV3eAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dims = np.sum(factors,axis=1)\n",
    "\n",
    "plt.hist(dims,30)\n",
    "plt.xlabel('effect A.U.')\n",
    "plt.ylabel('number of parts')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, depending upon which random numbers were generated, you should see what looks like a normal distribution. \n",
    "\n",
    "Normal distributions come from the assumption that we have a large (or infinite) number of uncontrollable factors that can change our desired result. In our case, ideally each factor would have an effect of 0, because then it is exactly as specified, but the reality is that we can't control most factors. As engineers, we always have to consider the uncertainty in our models and measurements. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We learned:\n",
    "\n",
    "* How to generate \"random\" numbers in Python$^+$\n",
    "* The definition of a Monte Carlo model\n",
    "* How to calculate $\\pi$ with Monte Carlo\n",
    "* How to take the integral of a function with Monte Carlo\n",
    "* How to propagate uncertainty in a model with Monte Carlo\n",
    "* **Bonus**: use Sympy to do calculus and algebra for us! *no need for Wolfram, sorry Stephen*\n",
    "* How to generate a normal distribution using uniformly random numbers\n",
    "\n",
    "$^+$ Remember, the computer only generates pseudo-random numbers. For further information **and** truly random numbers  check [www.random.org](https://www.random.org/randomness/) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# References\n",
    "\n",
    "1. [Why the Monte Carlo method is so important today\n",
    "Dirk P. Kroese, Tim Brereton *et al.* Wiley Interdisciplinary Reviews: Computational Statistics, 6, 6, 11 2014](https://onlinelibrary.wiley.com/doi/full/10.1002/wics.1314)\n",
    "\n",
    "2. [Wikipedia: Monte Carlo integration](https://en.wikipedia.org/wiki/Monte_Carlo_integration)\n",
    "\n",
    "3. [Weinzierl, S. (2000). \"Introduction to Monte Carlo methods\"](https://arxiv.org/abs/hep-ph/0006269)\n",
    "\n",
    "4. Meurer A, _et al._ (2017) SymPy: symbolic computing in Python. PeerJ Computer Science 3:e103 https://doi.org/10.7717/peerj-cs.103\n",
    "\n",
    "5. Whittaker, E. T. and Robinson, G. \"Normal Frequency Distribution.\" Ch. 8 in The Calculus of Observations: A Treatise on Numerical Mathematics, 4th ed. New York: Dover, p. 179, 1967."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problems\n",
    "\n",
    "1. Calculate the area of a unit circle using the Monte Carlo integration method that we devloped in Example 2. Calculate the area of a quarter of the circle, then multiply the final result by 4. \n",
    "\n",
    "    _Hint: You have to create a function that describes the perimeter of the circle._\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. 100 steel rods are going to be used to support a 1000 kg structure. The\n",
    "rods will buckle when the load in any rod exceeds the [critical buckling\n",
    "load](https://en.wikipedia.org/wiki/Euler%27s_critical_load)\n",
    "\n",
    "    $P_{cr}=\\frac{\\pi^3 Er^4}{16L^2}$\n",
    "\n",
    "    where E=200e9 Pa, r=0.01 m +/-0.001 m, and L is the \n",
    "    length of the rods supporting the structure. Create a Monte\n",
    "    Carlo model `montecarlo_buckle` that predicts \n",
    "    the mean and standard deviation of the buckling load for 100\n",
    "    samples with normally distributed dimensions r and L. \n",
    "\n",
    "    ```python\n",
    "    mean_buckle_load,std_buckle_load=\\\n",
    "    montecarlo_buckle(E,r_mean,r_std,L,N=100)\n",
    "    ```\n",
    "\n",
    "    a. What is the mean_buckle_load and std_buckle_load for L=5 m?\n",
    "\n",
    "    b. What length, L, should the beams be so that only 2.5% will \n",
    "    reach the critical buckling load?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def montecarlo_buckle(E,r_mean,r_std,L,N=100):\n",
    "    '''Generate N rods of length L with radii of r=r_mean+/-r_std\n",
    "    then calculate the mean and std of the buckling loads in for the\n",
    "    rod population holding a 1000-kg structure\n",
    "    Arguments\n",
    "    ---------\n",
    "    E: Young's modulus [note: keep units consistent]\n",
    "    r_mean: mean radius of the N rods holding the structure\n",
    "    r_std: standard deviation of the N rods holding the structure\n",
    "    L: length of the rods (or the height of the structure)\n",
    "    N: number of rods holding the structure, default is N=100 rods\n",
    "    Returns\n",
    "    -------\n",
    "    mean_buckle_load: mean buckling load of N rods under 1000*9.81/N-Newton load\n",
    "    std_buckle_load: std dev buckling load of N rods under 1000*9.81/N-Newton load\n",
    "    '''\n",
    "    \n",
    "    # your code here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Generate your own normal distribution using uniformly random numbers between -1/2 and 1/2. \n",
    "\n",
    "    a. What is the effect of changing the number of factors?\n",
    "    \n",
    "    b. What is the effect of changing the number of samples?\n",
    "    \n",
    "    *Hint: for a-b try plotting histograms of the results.*\n",
    "    \n",
    "    c. How would you change the mean in your generated distribution?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}