10_MonteCarlo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline\n",
    "%precision 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Psedorandom number generators\n",
    "\n",
    "While psuedorandom numbers are generated by a deterministic algorithm, we can mostly treat them as if they were true random numbers and we will drop the “pseudo” prefix. Fundamentally, the algorithm generates random integers which are then normalized to give a floating point number from the standard uniform distribution. Random numbers from other distributions are in turn generated using these uniform random deviates, see later.\n",
    "\n",
    "## Linear congruential generators (LCG)\n",
    "\n",
    "[LCG](https://en.wikipedia.org/wiki/Linear_congruential_generator) is among the simplest and most popular pseudo random number generators. It relies on the recursive and fully deterministic relation:\n",
    "\n",
    "$$\n",
    "z_{i+1}=(a z_i+c)\\mod{m}\n",
    "$$\n",
    "\n",
    "Hull-Dobell Theorem: The LCG will have a period $m$ for all seeds if and only if\n",
    "\n",
    "* $c$ and $m$ are relatively prime,\n",
    "* $a−1$ is divisible by all prime factors of $m$\n",
    "* $a−1$ is a multiple of 4 if $m$ is a multiple of $4$.\n",
    "\n",
    "The number $z_0$ is called the *seed*, and setting it allows us to have a reproducible sequence of (pseudo) random numbers. The LCG is typically coded to return $z/m$, a floating point number in $(0, 1)$. Obviosuly, this can be easily scaled to any other range $(a,b)$.\n",
    "Note that $z \\le m-1$ always holds, the yielded $z/m$ result is thus on purpose strictly smaller than 1. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lcg(m=2**32, a=1103515245, c=12345):\n",
    "    lcg.current = (a*lcg.current + c) % m\n",
    "    return lcg.current/m\n",
    "\n",
    "# setting the seed\n",
    "lcg.current = 12345"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[lcg() for i in range(10)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rn=[lcg() for i in range(1000)]\n",
    "print (np.mean(rn))\n",
    "print (np.std(rn),1/np.sqrt(12))\n",
    "plt.plot(rn,\"o\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LCG though is not sufficiently \"random\" for several complex modern applications. There are nowadays better performing algorithms, like [Mersenne twister](https://en.wikipedia.org/wiki/Mersenne_Twister), a generalized feedback shift-register generator, is used, in particular the numpy random package features it.\n",
    "\n",
    "Numpy uses as default [PCG-64](https://numpy.org/doc/stable/reference/random/bit_generators/index.html) of the [PCG family](https://www.pcg-random.org/), which are considered the ultimate random number generators.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python libraries for Random Number Generation \n",
    "\n",
    "There are two modules for (pseudo) random numbers that are commonly used. When all you need is to generate random numbers from some distribution, the `numpy.random` moodule is the simplest to use. When you need more information realted to a distribution such as quantiles or the PDF, you can use the `scipy.stats` module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.random as npr\n",
    "npr.seed(123) # fix seed for reproducible result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# random values from a uniform distrubtion in [0.0, 1.0) for a given shape\n",
    "print (\"uniform:\",'\\n', npr.rand(3,4),'\\n') # shape=(3,4)\n",
    " \n",
    "# random values from a standard normal distrubtion  for a given shape\n",
    "print (\"normal:\",'\\n', npr.randn(2,5),'\\n') # shape=(2,5)\n",
    "\n",
    "# random integers between \"low\" and \"high\" edges for a given shape (size)\n",
    "print (\"integers:\",'\\n', npr.randint(1,50,size=(3,6)),'\\n') # shape=(3,6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Operations on arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# randomly shuffling a vector\n",
    "x = np.arange(10)\n",
    "npr.shuffle(x)\n",
    "print (\"reshuffling:\",'\\n', x,'\\n')\n",
    "# npr.permutation works similarly\n",
    "\n",
    "# Generates a random sample from a given 1-D array (without replacement)\n",
    "x = np.arange(10,20)\n",
    "print (\"choice w/o replace\",'\\n',npr.choice(x, 10, replace=False),'\\n')\n",
    "\n",
    "# Generates a random sample of a given dimension from a given 1-D array\n",
    "print (\"choice w/ replace\",'\\n',npr.choice(x, (5, 10), replace=True),'\\n') # this is default"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probability Density Functions\n",
    "\n",
    "numbers can be drawn from basically any common PDF (normal, binomial, chisquare, etc..)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# exponential\n",
    "print (\"exponential:\",'\\n', npr.exponential(10, size=(3,2)),'\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# tossing a dice\n",
    "roll = 1./6\n",
    "x = npr.multinomial(20000, [roll]*6, size=2)\n",
    "print (\"frequencies of 6-faces dice tossing (two players):\",'\\n', x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's pretend we haven't seen any of the above functionalties available in numpy and scipy and walk through the basics.."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Non-uniform random numbers\n",
    "\n",
    "In several cases the actual random process occur with non-uniform probability, i.e. with a given probability density function (pdf), different from the uniform distribution. Several methods are available, we will see a few of them\n",
    "\n",
    "### Inverst transform method\n",
    "\n",
    "Let'start from a uniform distribution $u(z)$:\n",
    "\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{array}{ll}\n",
    "      1 & 0\\leq z\\leq 1 \\\\\n",
    "      0 & {\\rm elsewhere}\n",
    "\\end{array} \n",
    "\\right.\n",
    "$$\n",
    "\n",
    "and let's look for a function $x(z)$ such that $x$ is distributed accordingly to a given pdf $p(x)$. The probability to find $x$ between $x$ and $x+dx$ is equal to:\n",
    "\n",
    "$$\n",
    "p(x)dx = dz\n",
    "$$\n",
    "\n",
    "and thus:\n",
    "\n",
    "$$\n",
    "\\int_{-\\infty}^{x(z)} p(x') dx' = \\int_0^z dz'= z\n",
    "$$\n",
    "\n",
    "If (a) we could solve the integral and (b) solve for $x$, then we are done. For most of the pdf at least one of the two is not possible.. The typical solvable analitical example is:\n",
    "\n",
    "$$\n",
    "p(x) = \\mu e^{-\\mu x}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\int_{0}^{x(z)} p(x') dx' = 1 - e^{-\\mu x} = z\n",
    "$$\n",
    "\n",
    "and thus:\n",
    "\n",
    "$$\n",
    "x(z) = - \\frac{1}{\\mu}\\log{(1-z)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def expon_pdf(x, mu=1):\n",
    "    \"\"\"PDF of exponential distribution.\"\"\"\n",
    "    return mu*np.exp(-mu*x)\n",
    "\n",
    "def expon_cdf(x, mu=1):\n",
    "    \"\"\"CDF of exponetial distribution.\"\"\"\n",
    "    return 1 - np.exp(-mu*x)\n",
    "\n",
    "def expon_icdf(z, mu=1):\n",
    "    \"\"\"Inverse CDF of exponential distribution - i.e. quantile function.\"\"\"\n",
    "    return -np.log(1-z)/mu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist = stats.expon()\n",
    "x = np.linspace(0,4,100)\n",
    "y = np.linspace(0,1,100)\n",
    "\n",
    "plt.figure(figsize=(12,4))\n",
    "plt.subplot(121)\n",
    "plt.plot(x, expon_cdf(x))\n",
    "plt.axis([0, 4, 0, 1])\n",
    "for q in [0.5, 0.8]:\n",
    "    plt.arrow(0, q, expon_icdf(q)-0.1, 0, head_width=0.05, head_length=0.1, fc='b', ec='b')\n",
    "    plt.arrow(expon_icdf(q), q, 0, -q+0.1, head_width=0.1, head_length=0.05, fc='b', ec='b')\n",
    "plt.ylabel('1: Generate a (0,1) uniform PRNG')\n",
    "plt.xlabel('2: Find the inverse CDF')\n",
    "plt.title('Inverse transform method');\n",
    "\n",
    "plt.subplot(122)\n",
    "u = np.random.random(10000)\n",
    "v = expon_icdf(u)\n",
    "plt.hist(v, histtype='step', bins=100, density=True, linewidth=2)\n",
    "plt.plot(x, expon_pdf(x), linewidth=2)\n",
    "plt.axis([0,4,0,1])\n",
    "plt.title('Histogram of exponential PRNGs');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Box-Muller for generating normally distributed random numbers\n",
    "\n",
    "The inverse method is not applicable even for the gaussian distribution:\n",
    "\n",
    "$$\n",
    "\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\int_{-\\infty}^{x} \\exp{-\\frac{x'^2}{2\\sigma^2}} dx' = z\n",
    "$$\n",
    "\n",
    "is not solvable. The trick is to consider a two dimensional gaussian function with the same $\\sigma$ on both coordinates:\n",
    "\n",
    "$$\n",
    "p(x)dx \\times p(y)dy = \\frac{1}{\\sqrt{2\\pi\\sigma^2}} \\exp{-\\frac{x^2}{2\\sigma^2}} \\times \\frac{1}{\\sqrt{2\\pi\\sigma^2}} \\exp{-\\frac{y^2}{2\\sigma^2}} = \\frac{1}{2\\pi\\sigma^2} \\exp{-\\frac{(x^2+y^2)}{2\\sigma^2}}dxdy\n",
    "$$\n",
    "\n",
    "which written in radial coordinates:\n",
    "\n",
    "$$\n",
    "x=r\\cos{\\theta};\\,\\,\\,\n",
    "y=r\\sin{\\theta}\n",
    "$$\n",
    "\n",
    "$$\n",
    "p(r,\\theta)dr d\\theta = \\frac{r}{\\sigma^2} \\exp{-\\frac{r^2}{2\\sigma^2}} dr \\times \\frac{d\\theta}{2\\pi} = p(r)dr \\times p(\\theta)d\\theta\n",
    "$$\n",
    "\n",
    "with both $p(r)$ and $p(\\theta)$ normalized to 1. Now, the latter is a simple uniform distribution, whereas the former is solvable:\n",
    "\n",
    "$$\n",
    "\\frac{1}{\\sigma^2} \\int_{0}^{r} \\exp{-\\frac{r^2}{2\\sigma^2}} rdr = z\n",
    "$$\n",
    "\n",
    "which gives:\n",
    "\n",
    "$$\n",
    "r=\\sqrt{-2\\sigma^2\\log{1-z}}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10000\n",
    "z = np.random.random(n)\n",
    "theta = 2*np.pi*np.random.random(n)\n",
    "r_squared = -2*np.log(z)\n",
    "r = np.sqrt(r_squared)\n",
    "data = pd.DataFrame({'x':r*np.cos(theta), 'y':r*np.sin(theta)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.jointplot(data=data, x='x',y='y', kind=\"hist\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating a random number generator for arbitrary distributions\n",
    "\n",
    "Suppose we have some random samples with an unknown distribution. We can still use the inverse transform method to create a random number generator from a random sample, by estimating the inverse CDF function using interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import interp1d\n",
    "\n",
    "def extrap1d(interpolator):\n",
    "    \"\"\"From StackOverflow http://bit.ly/1BjyRfk\"\"\"\n",
    "    xs = interpolator.x\n",
    "    ys = interpolator.y\n",
    "\n",
    "    def pointwise(x):\n",
    "        if x < xs[0]:\n",
    "            return ys[0]+(x-xs[0])*(ys[1]-ys[0])/(xs[1]-xs[0])\n",
    "        elif x > xs[-1]:\n",
    "            return ys[-1]+(x-xs[-1])*(ys[-1]-ys[-2])/(xs[-1]-xs[-2])\n",
    "        else:\n",
    "            return interpolator(x)\n",
    "\n",
    "    def ufunclike(xs):\n",
    "        return np.array(list(map(pointwise, np.array(xs))))\n",
    "\n",
    "    return ufunclike"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.distributions.empirical_distribution import ECDF\n",
    "\n",
    "# Make up some random data\n",
    "x = np.concatenate([np.random.normal(0, 1, 10000),\n",
    "                    np.random.normal(4, 1, 10000)])\n",
    "\n",
    "ecdf = ECDF(x)\n",
    "inv_cdf = extrap1d(interp1d(ecdf.y, ecdf.x,\n",
    "                            bounds_error=False, assume_sorted=True))\n",
    "r = np.random.uniform(0, 1, 1000)\n",
    "ys = inv_cdf(r)\n",
    "\n",
    "plt.hist(x, 25, histtype='step', color='red', density=True, linewidth=1)\n",
    "plt.hist(ys, 25, histtype='step', color='blue', density=True, linewidth=1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rejection sampling (Accept-reject method)\n",
    "\n",
    "Suppose we want to sample from a (truncated) Cauchy distribution (a Student-t distribution with 1 degree of freedom).\n",
    "We use the uniform as a proposal distibution (highly inefficient)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-4, 4)\n",
    "\n",
    "dist = stats.cauchy()\n",
    "upper = dist.pdf(0)\n",
    "\n",
    "plt.figure(figsize=(12,4))\n",
    "plt.subplot(121)\n",
    "plt.plot(x, dist.pdf(x))\n",
    "plt.axhline(upper, color='grey')\n",
    "px = 1.0\n",
    "plt.arrow(px,0,0,dist.pdf(1.0)-0.01, linewidth=1,\n",
    "         head_width=0.2, head_length=0.01, fc='g', ec='g')\n",
    "plt.arrow(px,upper,0,-(upper-dist.pdf(px)-0.01), linewidth=1,\n",
    "         head_width=0.3, head_length=0.01, fc='r', ec='r')\n",
    "plt.text(px+.25, 0.2, 'Reject', fontsize=16)\n",
    "plt.text(px+.25, 0.01, 'Accept', fontsize=16)\n",
    "plt.axis([-4,4,0,0.4])\n",
    "plt.title('Rejection sampling concepts', fontsize=20)\n",
    "\n",
    "plt.subplot(122)\n",
    "n = 100000\n",
    "# generate from sampling distribution\n",
    "u = np.random.uniform(-4, 4, n)\n",
    "# accept-reject criterion for each point in sampling distribution\n",
    "r = np.random.uniform(0, upper, n)\n",
    "# accepted points will come from target (Cauchy) distribution\n",
    "v = u[r < dist.pdf(u)]\n",
    "\n",
    "plt.plot(x, dist.pdf(x), linewidth=2)\n",
    "\n",
    "# Plot scaled histogram\n",
    "factor = dist.cdf(4) - dist.cdf(-4)\n",
    "hist, bin_edges = np.histogram(v, bins=100, density=True)\n",
    "bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.\n",
    "plt.step(bin_centers, factor*hist, linewidth=2)\n",
    "\n",
    "plt.axis([-4,4,0,0.4])\n",
    "plt.title('Histogram of accepted samples', fontsize=20);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mixture representations\n",
    "Sometimee, the target distribution  which we need to generate random numbers from can be expressed as a mixture of “simpler” distributions that we already know how to sample from\n",
    "\n",
    "$$\n",
    "f(x)=\\int g(x|y)p(y)dy\n",
    "$$\n",
    "\n",
    "For example, if $y$ is drawn from the $\\chi^2_\\nu$ distrbution, then ${\\cal N}(0,\\nu/y)$ is a sample from the Student-t distribution with $\\nu$ degrees fo freedom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10000\n",
    "df = 5\n",
    "dist = stats.t(df=df)\n",
    "y = stats.chi2(df=df).rvs(size=n)\n",
    "r = stats.norm(0, df/y).rvs(n)\n",
    "\n",
    "\n",
    "plt.plot(x, dist.pdf(x), linewidth=2)\n",
    "\n",
    "# Plot scaled histogram\n",
    "factor = dist.cdf(4) - dist.cdf(-4)\n",
    "hist, bin_edges = np.histogram(v, bins=100, density=True)\n",
    "bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.\n",
    "plt.step(bin_centers, factor*hist, linewidth=2)\n",
    "\n",
    "plt.axis([-4,4,0,0.4])\n",
    "plt.title('Histogram of accepted samples', fontsize=20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Draw from an analytic pdf\n",
    "\n",
    "Obviously scipy stats module features all possible pdf that can come to your mind. You can draw random data from each of them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import gamma\n",
    "a = 1.99\n",
    "\n",
    "x = np.linspace(gamma.ppf(0.01, a),\n",
    "                gamma.ppf(0.99, a), 100)\n",
    "rv = gamma(a)\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(x, rv.pdf(x), 'k-', lw=2)\n",
    "r = gamma.rvs(a, size=10000)\n",
    "_ = ax.hist(r, density=True, histtype='stepfilled', alpha=0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monte Carlo integration\n",
    "\n",
    "Monte Carlo integration is typically less accurate than other integration methods, but very often is the only available tool, e.g. when the integrand has very rapid variations or singular points, or, most importantly, when dealing with high dimensional integrals.\n",
    "\n",
    "The idea is simple, let's the area under the function be $I$ whereas the all possible outcomes lay in a box of area $A$. The probability for a point to fall under the function is $p=I/A$. If we generate $N$ random points, the fraction $k$ which fall under the curve is $k/N$ and approximate $I/A$, thus:\n",
    "\n",
    "$$\n",
    "I\\simeq\\frac{k A}{N}\n",
    "$$\n",
    "\n",
    "Let's try this with the function $f(x) =\\sin^2{\\frac{1}{x(2-x)}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return (np.sin(1/(x*(2-x))))**2\n",
    "\n",
    "x=np.linspace(0.001,1.999,1000)\n",
    "plt.plot(x,f(x),'r-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monte Carlo integration\n",
    "N=100000\n",
    "count=0\n",
    "for i in range(N):\n",
    "    x=2*np.random.random()\n",
    "    y=np.random.random()\n",
    "    if y<f(x): count+=1\n",
    "I=2*count/N\n",
    "print(I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The mean value method\n",
    "\n",
    "Let's take the integral:\n",
    "\n",
    "$$\n",
    "I=\\int_a^b f(x) dx\n",
    "$$\n",
    "\n",
    "defining $\\langle f \\rangle$ as the mean of $f$:\n",
    "\n",
    "$$\n",
    "\\langle f \\rangle = \\frac{1}{b-a}\\int_a^b f(x) dx\n",
    "$$\n",
    "\n",
    "and estimating $\\langle f \\rangle$ by uniformely probing at random the function domain, such as\n",
    "\n",
    "$$\n",
    "\\langle f \\rangle = \\frac{1}{N} \\sum_{i=1}^{N} f(x_i)\n",
    "$$\n",
    "\n",
    "we get:\n",
    "\n",
    "$$\n",
    "I=\\frac{b-a}{N} \\sum_{i=1}^{N} f(x_i)\n",
    "$$\n",
    "\n",
    "this easily generalize to higher dimensions:\n",
    "\n",
    "$$\n",
    "I=\\frac{V}{N} \\sum_{i=1}^{N} f(\\vec{r}_i)\n",
    "$$\n",
    "\n",
    "where the sampling points $\\vec{r}_i$ are drawn uniformly at random from integration space of volume $V$.\n",
    "\n",
    "It can be proven that the standard deviation of the method scales $1/\\sqrt{N}$:\n",
    "\n",
    "$$\n",
    "\\sigma = V\\frac{\\sqrt{{\\rm var}\\, f}}{\\sqrt{N}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importance sampling\n",
    "\n",
    "There are several general techinques for variance reduction, sometimes known as Monte Carlo swindles since these methods improve the accuracy and convergence rate of Monte Carlo integration without increasing the number of Monte Carlo samples. *Importance sampling* is among the most commonly used.\n",
    "\n",
    "We can define a weighted average of a function $g(x)$:\n",
    "\n",
    "$$\n",
    "\\langle g_w \\rangle = \\frac{\\int_a^b w(x) g(x)dx}{\\int_a^b w(x)dx}\n",
    "$$\n",
    "\n",
    "Consider again the integral of $f(x)$:\n",
    "\n",
    "$$\n",
    "I=\\int_a^b f(x)dx\n",
    "$$\n",
    "\n",
    "Setting $g(x)=f(x)/w(x)$ we have:\n",
    "\n",
    "$$\n",
    "\\left\\langle \\frac{f(x)}{w(x)}\\right\\rangle = \\frac{\\int_a^b w(x)f(x)/w(x) )dx}{\\int_a^b w(x)dx} = \\frac{I}{\\int_a^b w(x)dx}\n",
    "$$\n",
    "\n",
    "and thus:\n",
    "\n",
    "$$\n",
    "I = \\left\\langle \\frac{f(x)}{w(x)}\\right\\rangle \\int_a^b w(x)dx \\simeq \\frac{1}{N}\\sum_{i=1}^N \\frac{f(x_i)}{w(x_i)}\\int_a^b w(x)dx\n",
    "$$\n",
    "\n",
    "which generalizes the mean value method if $w(x)$ is the uniform distribution between $a$ and $b$\n",
    "\n",
    "### Example\n",
    "\n",
    "Suppose we want to estimate the tail probability of ${\\cal N}(0,1)$ for $P(x>5)$. Regular MC integration using samples from ${\\cal N}(0,1)$ is hopeless since nearly all samples will be rejected. However, we can use the exponential density truncated at 5 as the importance function and use importance sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(4, 8, 100)\n",
    "plt.plot(x, stats.expon(5).pdf(x))\n",
    "plt.plot(x, stats.norm().pdf(x))\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%precision 10\n",
    "h_true =1 - stats.norm().cdf(5)\n",
    "h_true"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 1000000\n",
    "y = stats.norm().rvs(n)\n",
    "print (\"integral (sum) above 5:\", np.sum(y > 5))\n",
    "h_mc = 1.0/n * np.sum(y > 5)\n",
    "# estimate and relative error\n",
    "print (\"estimate:\", h_mc)\n",
    "print (\"relative error:\", np.abs(h_mc - h_true)/h_true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10000\n",
    "y = stats.expon(loc=5).rvs(n)\n",
    "h_is = 1.0/n * np.sum(stats.norm().pdf(y)/stats.expon(loc=5).pdf(y))\n",
    "# estimate and relative error\n",
    "print (\"estimate:\", h_is)\n",
    "print (\"relative error:\", np.abs(h_is - h_true)/h_true)"
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}