Lab_Lec2_unsolved.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The %... is an iPython thing, and is not part of the Python language.\n",
    "# In this case we're just telling the plotting library to draw things on\n",
    "# the notebook, instead of on a separate window.\n",
    "%matplotlib inline\n",
    "# See all the \"as ...\" contructs? They're just aliasing the package names.\n",
    "# That way we can call methods like plt.plot() instead of matplotlib.pyplot.plot().\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib as mpl\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import time\n",
    "pd.set_option('display.width', 500)\n",
    "pd.set_option('display.max_columns', 100)\n",
    "pd.set_option('display.notebook_repr_html', True)\n",
    "import seaborn as sns\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.set_context(\"poster\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The properties of PDFs\n",
    "Probability density functions (PDFs) are characterized by certain parameters, that give the probability as a function of some quantity. The normal distribution, for example, is characterized by a mean value $\\mu$ and and a standard deviation $\\sigma$. We can play with PDFs using Pyhton's scipy and its stats package. The \"standarized\" normal is always centerer at 0 and has standard deviation of 1:\n",
    "\n",
    "$$f(x) = \\frac{\\exp(-x^2/2)}{\\sqrt{2\\pi}}$$\n",
    "\n",
    "But we can always change the $\\mu$ and $\\sigma$ using the parameters *loc* and *scale*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm,gamma\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "x = np.linspace(norm.ppf(0.001),norm.ppf(0.999), 1000)\n",
    "ax.plot(x, norm.pdf(x), 'b-', lw=3, alpha=0.6, label='norm pdf')\n",
    "x = np.linspace(norm.ppf(0.001,loc=2.0,scale=1.5),norm.ppf(0.999,loc=2.0,scale=1.5), 1000)\n",
    "ax.plot(x, norm.pdf(x,loc=2.0,scale=1.5), 'r-', lw=3, alpha=0.6, label='norm pdf')\n",
    "ax.set_xlabel(r'$x$')\n",
    "ax.set_ylabel(r'PDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can obtain samples from the distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = norm.rvs(size=1000)\n",
    "x = np.linspace(norm.ppf(0.001),norm.ppf(0.999), 1000)\n",
    "plt.plot(x, norm.pdf(x), 'b-', lw=3, alpha=0.6, label='norm pdf')\n",
    "plt.hist(r, density=True, histtype='stepfilled', alpha=0.2)\n",
    "plt.legend(loc='best', frameon=False)\n",
    "plt.xlabel(r'$x$')\n",
    "plt.ylabel(r'N')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *cumulative distribution function* measures the cumulative probability of the distribution, namely, the accumulated probability up to certain value. It is useful, for example, to answer the question: what is the probability of obtaining a value that is less than certain $x_0$?:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "x = np.linspace(norm.ppf(0.001),norm.ppf(0.999), 1000)\n",
    "ax.plot(x, norm.cdf(x), 'b-', lw=3, alpha=0.6, label='norm pdf')\n",
    "ax.set_xlabel(r'$x$')\n",
    "ax.set_ylabel(r'CDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Poisson distribution\n",
    "\n",
    "The Poisson distribution is another discrete distribution, it expresses the probability of a given number of events occurring in a fixed interval of time (or space, volume, etc.). One assumption made is that these events occur with a known average rate and independently of each other. An example is the number of electrons detected by a sensor in an electron microscope during a time interval, or the number of soldiers in the Prussian army killed accidentally by horse kicks [(see here)](http://en.wikipedia.org/wiki/Poisson_distribution).\n",
    "\n",
    "The Poisson distribution is defined as:\n",
    "\n",
    "\n",
    "$$ f(k; \\mu)= \\frac{\\mu^k e^{-\\mu}}{k!}, $$\n",
    "\n",
    "where $k$ is the number of events, $\\mu$ is a positive real number, and $e$ is Euler's number ($e = 2.71828 \\ldots$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import poisson\n",
    "# generate samples for different values of mu\n",
    "kpts=np.arange(0,25)\n",
    "for mu, c in zip([1,2, 4, 6], sns.color_palette()[:4]):\n",
    "    randomVariates = poisson.rvs(mu, size=1000)\n",
    "    plt.hist(randomVariates, normed=True, color=c, alpha=0.2, bins=range(0,26), label='$M=' + np.str(mu) + '$')\n",
    "    plt.plot(kpts, poisson.pmf(kpts, mu), '.', color=c)\n",
    "\n",
    "plt.legend()\n",
    "plt.title(\"Poisson Distribution\")\n",
    "plt.xlabel(\"Number of Events\")\n",
    "plt.ylabel(\"Normed Counts\")\n",
    "plt.savefig('Posson_dist.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### X-ray counts in a Chandra Observation\n",
    "We now read an X-ray dataset from Chandra, and try to look at the distribution of events as a function of time. The dataset is the same as for the first exercise. We first need to read the event file, and then we look at the photons hotto"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "hdul = fits.open('/Users/jmartine/teaching/isya2018/lecture1/data/18344/primary/acisf18344N001_evt2.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us investigate how many photons hit the detector at a particular pixel position, per unit time. Is it always the same number? Let us plot the distribution of these numbers. Let us also try to guess a Poisson distribution that matches the observed values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a mask to select the desired pixel\n",
    "mask = ((hdul[1].data['ccd_id']==7) \n",
    "        & (hdul[1].data['dety']>=3961.)\n",
    "        & (hdul[1].data['dety']<3962.)\n",
    "        & (hdul[1].data['detx']>=4519.)\n",
    "        & (hdul[1].data['detx']<4920.))\n",
    "\n",
    "# Get the times of arrival of the photos for that pixel\n",
    "times = np.array(hdul[1].data['time'][mask])\n",
    "\n",
    "# Split the time into equal intervals\n",
    "t_intvls = np.linspace(np.min(times),np.max(times),300)\n",
    "\n",
    "# Count how many photos were detected in each equal time interval\n",
    "n_events = []\n",
    "for i in range(len(t_intvls[0:-1])):\n",
    "    n_events.append(len(times[(times>t_intvls[i]) & (times<t_intvls[i+1])]))\n",
    "\n",
    "#Plot the histogram\n",
    "plt.hist(n_events,alpha=0.4)\n",
    "plt.title('No. of photons detected on pixel (3961,4919) per unit time')\n",
    "plt.xlabel('No. of photons')\n",
    "plt.ylabel('N')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now try to see if we can match this with a Poisson distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kpts=np.arange(0,40)\n",
    "mu = 15\n",
    "plt.hist(n_events,alpha=0.4,normed=True,bins=15)\n",
    "plt.plot(kpts, poisson.pmf(kpts, mu), '.', color='r')\n",
    "plt.title('No. of photons detected on pixel (3961,4919) per unit time')\n",
    "plt.xlabel('No. of photons')\n",
    "plt.ylabel('N')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now look at the energy of each photon. What is the distribution of photon energies? Let's have a look. Can we adjust a normal distribution to the energies?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = (hdul[1].data['ccd_id']==7)\n",
    "plt.hist(np.log10(hdul[1].data['energy'][mask]),bins=15,alpha=0.4,normed=True)\n",
    "x = np.linspace(norm.ppf(0.0001,loc=3.25,scale=0.15),norm.ppf(0.9999,loc=3.25,scale=0.15), 1000)\n",
    "plt.plot(x, norm.pdf(x,loc=3.25,scale=0.15), 'b-', lw=3, alpha=0.6, label='norm pdf')\n",
    "plt.title('Distribution of photon energies')\n",
    "plt.xlabel('log Energy (eV)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next question is how many photons in a given energy range hit the detector per unit time, regardless of their position. Let us have a look."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = ((hdul[1].data['ccd_id']==7) \n",
    "        & (hdul[1].data['energy']>=10**3.264)\n",
    "        & (hdul[1].data['energy']<10**3.266))\n",
    "\n",
    "times = np.array(hdul[1].data['time'][mask])\n",
    "\n",
    "t_intvls = np.linspace(np.min(times),np.max(times),300)\n",
    "\n",
    "n_events = []\n",
    "for i in range(len(t_intvls[0:-1])):\n",
    "    n_events.append(len(times[(times>t_intvls[i]) & (times<t_intvls[i+1])]))\n",
    "    \n",
    "plt.hist(n_events,alpha=0.4)\n",
    "plt.title('No. of 3.26 keV photons detected per unit time')\n",
    "#plt.hist(np.log10(hdul[1].data['energy'][mask]),bins=15,alpha=0.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Joint, marginal, and conditional probabilities\n",
    "We now illustrate the concept of the different types of probabilities with an example. Suppose you roll two fair dice. What is the probability of getting 3 in one die and 2 in the other? Or 1 and 4, or 6 and 6? Also, what is the probability of getting 3 in one of them regardless of what you get in the other? And finally, what is the probability of getting 3 in one of them, given that you got exactly 3 in the other one? Let us perform the experiment by rolling our digital dice many times. First let us roll the dice a few times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "repeat = True\n",
    "while repeat:\n",
    "    print(\"You rolled\",np.random.randint(1,7))\n",
    "    print(\"Do you want to roll again?\")\n",
    "    repeat = (\"y\" or \"yes\") in input().lower()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now let us roll two dice 100000 times, and illustrate joint, marginal, and conditional probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_rolls = 100000\n",
    "d1 = np.random.randint(1,7,n_rolls)\n",
    "d2 = np.random.randint(1,7,n_rolls)\n",
    "\n",
    "# Here we count how many rolls of each pair we have\n",
    "freqs = np.zeros((6,6))\n",
    "for i in range(6):\n",
    "    for j in range(6):\n",
    "        for roll in list(zip(d1,d2)):\n",
    "            if (roll == (i+1,j+1)): \n",
    "                freqs[i,j]+= 1\n",
    "                \n",
    "freqs = freqs.astype(int)\n",
    "print(freqs)\n",
    "\n",
    "\n",
    "print(freqs[2,2]/100000.)\n",
    "print(sum(freqs[:,2])/100000.)\n",
    "print(freqs[:,1])\n",
    "\n",
    "#print(list(zip(d1,d2)))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us define $n_{i,j}$ as the number of times the die 1 ends in number $i$ and die 2 ends in number $j$. Then the **joint** probability that both dice end up in three is:\n",
    "\n",
    "$$P(d_1=3,d_2=3) = \\frac{n_{3,3}}{N}$$\n",
    "\n",
    "where $N$ is the total number of joint rolls:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P_33 = freqs[2,2]/n_rolls\n",
    "print(P_33)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **marginal** probability of die 1 ending up in three is:\n",
    "\n",
    "$$P(d_1=3) = \\frac{\\sum_j n_{3,j}}{N}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P_3j = sum(freqs[:,2])/n_rolls\n",
    "print(P_3j)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the **conditional** distribution that first die 1 ends in 3 given the second ones ends in 2 is:\n",
    "\n",
    "$$P(d_1=3|d_2=2) = \\frac{n_{3,2}}{\\sum_i n_{i,2}}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P_3_given_2 = freqs[2,1]/sum(freqs[:,1])\n",
    "print(P_3_given_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise:\n",
    "Consider the photons arriving to CCD 7 in the Cassiopea A dataset.\n",
    "* Using CIAO, download observation 12020 from the Chandra archive. This corresponds to the Cassiopeia A supernova remnant.\n",
    "* Produce an RGB image like last time, and identify regions dominated by photons of different energies. \n",
    "* Using DS9, load and save a few regions of interest.\n",
    "* What is the distribution of the logarithmic photon energies in the for Cas A? What kind of distribution do you think that is?\n",
    "* Find the mean and the standard deviations of the distribution, and record them\n",
    "* Now do the same, but now limit yourself to each of the regions of interest.\n",
    "* In each case, compare the histograms of energy with normal distributions with μ corresponding to the mean and σ corresponding to the standard deviation. What is expected value of the energy in each case?\n",
    "* What is the probability of detecting a photon in each of the three ands s,m, and h energy bands for each of the regions? Discuss this in terms of conditional probability distributions.\n",
    "* What is the probability of detecting hard or soft X-ray photons coming from a white dwarf?\n",
    "* What does this tell you about the physical processes taking place in the different regions?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this exercise we can either read the evt2 file using astropy, or just do the different searches in CIAO, save the tables, and then read the tables here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the evt2 file\n",
    "hdul = fits.open('/Users/jmartine/teaching/isya2018/lecture2/data/12020/repro/acisf12020_repro_evt2.fits')\n",
    "\n",
    "#maks for only CCD 7\n",
    "mask = ((hdul[1].data['ccd_id']==7))\n",
    "energies = hdul[1].data['energy'][mask]\n",
    "\n",
    "# Plot the histogram of energies\n",
    "plt.hist(np.log10(hdul[1].data['energy'][mask]),bins=15,alpha=0.4,normed=True)\n",
    "plt.title('Distribution of photon energies')\n",
    "plt.xlabel('log Energy (eV)')\n",
    "\n",
    "# Find the mean and standard deviation\n",
    "# YOUR CODE HERE\n",
    "\n",
    "# plot a normal distribution mathing the data\n",
    "loca = \n",
    "stda = \n",
    "x = np.linspace(norm.ppf(0.0001,loc=loca,scale=stda),norm.ppf(0.9999,loc=loca,scale=stda), 1000)\n",
    "plt.plot(x, norm.pdf(x,loc=loca,scale=stda), 'b-', lw=3, alpha=0.6, label='norm pdf')\n",
    "\n",
    "\n",
    "# Creat masks with the different ehergies\n",
    "h_band_mask = \n",
    "m_band_mask = \n",
    "s_band_mask = \n",
    "\n",
    "\n",
    "#The propbability of detecting an h band photon:\n",
    "\n",
    "\n",
    "#The propbability of detecting an h band photon:\n",
    "\n",
    "\n",
    "#The propbability of detecting an s band photon:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now using the Cumulative Density Function:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we the regions with tables generated with CIAO with dmlist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Region 1:\n",
    "from astropy.io import ascii\n",
    "data = ascii.read(\"reg1_energies.dat\")  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energies = data['energy'].data\n",
    "plt.hist(np.log10(energies),bins=10,alpha=0.4,normed=True)\n",
    "plt.title('Distribution of photon energies for Region 1')\n",
    "plt.xlabel('log Energy (eV)')\n",
    "\n",
    "# Find the mean and standard deviation\n",
    "\n",
    "\n",
    "# plot a normal distribution mathing the data\n",
    "\n",
    "\n",
    "# Create energy masks\n",
    "\n",
    "\n",
    "#The propbability of detecting an h band photon:\n",
    "\n",
    "\n",
    "#The propbability of detecting an h band photon:\n",
    "\n",
    "#The propbability of detecting an s band photon:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Region 2\n",
    "data = ascii.read(\"reg2_energies.dat\")  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energies = data['energy'].data\n",
    "plt.hist(np.log10(energies),bins=10,alpha=0.4,normed=True)\n",
    "plt.title('Distribution of photon energies for Region 1')\n",
    "plt.xlabel('log Energy (eV)')\n",
    "\n",
    "# Find the mean and standard deviation\n",
    "\n",
    "# plot a normal distribution mathing the data\n",
    "\n",
    "# Create energy bands\n",
    "\n",
    "\n",
    "#The propbability of detecting an h band photon:\n",
    "\n",
    "#The propbability of detecting an h band photon:\n",
    "\n",
    "\n",
    "#The propbability of detecting an s band photon:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Derivation of the Bayes' Theorem\n",
    "We start with the product rule of probabilities. We know that:\n",
    "$$P(A|B) = \\frac{P(A,B)}{P(B)}$$\n",
    "And also:\n",
    "$$P(B|A) = \\frac{P(A,B)}{P(A)}$$\n",
    "Rearranging and combining these two equations, we get:\n",
    "$$P(A|B)P(B) = P(A|B)P(A)$$\n",
    "Now we divide both sides by $P(B)$:\n",
    "$$P(A|B) = \\frac{P(A|B)P(A)}{P(B)}$$\n",
    "This is know as the Bayes' rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}