Week7/PiEstimate/PiEstimate_original.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "# Estimating the value of $\\pi$ using a Monte Carlo technique \n",
    "\n",
    "The program picks random points in the area $[-1,1] \\times [-1,1]$, and determines which fraction of these are within the unit circle. This in turn gives a measure of $\\pi$ with associated statistical uncertainty. Performing such \"experiments\" many times not only gives the value of $\\pi$, but also a feel for the use of Monte Carlo, and experience in calculating averages, RMSs, and the error on these.\n",
    "The first part of the program is mostly for illustration, while the exercise is question 4 on using this in higher dimension problems.\n",
    "\n",
    "For more information see:\n",
    "- P. R. Bevington: page 75-78\n",
    "\n",
    "***\n",
    "\n",
    "### Authors: \n",
    "- Troels C. Petersen (Niels Bohr Institute)\n",
    "\n",
    "### Date:    \n",
    "- 06-01-2023 (latest update)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                                     # Matlab like syntax for linear algebra and functions\n",
    "import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab\n",
    "from iminuit import Minuit                             # The actual fitting tool, better than scipy's\n",
    "import sys                                             # Modules to see files and folders in directories\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and define the parameters of the program:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_plots = True         # Determining if plots are saved or not\n",
    "r = np.random\n",
    "r.seed(42)\n",
    "\n",
    "# Set parameters:\n",
    "N_experiments = 100        # Number of \"experiments\" determining pi\n",
    "N_points      = 2000       # Number of points per experiment in determining pi\n",
    "\n",
    "pi_true = np.pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loop over process:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi_dist    = np.zeros(N_experiments)\n",
    "hit_dist_x = np.zeros(N_points)\n",
    "hit_dist_y = np.zeros(N_points)\n",
    "\n",
    "for iexp in range(N_experiments):\n",
    "\n",
    "    # Fill random numbers into square, and see how many are inside unit circle:\n",
    "    hit_dist_x = r.uniform(-1, 1, N_points)\n",
    "    hit_dist_y = r.uniform(-1, 1, N_points)\n",
    "    Nhit = np.sum(hit_dist_x**2 + hit_dist_y**2 < 1)\n",
    "    \n",
    "    # Calculate the fraction of points within the circle and its error:\n",
    "    f  = Nhit / N_points\n",
    "    sf = -1         # Remind yourself the uncertainty on a fraction n/N!!!\n",
    "\n",
    "    # From this we can get pi and its error, and save it:\n",
    "    pi_estm  = 4.0 * f\n",
    "    pi_error = 4.0 * sf\n",
    "    pi_dist[iexp] = pi_estm\n",
    "\n",
    "    # Print first couple of pi measurements:\n",
    "    if (iexp < 5) :\n",
    "        print(f\"{iexp:2d}. pi estimate:   {pi_estm:7.4f} +- {pi_error:6.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "if (N_experiments > 1):\n",
    "    pi_mean = pi_dist.mean()\n",
    "    pi_std = pi_dist.std(ddof=1)\n",
    "    pi_standard_error = pi_std / np.sqrt(N_experiments)\n",
    "    print(f\"The {N_experiments:d} experiments yield: pi = {pi_mean:7.5f} +- {pi_standard_error:7.5f}     (std = {pi_std:5.3f})\")\n",
    "    print(f\"and the agreement with the true value of pi is: {(pi_mean - pi_true) / pi_standard_error:6.2f} sigma\")\n",
    "else :\n",
    "    print(\"Since only one experiment was done, no average and RMS can be calculated.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the histograms:\n",
    "\n",
    "Distribution of points from one experiment:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig1, ax1 = plt.subplots(figsize=(8, 8))    \n",
    "\n",
    "# Make mask of points inside vs outside the circle:\n",
    "mask_inside = hit_dist_x**2 + hit_dist_y**2 < 1\n",
    "mask_outside = ~mask_inside\n",
    "\n",
    "# Make scatter plot of points:\n",
    "ax1.plot(hit_dist_x[mask_inside], hit_dist_y[mask_inside], '.')\n",
    "ax1.plot(hit_dist_x[mask_outside], hit_dist_y[mask_outside], '.')\n",
    "ax1.set(xlim=(-1, 1), ylim=(-1,1))\n",
    "\n",
    "# Add circle to plot:\n",
    "ax1.add_artist(plt.Circle((0, 0), 1, color='k', fill=False))\n",
    "\n",
    "if save_plots: \n",
    "    fig1.savefig(\"HitDist.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "First acquaint yourself with the program, and make sure that you understand what the parameters `N_experiment` and `N_points` refer to! Also, before running the program,\n",
    "calculate what precision you expect on pi in each experiment, when using the number of points chosen in the program (i.e. 2000 points).\n",
    "\n",
    "Then, run the program, and then take a look at the result... which requires that you fill in the calculations yourself!\n",
    "\n",
    "# Initial simple questions:\n",
    "\n",
    "0. Plot the distribution of $\\pi$ values obtained.\n",
    "\n",
    "1. Try to run 100 experiments with 2000 points in each. What is the expected approximate uncertainty on $\\pi$ in each experiment? Does the result agree with what you calculated before running the program? What is the uncertainty on the AVERAGE of all 100 experiments? \n",
    "\n",
    "2. How do you expect the values of $\\pi$ to distribute themselves? And is this the case here?\n",
    "\n",
    "3. Does it make any difference on the precision of the final $\\pi$ value, whether you make many experiments with few points, or one experiment with many points, as long as the product of `N_experiment` $\\times$ `N_points` remains constant?\n",
    "\n",
    "# The real exercise question:\n",
    "\n",
    "4. Now try to use this method in three dimensions to estimate the constant in front of the $r^3$ expression for the volume. Do you get $\\frac{4}{3} \\pi$?\n",
    "    Increase the dimensionality (say up to 10), and see if you can figure out the constants needed to calculate the hyper-volumes!\n",
    "\n",
    "    HINT: I'll reveal that for `N_dim` of 4 and 5, the constant contains $\\pi^2$ and some simple rational fraction, while for `N_dim` 6 and 7, it contains $\\pi^3$ and a\n",
    "      rational fraction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning points:\n",
    "\n",
    "This is an exercise in using simulation for numerical estimates, that are otherwise rather complicated. Especially in high dimensionality, it can be hard to produce good solutions, while simulation (such as this exercise) gives an easier alternative.\n",
    "\n",
    "From the exercise you should:\n",
    "1. Get experience with a simulation solution example.\n",
    "2. Train estimation of uncertainties on multiple estimates (here of pi).\n",
    "3. Exercise hypothesis testing, comparing simulated estimates to multiple hypothesis (rational fractions). \n",
    "\n",
    "Following this exercise, you should feel comfortable in working out a simulated answer to a problem/case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "executable": "/usr/bin/env python",
  "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.7.6"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}