In [22]:
import numpy as np
import scipy.integrate as integrate
def E(z,OmDE,OmM):
"""
This function computes the integrand for the computation of the luminosity distance for a flat universe
z -> float
OmDE -> float
OmM -> float
gives
E -> float
"""
Omk=1-OmDE-OmM
return 1/np.sqrt(OmM*(1+z)**3+OmDE+Omk*(1+z)**2)
def dl(z,OmDE,OmM,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,OmM))
# Velocidad del sonido en km/s
c = 299792.458
# Factor de Hubble
Ho = 100*h
Omk=1-OmDE-OmM
distance_factor = c*(1+z)/Ho
if Omk>1e-10:
omsqrt = np.sqrt(Omk)
return distance_factor / omsqrt * np.sinh(omsqrt * inte[0])
elif Omk<-1e-10:
omsqrt = np.sqrt(-Omk)
return distance_factor / omsqrt * np.sin(omsqrt * inte[0])
else:
return distance_factor * inte[0]
zandmu = np.loadtxt('../data/SCPUnion2.1_mu_vs_z.txt', skiprows=5,usecols=(1,2))
covariance = np.loadtxt('../data/SCPUnion2.1_covmat_sys.txt')
inv_cov = np.linalg.inv(covariance)
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, Omega Matter)
"""
OmDE = params[0]
OmM = params[1]
# Ahora quiero calcular la diferencia entre el valor reportado y el calculado
muteo = 5.*np.log10(dl(zandmu[:,0],OmDE,OmM,h))+25
print dl(zandmu[:,0],OmDE,OmM,h)
delta = muteo-zandmu[:,1]
chisquare=np.dot(delta,np.dot(inv_cov,delta))
return -chisquare/2
In [23]:
loglike([0.6,0.3])
Out[23]:
In [2]:
def markovchain(steps, step_width, pasoinicial):
chain=[pasoinicial]
likechain=[loglike(chain[0])]
accepted = 0
for i in range(steps):
rand = np.random.normal(0.,1.,len(pasoinicial))
newpoint = chain[i] + step_width*rand
liketry = loglike(newpoint)
if np.isnan(liketry) :
print 'Paso algo raro'
liketry = -1E50
accept_prob = 0
elif liketry > likechain[i]:
accept_prob = 1
else:
accept_prob = np.exp(liketry - likechain[i])
if accept_prob >= np.random.uniform(0.,1.):
chain.append(newpoint)
likechain.append(liketry)
accepted += 1
else:
chain.append(chain[i])
likechain.append(likechain[i])
chain = np.array(chain)
likechain = np.array(likechain)
print "Razon de aceptacion =",float(accepted)/float(steps)
return chain, likechain
In [3]:
chain1, likechain1 = markovchain(100,[0.1,0.1],[0.7,0.3])
In [4]:
import matplotlib.pyplot as plt
%matplotlib inline
In [5]:
plt.plot(chain1[:,0],'o')
plt.plot(chain1[:,1],'o')
Out[5]:
In [6]:
columna1=chain1[:,0]
np.sqrt(np.mean(columna1**2)-np.mean(columna1)**2)
Out[6]:
First just an histogram
In [7]:
import seaborn as sns
In [8]:
omegade=chain1[:,0]
omegadm=chain1[:,1]
In [9]:
sns.distplot(omegade)
Out[9]:
In [10]:
sns.distplot(omegadm)
Out[10]:
In [10]:
sns.jointplot(x=omegadm,y=omegade)
Out[10]:
In [11]:
import corner
In [12]:
corner.corner(chain1)
Out[12]:
How to obtain the confidence regions of a variable?
In [13]:
mean_de=np.mean(omegade)
print(mean_de)
You can use the "central credible interval"
In [14]:
dimde=len(omegade)
points_outside = np.int((1-0.68)*dimde/2)
de_sorted=np.sort(omegade)
print de_sorted[points_outside], de_sorted[dimde-points_outside]
In [28]:
dimde=len(omegade)
points_outside = np.int((1-0.95)*dimde/2)
de_sorted=np.sort(omegade)
print de_sorted[points_outside], de_sorted[dimde-points_outside]