<jupyter_start><jupyter_text>MCMC from scratch Here we will write a simple python program that will perform the Metropolis algorithm. In order to sample the posterior of the probability function given supernovae data.  Loglike computation First we need to be able to compute the likelihood of the data given some parameters. We will actually compute the logarithm of the likelihood as it is easier to handle for the computer.First let's import some basic libraries:<jupyter_code>import numpy as np
import scipy.integrate as integrate<jupyter_output><empty_output><jupyter_text>Now let's use a scipy integrator in order to obtain the luminosity distance<jupyter_code>def E(z,OmDE):
    """
    This function computes the integrand for the computation of the luminosity distance for a flat universe
    z -> float
    OmDE -> float
    gives
    E -> float
    """
    return 1/np.sqrt((1-OmDE)*(1+z)**3+OmDE)

def dl(z,OmDE,h=0.7):
    """
    This function computes the luminosity distance
    z -> float
    OmDE -> float
    h ->float
    returns
    dl -> float
    """
    inte=integrate.quad(E,0,z,args=(OmDE))
    # Velocidad del sonido en km/s
    c = 299792.458
    # Factor de Hubble
    Ho = 100*h
    return c*(1+z)/Ho * inte[0]<jupyter_output><empty_output><jupyter_text>We now have to load the supernovae data in order to compare it with the theoretical data<jupyter_code>zandmu = np.loadtxt('../data/SCPUnion2.1_mu_vs_z.txt', skiprows=5,usecols=(1,2))
covariance = np.loadtxt('../data/SCPUnion2.1_covmat_sys.txt')
type(zandmu)
np.shape(zandmu[1])<jupyter_output><empty_output><jupyter_text>How does this data look like?<jupyter_code>import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sb<jupyter_output><empty_output><jupyter_text>The supernovae<jupyter_code>plt.plot(zandmu[:,0],zandmu[:,1],'o')
heat = np.log(abs(covariance))
sb.set()
sb.heatmap(heat)
yerr = np.loadtxt('../data/SCPUnion2.1_mu_vs_z.txt', skiprows=5,usecols=[3])
plt.errorbar(zandmu[:,0],zandmu[:,1], yerr=yerr, fmt='o')<jupyter_output><empty_output><jupyter_text>Now we can compute the log like.<jupyter_code>dl = np.vectorize(dl)
def loglike(params,h=0.7):
    """
    This function computes the logarithm of the likelihood. It recieves a vector
    params-> vector with one component (Omega Dark Energy)
    """
    OmDE = params[0]
# Ahora quiero calcular la diferencia entre el valor reportado y el calculado
    delta = 5.*np.log10(dl(zandmu[:,0],OmDE,h))+25-zandmu[:,1]
    chisquare=np.dot(delta,np.dot(np.linalg.inv(covariance),delta))
    return -chisquare/2

table_omega = np.arange(0.,1.,0.01)
tableprob=[loglike([hola]) for hola in table_omega]

plt.plot(table_omega,tableprob)
loglike([0.6])<jupyter_output><empty_output><jupyter_text>All for today!  You can try jla data From http://supernovae.in2p3.fr/sdss_snls_jla/ReadMe.html download  jla_likelihood_v6.tgz and covmat_v6.tgz.- The supernovae info is contained in the file jla_lcparams.txt go that file and remove the "" from the first line- Load it using genfromtxt<jupyter_code>jla_dataset = np.genfromtxt('../data/jla/data/jla_lcparams.txt', names=True, dtype=None)<jupyter_output><empty_output><jupyter_text>This function is very clever and you can use the names of the variables instead of the number of the corresponding column, like that:<jupyter_code>jla_dataset['zcmb']<jupyter_output><empty_output><jupyter_text>You can ask which are the different columns or info in that file:<jupyter_code>jla_dataset.dtype.names
plt.plot(jla_dataset['zcmb'],jla_dataset['mb'],'o')<jupyter_output><empty_output><jupyter_text>From http://supernovae.in2p3.fr/sdss_snls_jla/ReadMe.html download  covmat_v6.tgz and use the example.py to compute the covariance matrix for a particular alpha, beta.- Copy the contents of example.py here- Change the names of the files according to your own file order- Change all the "numpy" for "np"<jupyter_code>import pyfits
import glob
def mu_cov(alpha, beta):
    """ Assemble the full covariance matrix of distance modulus

    See Betoule et al. (2014), Eq. 11-13 for reference
    """
    Ceta = sum([pyfits.getdata(mat) for mat in glob.glob('../data/jla/covmat/C*.fits')])

    Cmu = np.zeros_like(Ceta[::3,::3])
    for i, coef1 in enumerate([1., alpha, -beta]):
        for j, coef2 in enumerate([1., alpha, -beta]):
            Cmu += (coef1 * coef2) * Ceta[i::3,j::3]

    # Add diagonal term from Eq. 13
    sigma = np.loadtxt('../data/jla/covmat/sigma_mu.txt')
    sigma_pecvel = (5 * 150 / 3e5) / (np.log(10.) * sigma[:, 2])
    Cmu[np.diag_indices_from(Cmu)] += sigma[:, 0] ** 2 + sigma[:, 1] ** 2 + sigma_pecvel ** 2
    
    return Cmu<jupyter_output><empty_output><jupyter_text>Now you can compute the covariance matrix but it will depend on the nuisanse parammeters alfa and beta<jupyter_code>Cmu = mu_cov(0.13, 3.1)
np.shape(Cmu)<jupyter_output><empty_output><jupyter_text>Now I can compute the loglike in the same way as for the union2.1 data. Just copy the loglike function from before, add as params - OmDE - alpha - beta - M_bFor the covariance matrix use the function from before, and instead of the "zandmu" variable use jla_dataset<jupyter_code>def loglike_jla(params,h=0.7):
    """
    This function computes the logarithm of the likelihood. It recieves a vector
    params-> vector with three components (Omega Dark Energy, alpha, beta and M_b)
    """
    OmDE = params[0]
    alpha = params[1]
    beta = params[2]
    MB = params[3]
    covariance = mu_cov(alpha, beta)
    inv_covariance=np.linalg.inv(covariance)
# Ahora quiero calcular la diferencia entre el valor reportado y el calculado
    mu_obs = jla_dataset['mb']-(MB-alpha*jla_dataset['x1']+beta*jla_dataset['color'])
    mu_teo = 5.*np.log10(dl(jla_dataset['zcmb'],OmDE,h))+25
    delta = mu_teo - mu_obs
    chisquare = np.dot(delta,np.dot(inv_covariance,delta))
    return -chisquare/2
param = [0.65,0.13,3.1,-20]
loglike_jla(param)
param = [0.7,0.13,3.1,-19]
loglike_jla(param)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11974649766511007, 'fraction_numerical': 0.022515010006671114, 'mean_word_length': 4.609915809167446, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 58, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '46517', 'license': 'MIT', 'n_tokens_mistral': 2214, 'n_tokens_neox': 2113, 'n_words': 574, 'path': 'Projects/mcmc/first_day.ipynb', 'repo_name': 'cosmostatschool/MACSS2017', 'size': '5996'}