Quit[]

(* Monte Carlo integration *)

data=RandomReal[{-1,1},{10000,2}];
Ntot=Length[data];
datacircle=Select[data,#[[1]]^2+#[[2]]^2<1&];
ListPlot[data,Frame->True,FrameLabel->{"x","y"},AspectRatio->1,Epilog->{{Thickness[.01],Red,Circle[]},{Green,Point[datacircle]}},PlotRange->{{-1.1,1.1},{-1.1,1.1}}]
Nc=datacircle//Length;
Print["((Subscript[Area, cicle])/(Subscript[Area, square])) = ",(Nc/Ntot)//N,", Exact = (π/4) = ",N[Pi/4,10]]
Clear[data,Ntot,datacircle,Nc]

(* end *)

(* Monte Carlo sampling // Metropolis Hastings *)

(* 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 *)
chi2[a_,b_]:=Sum[((data[[i,2]]-model[data[[i,1]],a,b])/data[[i,3]])^2,{i,1,ndata}]

(* This is the jump function // Proposes the next point *)
jumpfunc[point_,s_]:=point+RandomVariate[NormalDistribution[0,s],Length[point]]

jumpfunc[{1,1},0.0001]
jumpfunc[{1,1},0.5]
jumpfunc[{1,1},10]

(* Do the MCMC *)
(* http://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm#Intuition *)
(* Choose the total number of steps *)
totalsteps=30000;
(* Choose the initial point *)
initpoint={2.5,-0.5};
(* Initialize the MCMC *)
step=0;
(* This is the current point *)
currentpoint=initpoint;
(* This is the (χ^2) at the current point // Notice the @@ trick (Apply) // f@@{a,b,c}=f[a,b,c] *)
chi2current=chi2@@currentpoint;
(* Initialize the chain. This is were the points are stored *)
chain={};
(* Create a progress indicator. This just shows the progress *)
ProgressIndicator[Dynamic[step],{0,totalsteps}]
(* Start the MCMC *)
While[step<totalsteps+1,
(
(* Propose a new point based on the jump function. In this case the standard dev is 0.01 *)
proposedpoint=jumpfunc[currentpoint,0.01];
(* This is the (χ^2) at the proposed point // Notice the @@ trick *)
chi2prop=chi2@@proposedpoint;
(* This is the probability to jump to the new point=((Exp[-(χ^2)(prop)])/(Exp[-(χ^2)(current)]))*)
jumpprob=Exp[-(chi2prop-chi2current)/2];
(* Decide whether to keep the new point. If so, store it in the chain, otherwise start over. *)
If[RandomReal[]≤Min[{1,jumpprob}],currentpoint=proposedpoint;chi2current=chi2prop;AppendTo[chain,{currentpoint,chi2current}];,currentpoint=currentpoint;];
);
step++]

(* Uncomment the following if you want to save the chain. Very helpful if the calculation takes for ever :) . Later, use Import[] to read it back in *)
(* Export["chain_1.m",chain] *)

(* Find the best-fit by minimizing the (χ^2) *)
fmin=FindMinimum[chi2[a,b],{a,1},{b,1}]
fmin[[1]]/(ndata-2)

(* MCMC plot *)
mcmcplot=ListPlot[chain[[All,1]],Frame->True,FrameLabel->{"a","b"},PlotStyle->{Red,PointSize[0.003]},PlotRange->{{2.35,3.5},{-1.35,-0.45}},Axes->False,Epilog->{{Blue,Point[initpoint],{Black,Point[fmin[[2]][[All,2]]]}}},AspectRatio->1]

(* Fisher Matrix estimation of errors *)
p={a,b};
(* Fisher matrix *)
Fij=Table[(1/2)D[chi2[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[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[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[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]*)

sh1=Show[mcmcplot,imp3s[[2]],imp2s[[2]],imp1s[[2]],PlotRange->{{2.45,3.45},{-1.35,-0.45}},AspectRatio->1,FrameLabel->{"a","b"}]

(* 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->2,SplineKnots->"Unclamped"]],Antialiasing->True}],Frame->True];qq]

(* Best-fit from the MCMC *)
chi2min=Min[chain[[All,2]]]
bf=chain[[Position[chain,chi2min][[1,1]]]][[1]]

(* Select the points in the 1,2,3 sigma region *)
data1s=Select[chain,#[[2]]≤(chi2min+dx2[[1]])&][[All,1]];
data2s=Select[chain,#[[2]]≤(chi2min+dx2[[2]])&][[All,1]];
data3s=Select[chain,#[[2]]≤(chi2min+dx2[[3]])&][[All,1]];

(* 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[bf]},FrameLabel->{"a","b"},AspectRatio->1]

(* Overlay the MCMC points *)
sh11=Show[mcmcplot,imp3s[[2]],imp2s[[2]],imp1s[[2]],PlotRange->{{2.75,3.45},{-1.35,-0.8}},AspectRatio->1,FrameLabel->{"a","b"},Epilog->{{Blue,Point[data3s]},{Green,Point[data2s]},{Red,Point[data1s]}}];

(* Compare the plots *)
GraphicsGrid[{{sh11,sh2}}]

(* Mean errors for a and b from contours *)
bd=FindRoot[{chi2[a,b]==fmin[[1]]+2.30,(D[chi2[a,b],a])==0},{a,3.2},{b,-1.1}][[All,2]][[2]];
bu=FindRoot[{chi2[a,b]==fmin[[1]]+2.30,(D[chi2[a,b],a])==0},{a,2.2},{b,-0.6}][[All,2]][[2]];
au=FindRoot[{chi2[a,b]==fmin[[1]]+2.30,(D[chi2[a,b],b])==0},{a,3.2},{b,-1.1}][[All,2]][[1]];
ad=FindRoot[{chi2[a,b]==fmin[[1]]+2.30,(D[chi2[a,b],b])==0},{a,2.2},{b,-0.6}][[All,2]][[1]];

(* Cramer-Rao bound works!!!!!! *)
{{au-fmin[[2,1,2]],fmin[[2,1,2]]-ad}//Mean,{bu-fmin[[2,2,2]],fmin[[2,2,2]]-bd}//Mean}
(Sqrt[Cij])//Diagonal

(* Evolution of (χ^2)(N), a, b

ListPlot[chain[[All,2]],Joined->True,PlotRange->{All,{50,100}},Frame->True,FrameLabel->{"N","(χ^2)"}]

ListPlot[chain[[All,2]],Joined->True,PlotRange->{{0,2000},{50,100}},Frame->True,FrameLabel->{"N","(χ^2)"}]

ListPlot[chain[[All,1,1]],Joined->True,PlotRange->{2,4},Frame->True,FrameLabel->{"N","a"}]
ListPlot[chain[[All,1,2]],Joined->True,PlotRange->{-2,0},Frame->True,FrameLabel->{"N","b"}]

(* Estimate of the Covariance matrix *)

(* Fisher matrix approach *)
Cij

(* Directly from the chain // Remove the "burnin" *)
burnin=500;
chlen=Length[chain];
c11=Mean[Table[chain[[i,1,1]]^2,{i,burnin,chlen}]]-Mean[Table[chain[[i,1,1]],{i,burnin,chlen}]]^2;
c22=Mean[Table[chain[[i,1,2]]^2,{i,burnin,chlen}]]-Mean[Table[chain[[i,1,2]],{i,burnin,chlen}]]^2;
c12=Mean[Table[chain[[i,1,1]]chain[[i,1,2]],{i,burnin,chlen}]]-Mean[Table[chain[[i,1,1]],{i,burnin,chlen}]]Mean[Table[chain[[i,1,2]],{i,burnin,chlen}]];
Cijchain={{c11,c12},{c12,c22}}

(Sqrt[Diagonal[Cij]])
(Sqrt[Diagonal[Cijchain]])

(* end *)