Astronomy 8824 - Problem Set 4 Start

Paul Martini, Autumn 2019

import numpy as np
from numpy import matrix
from numpy import linalg
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import corner

def gaussian(x, mu, sig):
    '''
    calculate a gaussian with mean mu and dispersion sig at input points x
    '''
    return np.exp( -0.5 * np.power( (x - mu)/sig, 2.)) / (np.sqrt(2*np.pi)*sig)

def PlotTwoDist(xy1, xy2, label1, label2, size, addgauss=False, gxsig=False, gysig=False, connect=20): 
    '''
    xy1, xy2: (x,y) points for the two distributions
    label1, label2: labels for the two distributions
    size: (float) half width of plot range
    addgauss: (bool) True to overplot Gaussian on each histogram
    gxsig, gysig: (float) sigma values for the two histograms
    '''
    bins = np.arange(-1.*size, size, 0.25)
    fig = plt.figure(figsize=(8,8))
    gs = GridSpec(4,4)
    ax_scatter = fig.add_subplot(gs[1:4, 0:3])
    ax_xhist = fig.add_subplot(gs[0,0:3])
    ax_yhist = fig.add_subplot(gs[1:4, 3])
    ax_scatter.scatter(xy1[0], xy1[1], color='k', s=1, label=label1)
    ax_scatter.scatter(xy2[0], xy2[1], color='r', s=1, label=label2)
    ax_scatter.plot(xy2[0][:connect], xy2[1][:connect], color='r', ls='-')
    ax_xhist.hist(xy1[0], bins=bins, histtype='step', color='k', density=True)
    ax_xhist.hist(xy2[0], bins=bins, histtype='step', color='r', density=True)
    ax_yhist.hist(xy1[1], bins=bins, histtype='step', color='k', density=True, orientation='horizontal')
    ax_yhist.hist(xy2[1], bins=bins, histtype='step', color='r', density=True, orientation='horizontal')

    if addgauss and gxsig and gysig: 
        gg = np.linspace(-1*size, size, 100)
        ggx = gaussian(gg, 0., gxsig)
        ggy = gaussian(gg, 0., gysig)
        ax_xhist.plot(gg, ggx)
        ax_yhist.plot(ggy, gg)
    
    plt.setp(ax_xhist.get_xticklabels(), visible=False)
    plt.setp(ax_xhist.get_yticklabels(), visible=False)
    plt.setp(ax_yhist.get_xticklabels(), visible=False)
    plt.setp(ax_yhist.get_yticklabels(), visible=False)
    ax_scatter.set_xlabel("X", fontsize=16)
    ax_scatter.set_ylabel("Y", fontsize=16)
    ax_scatter.set_xlim(-8, 8)
    ax_scatter.set_ylim(-8, 8)
    ax_xhist.set_ylabel("N", fontsize=16)
    ax_yhist.set_xlabel("N", fontsize=16)
    ax_scatter.legend(frameon=False, fontsize=16)

1. Bivariate Gaussian

# generate bivariate Gaussian distributions
# following Ivezic et al. sections 3.5.2 and 3.5.4
nelem = 5000
seed = 12         # starting seed for random number generator

# means and dispersions of initial Gaussians
mux = 0
muy = 0
sig1 = 2
sig2 = 0.5

# Generate (p1, p2)
np.random.seed(seed)
p1 = sig1*np.random.normal(size=nelem)
p2 = sig2*np.random.normal(size=nelem)

# Now compute (x, y) from (p1, p2) rotated by angle alpha
alpha = np.pi/6
cosa = np.cos(alpha)
sina = np.sin(alpha)
x = ...
y = ...

# Variances of x and y are
# (See Ivezic et al. eqns 3.85 and 3.86)
sigx2 = ...
sigy2 = ...

# Use python multivariate Gaussian call with expected covariance matrix
# Ivezic et al. eqs. 3.85-3.87
sigxy = ...

mu = ...
cov = ...
xy=np.random.multivariate_normal()

Compare (p1, p2) and the rotated (x, y) distributions

PlotTwoDist([x,y], [p1,p2], "(x,y)", "(p1,p2)", 8, addgauss=True, gxsig=, gysig=)

2. MCMC Realization of a 2-D Probability Distribution

# DHW mc2d.py, modified for notebook 

def MultiGaussProb(x,cinv,prefac):
    '''
    Return multivariate Gaussian probability
    x: vector of data values(matrix)
    cinv: inverse covariance matrix
    prefac:  prefactor for properly normalized Gaussian (float),
             taken as input so that it isn't computed every call
             should be [(2\pi)^{M/2} \sqrt{det(C)}]^{-1}
    '''
    arg = float(x * cinv * x.T)
    return (prefac * np.exp(-arg/2.))

