Test Notebook for NE506 Exercises

random-sampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In-Class Exercise #3: Random Sampling\n",
    "\n",
    "In this exercise, we will use a Jupyter notebook to explore discrete and direct continuous sampling. You can use these templates as a starting point. Some of the following steps have already been incorporated into the templates as denoted by a ✅. Study these portions of the files, and make your own additions where necessary. \n",
    "\n",
    "Download this file and submit it with your last name prepended to the filename.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### I. Set the Random Number Seed\n",
    "\n",
    "In order to ensure reproducibility, initialize the seed of NumPy's random number generator. It is good practice to have a random seed that stays constant for testing your results.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### II. Direct Discrete Sampling\n",
    "\n",
    "#### A. Load and Visualize PDF\n",
    "\n",
    "In the template, create arrays with the data in the table below. This is a small data set representing the pets owned by 30 different people. The data contain five possibilities for pets (denoted by a \"bin number\" in column 1 of the table and labeled in column two in the table below) and their corresponding (unnormalized) probabilities. \n",
    "\n",
    "  | Bin | Pet | Unnormalized PDF | Normalized PDF Value | CDF |\n",
    "  |---|---|---|---|---|\n",
    "  | 0   | Fish | 5 | | |\n",
    "  | 1 | Cat | 4 | | |\n",
    "  | 2 | Dog | 12 | | |\n",
    "  | 3 | Rabbit | 2 | | \n",
    "  | 4 | Bird | 7 | | \n",
    "\n",
    "\n",
    "\n",
    "Task: Starting on paper, normalize the PDF data and determine the CDF (i.e., fill in the empty columns in the table). Next, implement these calculations in the Jupyter notebook template. Plot the PDF and CDF of this distribution ✅."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "vscode": {
     "languageId": "plaintext"
    }
   },
   "source": [
    "### B. Create Discrete Sampling Function\n",
    "\n",
    "Write a function, samplePDF, to perform direct discrete sampling.  Your function will take three arguments:\n",
    "\n",
    "  - a uniform random variable,\n",
    "  - a vector representing the set of possible pets for the distribution being sampled,\n",
    "  - a vector of probabilities for each of those pets,\n",
    "\n",
    "and will return a single value:\n",
    "\n",
    "a sampled value (i.e., \"bin number\" corresponding to a particular pet) from the set of possibilities.\n",
    "\n",
    "The beginning of `samplePDF` in the cell below provides some checks on the arguments. The second, commented, half of this function is left for you to complete.\n",
    "\n",
    "\n",
    "1) determine the CDF (this should be identical to your solution in the previous part\n",
    "2) return a \"bin\" value (Hint: see the while loop). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def samplePDF(prn, vals, pdf):\n",
    "    ## some defensive programming\n",
    "    if prn < 0.0 or prn > 1.0:\n",
    "        raise ValueError(f'Random variable {prn} is out of [0, 1) bounds')\n",
    "\n",
    "    if len(vals) != len(pdf):\n",
    "        raise ValueError(f'Length of values ({len(vals)}) and probabilities ({len(pdf)}) are not the same')\n",
    "\n",
    "    if sum(pdf) != 1.0:\n",
    "        print(f'Original PDF was not normalized and has been normalized')\n",
    "        pdf /= sum(pdf)\n",
    "\n",
    "    # calculate the CDF\n",
    "    cdf = None # <-- fill this in\n",
    "\n",
    "    # search for the sample in the PDF\n",
    "    bin = 0\n",
    "    while False: # replace False with the correct expresion\n",
    "        bin += 1\n",
    "\n",
    "    return vals[bin]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### C. Test Discrete Sampling Function\n",
    "\n",
    "Use your function, `samplePDF`` to sample the PDF and\n",
    "plot a histogram of the results. ✅  \n",
    "\n",
    "Do this for 10 samples, 10,000 samples, and 1,000,000 samples. \n",
    "\n",
    "Some convenience functions have been provided below to help with plotting.\n",
    "\n",
    "Note that code is included in these functions to plot the true PDF from the table above and the sampled data on the same plot.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "pdf = np.array([5, 4, 12, 2, 7], dtype=float)\n",
    "vals = np.array([0, 1, 2, 3, 4], dtype=float)\n",
    "pets = ['fish', 'cat', 'dog', 'rabbit', 'bird']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "pdf /= pdf.sum()\n",
    "cdf = np.cumsum(pdf)\n",
    "# plot the CDF using the \"stairs\" plot type\n",
    "plt.figure()\n",
    "plt.title('Cumulative Distribution Function (CDF)')\n",
    "plt.xlabel('Pets')\n",
    "plt.stairs(vals[1:], cdf)\n",
    "plt.show()\n",
    "\n",
    "num_samples = 100\n",
    "results = np.array([samplePDF(prn, vals, pdf) for prn in np.random.rand(num_samples)])\n",
    "# plot the PDF from the imported data and the sampled data\n",
    "# using the \"stem\" plot type\n",
    "plt.figure()\n",
    "plt.title('Probability Distribution Function (PDF)')\n",
    "l1 = plt.stem(..., label='True PDF') # <-- complete this line\n",
    "hist = None # <-- fill this in\n",
    "\n",
    "plt.xlim([-1, 5])\n",
    "l2 = plt.stem(... label='Sampled PDF', markerfmt='ro') # <-- complete this line\n",
    "plt.xticks(range(len(vals)), labels=pets)\n",
    "plt.xlabel('Pets')\n",
    "plt.ylabel('p(pets) PDF')\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III Direct Continuous Sampling\n",
    "\n",
    "We will now sample the scattering conditions of a particle undergoing a collision with the following (fictional) scattering laws:\n",
    "\n",
    "- $f(\\mu) = \\mu^{2}$  for 0 ≤ $\\mu$ < 1, where μ is the cosine of the scattering angle\n",
    "- $f(\\epsilon) = \\large{\\frac{1-\\epsilon}{2-\\mu}}$ for  0 < $\\epsilon$ < 1, where ε is the ratio of the exit energy to the incoming energy\n",
    "\n",
    "### A. Pre-coding Preparation\n",
    "\n",
    "Answer the questions in the following cells or on paper. If done on paper, hand them in on the exercise due date.\n",
    "\n",
    "#### Scattering angle:\n",
    "\n",
    "- What is the normalized PDF for the scattering angle cosine, μ?\n",
    "\n",
    "- What is the CDF for the scattering angle cosine, μ?\n",
    "\n",
    "- What is the inverse of this CDF for μ?\n",
    "\n",
    "#### Exit energy:\n",
    "\n",
    "- What is the normalized PDF for the outgoing-to-incoming energy ratio, ε? Does it depend on μ?\n",
    "\n",
    "- What is the CDF for the outgoing-to-incoming energy ratio, ε?\n",
    "\n",
    "- What is the inverse of this CDF for ε?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### B. Sample the angle\n",
    "\n",
    "- Write a function, sampleAngle, to sample the scattering angle given a uniform random variable.\n",
    "- Plot a histogram with 50 bins from 50,000 samples found using your function (sampleAngle). ✅\n",
    "- Determine the mean, standard deviation and standard error of this result.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def sampleAngle(prn):\n",
    "    if (prn < 0 or prn >= 1):\n",
    "        raise ValueError(f'Random number {prn} is outside the valid interval [0, 1)')\n",
    "\n",
    "    return None # <-- fill this in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "def sampleEnergy(prn):\n",
    "    return None # <-- fill this in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### D. Sample a collision\n",
    "\n",
    "Consider a beam of particles starting with 5 MeV and having a collision at the origin. In the cells below, do the following:\n",
    "\n",
    "- Create a list of the energy and direction of the first 5 particles after the collision.\n",
    "- Determine and plot the location of these particles as they leave across a unit circle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E0 = 5\n",
    "n_particles = 5\n",
    "\n",
    "for i in range(n_particles):\n",
    "    mu = None # <-- fill this in\n",
    "    E = None # <-- fill this in\n",
    "\n",
    "    x = None # <-- fill this in\n",
    "    y = None # <-- fill this in\n",
    "\n",
    "    print(f'Particle energy: {E:.2f} MeV, Particle crossing: x: {x:.4f}, y: {y:.4f}')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}