{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Homework 7\n",
"\n",
"## MLE, MCMC, Interpolation, Expectation Maximization (EM), and Resampling Methods\n",
"
\n",
"This notebook is arranged in cells. Texts are usually written in the markdown cells, and here you can use html tags (make it bold, italic, colored, etc). You can double click on this cell to see the formatting.
\n",
"
\n",
"The ellipsis (...) are provided where you are expected to write your solution but feel free to change the template (not over much) in case this style is not to your taste.
\n",
"
\n",
"Hit \"Shift-Enter\" on a code cell to evaluate it. Double click a Markdown cell to edit.
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***\n",
"### Link Okpy"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from client.api.notebook import Notebook\n",
"ok = Notebook('hw7.ok')\n",
"_ = ok.auth(inline = True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Imports"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy.integrate import quad\n",
"#For plotting\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Supernova Cosmology Project\n",
"\n",
"In this homework, we use a compilation of supernovae data to show that the expansion of the universe is accelerating, and hence it contains dark energy. This is the Nobel prize winning research in 2011 (https://www.nobelprize.org/nobel_prizes/physics/laureates/2011/), and Saul Perlmutter, a professor of physics at Berkeley, shared a prize in 2011 for this discovery.\n",
"
\n",
"\"The expansion history of the universe can be determined quite easily, using as a “standard candle” any distinguishable class of astronomical objects of known intrinsic brightness that can be identified over a wide distance range. As the light from such beacons travels to Earth through an expanding universe, the cosmic expansion stretches not only the distances between galaxy clusters, but also the very wavelengths of the photons en route. By the time the light reaches us, the spectral wavelength $\\lambda$ has thus been redshifted by precisely the same incremental factor $z = \\Delta \\lambda/\\lambda$ by which the cosmos has been stretched in the time interval since the light left its source. The recorded redshift and brightness of each such object thus provide a measurement of the total integrated expansion of the universe since the time the light was emitted. A collection of such measurements, over a sufficient range of distances, would yield an entire historical record of the universe’s expansion.\" (Saul Perlmutter, http://supernova.lbl.gov/PhysicsTodayArticle.pdf).\n",
"
\n",
"Supernovae emerge as extremely promising candidates for measuring the cosmic expansion. Type I Supernovae arises from the collapse of white dwarf stars when the Chandrasekhar limit is reached. Such nuclear chain reaction occurs in the same way and at the same mass, the brightness of these supernovae are always the same. The relationship between the apparent brightness and distance of supernovae depend on the contents and curvature of the universe.\n",
"
\n",
"We can infer the \"luminosity distance\" $D_L$ from measuring the inferred brightness of a supernova of luminosity $L$. Assuming a naive Euclidean approach, if the supernova is observed to have flux $F$, then the area over which the flux is distributed is a sphere radius $D_L$, and hence
\n",
"$$F = \\frac{L}{4\\pi D_L^2}.$$\n",
"
\n",
"In Big Bang cosmology, $D_L$ is given by:\n",
"
\n",
"$$ D_L = \\frac{\\chi(a)}{a} $$\n",
"
\n",
"where $a$ is the scale factor ($\\frac{\\lambda_0}{\\lambda} = 1 + z = \\frac{a_0}{a}$, and the quantity with the subscript 0 means the value at present. Note that $a_0 = 1, z_0 = 0$.), and $\\chi$ is the comoving distance, the distance between two objects as would be measured instantaneously today. For a photon, $cdt = a(t)d\\chi$, so $\\chi(t) = c\\int_t^{t_0} \\frac{dt'}{a(t')}$. We can write this in terms of a Hubble factor ($H(t) = \\frac{1}{a}\\frac{da}{dt}$), which tells you the expansion rate: $\\chi(a) = c\\int_a^1 \\frac{da'}{a'^2H(a')} = c\\int_0^z \\frac{dz'}{H(z')}$. (change of variable using $a = \\frac{1}{1+z}$.)\n",
"
\n",
"Using the Friedmann equation (which basically solves Einstein's equations for a homogenous and isotropic universe), we can write $H^2$ in terms of the mass density $\\rho$ of the components in the universe: $H^2(z) = H_0^2[\\Omega_m(1+z)^3 + (1-\\Omega_m)(1+z)^2].$
\n",
"$\\Omega$ is the density parameter; it is the ratio of the observed density of matter and energy in the universe ($\\rho$) to the critical density $\\rho_c$ at which the universe would halt is expansion. So $\\Omega_0$ (again, the subscript 0 means the value at the present) is the total mass and energy density of the universe today, and consequently $\\Omega_0 = \\Omega_{m}$ (matter density parameter today; remember we obtained the best-fit value of this parameter in Project 1?) = $\\Omega_{\\mathrm{baryonoic\\ matter}}$ + $\\Omega_{\\mathrm{dark\\ matter}}$. If $\\Omega_0 < 1$, the universe will continue to expand forever. If $\\Omega_0 > 1$, the expansion will stop eventually and the universe will start to recollapse. If $\\Omega_0 = 1$, then the universe is flat and contains enough matter to halt the expansion but not enough to recollapse it. So it will continue expanding, but gradually slowing down all the time, finally running out of steam only in the infinite future. Even including dark matter in this calculation, cosmologists found that all the matters in the universe only amounts to about a quarter of the required critical mass, suggesting a continuously expanding universe with deceleration. Then, using all this, we can write the luminosity distance in terms of the density parameters:
\n",
"$$ D_L = \\frac{\\chi(a)}{a} = c(1+z)\\int_0^z \\frac{dz'}{H(z')} = c(1+z)\\int_0^z \\frac{dz'}{H_0[\\Omega_m(1+z')^3 + (1-\\Omega_m)(1+z')^2]^{1/2}} $$
\n",
"$$ = \\frac{2997.92458}{h} (1+z)\\int_0^z \\frac{dz'}{[\\Omega_m(1+z')^3 + (1-\\Omega_m)(1+z')^2]^{1/2}}\\ [unit\\ of\\ Mpc] $$\n",
"
\n",
"where $H_0 = 100\\cdot h\\ [km\\cdot s^{-1} Mpc^{-1}]$.\n",
"
\n",
"Fluxes can be expressed in magnitudes $m$, where $m = -2.5\\cdot\\mathrm{log}_{10}F$ + const. The distance modulus is $\\mu = m - M$ ($M$ is the absolute magnitude, the value of $m$ if the supernova is at a distance 10pc. Then, we have:\n",
"
\n",
"$$ \\mu = 25 + 5\\cdot \\mathrm{log}_{10}\\Big(D_L\\ [in\\ the\\ unit\\ of \\ Mpc]\\Big)$$\n",
"
\n",
"In this assignment, we use the SCP Union2.1 Supernova (SN) Ia compilation. (http://supernova.lbl.gov/union/)\n",
"
\n",
"First, load the measured data: $z$ (redshift), $\\mu$ (distance modulus), $\\sigma(\\mu)$ (error on distance modulus)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"data = np.loadtxt(\"sn_z_mu_dmu_plow_union2.1.txt\", usecols=range(1,5))\n",
"# z\n",
"z_data = data[:,0]\n",
"# mu\n",
"mu_data = data[:,1]\n",
"# error on mu (sigma(mu))\n",
"mu_err_data = data[:,2]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" 1. Plot the measured distance modulus as a function of redshift with errorbars. Then, assume three different scenarios: $\\Omega_m = 0, 0.3, 1.$
\n",
"Remember:\n",
"$$ D_L = \\frac{2997.92458}{h} (1+z)\\int_0^z \\frac{dz'}{[\\Omega_m(1+z')^3 + (1-\\Omega_m)(1+z')^2]^{1/2}}$$
\n",
"$$ \\mu = 25 + 5\\cdot \\mathrm{log}_{10}(D_L)$$
\n",
" Now, plot three curves of $\\mu$ as a function of $z$ for $\\Omega_m = 0, 0.3, 1$ on top of the measured data (Calculate $D_L$ using quad. For now, assume $h = 0.7$.) How do they fit?
\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import math\n",
"\n",
"..."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"plt.figure(figsize = (20,14))\n",
"\n",
"...\n",
"\n",
"plt.legend()\n",
"plt.xlim(0.01, 1.5)\n",
"plt.xlabel('$z$')\n",
"plt.ylabel('$\\mu$')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You should find that the measured data do not fit well to all three scenarios. \"The high-redshift supernovae are fainter than would be expected even for an empty cosmos (corresponding to $\\Omega_m = 0$).\" So what's wrong? \n",
"
\n",
"\"If these data are correct, the obvious implication is that the simplest cosmological model must be too simple. The next simplest model might be one that Einstein entertained for a time. Believing the universe to be static, he tentatively introduced into the equations of general relativity an expansionary term he called the “cosmological constant” ($\\Lambda$) that would compete against gravitational collapse. After Hubble’s discovery of the cosmic expansion, Einstein famously rejected $\\Lambda$ as his “greatest blunder.” In later years, $\\Lambda$ came to be identified with the zero-point vacuum energy of all quantum fields. It turns out that invoking a cosmological constant allows us to fit the supernova data quite well.\" (Saul Perlmutter, https://www.nobelprize.org/nobel_prizes/physics/laureates/2011/)\n",
"
\n",
"So in short, the data indicates that faint supernovae are further away from the earth than had been theoretically expected. The expansion rate of the universe is increasing indeed. It seems that some mysterious material (which we call \"dark energy\") is causing such antigravity effects. The cosmological constant, $\\Lambda$, the value of the energy density of the vacuum of space is widely accepted as a leading candidate of dark energy. \n",
"
\n",
"Now let us add a general form of dark energy to our model.\n",
"
\n",
"$$H^2(z) = H_0^2[\\Omega_m(1+z)^3 + \\Omega_{DE}(1+z)^{3(1+w)} + (1-\\Omega_m-\\Omega_{DE})(1+z)^2].$$
$w$ is the dark energy equation of state, which is the ratio of its pressure to its energy density. $w = -1$ for the cosmological constant $\\Lambda$.
\n",
"$\\Omega_0 = \\Omega_{m}$ (matter density parameter today) + $\\Omega_{DE}$ (dark energy density parameter today), and \n",
"
\n",
"$$ D_L = \\frac{\\chi(a)}{a} = c(1+z)\\int_0^z \\frac{dz'}{H(z')} = c(1+z)\\int_0^z \\frac{dz'}{H_0[\\Omega_m(1+z')^3 + \\Omega_{DE}(1+z')^{3(1+w)} + (1-\\Omega_m-\\Omega_{DE})(1+z')^2]^{1/2}} $$
\n",
"$$ = \\frac{2997.92458}{h} (1+z)\\int_0^z \\frac{dz'}{[\\Omega_m(1+z')^3 + \\Omega_{DE}(1+z')^{3(1+w)} + (1-\\Omega_m-\\Omega_{DE})(1+z')^2]^{1/2}}\\ [unit\\ of\\ Mpc] $$\n",
"
\n",
"where $H_0 = 100\\cdot h\\ [km\\cdot s^{-1} Mpc^{-1}]$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" 2. Now assume three different scenarios: ($\\Omega_m = 0.3, \\Omega_{DE} = 0$), ($\\Omega_m = 0, \\Omega_{DE} = 1, w = -1$), and ($\\Omega_m = 0.3, \\Omega_{DE} = 0.7, w = -1$). Again, plot three curves of $\\mu$ as a function of $z$ on top of data (assume $h = 0.7$)
\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"..."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"plt.figure(figsize = (20,14))\n",
"\n",
"...\n",
"\n",
"plt.legend()\n",
"plt.xlim(0.01, 1.5)\n",
"plt.xlabel('$z$')\n",
"plt.ylabel('$\\mu$')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You basically reproduced the below figure!\n",
"![alt text](SN.png \"Title\")\n",
"You should see that $\\Omega_m = 0.3$ and $\\Omega_m = 0.7$ fits the data best. In combination with the CMB data, this shows that about 70% of the total energy density is vacuum energy and 30% is mass."
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"***"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, with measurements of the distance modulus $\\mu$, use Bayesian analysis to estimate the cosmological parameters.\n",
"
\n",
"let us assume that the universe is flat (which is a fair assumption since the CMB measurements indicate that the universe has no large-scale curvature). $\\Omega_0 = \\Omega_m + \\Omega_{DE} = 1$. Then, we do not need to worry about the curvature term:
\n",
"$$ D_L = \\frac{\\chi(a)}{a} = c(1+z)\\int_0^z \\frac{dz'}{H(z')} = c(1+z)\\int_0^z \\frac{dz'}{H_0[\\Omega_m(1+z')^3 + (1-\\Omega_m)(1+z')^{3(1+w)}]^{1/2}} $$
\n",
"$$ = \\frac{2997.92458}{h} (1+z)\\int_0^z \\frac{dz'}{[\\Omega_m(1+z')^3 + (1-\\Omega_m)(1+z')^{3(1+w)}]^{1/2}}\\ [unit\\ of\\ Mpc] $$
\n",
"where $H_0 = 100\\cdot h\\ [km\\cdot s^{-1} Mpc^{-1}]$.
\n",
"Assuming that errors are Gaussian (can be justified by averaging over large numbers of SN; central limit theorem), we calculate the likelihood $L$ as:
\n",
"$$ L \\propto \\mathrm{exp}\\Big( -\\frac{1}{2} \\sum_{i = 1}^{N_{\\mathrm{SN}}} \\frac{[\\mu_{i,\\ data}(z_i) - \\mu_{i,\\ model}(z_i, \\Omega_m, w)]^2}{\\sigma(\\mu_i)^2} \\Big) $$\n",
"
\n",
"where $z_i, \\mu_i, \\sigma(\\mu_i)$ are from the measurements, and we compute $\\mu_{model}$ as a function of $z, \\Omega_m, w$.\n",
"\n",
"
\n",
"First, try the maximum likelihood estimation (MLE). \n",
"
\n",
"\n",
" 3. Assuming that $h$ = 0.7, find the maximum likelihood estimation of $\\Omega_m$ and $w$ (i.e. find $\\Omega_m$ and $w$ which maximizes the likelihood. \n",
"
\n",
"(Hint: This is very similar to Problem 2-1, HW6. Take the log of the likelihood and maximize it using scipy.optimize.fmin (https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.optimize.fmin.html). Note that you need to make initial guesses on the parameters in order to use fmin. You can set them to be 0. Caveat: \"fmin\" minimizes a given function, so you should multiply the log-likelihood by $-1$ in order to maximize it using fmin.)\n",
"
\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"..."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from scipy import optimize\n",
"\n",
"def minus_log_likelihood(param):\n",
" \n",
" Omegam, w = param\n",
" \n",
" if(Omegam<=0 or w>=0):\n",
" lnL = -1.e100\n",
" else:\n",
" \n",
" ...\n",
" lnL = ...\n",
"\n",
" return -lnL"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"Omegam, w = optimize.fmin(...)\n",
"print('MAP solution')\n",
"print('Omega_m = ', Omegam, ', w = ', w)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, write an MCMC code using the Metropolis algorithm. In this problem, assume that the universe is flat, and $w = -1$ (dark energy is $\\Lambda$.) We call this flat $\\Lambda$CDM cosmology. By fixing $h$ and $w$, $D_L = D_L(z, \\Omega_m)$.\n",
"
\n",
"First, I precalculated $D_L$ for $h = 0.7$ from 2 one-dimensional vectors giving the tabulated values of the parameters $z$ and $\\Omega_m$ in the range $0.01 < z < 1.5$ and $0.01 < \\Omega_m < 1$. We call the tabulated values of $z$ and $\\Omega_m$ as \"z_fit\" (length 200) and \"Om0_fit\" (length 100). Then, \"DL_fit\" is a 2-dimensional grid of tabulated values of $D_L$ (its dimension 200 $\\times$ 100. i.e. $D_L$[i,j] is $D_L$$\\big($$z$ = z_fit[i], $\\Omega_m$ = Om0_fit[j]$\\big)$\n",
"
\n",
"Now using a 2-D spline interpolation, estimate $D_L(z, \\Omega_m)$ for any $z$ and $\\Omega_m$.\n",
"
\n",
" 4. Using scipy.interpolate.RectBivariateSpline, estimate $D_L(z, \\Omega_m)$ for any $z$ and $\\Omega_m$. Plot $\\mu = 25 + 5\\cdot \\mathrm{log}_{10}(D_L(z, \\Omega_m = 0.3))$ as a function of $z$ on top of the measured data. How does it fit to the data? \n",
"
\n",
"(Hint: Let z_spline = RectBivariateSpline(x_fit,y_fit,z_fit). Then, z_spline.ev(x,y) will evaulate the spline at given positions x and y.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"DL_fit = np.loadtxt(\"DL_fit.txt\").T\n",
"Om0_fit = np.loadtxt(\"Om0_fit.txt\")\n",
"z_fit = np.loadtxt(\"z_fit.txt\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from scipy.interpolate import RectBivariateSpline"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"..."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"plt.figure(figsize = (20,14))\n",
"\n",
"...\n",
"\n",
"plt.legend()\n",
"plt.xlim(0.01, 1.5)\n",
"plt.xlabel('$z$')\n",
"plt.ylabel('$\\mu$')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"Now, run the MCMC code to estimate $\\Omega_m$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# The following is a modified version of the code written by Alan Heavens (http://www.imperial.ac.uk/media/imperial-college/research-centres-and-groups/astrophysics/public/icic/data-analysis-workshop/2016/SNcodePython.txt)\n",
"\n",
"# Import data\n",
"data = np.loadtxt(\"sn_z_mu_dmu_plow_union2.1.txt\", usecols=range(1,5))\n",
"# z\n",
"z_data = data[:,0]\n",
"# mu\n",
"mu_data = data[:,1]\n",
"# error on mu (sigma(mu))\n",
"mu_err_data = data[:,2]\n",
"\n",
"# length of MCMC chains\n",
"nsamples = 10000\n",
"# number of parameters\n",
"npars = 1\n",
"\n",
"# Define (gaussian) width of the proposal distribution, one for each parameter. This determines how far you propose jumps\n",
"Sigma = [0.01]\n",
"\n",
"# Number of supernova:\n",
"nSN = len(z_data)\n",
"\n",
"# Declare an empty array of the parameter values of each point. \n",
"# Theta[:,0] stores a trace of the parameter \\Omega_m \n",
"# Theta[:,1] stores log-likelihood values at each point\n",
"Theta = np.empty([nsamples,npars+1])\n",
"\n",
"# Dmu stores mu(data)-mu(theory), temporarily:\n",
"Dmu = np.empty(nSN)\n",
"\n",
"# Random starting point in parameter space\n",
"# Set initial likelihood to low value so next point is accepted (could compute it instead):\n",
"Theta[0,:] = [np.random.uniform(), -1.e100]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In Part 4, you calculated $\\mu$ from the $\\Lambda$CDM model given $z$ and $\\Omega_m$, using $D_L$ from the 2D spline interpolation. Using this result, define a function $\\mu_{model}$ which outputs $\\mu$ from the $\\Lambda$CDM theory model given $z$ and $\\Omega_m$. Then, Dmu[j] = mu_data[j]-mu_model(z_data[j],Omegam), and the log-likelihood is \n",
"
\n",
"$$ \\mathrm{ln}(L) \\approx -\\frac{1}{2} \\sum_{i = 1}^{N_{\\mathrm{SN}}} \\frac{[\\mu_{i,\\ data}(z_i) - \\mu_{i,\\ model}(z_i, \\Omega_m)]^2}{\\sigma(\\mu_i)^2} = -\\frac{1}{2} \\sum_{i = 1}^{N_{\\mathrm{SN}}} \\frac{Dmu_i^2}{\\sigma(\\mu_i)^2} $$\n",
"\n",
" 5. Define a function for mu_model (mu predicted from theory) and lnL (log-likelihood). Then, you can run the MCMC code and plot the posterior (the routine already given).\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def mu_model(z,Omegam):\n",
" ...\n",
" return ...\n",
"\n",
"# Define the likelihood function:\n",
"\n",
"def lnL(Omegam):\n",
"\n",
" # Treat unphysical regions by setting likelihood to (almost) zero: \n",
" if(Omegam<=0):\n",
" lnL = -1.e100\n",
" else:\n",
" \n",
" # Compute difference with theory mu at redshifts of the SN, for trial Omegam\n",
" ...\n",
"\n",
" # Compute ln(likelihood) assuming gaussian errors\n",
" lnL = ...\n",
" \n",
" return lnL\n",
"\n",
"\n",
"# Now run the MCMC code (all routine already given below)\n",
"\n",
"# Draw new proposed samples from a proposal distribution, centred on old values Omegam[i-1]\n",
"# Accept or reject, and colour points according to ln(likelihood):\n",
"\n",
"# Compute initial likelihood value:\n",
"Theta[0,npars] = lnL(Theta[0,0])\n",
"\n",
"progress = nsamples/10; val = 0\n",
"for i in range(1,nsamples): \n",
" \n",
" if i%progress == 0:\n",
" val = val + 10\n",
" print(\"%d percent done\" %val)\n",
" \n",
" lnLPrevious = Theta[i-1,npars]\n",
" OmegamProp = np.random.normal(Theta[i-1,0],Sigma[0])\n",
" \n",
" lnLProp = lnL(OmegamProp)\n",
"\n",
" # Metroplis-Hastings algorithm:\n",
"\n",
" if(lnLProp > lnLPrevious):\n",
" # Accept point if likelihood has gone up:\n",
" Theta[i,0] = OmegamProp\n",
" Theta[i,npars] = lnLProp\n",
" else:\n",
" # Otherwise accept it with probability given by ratio of likelihoods:\n",
" alpha = np.random.uniform()\n",
" \n",
" if(lnLProp - lnLPrevious > np.log(alpha)):\n",
" Theta[i,0] = OmegamProp\n",
" Theta[i,npars] = lnLProp\n",
" else:\n",
" # Reject; Repeat the previous point in the chain:\n",
" Theta[i,0] = Theta[i-1,0]\n",
" Theta[i,npars] = lnLPrevious\n",
"\n",
"# Remove a burn in period, arbitrarily chosen to be the first 10% of the chain:\n",
"nburn = math.floor(nsamples/10)\n",
" \n",
"\n",
"# Plot the histogram of Omegam after the burn-in phase:\n",
"plt.hist(Theta[nburn:,0],bins=30)\n",
"plt.xlabel(r'$\\Omega_m$')\n",
"plt.show()\n",
"\n",
"# Determine the best-fit value and constraint\n",
"print ('Omega_m = ',np.mean(Theta[nburn:nsamples,0]), '+/-' ,np.std(Theta[nburn:nsamples,0]))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
Now, assume a more general form of dark energy, as in Part 3. (Do not fix $w$ to -1; add $w$ as a parameter.)\n",
"\n",
"In the flat universe,
\n",
"$$ D_L = \\frac{\\chi(a)}{a} = c(1+z)\\int_0^z \\frac{dz'}{H(z')} = c(1+z)\\int_0^z \\frac{dz'}{H_0[\\Omega_m(1+z')^3 + (1-\\Omega_m)(1+z')^{3(1+w)}]^{1/2}} $$
\n",
"$$ = \\frac{2997.92458}{h} (1+z)\\int_0^z \\frac{dz'}{[\\Omega_m(1+z')^3 + (1-\\Omega_m)(1+z')^{3(1+w)}]^{1/2}}\\ [unit\\ of\\ Mpc] $$\n",
"
\n",
"where $H_0 = 100\\cdot h\\ [km\\cdot s^{-1} Mpc^{-1}]$. Here, we fix $h = 0.7$.
\n",
"We calculate the likelihood $L$ as:
\n",
"$$ \\mathrm{ln}(L) \\approx -\\frac{1}{2} \\sum_{i = 1}^{N_{\\mathrm{SN}}} \\frac{[\\mu_{i,\\ data}(z_i) - \\mu_{i,\\ model}(z_i, \\Omega_m, w)]^2}{\\sigma(\\mu_i)^2} = -\\frac{1}{2} \\sum_{i = 1}^{N_{\\mathrm{SN}}} \\frac{Dmu_i^2}{\\sigma(\\mu_i)^2} $$\n",
"
\n",
"where $$ \\mu_{i,\\ model}(z_i, \\Omega_m, w) = 25 + 5\\cdot \\mathrm{log}_{10}(D_{L,\\ i})$$
\n",
"$$ D_{L,\\ i} = \\frac{2997.92458}{0.7} (1+z_i)\\int_0^{z_i} \\frac{dz'}{[\\Omega_m(1+z')^3 + (1-\\Omega_m)(1+z')^{3(1+w)}]^{1/2}} $$\n",
"
\n",
" 6. Modify the code in Part 5 (most routine already given) and run the MCMC code to estimate $w$ and $\\Omega_m$. (npars = 2 in this case). Plot 1-d posterior of $w$ and $\\Omega_m$ as well as 2-d posterior (i.e. plot the chain in two-dimensional parameter space. Make sure that the chain has converged (you can change nsamples, nburn).
\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Import data\n",
"data = np.loadtxt(\"sn_z_mu_dmu_plow_union2.1.txt\", usecols=range(1,5))\n",
"# z\n",
"z_data = data[:,0]\n",
"# mu\n",
"mu_data = data[:,1]\n",
"# error on mu (sigma(mu))\n",
"mu_err_data = data[:,2]\n",
"\n",
"# length of MCMC chains\n",
"nsamples = 15000\n",
"# number of parameters\n",
"npars = 2\n",
"\n",
"# Define (gaussian) width of the proposal distribution, one for each parameter. This determines how far you propose jumps\n",
"Sigma = [0.01, 0.01]\n",
"\n",
"# Number of supernova:\n",
"nSN = len(z_data)\n",
"\n",
"# Declare an empty array of the parameter values of each point. \n",
"# Theta[:,0] stores a trace of the parameter \\Omega_m \n",
"# Theta[:,1] stores a trace of the parameter w \n",
"# Theta[:,2] stores log-likelihood values at each point\n",
"Theta = np.empty([nsamples,npars+1])\n",
"\n",
"# Dmu stores mu(data)-mu(theory), temporarily:\n",
"Dmu = np.empty(nSN)\n",
"\n",
"# Random starting point in parameter space\n",
"# Set initial likelihood to low value so next point is accepted (could compute it instead):\n",
"Theta[0,:] = [np.random.uniform(), -np.random.uniform(), -1.e100]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def mu_model(z,Omegam,w):\n",
" ...\n",
" return ..."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Define the likelihood function:\n",
"\n",
"def lnL(Omegam, w):\n",
"\n",
" # Treat unphysical regions by setting likelihood to (almost) zero: \n",
" if(Omegam<=0 or w>=0):\n",
" lnL = -1.e100\n",
" else:\n",
" \n",
" # Compute difference with theory mu at redshifts of the SN, for trial Omegam\n",
" ...\n",
"\n",
" # Compute ln(likelihood) assuming gaussian errors\n",
" lnL = ...\n",
" \n",
" return lnL\n",
"\n",
"\n",
"# Draw new proposed samples from a proposal distribution, centred on old values Omegam[i-1]\n",
"# Accept or reject, and colour points according to ln(likelihood):\n",
"\n",
"# Compute initial likelihood value:\n",
"Theta[0,npars] = lnL(Theta[0,0], Theta[0,1])\n",
"\n",
"progress = nsamples/10; val = 0\n",
"for i in range(1,nsamples): \n",
" \n",
" if i%progress == 0:\n",
" val = val + 10\n",
" print(\"%d percent done\" %val)\n",
" \n",
" lnLPrevious = Theta[i-1,npars]\n",
" OmegamProp = np.random.normal(Theta[i-1,0],Sigma[0])\n",
" wProp = np.random.normal(Theta[i-1,1],Sigma[1])\n",
" \n",
" lnLProp = lnL(OmegamProp, wProp)\n",
"\n",
" # Metroplis-Hastings algorithm:\n",
"\n",
" if(lnLProp > lnLPrevious):\n",
" # Accept point if likelihood has gone up:\n",
" Theta[i,0] = OmegamProp\n",
" Theta[i,1] = wProp\n",
" Theta[i,npars] = lnLProp\n",
" else:\n",
" # Otherwise accept it with probability given by ratio of likelihoods:\n",
" alpha = np.random.uniform()\n",
" \n",
" if(lnLProp - lnLPrevious > np.log(alpha)):\n",
" Theta[i,0] = OmegamProp\n",
" Theta[i,1] = wProp\n",
" Theta[i,npars] = lnLProp\n",
" else:\n",
" # Reject; Repeat the previous point in the chain:\n",
" Theta[i,0:2] = Theta[i-1,0:2]\n",
" Theta[i,npars] = lnLPrevious\n",
"\n",
"# Remove a burn in period, arbitrarily chosen to be the first 20% of the chain:\n",
"nburn = 2*math.floor(nsamples/10)\n",
" \n",
"\n",
"# Plot the histogram of Omegam after the burn-in phase:\n",
"plt.hist(Theta[nburn:,0],bins=30)\n",
"plt.xlabel(r'$\\Omega_m$')\n",
"plt.show()\n",
"\n",
"# Plot the histogram of w after the burn-in phase:\n",
"plt.hist(Theta[nburn:,1],bins=30)\n",
"plt.xlabel('w')\n",
"plt.show()\n",
"\n",
"# Scatter plot of the samples (2-d posterior):\n",
"plt.scatter(Theta[nburn:,0], Theta[nburn:,1], c = -Theta[nburn:,npars])\n",
"plt.xlabel(r'$\\Omega_m$')\n",
"plt.ylabel('w')\n",
"plt.show() \n",
"\n",
"# Print best-fit values and constraints\n",
"print ('Omega_m = ',np.mean(Theta[nburn:nsamples,0]), '+/-' ,np.std(Theta[nburn:nsamples,0]))\n",
"print ('w = ',np.mean(Theta[nburn:nsamples,1]), '+/-' ,np.std(Theta[nburn:nsamples,1]))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"Now, include the distance modulus of 12 additional supernovae, which are not-so-good standard candles. They are 3$\\sigma$ away from the best-fit mode. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"data = np.loadtxt(\"sn_z_mu_dmu_plow_union2.1_outlier.txt\", usecols=range(1,5))\n",
"# z\n",
"z_data = data[:,0]\n",
"# mu\n",
"mu_data = data[:,1]\n",
"# error on mu (sigma(mu))\n",
"mu_err_data = data[:,2]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So we have a total of 592 supernovae, and we can see that the last 12 supernovae seem to be outliers. (i.e. mu_data[580:] contains the distance modulus measurements of these 12 supernovae.)\n",
"
\n",
" 7. First run the MCMC code in Part 6 with the new data (total of 592 supernovae). Then, using the estimates of $\\Omega_m$ and $w$ from the MCMC chain, calculate the distance modulus from theory and plot the curve on top of the measured data. Plot the measurements of the last 12 supernovae with different color.
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"..."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"plt.figure(figsize = (20,14))\n",
"\n",
"...\n",
"\n",
"plt.legend()\n",
"plt.xlim(0.01, 1.5)\n",
"plt.xlabel('$z$')\n",
"plt.ylabel('$\\mu$')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"Remember that in HW6, we used the Gaussian mixture to better model the measurements with outliers. Let us apply the same technique in this case.\n",
"
\n",
"$$ L = \\prod_{i = 1}^{N_{\\mathrm{SN}}} \\Big[ \\frac{g}{\\sqrt{2\\pi\\sigma(\\mu_i)^2}}\\mathrm{exp}\\Big( -\\frac{1}{2} \\frac{[\\mu_{i,\\ data}(z_i) - \\mu_{i,\\ model}(z_i, \\Omega_m, w)]^2}{\\sigma(\\mu_i)^2} \\Big) + \\frac{1-g}{\\sqrt{2\\pi\\sigma_B^2}}\\mathrm{exp}\\Big( -\\frac{1}{2} \\frac{[\\mu_{i,\\ data}(z_i) - \\mu_{i,\\ model}(z_i, \\Omega_m, w) + \\Delta \\mu]^2}{\\sigma_B^2} \\Big) \\Big] $$\n",
"
\n",
"Here, we have 5 free parameters: $\\Omega_m, w, g, \\sigma_B, \\Delta \\mu$.\n",
"
\n",
"With outliers, we think there is something in the noise we do not really understand, which makes error distribution non-Gaussian. So we hope adding a second Gaussian to the model would better describe the pdf. $g$ determines weights on the two Gaussians. $\\sigma_B^2$ is the variance of the second Gaussian, which we assume to be larger than the variance of the first Gaussian. $\\Delta \\mu$ is the distance modulus offset in the second Gaussian.\n",
"
\n",
"\n",
" 8. Re-run the MCMC code with this new model. Plot 1-d and 2-d constrains of $\\Omega_m$ and $w$ as in Part 6 and 7.
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Setup (all routine already given)\n",
"\n",
"data = np.loadtxt(\"sn_z_mu_dmu_plow_union2.1_outlier.txt\", usecols=range(1,5))\n",
"# z\n",
"z_data = data[:,0]\n",
"# mu\n",
"mu_data = data[:,1]\n",
"# error on mu (sigma(mu))\n",
"mu_err_data = data[:,2]\n",
"\n",
"# length of MCMC chains\n",
"nsamples = 15000\n",
"# number of parameters\n",
"npars = 5\n",
"\n",
"# Define (gaussian) width of the proposal distribution, one for each parameter. This determines how far you propose jumps\n",
"Sigma = 0.01*np.ones(npars)\n",
"Sigma[2] = 0.03*np.ones(1)\n",
"Sigma[3] = 0.1*np.ones(1)\n",
"Sigma[4] = 0.01*np.ones(1)\n",
"\n",
"# Number of supernova:\n",
"nSN = len(z_data)\n",
"\n",
"# Declare an empty array of the parameter values of each point. \n",
"# Theta[:,0] stores a trace of the parameter \\Omega_m \n",
"# Theta[:,1] stores a trace of the parameter w \n",
"# Theta[:,2] stores a trace of the parameter g\n",
"# Theta[:,3] stores a trace of the parameter sigma_B \n",
"# Theta[:,4] stores a trace of the parameter delta mu\n",
"# Theta[:,5] stores log-likelihood values at each point\n",
"Theta = np.empty([nsamples,npars+1])\n",
"\n",
"# Dmu stores mu(data)-mu(theory), temporarily:\n",
"Dmu = np.empty(nSN)\n",
"\n",
"# Random starting point in parameter space\n",
"# Set initial likelihood to low value so next point is accepted (could compute it instead):\n",
"Theta[0,:2] = [np.random.uniform(), -np.random.uniform()]\n",
"Theta[0,2] = np.random.normal(0.5, 0.1)\n",
"Theta[0,3] = np.random.normal(10, 0.1)\n",
"Theta[0,4] = np.random.uniform()\n",
"Theta[0,npars] = -1.e100"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Define mu from theory\n",
"def mu_model(z,Omegam,w):\n",
" ...\n",
" return ...\n",
"\n",
"# Define the likelihood function:\n",
"\n",
"def lnL(Omegam, w, sigmaB, gs, errs, muoffset):\n",
"\n",
" # Treat unphysical regions by setting likelihood to (almost) zero: \n",
" if(Omegam<=0 or w>=0 or gs < 0 or gs > 1.):\n",
" lnL = -1.e100\n",
" else:\n",
" \n",
" # Compute difference with theory mu at redshifts of the SN, for trial Omegam\n",
" ...\n",
" lnL = ...\n",
" \n",
" return lnL\n",
"\n",
"\n",
"# Run the MCMC code\n",
"\n",
"# Draw new proposed samples from a proposal distribution, centred on old values Omegam[i-1]\n",
"# Accept or reject, and colour points according to ln(likelihood):\n",
"\n",
"# Compute initial likelihood value:\n",
"Theta[0,npars] = lnL(Theta[0,0], Theta[0,1], Theta[0,3], Theta[0,2], mu_err_data, Theta[0,4])\n",
"\n",
"progress = nsamples/10; val = 0\n",
"for i in range(1,nsamples): \n",
" \n",
" if i%progress == 0:\n",
" val = val + 10\n",
" print(\"%d percent done\" %val)\n",
" \n",
" lnLPrevious = Theta[i-1,npars]\n",
" OmegamProp = np.random.normal(Theta[i-1,0],Sigma[0])\n",
" wProp = np.random.normal(Theta[i-1,1],Sigma[1])\n",
" gvalProp = np.random.normal(Theta[i-1,2],Sigma[2])\n",
" sigmaProp = np.random.normal(Theta[i-1,3],Sigma[3])\n",
" offsetProp = np.random.normal(Theta[i-1,4],Sigma[4])\n",
" \n",
" lnLProp = lnL(OmegamProp,wProp,sigmaProp,gvalProp,mu_err_data,offsetProp)\n",
"\n",
" # Metroplis-Hastings algorithm:\n",
"\n",
" if(lnLProp > lnLPrevious):\n",
" # Accept point if likelihood has gone up:\n",
" Theta[i,0] = OmegamProp\n",
" Theta[i,1] = wProp\n",
" Theta[i,2] = gvalProp\n",
" Theta[i,3] = sigmaProp\n",
" Theta[i,4] = offsetProp\n",
" Theta[i,npars] = lnLProp\n",
" else:\n",
" # Otherwise accept it with probability given by ratio of likelihoods:\n",
" alpha = np.random.uniform()\n",
" \n",
" if(lnLProp - lnLPrevious > np.log(alpha)):\n",
" Theta[i,0] = OmegamProp\n",
" Theta[i,1] = wProp\n",
" Theta[i,2] = gvalProp\n",
" Theta[i,3] = sigmaProp\n",
" Theta[i,4] = offsetProp\n",
" Theta[i,npars] = lnLProp\n",
" else:\n",
" # Reject; Repeat the previous point in the chain:\n",
" Theta[i,0:5] = Theta[i-1,0:5]\n",
" Theta[i,npars] = lnLPrevious\n",
" \n",
"# Remove a burn in period, arbitrarily chosen to be the first 40% of the chain:\n",
"nburn = 4*math.floor(nsamples/10)\n",
" \n",
"\n",
"# Plot the histogram of Omegam after the burn-in phase:\n",
"plt.hist(Theta[nburn:,0],bins=30)\n",
"plt.xlabel(r'$\\Omega_m$')\n",
"plt.show()\n",
"\n",
"# Plot the histogram of w after the burn-in phase:\n",
"plt.hist(Theta[nburn:,1],bins=30)\n",
"plt.xlabel('w')\n",
"plt.show()\n",
"\n",
"# Scatter plot of the samples (2-d posterior):\n",
"plt.scatter(Theta[nburn:,0], Theta[nburn:,1], c = -Theta[nburn:,npars])\n",
"plt.xlabel(r'$\\Omega_m$')\n",
"plt.ylabel('w')\n",
"plt.show() \n",
"\n",
"# Print best-fit values and constraints\n",
"print ('Omega_m = ',np.mean(Theta[nburn:nsamples,0]), '+/-' ,np.std(Theta[nburn:nsamples,0]))\n",
"print ('w = ',np.mean(Theta[nburn:nsamples,1]), '+/-' ,np.std(Theta[nburn:nsamples,1]))\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Reference: See pg. 8-16 (https://lear.inrialpes.fr/~jegou/bishopreadinggroup/chap9.pdf)\n",
"
\n",
"For this Gaussian mixture model, we wish to maximize the likelihood function with respect to the parameters $g, \\sigma_B, \\Delta \\mu$ for $\\Omega_m = 0.3, w = -1$. In order to do this, we will apply the expectation-maximization (EM) algorithm. This is an iterative method to find maximum likelihood in the case where the model depends on the hidden/latent variable. Here, we call binary variable a as our latent variable such that $p(a_k = 1) = \\pi_k$\n",
"
\n",
"Re-write the likelihood as:\n",
"
\n",
"$$ L = \\prod_{i = 1}^{N_{\\mathrm{SN}}} \\Big[ \\frac{\\pi_1}{\\sqrt{2\\pi\\sigma(\\mu_i)^2}}\\mathrm{exp}\\Big( -\\frac{1}{2} \\frac{[\\mu_{i,\\ data}(z_i) - \\mu_{i,\\ model}(z_i, \\Omega_m = 0.3, w = -1)]^2}{\\sigma(\\mu_i)^2} \\Big) + \\frac{\\pi_2}{\\sqrt{2\\pi\\sigma_B^2}}\\mathrm{exp}\\Big( -\\frac{1}{2} \\frac{[\\mu_{i,\\ data}(z_i) - \\mu_{i,\\ model}(z_i, \\Omega_m = 0.3, w = -1) - \\mu_{\\mathrm{offset}}]^2}{\\sigma_B^2} \\Big) \\Big] $$\n",
"
\n",
"$$ = \\prod_{i = 1}^{N_{\\mathrm{SN}}} \\Big[ \\pi_1 \\cdot \\mathrm{Normal}\\big(\\Delta \\mu_i = \\mu_{i,\\ data}-\\mu_{i,\\ model}\\big|\\ \\overline{\\Delta \\mu}_{\\mathrm{class\\ 1}} = 0, \\sigma(\\mu_i)^2 \\big) + \\pi_2 \\cdot \\mathrm{Normal}\\big(\\Delta \\mu_i = \\mu_{i,\\ data}-\\mu_{i,\\ model}\\big|\\ \\overline{\\Delta \\mu}_{\\mathrm{class\\ 2}} = \\mu_{\\mathrm{offset}}, \\sigma_B^2 \\big) \\Big] $$\n",
"where $\\mu_{i,\\ model}$ assumes $\\Omega_m = 0.3, w = -1$. Suppose that we measure $\\Delta \\mu = \\mu_{i,\\ data}-\\mu_{i,\\ model}$. For the first Gaussian (expected to describe the distribution of 580 non-outlier, standard-candle supernovae), the mean value of $\\Delta \\mu$ is $0$, and its variance is the measurement noise $\\sigma(\\mu)^2$. For the second Gaussian which expects to describe the distribution of 12 outliers, we assume that there will be some offset in $\\mu$ ($\\mu_{\\mathrm{offset}}$), so the mean value of $\\Delta \\mu$ is $\\mu_{\\mathrm{offset}}$, and it has some unknown variance $\\sigma_B^2$.\n",
"
\n",
"Now apply the EM algorithm.\n",
"
\n",
"1. First, initialize: choose $\\pi_1 = 0.95$ and $\\pi_2 = 0.05$. Let $\\mu_{\\mathrm{offset}} = 0, \\sigma_B = 0.5$ initially. \n",
"
\n",
"2. Expectation (E) step: Evaluate the responsibilities using the current parameter values.\n",
"
\n",
"$$ \\gamma_{1,\\ i} = \\frac{\\pi_1 \\cdot \\mathrm{Normal}\\big(\\Delta \\mu_i \\big|\\ \\overline{\\Delta \\mu}_{\\mathrm{class\\ 1}}, \\sigma(\\mu_i)^2 \\big)}{\\pi_1 \\cdot \\mathrm{Normal}\\big(\\Delta \\mu_i \\big|\\ \\overline{\\Delta \\mu}_{\\mathrm{class\\ 1}}, \\sigma(\\mu_i)^2 \\big) + \\pi_2 \\cdot \\mathrm{Normal}\\big(\\Delta \\mu_i \\big|\\ \\overline{\\Delta \\mu}_{\\mathrm{class\\ 2}}, \\sigma_B^2 \\big)} $$\n",
"$$ \\gamma_{2,\\ i} = \\frac{\\pi_2 \\cdot \\mathrm{Normal}\\big(\\Delta \\mu_i \\big|\\ \\overline{\\Delta \\mu}_{\\mathrm{class\\ 2}}, \\sigma_B^2 \\big)}{\\pi_1 \\cdot \\mathrm{Normal}\\big(\\Delta \\mu_i \\big|\\ \\overline{\\Delta \\mu}_{\\mathrm{class\\ 1}}, \\sigma(\\mu_i)^2 \\big) + \\pi_2 \\cdot \\mathrm{Normal}\\big(\\Delta \\mu_i \\big|\\ \\overline{\\Delta \\mu}_{\\mathrm{class\\ 2}}, \\sigma_B^2 \\big)} $$\n",
"
\n",
"where $i = 1, ..., N_{SN}$ (number of measurements). Note that $\\gamma_{1}$ and $\\gamma_{2}$ are vectors of length $N_{SN}$. For a supernova $i$, $\\gamma_{1,\\ i}$ describes its probability of belonging to the first class (described by the first Gaussian). (Note: Therefore, in the end, we expect 12 outliers have much higher values of $\\gamma_{2}$ than normal 580 supernovae - i.e. they have much greater probability of belonging to the second class. This is a systematic way to identify an outlier.)\n",
"
\n",
"3. Maximization (M) step: Re-estimate the parameters using the current responsibilities\n",
"
\n",
"$$ \\mathrm{The\\ mean\\ (\\overline{\\Delta \\mu}_{\\mathrm{class\\ 1}} = 0)\\ and\\ variance,\\ \\sigma(\\mu_i)^2,\\ of\\ the\\ first\\ Gaussian\\ are\\ fixed\\ at\\ initial\\ values} $$\n",
"$$ N_1 = \\sum_{i=1}^{N_{SN}} \\gamma_{1,\\ i},\\ \\ N_2 = \\sum_{i=1}^{N_{SN}} \\gamma_{2,\\ i}$$\n",
"$$ \\overline{\\Delta \\mu}_{\\mathrm{class\\ 2}} = \\frac{1}{N_2} \\sum_{i=1}^{N_{SN}} \\gamma_{2,\\ i} \\cdot \\Delta \\mu_i $$\n",
"$$ \\sigma_B^2 = \\frac{1}{N_2} \\sum_{i=1}^{N_{SN}} \\gamma_{2,\\ i} \\cdot (\\Delta \\mu_i - \\overline{\\Delta \\mu}_{\\mathrm{class\\ 2}})^2$$\n",
"$$ \\pi_1 = \\frac{N_1}{{N_{SN}}},\\ \\ \\pi_2 = \\frac{N_2}{{N_{SN}}} $$\n",
"
\n",
"4. Evaluate the log-likelihood and check for convergence of either the parameters or the log likelihood. If the convergence criterion is not satisfied return to step 2.\n",
"
\n",
"\n",
" 9. Using EM, calculate the converged values of $\\pi_1$, $\\pi_2$, and $N_2$. $N_2$ is the total number of SN in the second class (can be identified as outliers). Iterate until you reach the convergence (parameters not changing). Then, print out the values of $\\gamma_{2}$ and show that 12 outliers have higher values of $\\gamma_{2}$ than other supernovae.
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"data = np.loadtxt(\"sn_z_mu_dmu_plow_union2.1_outlier.txt\", usecols=range(1,5))\n",
"# z\n",
"z_data = data[:,0]\n",
"# mu\n",
"mu_data = data[:,1]\n",
"# error on mu (sigma(mu))\n",
"mu_err_data = data[:,2]\n",
"\n",
"...\n",
"..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we use a Boostrap resampling method to estimate the posterior of $\\Omega_m$ and $w$.\n",
"
\n",
"Suppose that we have 10 measurements of $x$: [3.7, 3.2, 3.3, 3.1, 3.2, 3.5, 2.9, 3.4, 3.0, 3.1]. Now, randomly take 5 samples of 10 data measurements \"with replacement.\""
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"After bootstrap re-sampling\n",
"[[ 2.9 3.1 3. 3. 3.2 3.1 3.3 3.5 3.4 2.9]\n",
" [ 3.4 3.3 3. 3.2 3.1 3.1 3.5 3.1 3.1 3.4]\n",
" [ 3.2 3.2 3.2 3.2 3.7 3.1 3.2 3.2 3.3 3. ]\n",
" [ 2.9 2.9 2.9 3.2 3. 3.2 3.5 3.1 3.7 3.5]\n",
" [ 3.2 2.9 3.2 3.2 3.3 3.3 3.2 3.2 3.4 3.7]]\n"
]
}
],
"source": [
"x = np.array([3.7, 3.2, 3.3, 3.1, 3.2, 3.5, 2.9, 3.4, 3.0, 3.1])\n",
"\n",
"num_samples = 5\n",
"len_x = len(x)\n",
"idx = np.random.randint(0, len_x, (num_samples, len_x))\n",
"print(\"After bootstrap re-sampling\")\n",
"print(x[idx])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Say you wish to see the probability distribution of $\\bar{x}$. Then, take 100 samples using bootstrap and plot the histogram of $\\bar{x}$.\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEMCAYAAADK231MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAElFJREFUeJzt3X2QXXddx/H3xwaqULWt2YbYtG7UFKZFEFxrxYcpFmyl\nDCkjg0HRqJ3JqBUfBsRURzs+ZCY+jE8oOhmoxJFpJ4PVZqigNYLV0bZsoUDTB5uxLU1Nm0UUxIdi\nytc/7kF3lmz27j139+7++n7NZPac3znn3u9vTvLZX3733HNSVUiS2vUFky5AkrSyDHpJapxBL0mN\nM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4zZMugCAjRs31vT09KTLkKR15a677vp4VU0ttd+a\nCPrp6WlmZ2cnXYYkrStJHhlmP6duJKlxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z\n6CWpcWvim7HSUqZ339Lr+If3XjmmSqT1xxG9JDXOoJekxhn0ktS4JYM+yfVJjie5Z0H7G5Lcn+Rw\nkl+d135tkiNJHkhy+UoULUka3jAfxr4D+F3gjz7XkOSlwHbghVX1ZJJzuvYLgR3ARcCXA3+V5IKq\nemrchUuShrPkiL6qbgM+saD5h4G9VfVkt8/xrn07cGNVPVlVDwFHgIvHWK8kaZlGvbzyAuBbkuwB\n/ht4U1V9ADgXuH3efke7ts+TZBewC+D8888fsQytJ30vkZQ0mlE/jN0AnA1cAvwUcCBJlvMCVbWv\nqmaqamZqasknYUmSRjRq0B8FbqqBO4HPAhuBx4Dz5u23pWuTJE3IqEH/Z8BLAZJcADwT+DhwENiR\n5PQkW4FtwJ3jKFSSNJol5+iT3ABcCmxMchS4DrgeuL675PIzwM6qKuBwkgPAvcAJ4BqvuJGkyVoy\n6KvqdYtsev0i++8B9vQpSpI0Pn4zVpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4\ng16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYtGfRJrk9yvHvIyMJtb0xSSTbOa7s2yZEk\nDyS5fNwFS5KWZ5gR/TuAKxY2JjkP+HbgY/PaLgR2ABd1x7w1yWljqVSSNJIlg76qbgM+cZJNvwm8\nGah5bduBG6vqyap6CDgCXDyOQiVJoxlpjj7JduCxqvrwgk3nAo/OWz/atUmSJmTJZ8YulORZwM8w\nmLYZWZJdwC6A888/v89LaRVN775l0iVIWqZRRvRfBWwFPpzkYWAL8MEkzwEeA86bt++Wru3zVNW+\nqpqpqpmpqakRypAkDWPZQV9VH62qc6pquqqmGUzPvLiqHgcOAjuSnJ5kK7ANuHOsFUuSlmWYyytv\nAP4BeG6So0muXmzfqjoMHADuBd4LXFNVT42rWEnS8i05R19Vr1ti+/SC9T3Ann5lSZLGxW/GSlLj\nDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6g\nl6TGGfSS1DiDXpIaN8wTpq5PcjzJPfPafi3J/Uk+kuRPk5w5b9u1SY4keSDJ5StVuCRpOMOM6N8B\nXLGg7Vbg+VX1AuAfgWsBklwI7AAu6o55a5LTxlatJGnZlgz6qroN+MSCtr+sqhPd6u3Alm55O3Bj\nVT1ZVQ8BR4CLx1ivJGmZxjFH/4PAe7rlc4FH52072rVJkiakV9An+VngBPDOEY7dlWQ2yezc3Fyf\nMiRJpzBy0Cf5fuCVwPdUVXXNjwHnzdttS9f2eapqX1XNVNXM1NTUqGVIkpYwUtAnuQJ4M/CqqvrP\neZsOAjuSnJ5kK7ANuLN/mZKkUW1YaockNwCXAhuTHAWuY3CVzenArUkAbq+qH6qqw0kOAPcymNK5\npqqeWqniJUlLWzLoq+p1J2l++yn23wPs6VOUJGl8/GasJDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJ\napxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuCXvdSM93U3vvmXkYx/ee+UYK5FG44hekhpn0EtS4wx6\nSWqcQS9JjVsy6JNcn+R4knvmtZ2d5NYkD3Y/z5q37dokR5I8kOTylSpckjScYUb07wCuWNC2GzhU\nVduAQ906SS4EdgAXdce8NclpY6tWkrRsSwZ9Vd0GfGJB83Zgf7e8H7hqXvuNVfVkVT0EHAEuHlOt\nkqQRjDpHv6mqjnXLjwObuuVzgUfn7Xe0a5MkTUjvD2OrqoBa7nFJdiWZTTI7NzfXtwxJ0iJGDfon\nkmwG6H4e79ofA86bt9+Wru3zVNW+qpqpqpmpqakRy5AkLWXUoD8I7OyWdwI3z2vfkeT0JFuBbcCd\n/UqUJPWx5L1uktwAXApsTHIUuA7YCxxIcjXwCPBagKo6nOQAcC9wArimqp5aodolSUNYMuir6nWL\nbLpskf33AHv6FCVJGh+/GStJjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuN8OLieFvo84Fta\n7xzRS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMZ5eeXTkJcaSk8vjuglqXG9gj7JTyY5nOSeJDck\n+cIkZye5NcmD3c+zxlWsJGn5Rg76JOcCPwbMVNXzgdOAHcBu4FBVbQMOdeuSpAnpO3WzAfiiJBuA\nZwH/DGwH9nfb9wNX9XwPSVIPIwd9VT0G/DrwMeAY8Mmq+ktgU1Ud63Z7HNjUu0pJ0sj6TN2cxWD0\nvhX4cuDZSV4/f5+qKqAWOX5Xktkks3Nzc6OWIUlaQp/LK18GPFRVcwBJbgJeAjyRZHNVHUuyGTh+\nsoOrah+wD2BmZuakvwyk9a7PpawP771yjJXo6azPHP3HgEuSPCtJgMuA+4CDwM5un53Azf1KlCT1\nMfKIvqruSPIu4IPACeBDDEboZwAHklwNPAK8dhyFSpJG0+ubsVV1HXDdguYnGYzuJUlrgN+MlaTG\nGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxB\nL0mNM+glqXEGvSQ1rlfQJzkzybuS3J/kviTfmOTsJLcmebD7eda4ipUkLV/fEf1vA++tqucBL2Tw\nzNjdwKGq2gYc6tYlSRMyctAn+VLgW4G3A1TVZ6rq34DtwP5ut/3AVX2LlCSNrs+IfiswB/xhkg8l\neVuSZwObqupYt8/jwKa+RUqSRtcn6DcALwZ+v6peBPwHC6ZpqqqAOtnBSXYlmU0yOzc316MMSdKp\nbOhx7FHgaFXd0a2/i0HQP5Fkc1UdS7IZOH6yg6tqH7APYGZm5qS/DHRy07tvmXQJktaRkUf0VfU4\n8GiS53ZNlwH3AgeBnV3bTuDmXhVKknrpM6IHeAPwziTPBP4J+AEGvzwOJLkaeAR4bc/3kCT10Cvo\nq+puYOYkmy7r87qSpPHxm7GS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4\ng16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMb1DvokpyX5UJJ3d+tnJ7k1yYPdz7P6lylJ\nGtU4RvQ/Dtw3b303cKiqtgGHunVJ0oT0CvokW4ArgbfNa94O7O+W9wNX9XkPSVI/fR8O/lvAm4Ev\nnte2qaqOdcuPA5t6vof0tDS9+5aRj31475VjrETr3cgj+iSvBI5X1V2L7VNVBdQix+9KMptkdm5u\nbtQyJElL6DN1803Aq5I8DNwIfFuSPwaeSLIZoPt5/GQHV9W+qpqpqpmpqakeZUiSTmXkoK+qa6tq\nS1VNAzuAv66q1wMHgZ3dbjuBm3tXKUka2UpcR78XeHmSB4GXdeuSpAnp+2EsAFX1fuD93fK/AJeN\n43UlSf35zVhJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuLFcXqnl63MfE0laDkf0ktQ4g16SGmfQ\nS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuD4PBz8vyfuS3JvkcJIf79rPTnJrkge7n2eN\nr1xJ0nL1GdGfAN5YVRcClwDXJLkQ2A0cqqptwKFuXZI0IX0eDn6sqj7YLf87cB9wLrAd2N/tth+4\nqm+RkqTRjWWOPsk08CLgDmBTVR3rNj0ObBrHe0iSRtM76JOcAfwJ8BNV9an526qqgFrkuF1JZpPM\nzs3N9S1DkrSIXkGf5BkMQv6dVXVT1/xEks3d9s3A8ZMdW1X7qmqmqmampqb6lCFJOoWR70efJMDb\ngfuq6jfmbToI7AT2dj9v7lWhpGXr+7yDh/deOaZKtBb0efDINwHfC3w0yd1d288wCPgDSa4GHgFe\n269ESVIfIwd9Vf0dkEU2Xzbq60qSxstvxkpS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiD\nXpIaZ9BLUuMMeklqnEEvSY0z6CWpcX3uXvm01/dWsNJa1efvtrc4Xnsc0UtS4wx6SWqcUzeSxspp\nn7VnxYI+yRXAbwOnAW+rqr0r9V59OM8uqXUrEvRJTgN+D3g5cBT4QJKDVXXvSryfpDb4v4GVsVIj\n+ouBI1X1TwBJbgS2AysS9I7KJU3yl8Ra/wW1Uh/Gngs8Om/9aNcmSVplE/swNskuYFe3+ukkD/R4\nuY3Ax/tXtWbYn7WvtT611h9YRp/yKytcyXje+2T9+YphDlypoH8MOG/e+pau7f9U1T5g3zjeLMls\nVc2M47XWAvuz9rXWp9b6A+31qU9/Vmrq5gPAtiRbkzwT2AEcXKH3kiSdwoqM6KvqRJIfBf6CweWV\n11fV4ZV4L0nSqa3YHH1V/Tnw5yv1+guMZQpoDbE/a19rfWqtP9Ben0buT6pqnIVIktYY73UjSY1b\nF0Gf5AuT3Jnkw0kOJ/mFk+yTJL+T5EiSjyR58SRqHdaQfXpekn9I8mSSN02izmEN2Z/v6c7NR5P8\nfZIXTqLWYQ3Zp+1dn+5OMpvkmydR6zCG6c+8fb8+yYkkr1nNGpdryHN0aZJPdufo7iQ/P4lahzHs\nOer6dHe3z98s+cJVteb/AAHO6JafAdwBXLJgn1cA7+n2vQS4Y9J1j6FP5wBfD+wB3jTpmsfQn5cA\nZ3XL39HIOTqD/58CfQFw/6Tr7tOfbttpwF8z+IztNZOuewzn6FLg3ZOudYz9OZPBXQbO79bPWep1\n18WIvgY+3a0+o/uz8MOF7cAfdfveDpyZZPNq1rkcw/Spqo5X1QeA/1nt+pZryP78fVX9a7d6O4Pv\nV6xZQ/bp09X9awOevXD7WjLkvyOANwB/AhxfrdpGtYw+rQtD9ue7gZuq6mPdMUuep3UR9DC4UVqS\nuxn85bu1qu5YsMu6u+3CEH1aV5bZn6sZ/A9sTRumT0leneR+4BbgB1e7xuVYqj9JzgVeDfz+JOob\nxZB/717STbG9J8lFq1zisgzRnwuAs5K8P8ldSb5vqddcN0FfVU9V1dcyGAVenOT5k66pr9b6NGx/\nkryUQdD/9GrWN4ph+lRVf1pVzwOuAn5ptWtcjiH681vAT1fVZ1e/utEM0acPMpjmeAHwFuDPVrvG\n5RiiPxuArwOuBC4Hfi7JBad6zXUT9J9TVf8GvA+4YsGmJW+7sFadok/r0qn6k+QFwNuA7VX1L6td\n26iGOUdVdRvwlUk2rlphIzpFf2aAG5M8DLwGeGuSq1a5vJEs1qeq+tTnpkNq8P2eZ6zzc3QU+Iuq\n+o+q+jhwG3DKCxvWRdAnmUpyZrf8RQzuc3//gt0OAt/XXX1zCfDJqjq2yqUObcg+rRvD9CfJ+cBN\nwPdW1T+ufpXLM2SfvjpJuuUXA6cDa/IX2DD9qaqtVTVdVdPAu4Afqao1OwIe8hw9Z945uphB7q3b\ncwTcDHxzkg1JngV8A3DfqV53vTxKcDOwP4MHmnwBcKCq3p3khwCq6g8YXCHwCuAI8J/AD0yq2CEt\n2ackzwFmgS8BPpvkJ4ALq+pTE6t6ccOco58HvozBKBHgRK3tm04N06fvZDDA+B/gv4Dvmvfh7Foz\nTH/Wm2H69Brgh5OcYHCOdqznc1RV9yV5L/AR4LMMnuB3z6le1G/GSlLj1sXUjSRpdAa9JDXOoJek\nxhn0ktQ4g16SGmfQS1Lj1st19NJEJPlb4ItZ+9f8S4vyOnpJapxTN9Iikrwvycu75V9O8pZJ1ySN\nwqkbaXHXAb+Y5BzgRcCrJlyPNBKnbqRT6B7TdgZwaVX9+6TrkUbh1I20iCRfw+AmU58x5LWeGfTS\nSXSPoXwng0dUfjpJE88K0NOTQS8t0N3j+ybgjVV1H4OnRl032aqk0TlHL0mNc0QvSY0z6CWpcQa9\nJDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mN+1+wzzpHS4vsSwAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"num_samples = 1000\n",
"len_x = len(x)\n",
"idx = np.random.randint(0, len_x, (num_samples, len_x))\n",
"x_bar = np.mean(x[idx], axis = 1)\n",
"plt.hist(x_bar,bins=20)\n",
"plt.xlabel(r'$\\bar{x}$')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now use bootstrap resampling technique to estimate the posterior of $\\Omega_m$ and $w$.\n",
"
\n",
" 10. Take 200 (or more) samples of 580 supernova distance modulus measurements and estimate $\\Omega_m$ and $w$ using maximum likelihood estimation, as in Part 3. Plot the 1-d posteriors (histogram).
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"data = np.loadtxt(\"sn_z_mu_dmu_plow_union2.1.txt\", usecols=range(1,5))\n",
"# z\n",
"z_data = data[:,0]\n",
"# mu\n",
"mu_data = data[:,1]\n",
"# error on mu (sigma(mu))\n",
"mu_err_data = data[:,2]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"***"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## To Submit\n",
"Execute the following cell to submit.\n",
"If you make changes, execute the cell again to resubmit the final copy of the notebook, they do not get updated automatically.
\n",
"__We recommend that all the above cells should be executed (their output visible) in the notebook at the time of submission.__
\n",
"Only the final submission before the deadline will be graded. \n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"_ = ok.submit()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.13"
}
},
"nbformat": 4,
"nbformat_minor": 2
}