def mcmc(xy_start, xy_step, nchain, nthin, seed, sig1=2, sig2=0.5):
    '''
    xy_start: [xstart, ystart]
    xy_step: [xstep, ystep]
    nchain: output chain length (after thinning)
    nthin: thinning factor
    '''
    np.random.random(seed)
    
    # generate a chain nthin times longer than the final chain
    nchain = nthin * nchain    
    
    mux = 0
    muy = 0
    alpha= np.pi/6
    cosa = np.cos(alpha)
    sina = np.sin(alpha)
    
    '''
    your code here
    '''

1. Start at xy_start = (1,1) and use xy_step = (1,1), compare to multivariate Gaussian

3. Cosmo MCMC: Parameters of the Universe

# Import DHW routines
import cosmodist_subs as cs

h0 = 70. 
om = 0.3
ok = 0. 
w = -1

### DHW cosmomc.py modified for notebook

def prob(model,data,errors):
    """
    Return e^{-chi^2/2}, for diagonal covariance matrix
    model = vector of model predictions
    data = vector of data values
    errors = vector of errors on data
    """
    dev = (model - data)/errors
    return (np.exp(-1.*sum(dev**2)/2.))

def modelvals(om,h,w,ok):
    """
    Return vector of the model predicted data values
    """
    omhsqr=om*(h**2)
    h0=100.*h                                   # value in km/s/Mpc
    [dmzrec,hzrec] = cs.cosmodist(1090, h0, om, ok, w)
    [dmz234,hz234] = cs.cosmodist(2.34, h0, om, ok, w)
    [dmz057,hz057] = cs.cosmodist(0.57, h0, om, ok, w)
    [dmz032,hz032] = cs.cosmodist(0.32, h0, om, ok, w)

    return (np.array([omhsqr,dmzrec,dmz234,hz234*h0,dmz057,hz057*h0,dmz032]))

def cosmomcmc(startvals, stepvals, nchain, nthin, seed):
    
    omstart = startvals['omega_m']
    hstart = startvals['H0']
    okstart = startvals['omega_k']
    wstart = startvals['w']
    
    omstep = stepvals['omega_m']
    hstep = stepvals['H0']
    okstep = stepvals['omega_k']
    wstep = stepvals['w']
    
    np.random.seed(seed)

    nchain = nchain * nthin     # number of elements needed for nchain outputs
    
    # data values and errors being used
    omhsqr=0.1386
    omhsqrerr=0.0027
    dmzrec=13962
    dmzrecerr=10
    dmz234=5381
    dmz234err=170
    hz234=222
    hz234err=5
    dmz057=2204
    dmz057err=31
    hz057=98
    hz057err=3
    dmz032=1249
    dmz032err=25

    data = np.array([omhsqr,dmzrec,dmz234,hz234,dmz057,hz057,dmz032])
    errors = np.array([omhsqrerr,dmzrecerr,dmz234err,hz234err,dmz057err,hz057err,dmz032err])

    # Save all elements of the chain in memory to make the i/o more efficient, 
    # but if memory were an issue we could write out steps one at a time (or in chunks)
    chain = np.zeros((nchain,5))
    chain[0][0] = om1 = omstart
    chain[0][1] = h1 = hstart
    chain[0][2] = w1 = wstart
    chain[0][3] = ok1 = okstart

    '''
    your code here
    '''

Flat Universe

Fix Omega_k = 0 and start with Omega_m = 0.3, h = 0.68, and w = -1. Use a step size Delta = 0.03 in each parameter

startvals = dict()
startvals['omega_m'] = 0.3
startvals['H0'] = 0.68
startvals['omega_k'] = 0.
startvals['w'] = -1.

stepvals = dict()
stepvals['omega_m'] = 0.03
stepvals['H0'] = 0.03
stepvals['omega_k'] = 0.0 # hold fixed
stepvals['w'] = 0.03

# Compute the chain
c_chain1 = cosmomcmc(startvals, stepvals, nchain=4000, nthin=50, seed=12)

# User corner to make a nice plot: 
corner_input1 = c_chain1.T[0:3].T
figure = corner.corner(corner_input1, labels=[r"$\Omega_m$", r"$h$", r"$w$"],
                      quantiles=[0.5], 
                      show_titles=True, title_kwargs={"fontsize": 14},
                      label_kwargs={"fontsize": 14},
                      hist2d_kwargs={"plot_contours": False, "plot_density": False}
                      )