{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Computing $\\pi$ using Sampling"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "nsamples = 3000"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start with a uniform distribution on the unit square $[0,1]^2$.\n",
        "\n",
        "Create a 2D array `samples` of shape `(2, nsamples)`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "samples = np.random.rand(2, nsamples)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot them:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.plot(samples[0], samples[1], \"o\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we would like to find $\\int_0^1\\int_0^1 p(x,y) dx dy$ where $p(x,y)$ is 1 if $(x,y)$ is in the unit circle around the origin, otherwise zero.\n",
        "\n",
        "$p(x)$ is not a probability distribution. But with a scaling factor $\\alpha$ it can be one, so that\n",
        "\n",
        "$$\\frac1{\\alpha} \\int_0^1\\int_0^1 p(x,y) dx dy=1.$$\n",
        "\n",
        "So it is $\\alpha=\\pi/4$ we're looking for. We can compute it by\n",
        "\n",
        "$$\\alpha=\\int_{0}^1\\int_{0}^1 p(x,y) dx dy\\approx \\sum_{i=1}^N p(x_i)$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "----------------------\n",
        "Next, evaluate $p$ for each sample:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "p = np.zeros(nsamples)\n",
        "\n",
        "r = samples[0]**2 + samples[1]**2\n",
        "p[r<=1] = 1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the value of $p$ for the samples:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.plot(samples[0, p==1], samples[1, p==1], \"o\", color=\"red\")\n",
        "plt.plot(samples[0, p==0], samples[1, p==0], \"o\", color=\"blue\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Evaluate $\\alpha$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "\n",
        "alpha = np.sum(p) / nsamples\n",
        "alpha"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "my_pi = 4*alpha\n",
        "my_pi"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": true
      },
      "outputs": [],
      "source": []
    }
  ],
  "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.5.2"
    },
    "widgets": {
      "state": {},
      "version": "1.1.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}