chap3_ex.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 3: Statistical models and inference (exercises)\n",
    "\n",
    "All the exercises are on distance inference from parallaxes, and build on one another, so **you need to do them in order**.\n",
    "\n",
    "The parallax of an object is its observed angular displacement with respect to a reference frame due to the movement of the observer over a baseline. Stellar parallaxes are measured using twice the Earth-Sun separation as a baseline. From geometry and the definition of the parallax, $\\omega$, the distance of a star from the Sun, $r$ pc, is $1/\\omega$ arcsec to a very good approximation.\n",
    "\n",
    "But inverting a parallax to give a distance is only appropriate when we have no measurement errors. As we always have measurement errors, determining the distance given a parallax becomes an inference problem. Let us investigate how to do this"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Excercise 1: Truncated uniform prior (30 points)\n",
    "\n",
    "Write computer programs to\n",
    "\n",
    "* compute the posterior PDF (which is a function of distance) given the parallax, $\\omega$, and its standard deviation, $\\sigma_{\\omega}$, using the Gaussian likelihood (equation 3.19) and an arbitrary prior (i.e., the program should be flexible to make use of different priors);\n",
    "$$ P(\\omega \\mid r, \\sigma_\\omega) = \\frac{1}{\\sqrt{2 \\pi}\\sigma_\\omega} \n",
    "\\exp{ \\left[ -\\frac{1}{2\\sigma_\\omega^2} \\left(\\omega-\\frac{1}{r} \\right)^2 \\right]} $$\n",
    "\n",
    "* perform the empirical test described in section 3.5.4 to determine the bias and standard deviation of an arbitrary distance estimator (in particular the mode) as a function of $f_{true}$.  \n",
    "\n",
    "Use this second program to investigate the performance of the two distance estimators given in\n",
    "section 3.5.3 (i.e. the mode of the improper uniform prior) and section 3.5.5 (i.e., the mode of the\n",
    "truncated uniform prior, equation 3.30). In both cases: draw simulated data from a distribution\n",
    "which is uniform in the range $(0, 10^3)$ pc and zero outside this; omit from your samples simulated\n",
    "parallax measurements which are negative. Specifically,\n",
    "\n",
    "1. plot the posterior PDF for the truncated uniform prior (i.e. equation 3.29) for $f = (0.1, 0.2, 0.5, 1.0)$\n",
    "and $r_{lim}=10^3$ pc; (i.e. reproduce at least the left panel of Figure 3.7). **posterior normalization can be numerically computed**.\n",
    "\n",
    "2. do the tests which led to Figures 3.6 and 3.8 and plot your results (i.e. reproduce these figures;\n",
    "use the same axes and ranges). Note: _the distance estimator is $1/\\omega$._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 2: Truncated uniform space density prior (40 points)\n",
    "\n",
    "Use your programs from exercise 1 to investigate the use of the truncated uniform space density prior (equation 3.26), to define the posterior PDF, and the bias and standard deviation of the mode of this as a distance estimator. Specifically, \n",
    "\n",
    "1. determine the posterior PDF, and then use your program to plot it with $\\omega = 1/100$ arcsec, $r_{lim} = 10^3$ pc for $f = (0.1, 0.2, 0.5, 1.0)$ (overplot these on the same plot so you can best see the differences). **Recall that a posterior should be normalized**: you will have to do this numerically here (it is sufficient to approximate the area under the curve using a sequence of narrow rectangles; section 3.4.3). Plot the **normalized** posteriors on one plot, **and** the unnormalized posteriors scaled to have their modes equal to unity on another plot. Comment on how the shape of the posterior in general, and the location of the mode in particular, changes with a function of $f$ . (You may want to plot at other values of $f$ to see what is going on.) There is of course another mode at $r = r_{lim}$ imposed by the prior, but this is independent of $f$. (2 plots to produce)\n",
    "\n",
    "2.  calculate **analytically** the value of the mode, $r_{mode}$, of this posterior (show your working) as a function of $f$ (this is not referring to the fixed mode at $r=r_{lim}$!). Plot $\\omega \\cdot r_{mode}$ vs. $f$. What do you notice about the solution? For what values of $f$, if any, is the mode not defined? What value of $r_{mode}$ would be sensible to use in this case (i.e., to be consistent with our prior assumptions)? Likewise, what could we/should we do if the parallax is negative?\n",
    "\n",
    "3. use your program from exercise 1 to calculate and plot the bias and standard deviation of $r_{mode}$ as a function of $f_{true}$ over the range $f_{true}=(0,1)$, i.e.\\ produce a plot similar to Figure 3.5.  To do this you need to draw random numbers from the uniform space density prior (equation 3.26). We will look at methods for drawing random numbers from arbitrary probability distributions in a later chapter. For now you can just use the following implemention of rejection sampling, given below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below you'll find a R and python function that will return _approximately_ `Nsamp/3` samples for the distribution truncated at `rlim`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting rsnip.r\n"
     ]
    }
   ],
   "source": [
    "%%file rsnip.r\n",
    "r.prior1 <- function(Nsamp, rlim) { # Vectorized in Nsamp\n",
    "    # Use rejection sampling. Expect 1/3 * Nsamp samples to be returned\n",
    "    r1 <- runif(Nsamp, min=0, max=1)\n",
    "    r2 <- runif(Nsamp, min=0, max=1)\n",
    "    r <- rlim*r1[which(r2<=r1^2)]\n",
    "    return(r)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def prior1(Nsamp, rlim):\n",
    "    # Use rejection sampling. Expect 1/3 * Nsamp samples to be returned\n",
    "    r1 = np.random.uniform(0, 1, Nsamp)\n",
    "    r2 = np.random.uniform(0, 1, Nsamp)\n",
    "    r = rlim * r1[r2 <= r1 ** 2]\n",
    "    return r    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To keep things simple when calculating the bias and standard deviation, reject samples which\n",
    "have negative parallaxes or which give undefined distances. For those cases where the mode\n",
    "given by your estimator (solution to part (b)) is undefined or exceeds $r_{lim}$ , what should you\n",
    "use as a distance estimator? _Hint: Think about what the posterior looks like in those cases._\n",
    "\n",
    "4. compare the performance of this distance estimator with the mode of the posterior from the truncated uniform prior (exercise 2). Which gives lower bias and standard deviation for low and high values of $f_{true}$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 3: Hipparcos data (25 points)\n",
    "\n",
    "Look at the data set `hipparcos.dat`. This contains the parallax (in mas), parallax standard deviations (`errors`) $\\sigma_\\omega$ (in mas), and V-band apparent magnitudes of all stars in the Hipparcos catalogue. The noise model for the parallaxes is a Gaussian. We want to derive distances to these stars using the mode estimators of the two posteriors (corresponding to the two priors) we studied in the previous two questions. Let's call\n",
    "them $r_{\\rm uni}$ for the mode of the truncated uniform prior in $r$, and $r_{\\rm usd}$ for the \n",
    "mode of the truncated uniform space density prior.\n",
    "\n",
    "1. Read in and inspect the data: look at the ranges, plot histograms of **each** the measurements and the fractional parallax error (as in exercice 2), plot the measurements against each other. Plot $\\sigma_\\omega$ vs. $\\omega$ and $\\sigma_\\omega$ $ vs. magnitudes. Think about what these tell you. Remove stars which have any missing data (indicated with `NA`). Find out how many of each of the three variables are negative, zero, infinite. Are such values valid for these variables? Be aware of the units!\n",
    "\n",
    "2. Calculate the distance estimates $r_{\\rm uni}$ and $r_{\\rm usd}$ with $r_{lim}=1000$ pc. Rather than rejecting stars with negative parallaxes or which give undefined solutions of the estimator, assign them a value which is consistent with the corresponding prior (justify your choice).  Plot histograms of $r_{\\rm uni}$, $r_{\\rm usd}$, and $r_{\\rm uni} - r_{\\rm usd}$, choosing sensible bins in each case. Explain the differences that you see.  Plot also $r_{\\rm uni}$ vs. $r_{\\rm usd}$ and overplot the line $r_{\\rm uni} = r_{\\rm usd}$. Do you see systematic differences? If so, explain why these occur.\n",
    "\n",
    "3. The apparent magnitude, $m$, is related to the absolute magnitude, $M$, by the relation $m - M = 5\\log_{10}r - 5$.  For a star of given $M$, how, qualitatively, does $m$ vary with $r$? (This is important for understanding the magnitude scale.) Compute $M$ for all stars using each of the distance estimators, and plot the histogram over $M$ in each case, as well as a histogram of the differences. What does this tell you? Plot $M$ vs.\\ $r$ for your two estimators. Decsribe the result. Given that $M$ is an _intrinsic_ property of a star, what does this plot tell you about the Hipparcos survey?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 4: A better prior (5 points) \n",
    "\n",
    "1. What is a significant drawback of the two priors we have considered so far (uniform in $r$ and uniform space density)?\n",
    "2. Can you think of a sensible prior which might avoid this? (explain)"
   ]
  },
  {
   "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}