Quit[]

(* Bootstrap Monte Carlo *)

(* Various options *)
SetDirectory[NotebookDirectory[]];
Needs["ErrorBarPlots`"];
Needs["ComputationalGeometry`"];
Off[CompiledFunction::cflist];

(* Some fake data with noise *)
SeedRandom[1264645];
(* A model *)
model[x_,a_,b_]:=a x^2+b x
(* Mock data with noise *)
data=Table[{x,model[x,3,-1]+RandomVariate[NormalDistribution[0,.1]],RandomReal[{.05,.1}]},{x,0,1,0.02}];
ndata=Length[data];
(* Plot of the data and model *)
Show[ErrorListPlot[data,PlotRange->All],Plot[model[x,3,-1],{x,0,1},PlotStyle->Red],Frame->True,FrameLabel->{"x","y"},Axes->False]

(* Define the chi^2 as a function of the data and the parameters a,b *)
chi2[dat_,a_,b_]:=Sum[((dat[[i,2]]-model[dat[[i,1]],a,b])/dat[[i,3]])^2,{i,1,Length[dat]}]

(* Highlight the two different ways of sampling *)

Range[30]

(* Unique sample *)
RandomSample[Range[30],20]//Sort
%//Length

(* Sample with replacement *)
RandomChoice[Range[30],20]//Sort
%//Length

(* Launch some Kernels to speed up the work *)
LaunchKernels[7];

