Base Code 10.26 in the book of Bayesian Models for Astrophysical Data Using R, JAGS, Python, and Stan for jla data :
stan_model="
functions {
real[] Ez(real z,
real[] H,
real[] params,
real[] x_r,
int[] x_i) {
real dEdz[1];
dEdz[1] = 1.0/sqrt(params[1]*(1+z)^3
+(1-params[1])*(1+z)^(3*(1+params[2])));
return dEdz;
}
}
data {
int<lower=1> nobs; // number of data points
real E0[1]; // integral(1/H) at z=0
real z0; // initial redshift, 0
real c; // speed of light
real H0; // Hubble parameter
vector[nobs] obs_mag; // observed magnitude at B max
real x1[nobs]; // stretch
real color[nobs]; // color
real redshift[nobs]; // redshift
real hmass[nobs]; // host mass
}
transformed data {
real x_r[0]; // required by ODE (empty)
int x_i[0];
}
parameters{
real<lower=0, upper=1> om; // dark matter energy density
real alpha; // stretch coefficient
real beta; // color coefficient
real Mint; // intrinsic magnitude
real deltaM; // shift due to host galaxy mass
real<lower=0> sigint; // magnitude dispersion
real<lower=-2, upper=0> w; // dark matter equation of state
parameter
}
transformed parameters{
real DC[nobs,1]; // co-moving distance
real pars[2]; // ODE input = (om, w)
vector[nobs] mag; // apparent magnitude
real dl[nobs]; // luminosity distance
real DH; // Hubble distance = c/H0
pars[1] = om;
pars[2] = w;
DH = (c/H0);
# Integral of 1/E(z)
DC = integrate_ode_rk45(Ez, E0, z0, redshift, pars, x_r, x_i);
for (i in 1:nobs) {
dl[i] = DH * (1 + redshift[i]) * DC[i, 1];
if (hmass[i] < 10) mag[i] = 25 + 5 * log10(dl[i])
+ Mint - alpha * x1[i] + beta
* color[i];
else
mag[i] = 25 + 5 * log10(dl[i])
+ Mint + deltaM - alpha * x1[i] + beta * color[i];
}
}
model {
# Priors and likelihood
sigint ~ gamma(0.001, 0.001);
Mint ~ normal(-20, 5.);
beta ~ normal(0, 10);
alpha ~ normal(0, 1);
deltaM ~ normal(0, 1);
obs_mag ~ normal(mag, sigint);
}
"
We can use this code when we have Analytical solutions and put these solutions in Function block
such as:\Omega_{D}(z)= (1-params[1])(1+z)^{(3(1+params[2]))} or params[1](1+z)^{3}. My question is
what if this parameter \Omega_{D} is solution of ode such as
\Omega_{D}'=\frac{d\Omega_{D}}{d z}=-3(\delta-1)\Omega_{D}\frac{1-\Omega_{D}}{1-(2-\delta)\Omega_{D}} that \delta
and initial value for \Omega_{D} are our parameters .
For a better understanding of problem and what do I want Solve problem with python.
def ode_func(x, t,a):
return -(3*pow(1+t,-1)*(a-1))*x*((1-x)/(1-2*x+a*x))
sol=lambda a,inits: odeint(ode_func,inits,t,args=(a,))
def model(a,inits):
y = np.ravel(sol(a,inits))
return interp1d(t, y, kind='cubic')
def comoving_integrand(z, omm,a):
inits=1-omm
f2= model(a,inits)
Ez = 1.0/np.sqrt((omm*(1.0+z)**3 + f2(z)))
return Ez
def luminosityDistance(z, omm,a):
integral = integrate.quad(comoving_integrand, 0, z, args=(omm,a))[0]
return 300000*(1+z)*integral
My goal is achieving luminosityDistance that is similar to DC = integrate_ode_rk45(Ez, E0, z0, redshift, pars, x_r, x_i);
in transformed parameters block of Stan code that I mentioned above .
Any suggestions would be greatly appreciated.