Overthinking a Monte Carlo Estimate of π

overthinking-mc-pi.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "602aa774-ddf6-4bcd-b219-f3734132cfe5",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "While preparing the first in an upcoming series of posts on [multi-armed bandits](https://en.wikipedia.org/wiki/Multi-armed_bandit), I realized that a post diving deep on a simple Monte Carlo estimate of $\\pi$ would be a useful companion, so here it is."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e6515b9-8be0-4036-9d90-b1db4bf21de5",
   "metadata": {},
   "source": [
    "## Monte Carlo Methods: the Standard Introductory Example\n",
    "\n",
    "[Monte Carlo methods](https://en.wikipedia.org/wiki/Monte_Carlo_method) are a powerful tool for approximating hard- or impossible-to-calculate quantities using randomness and simulation.  In my [talks](https://austinrochford.com/talks.html) on Bayesian computation with [PyMC](https://www.pymc.io/welcome.html), I like to follow convention by introducing Monte Carlo methods by approximating $\\pi$ as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "153b106e-9f50-4f83-832c-2bbbc5bba6f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = \"retina\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bf8e48cd-cf59-4f2f-9c93-8f171efebdd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt, ticker\n",
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cbad033c-b589-4d32-98f4-beaa7d501856",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(color_codes=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "963f4f1b-9b8f-450b-83e3-e1f7e714ae17",
   "metadata": {},
   "source": [
    "We begin by generating 5,000 uniformly random points in the unit square $0 \\leq x, y \\leq 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "16dafed0-be07-4bd6-a663-44b2775263eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for reproducibility\n",
    "SEED = 123456789\n",
    "\n",
    "rng = np.random.default_rng(SEED)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "100ea4ce-5a5c-4e28-9dc9-d2de709908a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 5_000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "564bd8cf-b91e-4916-bbb1-9213f2ae66db",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = rng.uniform(size=(2, N))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dc08bac-ca50-4633-b8f4-16e0f4f57a24",
   "metadata": {},
   "source": [
    "Here we visualize these samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "61bfc452-0ff6-4efe-8ab9-120f60133e2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_axes():\n",
    "    _, ax = plt.subplots()\n",
    "    ax.set_aspect(\"equal\")\n",
    "\n",
    "    ax.xaxis.set_major_locator(ticker.MultipleLocator(1))\n",
    "    ax.set_xlim(0, 1.01)\n",
    "    \n",
    "    ax.yaxis.set_major_locator(ticker.MultipleLocator(1))\n",
    "    ax.set_ylim(0, 1.01)\n",
    "\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "899f1956-5e66-4bf9-a64c-c6bfa0212074",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {