04-monte-carlo-simulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import library_data_science as lds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inferential Statictics\n",
    "\n",
    "A branch of statistics dealing with the problems of generalizing the results of a random sample study to the entire population and estimating the errors arising from such generalization.\n",
    "\n",
    "![Illustration of Inferential Statistics](https://datatab.net/assets/tutorial/Descriptive_statistics_and_inferential_statistics.png)\n",
    "\n",
    "**Population**: A set of examples.\n",
    "\n",
    "**Sample**: A proper subset of a population."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monte Carlo simulation\n",
    "\n",
    "A method of estimating the value of an unknown quantity using the principles of inferential statistics.\n",
    "\n",
    "1. Set number of **tests**.\n",
    "\n",
    "2. **Conduct** the set amount of tests.\n",
    "\n",
    "3. The result is the **average value** of all tests."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, in virtually every notebook, I have come across **Monte Carlo** simulations.\n",
    "\n",
    "* Probability that rolling five times we got combination `11111`.\n",
    "\n",
    "* The birthday problem.\n",
    "\n",
    "* Predicting the distance from the endpoint to the starting point using Random Walks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Law of Large Numbers\n",
    "\n",
    "It states that if you repeat an experiment independently a large number of times and average the result, what you obtain should be close to expected value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, we know that the mean of rolling a dice $n \\geq 2$ times equals $3.5$. We can check the truth of the **Law of Large Numbers** using Monte Carlo simulation. I will repeat the test four times for bigger amount of trials and sum up the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "dice = [1, 2, 3, 4, 5, 6]\n",
    "trials = [10, 100, 1000, 10000000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[4.9, 3.54, 3.417, 3.49947]\n"
     ]
    }
   ],
   "source": [
    "means = []\n",
    "\n",
    "for num_trials in trials:\n",
    "    results = []\n",
    "\n",
    "    for attempt in range(num_trials):\n",
    "        experiment = random.choice(dice)\n",
    "        results.append(experiment)\n",
    "    \n",
    "    current_mean = lds.mean(results)\n",
    "    means.append(current_mean)\n",
    "\n",
    "print(means)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|                 | Mean       |\n",
    "|-----------------|------------|\n",
    "| 10 Trials       | 4.9        |\n",
    "| 100 Trials      | 3.54       |\n",
    "| 1000 Trials     | 3.417      |\n",
    "| 10000000 Trials | 3.49947    |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculus methods using Monte Carlo simulation\n",
    "\n",
    "Instead, I will focus on a more practical application of the Monte Carlo simulation method. It is very useful for scientists who require computational methods for differential calculus in their calculations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Derivative\n",
    "\n",
    "The first function is the method for calculating the derivative at a given point. According to the mathematical formula, the derivative can be computed using the equation below.\n",
    "\n",
    "$$\\frac{d}{dx}\\bigg[ f(x_0) \\bigg] = \\lim_{\\Delta x \\to 0}{\\frac{f(x_0 + \\Delta x) - f(x_0)}{\\Delta x}}$$\n",
    "\n",
    "As we know, there are many more practical methods, but for computational purposes, it serves its role entirely. Writing a function to calculate the derivative at a point is not difficult."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def derivative(f, x_0: float, delta_x = 0.0001) -> float:\n",
    "    \"\"\"Function returns the derivative of the given function at specified point.\n",
    "    \n",
    "       Function must be continuous at specified point.\"\"\"\n",
    "    \n",
    "    return (f(x_0 + delta_x) - f(x_0)) / delta_x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integral\n",
    "\n",
    "The second function is the method for calculating the integral over a given interval. According to the mathematical formula, the integral can be computed using the riemann sum.\n",
    "\n",
    "$$\\int_{a}^{b}{f(x)}dx = \\lim_{n \\to +\\infty}{\\sum_{i=1}^{n}{f(c_i)\\Delta x}}$$\n",
    "\n",
    "Similarly, writing such a function in Python is not a significant challenge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integral(f, a: float, b: float, delta_x = 0.0001) -> float:\n",
    "    \"\"\"Function returns the integral of the given function at specified range.\n",
    "    \n",
    "       Function must be continuous at specified range.\"\"\"\n",
    "    \n",
    "    counter = 0\n",
    "\n",
    "    while a <= b:\n",
    "        counter = counter + (f(a) * delta_x)\n",
    "        a = a + delta_x\n",
    "    \n",
    "    return counter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, this approach has many limitations. It assumes that the function is continuous over the specified interval. Accounting for discontinuities can be very challenging for beginner programmers and even for brilliant mathematicians. Of course, there are many libraries that account for such function behaviors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Area using Monte Carlo simulation\n",
    "\n",
    "To estimate the area of some region $R$ follow the instruction.\n",
    "\n",
    "1. Pick an enclosing region $E$, such that the area of $E$ is easy to calculate and $R$ lies completely within $E$.\n",
    "\n",
    "2. Pick a set of random points that lie within $E$.\n",
    "\n",
    "3. Let $F$ be the fraction of the points that fall within $R$. $$F = \\frac{lying\\_points}{all\\_points}$$\n",
    "\n",
    "4. Multiply the area of $E$ by $F$. That is the result. $$R \\approx E \\times F$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def area_monte_carlo(f, a: float, b: float, c: float, d: float, num_samples = 1000000) -> float:\n",
    "    \"\"\"Function returns the integral of the given function at specified range.\n",
    "\n",
    "       Method used in the process is Monte Carlo Simulation.\n",
    "\n",
    "       Factors a, b are the range of x-axis, and factors c, d are the range of y-axis.\n",
    "    \n",
    "       Function must be continuous at specified range.\"\"\"\n",
    "\n",
    "    within_the_area = 0\n",
    "    for sample in range(num_samples):\n",
    "        x = random.uniform(a, b)\n",
    "        y = random.uniform(c, d)\n",
    "\n",
    "        try:\n",
    "            if lds.distance((x, 0), (x, y)) <= f(x):\n",
    "                within_the_area = within_the_area + 1\n",
    "        except ZeroDivisionError:\n",
    "            pass\n",
    "\n",
    "    estimated_area = (np.abs(a - b) * np.abs(c - d)) * (within_the_area / num_samples)\n",
    "\n",
    "    return estimated_area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example for $\\displaystyle\\int_{0}^{\\pi}{\\sin{(x)}}dx$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated area of sin(x) on [0, pi]: 1.999595449676021\n",
      "Real area of sin(x) on [0, pi]: 1.9999999986725736\n"
     ]
    }
   ],
   "source": [
    "f = lambda x: np.sin(x)\n",
    "num_samples = 100000\n",
    "a = 0\n",
    "b = np.pi\n",
    "c = 0\n",
    "d = 1\n",
    "\n",
    "print('Estimated area of sin(x) on [0, pi]:', area_monte_carlo(f, a, b, c, d))\n",
    "print('Real area of sin(x) on [0, pi]:', integral(f, a, b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also easily calculate the value of $\\pi$. To do this, we will compute the area of the upper part of a circle with a radius of 1. For this purpose, we use the function describing the upper half of the circle, $\\sqrt{1 - x^2}$ on range $[0, 1]$ and then multiply the result by $4$. This approach is called **Buffon-Laplace Method**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated value of pi: 3.141824\n",
      "Real area of pi: 3.141592653589793\n"
     ]
    }
   ],
   "source": [
    "f = lambda x: (1 - x ** 2) ** 0.5\n",
    "num_samples = 100000\n",
    "a = 0\n",
    "b = 1\n",
    "c = 0\n",
    "d = 1\n",
    "\n",
    "print('Estimated value of pi:', 4 * area_monte_carlo(f, a, b, c, d))\n",
    "print('Real area of pi:', np.pi)"
   ]
  }
 ],
 "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}