This is the notebook that computes the likelyhood functions for the DE/DM interaction parameter zeta.

Any messages printed when you run the notebook should be irrelevant (warnings about names and such crap).

					Last modified Aug. 23, 14:10

This is the data.
I have written it as a vector (lazy way...)
Below I will separate this vector into 34 rows and 9 columns.

data={661. ,774. ,549. ,0. ,0. ,0. ,693, 944, 442,915., 1000., 831., 1010., 1023., 992., 1097, 1173, 1034,861., 1010. ,713. ,913., 927., 899., 695 ,895, 587,922., 1107., 737., 1199., 1263., 1156., 0 ,0 ,0,875., 1008., 742., 944., 968., 921., 770, 857, 683,210., 706., 69., 988., 1023., 955., 695 ,895 ,587,1037., 1101. ,972., 1115., 1142., 1090., 1164, 1262, 1086,1037., 1145., 930., 1197., 1241., 1167., 1201, 1401, 1001,1027., 1135., 923., 1000., 1018., 982., 0, 0 ,0,929., 1045., 812., 1023., 1221., 857., 779, 846, 730,674., 833., 516., 0., 0., 0., 846, 946, 746,999., 1100., 900., 0., 0., 0., 729, 871, 587,980., 1105., 855., 1152., 1375., 964., 1102, 1239, 995,775., 876., 674., 831., 868., 762., 597, 714 ,507,705., 809., 601., 0., 0., 0., 625, 752, 530,1243, 1301, 1185, 1174, 1304, 1044, 1161, 1357, 965,1117., 1193., 1035., 1201., 1240., 1161., 1262, 1351, 1194,1164., 1315., 991., 1168., 1221., 1120., 1127, 1293, 1015,880., 979., 769., 978., 1010., 945., 886, 964, 818,837., 970., 679., 868., 970., 772., 831, 960, 774,566., 711., 371., 0., 0., 0., 686, 751, 636,1048., 1150., 935., 1025., 1070., 986., 1003, 1064, 951,964., 1051., 869., 840., 859., 821., 1032, 1162, 937,722., 867., 541., 751., 840., 674., 575, 713, 485,998., 1126., 852., 1010., 1106., 929., 839, 906, 786,1036., 1125., 939., 1120., 1175., 1077., 0, 0, 0,1008., 1098., 909., 962., 1002., 929., 0, 0, 0,844., 943., 732., 1017., 1048., 986., 0, 0, 0,1370., 1435., 1302., 1201., 1227., 1181., 1172, 1295, 1073,1060., 1147., 965., 1099., 1134., 1070., 0, 0, 0,1042., 1129., 948., 1048., 1077., 1025., 1222, 1369, 1113,1074., 1156., 985., 1240., 1283., 1201., 0, 0, 0,1359., 1449., 1263., 1063., 1120., 1002., 859, 977, 747,1140., 1251., 1017., 903., 1063., 782., 911, 992, 804};

To this data I will add the clusters with X-ray and Virial data (but no weak lensing data) which Laerte sent me on Aug. 24. From his data I have chosen:

Name           Ref. RX(h-1 Mpc) MX(1014 h-1 Msun)     Mopt(1014 h-1 Msun
A119...		6	0.75	2.65 +0.27 -0.27	1.63 +0.58 -0.48
A194...		5	0.22	0.25 +0.13 -0.13	0.25 +0.09 -0.06
A262...		1	1.50	1.20 +0.24 -0.23	1.78 +0.41 -0.34
A399...		5	0.73	4.28 +0.59 -0.51	6.89 +1.40 -1.35
A400...		5	0.34	0.63 +0.32 -0.15	1.14 +0.33 -0.28
A401...		7	0.80   15.64 +4.43 -4.36		8.25 +1.67 -1.51
A426...		8	0.65	2.50 +1.00 -0.75	3.98 +0.93 -0.66
A458...		5	0.37	1.13 +0.56 -0.56	1.69 +0.46 -0.35
A496...		5	0.61	1.88 +0.32 -0.28	2.11 +0.69 -0.63
A539...		1	1.50	1.00 +0.28 -0.31	2.21 +0.57 -0.46
A576...		5	0.30	1.07 +0.20 -0.14	1.45 +0.24 -0.22
A754...		5	0.61	5.00 +0.86 -0.75	2.61 +0.64 -0.45
A999...		5	0.13	0.19 +0.10 -0.10	0.09 +0.07 -0.04
A1060...	12	0.25	0.50 +0.05 -0.10	0.70 +0.14 -0.12
A1142...	5	0.21	0.40 +0.18 -0.18	0.38 +0.13 -0.07
A1146...	5	0.35	1.95 +0.97 -0.97	2.58 +0.73 -0.65
A1185...	5	0.22	0.93 +0.47 -0.26	0.50 +0.22 -0.14
A1314...	5	0.20	0.36 +0.32 -0.13	0.13 +0.09 -0.06
A1631...	5	0.22	0.63 +0.32 -0.32	0.87 +0.15 -0.14
A1644...	5	0.28	1.25 +0.14 -0.14	0.41 +0.07 -0.07
A1656...	14	2.50	9.25 +1.20 -1.20	6.43 +0.88 -0.74
A1795...	1	1.50	6.59 +0.58 -0.61	5.52 +1.25 -1.15
A1809...	5	0.58	2.02 +0.55 -0.55	2.80 +0.65 -0.57
A1983...	5	0.17	0.41 +0.20 -0.20	0.38 +0.08 -0.07
A1991...	5	0.33	1.90 +2.08 -0.77	1.02 +0.57 -0.55
A2029...	1	1.50	8.59 +1.69 -1.82	6.21 +2.19 -2.10
A2040...	5	0.31	0.54 +0.27 -0.27	0.65 +0.43 -0.33
A2063...	1	1.50	4.54 +0.89 -0.86	3.16 +0.61 -0.50
A2079...	5	0.37	0.87 +0.43 -0.43	1.49 +0.56 -0.39
A2092...	5	0.31	0.64 +0.32 -0.32	0.74 +0.37 -0.23
A2107...	5	0.33	1.43 +0.64 -0.38	1.08 +0.33 -0.32
A2124...	5	0.25	0.76 +0.38 -0.38	1.09 +0.29 -0.26
A2142...	5	1.08	10.70 +1.65 -1.46	11.69 +2.61 -2.29
A2151...	17	0.34	0.23 +0.00 -0.00	1.13 +0.29 -0.23
A2197...	5	0.27	0.39 +0.20 -0.20	0.98 +0.22 -0.22
A2199...	6	0.75	2.31 +0.10 -0.05	3.58 +0.98 -0.76
A2256...	18	1.40	9.50 +0.30 -0.30	16.77 +2.87 -2.49
A2319...	5	0.80	7.85 +0.63 -0.56	15.67 +2.32 -2.03
A2554...	5	0.48	2.10 +1.05 -1.05	2.98 +1.08 -0.73
A2589...	1	1.50	2.64 +0.52 -0.50	0.81 +0.56 -0.47
A2593...	5	0.55	1.63 +0.79 -0.47	1.31 +0.46 -0.30
A2634...	6	0.75	1.62 +0.09 -0.14	2.60 +0.82 -0.59
A2670...	6	0.75	2.33 +0.95 -0.54	3.43 +0.45 -0.36
A2717...	22	0.65	0.55 +0.10 -0.05	1.00 +0.26 -0.19
A2877...	5	0.23	0.86 +0.27 -0.19	1.53 +0.50 -0.45
A3158...	10	0.52	0.60 +0.30 -0.30	4.06 +0.82 -0.75
A3266...	5	0.73	4.26 +0.34 -0.28	3.64 +0.61 -0.51
A3376...	5	0.24	0.65 +0.32 -0.32	0.87 +0.19 -0.16
A3391...	5	0.41	2.00 +0.50 -0.34	1.42 +0.87 -0.53
A3528N...	23	0.50	0.16 +0.14 -0.05	0.72 +0.28 -0.16
A3532...	5	0.45	1.67 +0.57 -0.57	1.08 +0.35 -0.28
A3558...	24	0.75	1.50 +0.28 -0.28	5.49 +0.49 -0.43
A3667...	5	0.88	5.65 +0.87 -0.87	5.41 +0.74 -0.58
A3744...	5	0.23	0.81 +0.41 -0.41	0.36 +0.11 -0.08
S805...		5	0.13	0.15 +0.03 -0.03	0.11 +0.03 -0.02
S373...		26	0.22	0.15 +0.07 -0.07	0.12 +0.03 -0.03
AWM 4...	5	0.55	2.18 +1.18 -0.59	0.04 +0.07 -0.05
AWM 7...	5	0.31	1.26 +0.07 -0.07	2.07 +0.65 -0.54
MKW 3S...	5	0.38	1.24 +0.12 -0.12	0.81 +0.27 -0.24
MKW 4...	27	0.25	0.35 +0.18 -0.18	0.56 +0.17 -0.13
Virgo...		28	0.15	0.19 +0.66 -0.12	0.10 +0.01 -0.01

a) In the above, masses are estimated using:

	M(<r) = 3 beta sig_p^2 r /G,	with beta=0.8