(* Sample with replacement and make the new datasets *)
mockdata=Table[RandomChoice[data,Length[data]]//Sort,{1000}];

(* The best-fit of the real data *)
fmin0=FindMinimum[chi2[data,a,b],{a,1},{b,1}]

(* The best-fit of all of the mock data // This in practice will take the most time *)
fminmock=ParallelTable[FindMinimum[chi2[mockdata[[j]],a,b],{a,1},{b,1}],{j,1,Length[mockdata]},DistributedContexts:>Automatic];

(* Gather and plot the parameters (a,b) *)
mock2D=fminmock[[All,2,{1,2},2]];
amock=fminmock[[All,2,1,2]];
bmock=fminmock[[All,2,2,2]];
ListPlot[mock2D,Frame->True,Epilog->{PointSize[0.02],Red,Point[fmin0[[2,All,2]]]},FrameLabel->{"a","b"},AspectRatio->1]

fmin0

(* The mean values *)
{Mean[fminmock[[All,1]]],Mean[amock],Mean[bmock]}

(* ... and the errors *)
(Sqrt[Mean[Table[amock[[i]]^2,{i,1,Length[amock]}]]-Mean[amock]^2])
(Sqrt[Mean[Table[bmock[[i]]^2,{i,1,Length[bmock]}]]-Mean[bmock]^2])

(* Fisher Matrix estimation of errors *)
p={a,b};
(* Find the best-fit by minimizing the (χ^2) *)
fmin=FindMinimum[chi2[data,a,b],{a,1},{b,1}]
fmin[[1]]/(ndata-2)
(* Fisher matrix *)
Fij=Table[(1/2)D[chi2[data,a,b],p[[i]],p[[j]]]/.fmin[[2]],{i,1,Length[p]},{j,1,Length[p]}];
Print["(Subscript[F, ij])= ",Fij//MatrixForm]
(* Covariance matrix *)
Cij=Inverse[Fij];
Print["(Subscript[C, ij])= ",Cij//MatrixForm]
(* Parameter errors *)
Print["1σ errors = ",(Sqrt[Cij])//Diagonal]

(* Traditional contours *)

(* The values of ((Δχ)^2) // See Sigmas.nb on the site *)
Table[{"((Δχ)^2)("<>ToString[xσ]<>"σ,n="<>ToString[n]<>") =",N[2 InverseGammaRegularized[(n/2),1-Erf[((xσ)/(Sqrt[2]))]],6]},{n,1,4},{xσ,1,4}]//TableForm

(* The values of ((Δχ)^2) for 1,2,3σ *)
dx2=Table[2 InverseGammaRegularized[(Nn/2),1-Erf[((xσ)/(Sqrt[2]))]]/.Nn->2//N,{xσ,1,3}]

(* Find the contours by solving (χ^2)(a,b)=((Subscript[χ, min])^2)+((Δχ)^2)*)
(* 1sigma contour *)
imp1s=ContourPlot[chi2[data,a,b]==fmin[[1]]+dx2[[1]],{a,2.9,3.3},{b,-1.2,-0.95},AspectRatio->1,Frame->True,Axes->False]//AbsoluteTiming;
(* 2sigma contour *)
imp2s=ContourPlot[chi2[data,a,b]==fmin[[1]]+dx2[[2]],{a,2.5,3.8},{b,-1.5,-0.5},AspectRatio->1,Frame->True,Axes->False]//AbsoluteTiming;
(* 3sigma contour *)
imp3s=ContourPlot[chi2[data,a,b]==fmin[[1]]+dx2[[3]],{a,2.0,4.3},{b,-1.9,-0.0},AspectRatio->1,Frame->True,Axes->False]//AbsoluteTiming;
Show[imp1s[[2]],imp2s[[2]],imp3s[[2]],PlotRange->All]
(* Uncomment the following if you want to save the data of the contours. Very helpful if the calculation takes for ever :) . Later, use Import[] to read it back in *)
(*Export["contour_1s_model_1.m",imp1s]*)

(* Contours+color from the MCMC // Uses ConvexHull // Advanced++ *)
Confid2D[data_List,color_,inter_]:=Module[{cv,cv1,cv2,qq},cv=ConvexHull[data,AllPoints->True];cv1=data[[#]]&/@cv;cv2=Join[cv1,{cv1[[1]],cv1[[2]]}];qq=Show[Graphics[{color,FilledCurve[BSplineCurve[cv2,SplineDegree->inter,SplineKnots->"Unclamped"]],Antialiasing->True}],Frame->True];qq]

(* Cavemans approach to finding 68% region around best fit // Needs improvement *)

datatest=RandomVariate[NormalDistribution[0,1],{10000}];
Histogram[datatest]

stadchi2=StandardDeviation[datatest];
meanchi2=Mean[datatest];
data1s=Select[datatest,(meanchi2-stadchi2≤#≤meanchi2+stadchi2)&];
data2s=Select[datatest,(meanchi2-2stadchi2≤#≤meanchi2+2stadchi2)&];
data3s=Select[datatest,(meanchi2-3stadchi2≤#≤meanchi2+3stadchi2)&];

Histogram[{data3s,data2s,data1s},{0.5},ChartStyle->{Blue,Green,Red}]

(* Gather the results in the form ((χ^2),a,b)  *)
datachi2=SortBy[Table[{fminmock[[i,1]],a,b}/.fminmock[[i,2]],{i,1,Length[fminmock]}],First];

Histogram[datachi2[[All,1]],{5}]

(* Select the points in the 1,2,3 sigma region *)
stadchi2=StandardDeviation[datachi2[[All,1]]];
meanchi2=Mean[datachi2[[All,1]]];
data1s=Select[datachi2,(meanchi2-stadchi2≤#[[1]]≤meanchi2+stadchi2)&][[All,{2,3}]];
data2s=Select[datachi2,(meanchi2-2stadchi2≤#[[1]]≤meanchi2+2stadchi2)&][[All,{2,3}]];
data3s=Select[datachi2,(meanchi2-3stadchi2≤#[[1]]≤meanchi2+3stadchi2)&][[All,{2,3}]];

(* Make the plot *)
pl1=Confid2D[data1s,Red,2];
pl2=Confid2D[data2s,Green,2];
pl3=Confid2D[data3s,Blue,2];
sh2=Show[pl3,pl2,pl1,PlotRange->{{2.75,3.45},{-1.35,-0.8}},Epilog->{Black,Point[fmin0[[2,All,2]]]},FrameLabel->{"a","b"},AspectRatio->1]

(* Overlay the real contours *)
sh11=Show[sh2,imp3s[[2]],imp2s[[2]],imp1s[[2]],PlotRange->{{2.75,3.45},{-1.35,-0.8}},AspectRatio->1,FrameLabel->{"a","b"}]

(* The bootstrap is not very good at making contours!!!! Use MCMC instead!!! *)

(* Estimate of the Covariance matrix *)

(* Fisher matrix approach *)
Cij

(* Directly from the Bootstrap *)
c11=Mean[Table[amock[[i]]^2,{i,1,Length[amock]}]]-Mean[Table[amock[[i]],{i,1,Length[amock]}]]^2;
c22=Mean[Table[bmock[[i]]^2,{i,1,Length[bmock]}]]-Mean[Table[bmock[[i]],{i,1,Length[bmock]}]]^2;
c12=Mean[Table[amock[[i]]bmock[[i]],{i,1,Length[amock]}]]-Mean[Table[amock[[i]],{i,1,Length[amock]}]]Mean[Table[bmock[[i]],{i,1,Length[bmock]}]];
CijBoot={{c11,c12},{c12,c22}}

(Sqrt[Diagonal[Cij]])
(Sqrt[Diagonal[CijBoot]])

(* end *)