Monte_Carlo_integration_extremely_simple_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximating Pi by Monte Carlo Integration\n",
    "\n",
    "The area $ A $ of a circle is $ A = \\pi r^{2} $, where $ r > 0 $ is the radius. When $ r = 1 $, the area is $ \\pi $.\n",
    "\n",
    "The area of a $ \\frac{1}{4}$-th circle inside a $ 1 x 1 $ square is $ \\frac{\\pi}{4} $.\n",
    "\n",
    "We generate points distributed uniformly over the $ 1 x 1 $ square, and then count the percentage of the points which lie inside the quarter circle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.random import RandomState\n",
    "import matplotlib.pyplot as plt\n",
    "import timeit\n",
    "\n",
    "n = 1_000_000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## With NumPy's Vectorized Operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Approximate pi: 3.141372\n",
      "NumPy's value for pi: 3.141592653589793\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAUQ0lEQVR4nO3df4xddZ3G8eehFcgiithRkFYLbhG66gq5qYDLAgGl00i7BlansSsm4ARdyG4kC92wsab+RcmyyabdXVsx/kiVIqg7Ir9UICRC2Q7hhxRaGIvKFBYGZJtIlVL47B/3FG6ndzpnpufe7znnvl/JhHvO+TL3YYZ55sz58T2OCAEAqu+g1AEAAMWg0AGgJih0AKgJCh0AaoJCB4CamJnqjWfNmhVz585N9fYAUEkPPPDACxHR125bskKfO3euhoeHU709AFSS7d9OtI1DLgBQExQ6ANQEhQ4ANUGhA0BNUOgAUBMUOgDUxKSFbvubtp+3/egE2237322P2H7E9snFxwQATCbPHvq3JC3cz/Z+SfOyj0FJ/3ngsQAAUzVpoUfEPZJ+v58hSyR9J5o2SjrC9tFFBewZGzZIM2ZIdv6PPeMPOkg67DDp/POlM86QrrhCeuGF1P9FALqsiGPox0h6umV5NFu3D9uDtodtD4+NjRXw1jXy2c9Kr78+tX9nz/gIaedO6Yc/lO65R7rmGqmvr/0vgZkzpUMPbf4CAVArXT0pGhFrI6IREY2+vrZTEfSu117r3vu88oo0MLBv2c+Zw949UGFFFPp2SXNalmdn61A1o6Pt9+6POELauDF1OgCTKKLQhyR9Lrva5RRJOyLi2QI+L8pixw7p1FP3Lvkjj6TkgZLJc9ni9yXdJ+kDtkdtX2T7EtuXZENukbRN0oikdZK+1LG0KI+XXtq75A87TLrtttSpgJ426fS5EbF0ku0h6e8LS4Rq2rlT6u9/c/k975HuvFP6wAfSZQJ6DHeKojOeeUY64YQ39+Avvjh1IqD2KHR0x3XX7X0Mft261ImA2qHQy+Kww1In6K7BwWaxv/WtnFwFCkKhl8XLL6dOkMbLL795cnXmTE6sAgeAQkd5vPZa88QqV80A00Kho5z2XDVjS/PmSVu3pk4ElB6FjvIbGXnzipnLLkudBigtCh3Vsnp1s9iPOoqTqcA4FDqq6bnnmidTZ8zgWDuQodBRba+//uax9vPOS50GSIpCR33cfHOz2JctS50ESIJCR/2sX98s9k98grnd0VModNTXz37WnNv9xBO57BE9gUJH/W3Z0rzskWkGUHMUOnrHnmkG5sxhjx21RKGj94yONvfYDzmEPXbUCoWO3rVrV3OP/YILUicBCkGhAzfd1Lwq5txzUycBDgiFDuxxxx3NYr/yytRJgGmh0IHxVq2SDjpI2rAhdRJgSih0oJ0IaWBA+tCHuDkJlUGhA/vz6KPNm5N4yDUqgEIH8tjzkGtmdkSJUejAVPT3N5+gxGEYlBCFDkzVyEjzMMzKlamTAHuh0IHpWrFCOumk1CmAN1DowIF46CGuXUdpUOhAEVatkk4+mWPrSIpCB4ry4IMcW0dSFDpQtBUrpI99LHUK9KBchW57oe2ttkdsL2+z/b2277L9oO1HbC8qPipQIffeK73lLUzPi66atNBtz5C0RlK/pPmSltqeP27Yv0i6ISJOkjQg6T+KDgpUzu7dzel5uRkJXZJnD32BpJGI2BYRuyRdL2nJuDEh6W3Z67dLeqa4iEDF9fdLZ53FCVN0XJ5CP0bS0y3Lo9m6Vl+VtMz2qKRbJF3W7hPZHrQ9bHt4bGxsGnGBirr77uYJU2ZwRAcVdVJ0qaRvRcRsSYskfdf2Pp87ItZGRCMiGn19fQW9NVAhAwPSunWpU6Cm8hT6dklzWpZnZ+taXSTpBkmKiPskHSppVhEBgdoZHJSWLUudAjWUp9A3SZpn+1jbB6t50nNo3JjfSTpbkmyfqGahc0wFmMj69dLs2RxXR6EmLfSI2C3pUkm3S3pczatZNtteaXtxNuxySV+w/bCk70v6fEREp0IDtbB9u/T+91PqKIxT9W6j0Yjh4eEk711KduoESOmKK6Srr06dAhVg+4GIaLTbxp2iQBmsWsUEXzhgFDpQFqtWSQsWpE6BCqPQgTLZtEm6rO1tHMCkKHSgbFavlmbOlLZuTZ0EFUOhA2X02mvSCSdwZymmhEIHymxggFJHbhQ6UHaUOnKi0IEqYA4Y5EChA1XBHDCYBIUOVMn69dK116ZOgZKi0IGqufxyHkSNtih0oIpWrJDOOy91CpQMhQ5U1c03c1cp9kKhA1W2ejWTeuENFDpQdatWcfULJFHoQD2sX8/hF1DoQG2sXi1t3Jg6BRKi0IE6OfVU6bbbUqdAIhQ6UDf9/ZR6j6LQgTpatIiHT/cgCh2oowipr4+HZPQYCh2os9NPT50AXUShA3U2NiYdd1zqFOgSCh2ou6eeYt6XHkGhA73g5pu5m7QHUOhAr1i/nqce1RyFDvSSwcHUCdBBFDrQa971Li5nrCkKHeg1Y2PSSSelToEOoNCBXvTHP0rnnps6BQqWq9BtL7S91faI7eUTjPm07cdsb7b9vWJjAijcHXcw50vNTFrotmdIWiOpX9J8SUttzx83Zp6kf5b0sYj4C0n/2IGsAIrW38+UuzWSZw99gaSRiNgWEbskXS9pybgxX5C0JiJekqSIeL7YmAA65rTTUidAQfIU+jGSnm5ZHs3WtTpe0vG2f2l7o+2F7T6R7UHbw7aHx8bGppcYQLEiuD69Joo6KTpT0jxJZ0paKmmd7SPGD4qItRHRiIhGX19fQW8N4IANDkrXXps6BQ5QnkLfLmlOy/LsbF2rUUlDEfFqRDwl6Qk1Cx5AVVx+OXOoV1yeQt8kaZ7tY20fLGlA0tC4MT9Wc+9ctmepeQhmW4E5AXQD0+1W2qSFHhG7JV0q6XZJj0u6ISI2215pe3E27HZJL9p+TNJdkv4pIl7sVGgAHbJli3TllalTYJocEUneuNFoxPDwcJL3LiU7dQLgTbfeKi1se20DErP9QEQ02m3jTlEA+1q0KHUCTAOFDmBfEcyfXkEUOoD21q+XNmxInQJTQKEDmNjAQOoEmAIKHcD+XXxx6gTIiUIHsH/XXcesjBVBoQOYXH8/d5FWAIUOIB8uZSw9Ch1APps2sZdechQ6gPw++tHUCbAfFDqA/LZtk1auTJ0CE6DQAUzNihWpE2ACFDqAqWNGxlKi0AFM3apVqROgDQodwPQweVfpUOgApofJu0qHQgcwfUuXpk6AFhQ6gOmLkNatS50CGQodwIEZHEydABkKHcCBu+CC1AkgCh1AEW66KXUCiEIvDzt1AuDAMCVAchQ6gGKsWMFsjIlR6GURkToBcOBOPz11gp5GoQMozpYtqRP0NAodQLG4Lj0ZCh1AsbguPRkKHUDxuOIlCQodQPF4CEYSFDqAzuAhGF1HoQPoDB6C0XW5Ct32QttbbY/YXr6fcefbDtuN4iICqKyNG1Mn6CmTFrrtGZLWSOqXNF/SUtvz24w7XNI/SLq/6JAAKurMM1Mn6Cl59tAXSBqJiG0RsUvS9ZKWtBn3NUlXS/pTgfkAVNkrrzAdQBflKfRjJD3dsjyarXuD7ZMlzYmIn+7vE9ketD1se3hsbGzKYQFU0GmnpU7QMw74pKjtgyRdK+nyycZGxNqIaEREo6+v70DfGkAVPPlk6gQ9I0+hb5c0p2V5drZuj8MlfVDS3bZ/I+kUSUOcGAXwhosvTp2gJ+Qp9E2S5tk+1vbBkgYkDe3ZGBE7ImJWRMyNiLmSNkpaHBHDHUkMoHquuy51gp4waaFHxG5Jl0q6XdLjkm6IiM22V9pe3OmAPWPGjNQJgM7aujV1gtpzJJqHu9FoxPAwO/Fv4IlFqLujj5aeeSZ1isqz/UBEtD2kzZ2iALrj2WdTJ6g9Ch1A9zBXekdR6AC6h7nSO4pCB9BdnBztGAodQHcNDKROUFsUOoDuevjh1Alqi0IH0F0RTKvbIRQ6gO47++zUCWqJQi8LbixCL9m5M3WCWqLQyyLRHbtAMlyTXjgKHUAaXJNeOAodAGqCQgeQzm23pU5QKxQ6gHSWtHs8MaaLQgeQzq5dqRPUCoUOIK0NG1InqA0KHUBay5alTlAbFDqAtHbvTp2gNih0AOkxpW4hKPSy4NZ/9LJzzkmdoBYo9LLg1n/0stHR1AlqgUIHgJqg0AGgJih0AOVw2WWpE1QehQ6gHFavTp2g8ih0AKgJCh1AefCs0QNCoQMoj/POS52g0ih0AOXxwgupE1QahQ4ANZGr0G0vtL3V9ojt5W22f9n2Y7Yfsf0L2+8rPiqAnsBe+rRNWui2Z0haI6lf0nxJS23PHzfsQUmNiPiwpBslrSo6KIAeweWL05ZnD32BpJGI2BYRuyRdL2mv50ZFxF0RsTNb3ChpdrExAfSM++5LnaCy8hT6MZKeblkezdZN5CJJt7bbYHvQ9rDt4bGxsfwpAfSOO+5InaCyCj0panuZpIaka9ptj4i1EdGIiEZfX1+Rbw0APW9mjjHbJc1pWZ6drduL7XMkXSXpjIh4pZh4AIC88uyhb5I0z/axtg+WNCBpqHWA7ZMkfV3S4oh4vviYAHoKd4xOy6SFHhG7JV0q6XZJj0u6ISI2215pe3E27BpJb5X0A9sP2R6a4NMBwOQ+/vHUCSopzyEXRcQtkm4Zt+4rLa95fhSA4vzhD6kTVBJ3igJATVDoAFATFDoA1ASFDqCctm5NnaByKHQA5fSpT6VOUDkUOoByevzx1Akqh0IHgJqg0AGgJih0AKgJCh0AaoJCB1BePI5uSih0AOW1ZMnkY/AGCh1Aed17b+oElUKhA0BNUOgAUBMUOgDUBIUOADVBoQNATVDoAFATFDoA1ASFDgA1QaEDQE1Q6ADKjUfR5UahAyi3z3wmdYLKoNABlNuWLakTVAaFDqDc+vpSJ6gMCh1AuW3fnjpBZVDoAMotInWCyqDQAZSbnTpBZVDoAMqNPfTcchW67YW2t9oesb28zfZDbG/Itt9ve27RQQEA+zdpodueIWmNpH5J8yUttT1/3LCLJL0UEX8u6d8kXV10UADA/uXZQ18gaSQitkXELknXSxr/5NYlkr6dvb5R0tk2B74AoJvyFPoxkp5uWR7N1rUdExG7Je2Q9M7xn8j2oO1h28NjY2PTSwygtxx6aOoEldHVk6IRsTYiGhHR6ONmgb198pOpEwDl9KMfpU5QGTNzjNkuaU7L8uxsXbsxo7ZnSnq7pBcLSdgrfvKT1AkAVFyePfRNkubZPtb2wZIGJA2NGzMk6cLs9QWS7ozgWiMA6KZJ99AjYrftSyXdLmmGpG9GxGbbKyUNR8SQpOskfdf2iKTfq1n6AIAuynPIRRFxi6Rbxq37SsvrP0n622KjAQCmgjtFAaAmKHQAqAkKHQBqgkIHgJqg0AGgJpzqcnHbY5J+m+CtZ0l6IcH75kG26StzvjJnk8qdj2z7el9EtL3VPlmhp2J7OCIaqXO0Q7bpK3O+MmeTyp2PbFPDIRcAqAkKHQBqohcLfW3qAPtBtukrc74yZ5PKnY9sU9Bzx9ABoK56cQ8dAGqJQgeAmqh9ods+0vbPbD+Z/fMdbcZ8xPZ9tjfbfsT2ZzqcaaHtrbZHbC9vs/0Q2xuy7ffbntvJPFPM9mXbj2Vfp1/Yfl+3suXJ1zLufNthu2uXleXJZvvT2ddvs+3vlSWb7ffavsv2g9n3dlEXs33T9vO2H51gu23/e5b9EdsnlyjbZ7NMv7J9r+2/7Fa2tiKi1h+SVklanr1eLunqNmOOlzQve/0eSc9KOqJDeWZI+rWk4yQdLOlhSfPHjfmSpP/KXg9I2tClr1WebGdJ+rPs9Re7lS1vvmzc4ZLukbRRUqMs2STNk/SgpHdky+8qUba1kr6YvZ4v6Tdd/L7+taSTJT06wfZFkm6VZEmnSLq/RNlOa/l+9nczW7uP2u+hS1oi6dvZ629L+pvxAyLiiYh4Mnv9jKTnJXXqoacLJI1ExLaI2CXp+izjRJlvlHS2bXcoz5SyRcRdEbEzW9yo5iMJuyXP106Svibpakl/Klm2L0haExEvSVJEPF+ibCHpbdnrt0t6pkvZFBH3qPlgnIkskfSdaNoo6QjbR5chW0Tcu+f7qe7/POyjFwr93RHxbPb6fyW9e3+DbS9Qcy/m1x3Kc4ykp1uWR7N1bcdExG5JOyS9s0N5ppqt1UVq7jl1y6T5sj/H50TET7uYS8r3tTte0vG2f2l7o+2FJcr2VUnLbI+q+TCby7oTLZep/n+ZSrd/HvaR64lFZWf755KOarPpqtaFiAjbE16nmf3W/66kCyPi9WJT1ovtZZIaks5InWUP2wdJulbS5xNHmchMNQ+7nKnmntw9tj8UEf+XNFXTUknfioh/tX2qmo+U/CA/B/nYPkvNQv+rlDlqUegRcc5E22w/Z/voiHg2K+y2f+bafpukn0q6KvuzrlO2S5rTsjw7W9duzKjtmWr+CfxiBzNNJZtsn6PmL8szIuKVLuTaY7J8h0v6oKS7syNUR0kasr04IoYTZ5Oae5b3R8Srkp6y/YSaBb+pBNkukrRQkiLiPtuHqjn5VLcOC+1Prv8vU7H9YUnfkNQfEd34OZ1QLxxyGZJ0Yfb6Qkn/PX6A7YMl/UjN43Q3djjPJknzbB+bve9AlrFVa+YLJN0Z2VmX1NlsnyTp65IWd/EYcK58EbEjImZFxNyImKvmMc1ulPmk2TI/VnPvXLZnqXkIZltJsv1O0tlZthMlHSpprAvZ8hiS9LnsapdTJO1oOYyalO33SvqhpL+LiCdS50l2NrZbH2oee/6FpCcl/VzSkdn6hqRvZK+XSXpV0kMtHx/pYKZFkp5Q8zj9Vdm6lWqWj9T8YfqBpBFJ/yPpuC5+vSbL9nNJz7V8nYa6/P3cb75xY+9Wl65yyfm1s5qHhB6T9CtJAyXKNl/SL9W8AuYhSZ/oYrbvq3ll2atq/hVzkaRLJF3S8nVbk2X/VZe/p5Nl+4akl1p+Hoa7la3dB7f+A0BN9MIhFwDoCRQ6ANQEhQ4ANUGhA0BNUOgAUBMUOgDUBIUOADXx/xwaN++MpdjHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "rng = RandomState(1729)\n",
    "\n",
    "x = rng.uniform(low=0, high=1, size=n)\n",
    "y = rng.uniform(low=0, high=1, size=n)\n",
    "sum_of_squares = x * x + y * y\n",
    "mask = sum_of_squares < 1\n",
    "\n",
    "pi_estimate = (np.sum(x * x + y * y < 1) / n) * 4\n",
    "print(f\"Approximate pi: {pi_estimate}\")\n",
    "print(f\"NumPy's value for pi: {np.pi}\")\n",
    "\n",
    "plt.scatter(x[mask], y[mask], s=1, c=\"red\")\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## With For-Loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Approximate pi: 3.139312\n",
      "NumPy's value for pi: 3.141592653589793\n"
     ]
    }
   ],
   "source": [
    "rando = RandomState(1729)\n",
    "count = 0\n",
    "\n",
    "for i in range(n):\n",
    "    u = rando.uniform(low=0, high=1)\n",
    "    v = rando.uniform(low=0, high=1)\n",
    "    if (u * u + v * v) < 1:\n",
    "        count += 1\n",
    "\n",
    "pi_est = (count / n) * 4\n",
    "print(f\"Approximate pi: {pi_est}\")\n",
    "print(f\"NumPy's value for pi: {np.pi}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of Vectorization Versus For-Loop by Using Timeit Library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Vectorized average time over 100 repeats: 0.031701\n",
      "For-loop average time over 100 repeats: 1.650312\n"
     ]
    }
   ],
   "source": [
    "repeats = 100\n",
    "\n",
    "# Setup statement executed once\n",
    "setup = \"\"\"\\\n",
    "import numpy as np\n",
    "from numpy.random import RandomState\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "import timeit\n",
    "\n",
    "n = 1_000_000\n",
    "rng_x = RandomState(1729)\n",
    "rng_y = RandomState(88)\"\"\"\n",
    "\n",
    "# Actual statement executed for repeats times\n",
    "vectorized_stmt = \"\"\"\\\n",
    "x = rng_x.uniform(low=0, high=1, size=n)\n",
    "y = rng_y.uniform(low=0, high=1, size=n)\n",
    "sum_of_squares = x * x + y * y\n",
    "mask = sum_of_squares < 1\"\"\"\n",
    "\n",
    "vectorized_time = timeit.Timer(stmt=vectorized_stmt, setup=setup).timeit(number=repeats)\n",
    "print(f\"Vectorized average time over {repeats} repeats: {vectorized_time/repeats:.6f}\")\n",
    "\n",
    "# Actual statement executed for repeats times\n",
    "for_stmt = \"\"\"\\\n",
    "count = 0\n",
    "for i in range(n):\n",
    "    u = rng_x.uniform(low=0, high=1)\n",
    "    v = rng_y.uniform(low=0, high=1)\n",
    "    if (u * u + v * v) < 1:\n",
    "        count += 1\n",
    "\"\"\"\n",
    "\n",
    "for_time = timeit.Timer(stmt=for_stmt, setup=setup).timeit(number=repeats)\n",
    "print(f\"For-loop average time over {repeats} repeats: {for_time/repeats:.6f}\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}