Quit[]

(* Monte Carlo sampling // Metropolis Hastings *)

(* Various options and funcs *)

SetDirectory[NotebookDirectory[]];
Needs["ErrorBarPlots`"];
Needs["ComputationalGeometry`"];
Off[CompiledFunction::cflist];
(* This is the jump function // Proposes the next point *)
jumpfuncmat[point_,cijtest_]:=RandomVariate[MultinormalDistribution[point,cijtest]]
(* 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}];
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]
normalize[bins_,counts_]:=(counts/(Max[counts]));

(* end *)

(* 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}]

(* end *)

(* Do a small MCMC to get a test cov. matrix *)

(* Do a small MCMC to get a test cov. matrix *)
priors={{0,5},{-3,3}};
nparams=priors//Length;
step=0;
totalsteps=5000;
initpoint={1,1};
jumpsize={0.01,0.01};
currentpoint=initpoint;
chi2current=chi2@@currentpoint;
chain0={};
ProgressIndicator[Dynamic[step],{0,totalsteps}]
While[step<totalsteps+1,
(
proposedpoint=jumpfuncmat[currentpoint,DiagonalMatrix[jumpsize]^2];
chi2prop=chi2@@proposedpoint;
jumpprob=Exp[-(chi2prop-chi2current)/2];
If[(RandomReal[]≤Min[{1,jumpprob}])&&(And@@Table[IntervalMemberQ[Interval@priors[[i]],proposedpoint[[i]]],{i,1,Length[priors]}]),currentpoint=proposedpoint;chi2current=chi2prop;AppendTo[chain0,{currentpoint,chi2current}];,currentpoint=currentpoint;];
);
step++]

(* Save the chain *)
Export["chain_0.m",chain0]

(* Read in the chain // for postprocessing etc *)
chain0=Import["chain_0.m"];

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

(* (χ^2) plot *)
ListPlot[chain0[[All,2]],Joined->True,PlotRange->{{0,1.02Length[chain0[[All,2]]]},All},Frame->True,FrameLabel->{"N","(χ^2)"}]

(* MCMC plot *)
pp={1,2};
mcmcplot=ListPlot[chain0[[All,1]][[All,pp]],Frame->True,PlotStyle->{Red,PointSize[0.01]},PlotRange->All,Axes->False,Epilog->{{Blue,Point[initpoint[[pp]]],{Black,Point[bf[[pp]]]}}},AspectRatio->1]

(* Cov. Matrix directly from the chain *)
burnin=1000;
chlen=Length[chain0];
Cijchain=Covariance[chain0[[All,1,All]][[burnin;;-1]]]
%//Diagonal//Sqrt

(* end *)

(* This is the actual MCMC *)

priors={{0,5},{-3,3}};
nparams=priors//Length;
(* Import a trial chain to get an estimate for Cij *)
chain0=Import["chain_0.m"];
(* Cov. Matrix directly from the chain *)
burnin=1000;
chlen=Length[chain0];
Cijchain=Covariance[chain0[[All,1,All]][[burnin;;-1]]]

(* This is the main MCMC // Vanilla Metropolis for now *)
step=0;
totalsteps=100000;
initpoint={1,1};
currentpoint=initpoint;
chi2current=chi2@@currentpoint;
MCMC={};
ProgressIndicator[Dynamic[step],{0,totalsteps}]
While[step<totalsteps+1,
(
proposedpoint=jumpfuncmat[currentpoint,Cijchain];
chi2prop=chi2@@proposedpoint;
jumpprob=Exp[-(chi2prop-chi2current)/2];
If[(RandomReal[]≤Min[{1,jumpprob}])&&(And@@Table[IntervalMemberQ[Interval@priors[[i]],proposedpoint[[i]]],{i,1,Length[priors]}]),currentpoint=proposedpoint;chi2current=chi2prop;AppendTo[MCMC,{currentpoint,chi2current}];,currentpoint=currentpoint;];
);
step++]

(* Export the chain to save the result *)
Export["chain_big_1.m",MCMC]

(* Import the chain if starting from a fresh kernel and don't want to run the MCMC again (optional) *)
MCMC=Import["chain_big_1.m"];

(* Check the chains (for burnin etc) *)
ListPlot[MCMC[[All,2]],Joined->True,PlotRange->{{0,1.02Length[MCMC]},{40,150}},Frame->True,FrameLabel->{"N","(χ^2)"}]

(* Statistics *)
(* The best-fit *)
chi2min=Min[MCMC[[All,2]]];
bf=MCMC[[Position[MCMC,chi2min][[1,1]]]][[1]];
mn=Table[Mean[MCMC[[All,1,i]]],{i,1,nparams}];
(* Cov. Matrix directly from the chain *)
burnin=1000;
chlen=Length[MCMC];
Cijchain=Covariance[MCMC[[All,1,All]][[burnin;;-1]]]
Print["The Covariance Matrix is (Subscript[C, ij])=",Cijchain//MatrixForm]
errs=Cijchain//Diagonal//Sqrt;
Print["The bestfit params are ", Table[bf[[i]]±errs[[i]],{i,1,nparams}]]
Print["The  mean  params are  ", Table[mn[[i]]±StandardDeviation[MCMC[[All,1,i]]],{i,1,nparams}]]

datapl=MCMC[[burnin;;-1,1]];
skd=SmoothKernelDistribution[datapl];
pdf=PDF[skd,datapl];
cons=Quantile[pdf,1-Erf[(n/(Sqrt[2]))]/.n->{1,2}//N];

con=ContourPlot[PDF[skd,{a,b}],{a, 2.85,3.35},{b, -1.26, -0.9},Contours->cons,ContourShading->{White,Hue[0.5,0.9,0.9],Hue[0.55,0.8,0.8]},FrameLabel->{"a","b"},BaseStyle->FontSize->18,PlotRange->All,ClippingStyle->None]

hista=SmoothHistogram[datapl[[All,1]],Frame->True,FrameLabel->{"",""},FrameTicks->{Automatic,None},BaseStyle->FontSize->18,AspectRatio->1,Axes->False];
histb=SmoothHistogram[datapl[[All,2]],Frame->True,FrameLabel->{"b",""},FrameTicks->{Automatic,None},BaseStyle->FontSize->18,AspectRatio->1,Axes->False];

plots=GraphicsGrid[{{hista},{con,histb}},Spacings->Scaled[-0.2]]

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

(* Make the plot *)
pl1=Confid2D[data1s,Hue[0.5,0.9,0.9],2];
pl2=Confid2D[data2s,Hue[0.55,0.8,0.8],2];
sh2=Show[pl2,pl1,PlotRange->{{2.75,3.45},{-1.35,-0.8}},Epilog->{Black,Point[bf]},FrameLabel->{"a","b"},AspectRatio->1]

(* end *)

(* end *)