13-Random-Numbers-and-Stochastic-Simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 13.01 Stochastic Simulation and Randomness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stochastic Simulation** (def.) Mimic the behavior of a system by exploiting randomness to obtain a statistical sample of possible outcomes.\n",
    "\n",
    "Simulation Use Cases\n",
    "1. Nondeterministic systems\n",
    "2. Systems too complicated to model analytically\n",
    "3. Deterministic systems with high dimensionality.\n",
    "  * Example: Monte Carlo integration\n",
    "\n",
    "**Random** (def.) A sequence is random if it has no shorter description than itself.\n",
    "\n",
    "**Quasi-Random** (def.) A sequence which provides uniform coverage while maintaining reasonablly random appearance.\n",
    "  * Increasingly favored over uniform random number generators for Monte Carlo methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use Monte Carlo integration to estimate the definite integral $\\int_{-2}^2 (e^{-x^2} + 1) dx \\approx 5.76416$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ifx:  5.774353087151933\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import scipy.integrate as integrate\n",
    "\n",
    "def monte_carlo_integration(fx, a, b, n):\n",
    "    \"\"\"\n",
    "    Return definite integral of fx from [a, b] using Monte Carlo.\n",
    "\n",
    "    fx is the function to integrate. \n",
    "    a, b are the limits of integration.\n",
    "    n is the number of random samples.\n",
    "    \"\"\"\n",
    "    x = a + (b-a)*np.random.random_sample(n)\n",
    "    vfx = np.vectorize(fx)\n",
    "    return (b-a)*np.average(vfx(x))\n",
    "\n",
    "\n",
    "fx = lambda x: math.exp(-1.*x*x) + 1.\n",
    "a, b, n = -2., 2., np.power(10, 5)\n",
    "Ifx = monte_carlo_integration(fx, a, b, n)\n",
    "print(\"Ifx: \", Ifx)\n",
    "\n",
    "# Compare to numeric quadrature.\n",
    "expected, _  = integrate.quad(fx, a, b)\n",
    "np.testing.assert_almost_equal(Ifx, expected, decimal=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 13.02 Random Number Generators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Properties of a Good Random Number Generator\n",
    "1. Pass statistical tests of randomness.\n",
    "2. Long period.\n",
    "    * Period (def.) Length of sequence before starts repeating.\n",
    "3. Efficient.\n",
    "    * Execute quickly and hold little state.\n",
    "4. Repeatability.\n",
    "    * Produces repeatable sequence starting from same seed.\n",
    "5. Portability.\n",
    "    * Produce the same sequence from same seed on different computers.\n",
    "    \n",
    "#### Congruential Generators\n",
    "A congruential generator produces a sequence of integers in interval $[0, M-1]$:\n",
    "$$\n",
    "x_k = (a x_{k-1} + b) (\\bmod M)\n",
    "$$\n",
    "where\n",
    "* $a, b$ are integer constants associated with the generator.\n",
    "* $x_{k-1}$ is the previous value in the sequence or a seed if new sequence.\n",
    "* $M$ is an integer constant equal to the largest representable integer.\n",
    "\n",
    "Period of congruential generator cannot exceed $M$.\n",
    "\n",
    "#### Fibonacci Generators\n",
    "Produce floating point random numbers on interval $[0, 1]$ as difference or sum of previous values.\n",
    "$$\n",
    "x_k = x_{k-M} - x_{k-P}\n",
    "$$\n",
    "where\n",
    "* $M$ and $P$ are referred to as **lag**\n",
    "\n",
    "Fibonacci generators require special initialization in order to compute first terms without lag.\n",
    "\n",
    "Fibonacci generators have **longer** period than congruential generators since repitition of one sequence does not entail all subsequent members will repeat in same order.\n",
    "\n",
    "#### Nonuniform Distributions\n",
    "A nonuniform generator can be built from a uniform generator:\n",
    "1. Compute the inverse of the PDF of the distribution.\n",
    "2. Pass the uniform random deviate to the inverse."
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}