Notice that we only employ fractions of the sigmas in our analysis, so this is equivalent to just using the masses above (which are proportional to the same factors times the velocity dispersions sigma).

Directory[]

dat2=Import["/Users/lrwa/trabalho/Virial_Dark_Interactions/girardi.dat"];
TableForm[dat2]

l2=Length[dat2]

dat2[[61]][[8]]

Length[data]
l1=Length[data]/9

Now I generate the vectors for each individual set of data (corresponding to each row of the data that Laerte gave to us: 

	weak lensing (central value) - wl0
	weak lensing (superior value) - wls
	weak lensing (inferior value) - wli

	X-ray (central value) - rx0
	X-ray (superior value) - rxs
	X-ray (inferior value) - rxi

	virial (central value) - vi0
	virial (superior value) - vis
	virial (inferior value) - vii

wl0=Table[0,{Length[data]/9}];
wls=Table[0,{Length[data]/9}];
wli=Table[0,{Length[data]/9}];

rx0=Table[0,{Length[data]/9}];
rxs=Table[0,{Length[data]/9}];
rxi=Table[0,{Length[data]/9}];

vi0=Table[0,{Length[data]/9}];
vis=Table[0,{Length[data]/9}];
vii=Table[0,{Length[data]/9}];

Now I fill the vectors above with their actual values:

For[i=1,i<(l1+1),
wl0[[i]]=data[[9*(i-1)+1]];
wls[[i]]=data[[9*(i-1)+2]];
wli[[i]]=data[[9*(i-1)+3]];
rx0[[i]]=data[[9*(i-1)+4]];
rxs[[i]]=data[[9*(i-1)+5]];
rxi[[i]]=data[[9*(i-1)+6]];
vi0[[i]]=data[[9*(i-1)+7]];
vis[[i]]=data[[9*(i-1)+8]];
vii[[i]]=data[[9*(i-1)+9]];
i++]

For[j=1,j<(l2+1),
wl0=Append[wl0,0.];
wli=Append[wli,0.];
wls=Append[wls,0.];
rx0=Append[rx0,dat2[[j]][[3]]];
rxi=Append[rxi,dat2[[j]][[3]]+dat2[[j]][[4]]];
rxs=Append[rxs,dat2[[j]][[3]]+dat2[[j]][[5]]];
vi0=Append[vi0,dat2[[j]][[6]]];
vii=Append[vii,dat2[[j]][[6]]+dat2[[j]][[7]]];
vis=Append[vis,dat2[[j]][[6]]+dat2[[j]][[8]]];
j++]

Length[wl0]
Length[vis]

These are the central values of the X-ray data, for instance:

rx0

Now I will order the data using the X-ray data: first come the data for clusters with lowest X-ray velocity dispersions; last comes the data for clusters with highest X-ray velocity dispersions. This would enable us to see correlations between size (which is roughly proportional to the velocity dispersions) and our parameter (if there are any such correlations).

These are the permutations that would put the X-ray data in order:

ordrx0=Ordering[rx0]

Now I apply this ordering to all the data:

rx0

rx0r=rx0[[ordrx0]]
rxsr=rxs[[ordrx0]]
rxir=rxi[[ordrx0]]

wl0r=wl0[[ordrx0]]
wlsr=wls[[ordrx0]]
wlir=wli[[ordrx0]]

vi0r=vi0[[ordrx0]]
visr=vis[[ordrx0]]
viir=vii[[ordrx0]]

I will use multiple list-plots, so I need this...

Here are the (ordered) data for weak lensing velocity dispersions. 
Boxes are central values, stars are upper and lower limits (at 66% C.L.)

ListPlot[{wl0,wls,wli},Frame->True,Filling->{1->{2},1->{3}},FillingStyle->{Red,Blue}]

Here is the data for the X-ray velocities:

ListPlot[{rx0r,rxsr,rxir},Frame->True,Filling->{1->{2},1->{3}},FillingStyle->{Red,Blue},PlotRange->{{0,l2},{-0.1,10}}]

ListPlot[{rx0r,rxsr,rxir},Frame->True,Filling->{1->{2},1->{3}},FillingStyle->{Red,Blue},PlotRange->{{l2,l2+l1},{0,1500}}]

And here is the data for the virial velocities:

ListPlot[{vi0r,visr,viir},Frame->True,Filling->{1->{2},1->{3}},FillingStyle->{Red,Blue},PlotRange->{{0,66},{-0.1,5}}]

Now, let me define the fractional value of these velocity dispersions.
To estimate error bars, I will use here the upper value/smaller value and small value/upper value.
Notice that this is just an estimate, the error bars are very roughly estimated and have no real statistical value.

First, the virial/weak lensing fractions:

fracviwl=1.0*vi0r/wl0r;
fracviwls=1.0*visr/wlir;
fracviwli=1.0*viir/wlsr;

In our model, there should be a bias between the virial/weak lensing AS WELL AS between the virial/X-ray velocity dispersions, but there should be NO bias between the weak lensing/X-ray velocity dispersions.
(You can ignore the error messages: they just correspond to the missing data points, which are given the value 0 and which give an error when they appear in the denominator of a fraction.)

ListPlot[{fracviwl-1,fracviwls-1,fracviwli-1},Frame->True,Filling->{1->{2},1->{3}},FillingStyle->{Red,Blue}]

fracvirx=1.0*vi0r/rx0r;
fracvirxs=1.0*visr/rxir;
fracvirxi=1.0*viir/rxsr;

ListPlot[{fracvirx-1,fracvirxs-1,fracvirxi-1},Frame->True,Filling->{1->{2},1->{3}},FillingStyle->{Red,Blue}]

fracwlrx=1.0*wl0r/rx0r;
fracwlrxs=1.0*wlsr/rxir;
fracwlrxi=1.0*wlir/rxsr;

ListPlot[{fracwlrx-1,fracwlrxs-1,fracwlrxi-1},Frame->True,Filling->{1->{2},1->{3}},FillingStyle->{Red,Blue}]

Now I will begin to construct the Likelihood Function. 

The idea is that, having that likelihood, we can compute the probability distribution for the coupling ζ.
Usually, the Likelihood function is approximately Gaussian. This pressupposes that the errors are symmetric.

However, our data has asymmetric error bars, which means that our likelihood cannot be Gaussian. 
In fact, since our data is for some (σ^2)>0, the Likelihood function should be better approximated by a Poisson distribution, or a Maxwell-Boltzmann.

I will use the method described in arXiv: physics/0406120 (R. Barlow), in Sec. 3.6. This method uses the asymmetric errors to reconstruct a Likelihood function, and it does an optimal job for a Poisson-like distribution.

The method uses (SubPlus[σ]) and (SubMinus[σ])(the errors above and below the mean). The Likelihood is then given by:

	ln L((OverHat[x]);x)=-(1/2)(((((OverHat[x])-x))^2)/(U+V(x-(OverHat[x]))))  ,

where:

	U = (SubPlus[σ])(SubMinus[σ])

	V=(SubPlus[σ]) - (SubMinus[σ])

First, I will construct the data vectors for the (SubPlus[σ]) and (SubMinus[σ]) (they should all be positive, of course):

sigwlplus=wlsr-wl0r;
sigwlminus=wl0r-wlir;

sigrxplus=rxsr-rx0r;
sigrxminus=rx0r-rxir;

sigviplus=visr-vi0r;
sigviminus=vi0r-viir;

Now, compute the U and V for each dataset:

Uwl=sigwlplus*sigwlminus;
Vwl=sigwlplus-sigwlminus;

Urx=sigrxplus*sigrxminus;
Vrx=sigrxplus-sigrxminus;

Uvi=sigviplus*sigviminus;
Vvi=sigviplus-sigviminus;

Now we can construct the Likelihood for the unknown extra parameter ζ, given the data:

	L( (Subscript[(((Subscript[σ, a])^2)/((Subscript[σ, b])^2)), Theory]) | data ) = (∏_(data a/b)) Exp [ - (1/2)(((( (Subscript[(((Subscript[σ, a])^2)/((Subscript[σ, b])^2)), Theory])-(Subscript[(((Subscript[σ, a])^2)/((Subscript[σ, b])^2)), Obs])))^2)/(U+V((Subscript[(((Subscript[σ, a])^2)/((Subscript[σ, b])^2)), Obs])-(Subscript[(((Subscript[σ, a])^2)/((Subscript[σ, b])^2)), Theory]))))]

where the only way in which our extra parameter comes into the Theory is through the virial mass:

	(Subscript[(((Subscript[σ, vi])^2)), Corr])  =  (((1 - ζ/2)/(1-2ζ))) (Subscript[(((Subscript[σ, vi])^2)), 0])

This corrected virial mass should be equal to the other masses. 

Hence, our theory says that, if we use the uncorrected virial mass (Subscript[(((Subscript[σ, vi])^2)), 0]) to compute the ((Subscript[σ, vi])^2), then we should have:

	f1 := ((Subscript[σ, vi])^2)/((Subscript[σ, rx])^2)= (((1 - 2ζ)/(1-ζ/2)))
	
and

	f2 := ((Subscript[σ, vi])^2)/((Subscript[σ, wl])^2)= (((1 - 2ζ)/(1-ζ/2))).
	
As a consistency condition, we may also assume that the fraction:

	f3 := ((Subscript[σ, rx])^2)/((Subscript[σ, wl])^2)= 1 .


In order to construct the Likelihood function for the observables f1 and f2 we need to find the U and V for each one of these datasets.

I will use the fact that, for f=x/y, the asymmetric errors can be approximately given by:

	σf+ = f (Sqrt[((((σx+))^2)/x)+((((σy-))^2)/y)])
and
	σf- = f (Sqrt[((((σx-))^2)/x)+((((σy+))^2)/y)])

These are our basic datasets:

vi0r
rx0r
wl0r

Because every time there is a zero value we will get nonsense in our fractions f1, f2 and g3, I will now prune these datasets.

For f1, we eliminate all data points for which virial OR X-ray have a zero.
For f2, we eliminate all data points for which virial OR weak lensing have a zero.
For f3, we eliminate all data points for X-ray OR weak lensing have a zero.

The new datasets will be denoted by suffixes 1, 2 and 3.

First, I will construct the dataset #1.
The first step is find the zeros in the X-ray data, and delete those clusters both from X-ray and from Virial data:

list1a=Position[rx0r,_?(-0.1<#<0.1&)];

rx0r1a=Delete[rx0r,list1a];
rxsr1a=Delete[rxsr,list1a];
rxir1a=Delete[rxir,list1a];

vi0r1a=Delete[vi0r,list1a];
visr1a=Delete[visr,list1a];
viir1a=Delete[viir,list1a];

Next step is to delete the zeros that still remain in the Virial data from the X-ray and Virial data:

list1b=Position[vi0r1a,_?(-0.1<#<0.1&)];

rx0r1=Delete[rx0r1a,list1b]
rxsr1=Delete[rxsr1a,list1b]
rxir1=Delete[rxir1a,list1b]

vi0r1=Delete[vi0r1a,list1b]
visr1=Delete[visr1a,list1b]
viir1=Delete[viir1a,list1b]

As we can see, for the dataset #1 we remain with the following number of useful clusters:

Length[vi0r1]

Now, I will do the same for the dataset #2, which combines Weak Lensing with Virial data:

list2a=Position[wl0r,_?(-0.1<#<0.1&)];

wl0r2a=Delete[wl0r,list2a];
wlsr2a=Delete[wlsr,list2a];
wlir2a=Delete[wlir,list2a];

vi0r2a=Delete[vi0r,list2a];
visr2a=Delete[visr,list2a];
viir2a=Delete[viir,list2a];

list2b=Position[vi0r2a,_?(-0.1<#<0.1&)];

wl0r2=Delete[wl0r2a,list2b]
wlsr2=Delete[wlsr2a,list2b]
wlir2=Delete[wlir2a,list2b]

vi0r2=Delete[vi0r2a,list2b]
visr2=Delete[visr2a,list2b]
viir2=Delete[viir2a,list2b]

I will also construct a dataset #3, which will serve as a control. It combines Weak Lensing with X-ray data:

list3a=Position[wl0r,_?(-0.1<#<0.1&)];

wl0r3a=Delete[wl0r,list3a];
wlsr3a=Delete[wlsr,list3a];
wlir3a=Delete[wlir,list3a];

rx0r3a=Delete[rx0r,list3a];
rxsr3a=Delete[rxsr,list3a];
rxir3a=Delete[rxir,list3a];

list3b=Position[rx0r3a,_?(-0.1<#<0.1&)];

wl0r3=Delete[wl0r3a,list3b];
wlsr3=Delete[wlsr3a,list3b];
wlir3=Delete[wlir3a,list3b];

rx0r3=Delete[rx0r3a,list3b];
rxsr3=Delete[rxsr3a,list3b];
rxir3=Delete[rxir3a,list3b];

The final result is that we have the following number of useful clusters for the test f1:

Length[vi0r1]

and the following number of useful clusters for the test f2:

Length[vi0r2]

Now, compute the U and V for each dataset:

f1=1.0*vi0r1/rx0r1;
f2=1.0*vi0r2/wl0r2;
f3=1.0*rx0r3/wl0r3;

sigf1plus=f1*Sqrt[((((visr1-vi0r1)/vi0r1))^2)+((((rxir1-rx0r1)/rx0r1))^2)];
sigf1minus=f1*Sqrt[((((viir1-vi0r1)/vi0r1))^2)+((((rxsr1-rx0r1)/rx0r1))^2)];

sigf2plus=f2*Sqrt[((((visr2-vi0r2)/vi0r2))^2)+((((wlir2-wl0r2)/wl0r2))^2)];
sigf2minus=f2*Sqrt[((((viir2-vi0r2)/vi0r2))^2)+((((wlsr2-wl0r2)/wl0r2))^2)];

sigf3plus=f3*Sqrt[((((rxsr3-rx0r3)/rx0r3))^2)+((((wlir3-wl0r3)/wl0r3))^2)];
sigf3minus=f3*Sqrt[((((rxir3-rx0r3)/rx0r3))^2)+((((wlsr3-wl0r3)/wl0r3))^2)];

U1=sigf1plus*sigf1minus
V1=sigf1plus-sigf1minus

U2=sigf2plus*sigf2minus
V2=sigf2plus-sigf2minus

U3=sigf3plus*sigf3minus
V3=sigf3plus-sigf3minus

We need in fact values for U and V which are small, in order for our fit to the likelihood function to work. However, as one can see from the data above:

ListPlot[{U1,V1,U2,V2,U3,V3},Frame->True,PlotRange->All]

U1[[3]]
V1[[3]]
U2[[3]]
V2[[3]]
U3[[3]]
V3[[3]]

U1[[4]]
V1[[4]]
U2[[4]]
V2[[4]]
U3[[4]]
V3[[4]]

U1[[15]]
V1[[15]]
U2[[15]]
V2[[15]]
U3[[15]]
V3[[15]]

U1[[14]]
V1[[14]]
U2[[14]]
V2[[14]]
U3[[14]]
V3[[14]]

U1[[10]]
V1[[10]]
U2[[10]]
V2[[10]]
U3[[10]]
V3[[10]]

This is being caused mainly by the cluster A1651, for which the weak lensing data is extremely noisy, and by a couple of other clusters which have very noisy data. Below I will drop the data for the extremely problematic cases

f1=Delete[f1,3];
U1=Delete[U1,3];
V1=Delete[V1,3];
sigf1plus=Delete[sigf1plus,3];
sigf1minus=Delete[sigf1minus,3];

f1=Delete[f1,4];
U1=Delete[U1,4];
V1=Delete[V1,4];
sigf1plus=Delete[sigf1plus,4];
sigf1minus=Delete[sigf1minus,4];

f1=Delete[f1,15];
U1=Delete[U1,15];
V1=Delete[V1,15];
sigf1plus=Delete[sigf1plus,15];
sigf1minus=Delete[sigf1minus,15];

f2=Delete[f2,14];
U2=Delete[U2,14];
V2=Delete[V2,14];
sigf2plus=Delete[sigf2plus,14];
sigf2minus=Delete[sigf2minus,14];

f3=Delete[f3,10];
U3=Delete[U3,10];
V3=Delete[V3,10];
sigf3plus=Delete[sigf3plus,10];
sigf3minus=Delete[sigf3minus,10];

The total number of clusters being used in each case is now:

n1=Length[f1]
n2=Length[f2]
n3=Length[f3]

ListPlot[{f1,f1+sigf1plus,f1-sigf1minus},Frame->True,Filling->{1->{2},1->{3}},FillingStyle->{Red,Blue},PlotRange->All]

ListPlot[{f2,f2+sigf2plus,f2-sigf2minus},Frame->True,Filling->{1->{2},1->{3}},FillingStyle->{Red,Blue}]

ListPlot[{f3,f3+sigf3plus,f3-sigf3minus},Frame->True,Filling->{1->{2},1->{3}},FillingStyle->{Red,Blue}]

Now we can finally compute the likelihood functions, remembering that f1 and f2 are equal (in theory) to (1 - 2 zeta)/(1 - zeta/2).

I will also compute our control, the fraction f3.

fz[zeta_]:=(1-2 zeta)/(1-zeta/2);

n1=Length[f1];
n2=Length[f2];
n3=Length[f3];

lf1[zeta_,fudge_]:=(∏_(i=1)^n1)Exp[-0.5*(((f1[[i]]-fz[zeta]))^2)/(U1[[i]]+fudge*V1[[i]]*(fz[zeta]-f1[[i]]))];

lf2[zeta_,fudge_]:=(∏_(i=1)^n2)Exp[-0.5*(((f2[[i]]-fz[zeta]))^2)/(U2[[i]]+fudge*V2[[i]]*(fz[zeta]-f2[[i]]))];

lf3[zc_,fudge_]:=(∏_(i=1)^n3)Exp[-0.5*(((f3[[i]]-fz[zc]))^2)/(U3[[i]]+fudge*V3[[i]]*(fz[zc]-f3[[i]]))];

Notice that I have used a fudge factor (which is, in fact, a nuisance parameter) to tune the extent to which we take into account the asymmetry in the errors. A fudge=1 means we take full account of the asymmetries; a fudge=0 means we discard the asymmetries and assume that the likelihood is symmetric around the peak.

Below I plot the likelihood functions for zeta, for: 

	f1 test (virial/X-ray - black)
	f2 test (virial/weak lensing - green)
	f3 control (X-ray/weak lensing - red)

Plot[{Log[lf1[zeta,1]],Log[lf2[zeta,1]],Log[lf3[zeta,1]]},{zeta,-0.2,0.12},Frame->True,PlotStyle->{RGBColor[0,0,0],RGBColor[0,1,0],RGBColor[1,0,0]}]

Plot[{Log[lf1[zeta,0]],Log[lf2[zeta,0]],Log[lf3[zeta,0]]},{zeta,-0.3,0.2},Frame->True,PlotStyle->{RGBColor[0,0,0],RGBColor[0,1,0],RGBColor[1,0,0]}]

Both f1 and f2 seem to favor a slightly positive zeta, and f2 still is consistent with zero (as it should).

The delta-chi squared per degree of freedom is the best estimator of how good a job our extra parameter makes in fitting the data. For the test f2, this criterium gives:

dchi2f1a=2*(Log[lf1[0.05,1]]-Log[lf1[0.0,1]])/(n1+1)
dchi2f1b=2*(Log[lf1[0.09,0]]-Log[lf1[0.0,0]])/(n1+1)

dchi2f2a=2*(Log[lf2[0.023,1]]-Log[lf2[0.0,1]])/(n2+1)
dchi2f2b=2*(Log[lf2[0.0435,0]]-Log[lf2[0.0,0]])/(n2+1)

We can put a prior on the value of zeta, and then the likelihoods can be interpreted as probability distributions for zeta. 

I will assume that zeta can vary in the interval [-1,1].

First we will normalize the p.d.f.s on some interval roughly proportional to that:

Plot[{Log[lf1[zeta,1]]},{zeta,-0.2,0.15},Frame->True]

cn1a=NIntegrate[lf1[zeta,1],{zeta,-0.1,0.1}];
pf1a[zeta_]:=(1/cn1a)*lf1[zeta,1];
Plot[pf1a[zeta],{zeta,-0.1,0.1},Frame->True,PlotRange->All]

Plot[{Log[lf1[zeta,0]]},{zeta,-0.2,0.2},Frame->True]

cn1b=NIntegrate[lf1[zeta,0],{zeta,-0.2,0.2}];
pf1b[zeta_]:=(1/cn1b)*lf1[zeta,0];
Plot[pf1b[zeta],{zeta,-0.05,0.15},Frame->True,PlotRange->All]

Plot[{Log[lf2[zeta,1]]},{zeta,-0.2,0.2},Frame->True]

cn2a=NIntegrate[lf2[zeta,1],{zeta,-0.2,0.2}];
pf2a[zeta_]:=(1/cn2a)*lf2[zeta,1];
Plot[pf2a[zeta],{zeta,-0.05,0.15},Frame->True,PlotRange->All]

Plot[{Log[lf2[zeta,0]]},{zeta,-0.2,0.2},Frame->True]

cn2b=NIntegrate[lf2[zeta,0],{zeta,-0.2,0.2}];
pf2b[zeta_]:=(1/cn2b)*lf2[zeta,0];
Plot[pf2b[zeta],{zeta,-0.05,0.15},Frame->True,PlotRange->All]

Plot[{Log[lf3[zeta,1]]},{zeta,-0.2,0.2},Frame->True]

cn3a=NIntegrate[lf3[zeta,1],{zeta,-0.2,0.2}];
pf3a[zeta_]:=(1/cn3a)*lf3[zeta,1];
Plot[pf3a[zeta],{zeta,-0.1,0.1},Frame->True,PlotRange->All]

Plot[{Log[lf3[zeta,0]]},{zeta,-0.2,0.2},Frame->True]

cn3b=NIntegrate[lf3[zeta,0],{zeta,-0.2,0.2}];
pf3b[zeta_]:=(1/cn3b)*lf3[zeta,0];
Plot[pf3b[zeta],{zeta,-0.1,0.1},Frame->True,PlotRange->All]

Now here are the main results...

Plot[{pf1a[zeta],pf2a[zeta],pf3a[zeta]},{zeta,-0.05,0.08},Frame->True,Filling->Axis,PlotRange->All]

Plot[{pf1b[zeta],pf2b[zeta],pf3b[zeta]},{zeta,-0.05,0.15},Frame->True,Filling->Axis,PlotRange->All]

Plot[{pf1a[zeta],pf1b[zeta]},{zeta,-0.01,0.12},Frame->True,Filling->Axis,PlotRange->All]

Plot[{pf2a[zeta],pf2b[zeta]},{zeta,-0.05,0.12},Frame->True,Filling->Axis,PlotRange->All]

Plot[{pf3a[zeta],pf3b[zeta]},{zeta,-0.05,0.05},Frame->True,Filling->Axis,PlotRange->All]

The analysis of clusters was done above.

Below I will now consider the constraints imposed by the CMB.
For that, I will use the SHIFT PARAMETER (see Bond, Efsthatiou & Tegmark, 1997):

	R=(Sqrt[(Subscript[Ω, m0])])(∫_0^zrec)(1/(E(z')))dz'

where

	(E^2)( z, (Subscript[Ω, m0]), ω , ζ) = (Subscript[Ω, m0])(1+z)^( 3(1+ω))+((Subscript[Ω, m0])/(ζ+3ω))(((ζ(1+z))^(3(1+ω)))+3((ω(1+z))^(3-ζ)))
	
and I have assumed flat spatial section.

The shift parameter was measured by WMAP, and it is currently set at:

	R=1.70 +- 0.03		with	zrec=1089.

Eshift[z_,omm_,w_,zeta_]:=(1-omm)*(1+z)^(3*(1+w))+(omm/(zeta+3*w))(zeta*(1+z)^(3*(1+w))+3*w*(1+z)^(3-zeta));

Rshift[omm_,w_,zeta_]:=(Sqrt[omm])*NIntegrate[(1/(Sqrt[Eshift[x,omm,w,zeta]])),{x,0.,1089.}]

Rshift[0.3,-0.7,0]

Plot3D[Exp[-0.5*(Rshift[omm,w,0]-1.7)^2/(0.03)^2],{omm,0.1,0.4},{w,-0.5,-1.5}]

DensityPlot[Exp[-0.5*(Rshift[omm,w,0]-1.7)^2/(0.03)^2],{omm,0.1,0.4},{w,-0.5,-1.5}]

Plot3D[Exp[-0.5*(Rshift[omm,-1,zeta]-1.7)^2/(0.03)^2],{omm,0.1,0.4},{zeta,-0.15,0.15}]

DensityPlot[Exp[-0.5*(Rshift[omm,-1,zeta]-1.7)^2/(0.03)^2],{omm,0.1,0.4},{zeta,-0.15,0.15}]

Now let's marginalize over the equation of state w (which means, integrate the likelihood over w).

likcmb[omm_,w_,zeta_]:=Exp[-0.5*(Rshift[omm,w,zeta]-1.7)^2/(0.03)^2];

likcmbmarg[omm_,zeta_]:=NIntegrate[likcmb[omm,w,zeta],{w,-2.5,1}];

likcmbmarg[0.25,0.15]

DensityPlot[likcmbmarg[omm,zeta],{omm,0.1,0.4},{zeta,-0.1,0.1}]

DensityPlot[pf1a[zeta]*likcmbmarg[omm,zeta],{omm,0.1,0.4},{zeta,-0.1,0.1}]

DensityPlot[pf1b[zeta]*likcmbmarg[omm,zeta],{omm,0.1,0.4},{zeta,-0.0,0.15}]