fundamental-theorem-of-simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fundamental Theorem of Simulation\n",
    "The Fundamental Theorem of Simulation aka. Universality of the Uniform tells us that we can sample random-variables of any distribution, given its inverse CDF (or approximation) and a standard uniform rv. Basically $F^{-1}(U)=X$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sample binomial random-variable\n",
    "Here we have to use an approximation of the inverse CDF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.special import comb\n",
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def binom(n, k, p):\n",
    "    return comb(n, k) * np.power(p, k) * np.power(1-p, n-k)\n",
    "  \n",
    "def sample_binom(n, p):\n",
    "    k_prob = []\n",
    "    \n",
    "    for k in range(0, n+1):\n",
    "        k_prob.append(binom(n, k, p))\n",
    "\n",
    "    k_prob_cdf = 0\n",
    "    U = np.random.rand() # uniform rv.\n",
    "    for X, k_p in enumerate(k_prob):\n",
    "        # k_prob_cdf is the value of the binomial CDF at X (so P(X < X))\n",
    "        k_prob_cdf += k_p\n",
    "        if U <= k_prob_cdf:\n",
    "            return X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean of a rv $X \\sim{} Bin(n,p)$ is $np$ and its variance is $np(1-p)$. So lets check with a simple simulation if our sampling algorithm produces valid rvs. according to $Bin(n,p)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulated mean: 5.0037 mathematical mean: 5.0\n",
      "Simulated variance: 3.715 mathematical variance: 3.75\n"
     ]
    }
   ],
   "source": [
    "n_simulations = 10000\n",
    "\n",
    "n = 20\n",
    "p = 0.25\n",
    "mean = n*p\n",
    "variance = n*p*(1-p)\n",
    "\n",
    "x_sum = 0\n",
    "sim_var = 0\n",
    "for i in range(n_simulations):\n",
    "    x = sample_binom(n, p)\n",
    "    x_sum += x\n",
    "    sim_var += np.power(mean-x, 2) / n_simulations\n",
    "    \n",
    "sim_mean = x_sum / n_simulations\n",
    "\n",
    "print(\"Simulated mean:\", sim_mean, \"mathematical mean:\", mean)\n",
    "print(\"Simulated variance:\", np.around(sim_var,3), \"mathematical variance:\", variance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sample standard normal random-variable\n",
    "Since the inverse of $\\Phi(Z)$ is not defined, we will again have to use approximations. Just like students use the z-score table. So for us $\\Phi^{-1}(U)=Z$, where $Z$ is a standard normal rv.\n",
    "\n",
    "First we create the inverse table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 310,
   "metadata": {},
   "outputs": [],
   "source": [
    "inv_table = {}\n",
    "X_max = 10.0\n",
    "X_step = 0.001\n",
    "\n",
    "X_range = np.linspace(-X_max, X_max, int(1 / X_step))\n",
    "X_CDF_prob = norm.cdf(X_range)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 311,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create inverse-phi-table\n",
    "for i in range(len(X_range)):\n",
    "    inv_table[X_CDF_prob[i]] = X_range[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 312,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_normal(mu=0, sigma=1):\n",
    "    U = np.random.rand()\n",
    "    \n",
    "    for i in range(len(X_range)-1):\n",
    "        if X_CDF_prob[i] < U < X_CDF_prob[i+1]:\n",
    "            if np.abs(U - X_CDF_prob[i]) < np.abs(U - X_CDF_prob[i+1]):\n",
    "                return mu + sigma * X_range[i]\n",
    "            return mu + sigma * X_range[i+1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 313,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6306306306306304"
      ]
     },
     "execution_count": 313,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_normal()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simulate and check if rvs. are valid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 318,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_simulations = 10000\n",
    "\n",
    "X_mean = 0\n",
    "X_var = 0\n",
    "mu = 3\n",
    "sigma = 4\n",
    "\n",
    "for i in range(n_simulations):\n",
    "    z = sample_normal(mu=mu, sigma=sigma)\n",
    "    X_mean += z\n",
    "    X_var += np.power(3-z, 2) / n_simulations\n",
    "\n",
    "X_mean /= n_simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 320,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulated mean: 2.9991991991992006 mathematical mean: 3\n",
      "Simulated variance: 16.033124135145965 mathematical variance 16\n"
     ]
    }
   ],
   "source": [
    "print(\"Simulated mean:\", X_mean, \"mathematical mean:\", mu)\n",
    "print(\"Simulated variance:\", X_var, \"mathematical variance\", sigma**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Check ✓"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}