Analysis of  growth rate data

read current data

SetDirectory[NotebookDirectory[]]

We wish to analyse a dataset of growth rate, the so - called fs8.
We have a theory that depends on Omega_m and on the power spectrum normalization s8.
We wish to find the best fit of the parameters and their confidence regions.
To do so, we need to build the likelihood of the data.

(* Data vector: first 9 entries from Macaulay 2013;  11th - Chang & Wang 2013 (wrong!); 12th - Chuang 2013; 13th - Samushia 2013 *)
(* d = growth rate f*s8 data *)
zval={0.067,0.25,0.37,0.30,0.60,0.44,0.60,0.73,0.8,0.13,0.35,0.32,0.57};
nval={-0.064851,-0.223144,-0.314811,-0.262364,-0.470004,-0.364643,-0.470004,-0.548121,-0.587787,-0.122218,-0.300105,-0.277632,-0.451076};
d = {0.423,0.3512,0.4602,0.408,0.433,0.413,0.390,0.437,0.47,0.46,0.445,0.384,0.441}; 
errors={0.055,0.0583,0.0378,0.0552,0.0662,0.080,0.063,0.072,0.08,0.06,0.097,0.095,0.043};
len=Length[d];
data=Table[{nval[[i]],d[[i]]},{i,len}];
cov=DiagonalMatrix[errors^2];
invcov=Inverse[cov];

Here  (Subscript[fσ, 8])(z)=(Subscript[σ, 80])δ’/(Subscript[δ, 0])

data

ListPlot[Table[{zval[[i]],Around[d[[i]],errors[[i]]]},{i,len}],Frame->True,FrameLabel->{"z","growth ((Subscript[fσ, 8]))"}]

Likelihood function with parameters Omega_m and present amplitude sigma80

Values of the grid

xmin=0.02;
xmax=0.8;
xstep=0.01;
omegam=Table[i,{i,xmin,xmax,xstep}];
smin=.5;smax=1.3;sstep=0.01;
sigma=Table[i,{i,smin,smax,sstep}];
Length[omegam]*Length[sigma]
chi2r=Table[0,{iom,Length[omegam]},{sig,Length[sigma]}];
count=0;

perturbation differential equation

Clear[delta];
ni=-1.5;
nf=0;

f1[omegamx_,w0x_,wax_,n_]=Simplify[2-(3/2)((omegamx*Exp[-3n]+(1-omegamx)*(1+w0x+wax-Exp[n]wax)Exp[-3(wax-wax*Exp[n]+n(1+w0x+wax))])/(omegamx*Exp[-3n]+(1-omegamx)*Exp[Integrate[-3(1+w0x+wax(1-Exp[x])),{x,0,n}]]))];
f2[omegamx_,w0x_,wax_,n_]=Simplify[((3 )/2)((Exp[-3n]*omegamx)/(omegamx*Exp[-3n]+(1-omegamx)*Exp[Integrate[-3(1+w0x+wax(1-Exp[x])),{x,0,n}]]))];

delta[omegam_,w0_,wa_]:=NDSolve[{δ''[n]==-f1[omegam,w0,wa,n]*δ'[n]+f2[omegam,w0,wa,n]*δ[n],δ[ni]==Exp[ni],δ'[ni]==Exp[ni]},δ,{n,ni,nf}]

del=delta[.3,-1,0][[1,1,2]]
Plot[del'[x]/del[x],{x,ni,nf},Frame->True]

(* build the log-likelihood:  (data-theory)*InvCov*(data-theory) *)

count=0;
Timing[Do[{
len=Length[d];
del=delta[omegam[[iom]],-1,0][[1,1,2]];
t=Flatten[Table[(D[del[y],y]/del[nf])/.y->nval[[i]],{i,len}]];
chi2r[[iom,isig]]=Sum[(d[[i]]-sigma[[isig]]*t[[i]])*invcov[[i,j]]*(d[[j]]-sigma[[isig]]*t[[j]]),{i,len},{j,len}];
count++;
If[Mod[count,1000]==0 ,Print[count],0]

},{iom,Length[omegam]},{isig,Length[sigma]}]];

chi2=Table[chi2r[[iom,isig]],{iom,Length[omegam]},{isig,Length[sigma]}];

chi2tab=Table[{omegam[[iom]],sigma[[sig]],If[chi2r[[iom,sig]]<100,Exp[-chi2r[[iom,sig]]/2],0]},{iom,Length[omegam]},{sig,Length[sigma]}];
Export["chi-table.dat",chi2tab,"Table"]

lik=Interpolation[Flatten[chi2tab,1]]

Confidence Regions:

(* normalization of the likelihood *)
normint=NIntegrate[lik[x,y],{x,omegam[[1]],omegam[[Length[omegam]]]},{y,sigma[[1]],sigma[[-1]]},WorkingPrecision->10];
likn[x_,y_]=lik[x,y]/normint;
(* find the maximum *)
maxlin=Max[Table[likn[omegam[[i]],sigma[[j]]],{i,Length[omegam]},{j,Length[sigma]}]]

Plot3D[likn[x,y],{x,omegam[[1]],omegam[[-1]]},{y,sigma[[1]],sigma[[-1]]},PlotRange->All,AxesLabel->{"omega","sigma"}]

Cutting the likelihood with horizontal planes and finding the integral below

Plot[{Exp[-1/2*x^2],Exp[-1/2],Exp[-2],Exp[-9/2]},{x,-3.5,3.5},Frame->True]

line=ParallelTable[{NIntegrate[If[likn[x,y]>i,likn[x,y],0],{x,omegam[[1]],omegam[[-1]]},{y,sigma[[1]],sigma[[-1]]},Method->{"GlobalAdaptive","MaxErrorIncreases"->10000,Method->"GaussKronrodRule"},MaxRecursion->10000,WorkingPrecision->14,PrecisionGoal->4,AccuracyGoal->4],i},{i,0,maxlin,maxlin/15}]

j=0;eps=0.01;Do[If[line[[i,1]]>eps,j++,0],{i,Length[line]}];line2=Table[line[[i]],{i,j}];
regions=Interpolation[line2,InterpolationOrder->1]

ContourPlot[{likn[x,y]==regions[0.997],likn[x,y]==regions[0.95],likn[x,y]==regions[0.68]},{x,omegam[[1]],omegam[[-1]]},{y,sigma[[1]],sigma[[-1]]},FrameLabel->{"(Subscript[Ω, m])","(Subscript[σ, 8])"},Frame->{True},LabelStyle->{Large},LabelStyle->Thick,FrameStyle->Directive[Thick],PlotPoints->200]