MCMC based on Likelihood

data=Import["http://hep.fuw.edu.pl/u/zarnecki/saed23/11_homework_data.dat","Data"][[1]];

data[[1;;22]]

Options[myExponentialPlusConts]={p0->{1,1,1}}
myExponentialPlusConts[x_,opt:OptionsPattern[]]:=Module[{tdec=OptionValue[p0][[1]],ndec=OptionValue[p0][[2]],nb=OptionValue[p0][[3]]},
ndec/tdec*Exp[-x/tdec]+nb
]

ChiSquareTestStatistics[x_List,y_List,err_List,f_,parameters_List]:=Module[{vals},
vals=Map[f[#,p0->parameters]&,x];
Total[(((vals-y)/err))^2]
]

ChiSquareTestStatistics[{x1,x2},{y1,y2},{dy1,dy2},myExponentialPlusConts,{a,b,c}]

{bins,counts}=HistogramList[data,30,"Count"];
bincenters=((bins[[1;;-2]]+bins[[2;;-1]])/2);
Map[{#,myExponentialPlusConts[#,p0->{2.49,455,2.8}]}&,bincenters];
Show[Histogram[data,30,"LogCount",ChartStyle->Blue],ListLogPlot[%,Joined->True,PlotStyle->Red]]
ChiSquareTestStatistics[bincenters,counts,Sqrt[counts],myExponentialPlusConts,{2.49,455,2.8}]

?AppendTo

NMC=1000;
parsHistory={};
chiHistory={};
step={0.2,50,1};
parametersInitialValues={7.2,722,7.12};
par=parametersInitialValues;
chi2old=ChiSquareTestStatistics[bincenters,counts,Sqrt[counts],myExponentialPlusConts,par];
chiHistory={chi2old};

For[i=1,i<=NMC,i++,
u=RandomReal[{-1,1},3];
dpar=u*step;
chi2new=ChiSquareTestStatistics[bincenters,counts,Sqrt[counts],myExponentialPlusConts,par+dpar];
prob=Exp[-0.5*(chi2new-chi2old)];

r=RandomReal[{0,1},1][[1]];
If[r<prob,
 par=par+dpar;
chi2old=chi2new;
];
AppendTo[parsHistory,par];
AppendTo[chiHistory,chi2new];

]
Length[parsHistory]

Options[infoHistogram]={start->1,end->-1,Nbins->30,norm->"PDF",style->Opacity[0.33,Blue]}
infoHistogram[dataset_,opt:OptionsPattern[]]:=Module[{$start=OptionValue[start],$end=OptionValue[end],$nbins=OptionValue[Nbins],$data},
$data=dataset[[$start;;$end]];
Histogram[
$data,$nbins,OptionValue[norm],ChartStyle->OptionValue[style],Frame->True,ChartLegends->{Around[Mean[$data],StandardDeviation[$data]]}
]
]

Show[infoHistogram[parsHistory[[All,1]],Nbins->30,style->Opacity[0.5,Blue]],infoHistogram[parsHistory[[All,1]],Nbins->30,Nbins->30,start->-200,style->Opacity[0.5,Red]],FrameLabel->{"τ","PDF"}]

Show[infoHistogram[parsHistory[[All,2]],Nbins->30,style->Opacity[0.5,Blue]],infoHistogram[parsHistory[[All,2]],Nbins->30,Nbins->30,start->-200,style->Opacity[0.5,Red]],FrameLabel->{"(Subscript[N, μ])","PDF"}]

Show[infoHistogram[parsHistory[[All,3]],Nbins->30,style->Opacity[0.5,Blue]],infoHistogram[parsHistory[[All,3]],Nbins->30,Nbins->30,start->-200,style->Opacity[0.5,Red]],FrameLabel->{"(Subscript[N, B])","PDF"}]

Options[infoEvolution]={start->1,end->-1,offset->0,style->Opacity[0.33,Blue]}
infoEvolution[dataset_,opt:OptionsPattern[]]:=Module[{$start=OptionValue[start],$offset=OptionValue[offset],$end=OptionValue[end],$data},
$data=dataset[[$start;;$end]];
ListPlot[
Transpose[
{Table[t+$offset,{t,1,Length[$data]}],$data}
],Joined->True,PlotStyle->OptionValue[style],PlotLegends->{Around[Mean[$data],StandardDeviation[$data]]},PlotRange->All,Frame->True
]
]

Show[infoEvolution[parsHistory[[All,1]],style->Opacity[0.5,Blue]],infoEvolution[parsHistory[[All,1]],start->800,offset->800,style->Opacity[0.995,Red]],FrameLabel->{"time step","τ"}]

Show[infoEvolution[parsHistory[[All,2]],style->Opacity[0.5,Blue]],infoEvolution[parsHistory[[All,2]],start->800,offset->800,style->Opacity[0.99,Red]],FrameLabel->{"time step","(Subscript[N, μ])"}]

Show[infoEvolution[parsHistory[[All,3]],style->Opacity[0.5,Blue]],infoEvolution[parsHistory[[All,3]],start->800,offset->800,style->Opacity[0.995,Red]],FrameLabel->{"time step","(Subscript[N, B])"}]

ListPlot[parsHistory[[All,1]],PlotRange->All,Frame->True,FrameLabel->{"step","τ"}]
ListPlot[parsHistory[[All,2]],PlotRange->All,Frame->True,FrameLabel->{"step","(Subscript[N, μ])"}]
ListPlot[parsHistory[[All,3]],PlotRange->All,Frame->True,FrameLabel->{"step","(Subscript[N, b])"}]
ListPlot[chiHistory,PlotRange->All,Frame->True,FrameLabel->{"step","(Subscript[((χ^2)), proposed])"}]

prior[pars_]:=Tr[Map[If[(#>0),1,0]&,pars],Times]
prior[{1,1,-1}]
prior[{1,1,21}]
prior[{1,1,a}]

NMC=1000;
parsHistory={};
chiHistory={};
step={0.2,50,1};
parametersInitialValues={7.2,722,7.12};
par=parametersInitialValues;
chi2old=ChiSquareTestStatistics[bincenters,counts,Sqrt[counts],myExponentialPlusConts,par];
chiHistory={chi2old};

For[i=1,i<=NMC,i++,
u=RandomReal[{-1,1},3];
dpar=u*step;
chi2new=ChiSquareTestStatistics[bincenters,counts,Sqrt[counts],myExponentialPlusConts,par+dpar];

LLR=Exp[-0.5*(chi2new-chi2old)];
priorRatio=((prior[par+dpar])/(prior[par]));
prob=LLR*priorRatio;

r=RandomReal[{0,1},1][[1]];
If[r<prob,
 par=par+dpar;
chi2old=chi2new;
];
AppendTo[parsHistory,par];
AppendTo[chiHistory,chi2new];

]
Length[parsHistory]

Show[infoEvolution[parsHistory[[All,1]],style->Opacity[0.5,Blue]],infoEvolution[parsHistory[[All,1]],start->800,offset->800,style->Opacity[0.995,Red]],FrameLabel->{"time step","τ"}]

Show[infoEvolution[parsHistory[[All,2]],style->Opacity[0.5,Blue]],infoEvolution[parsHistory[[All,2]],start->800,offset->800,style->Opacity[0.99,Red]],FrameLabel->{"time step","(Subscript[N, μ])"}]

Show[infoEvolution[parsHistory[[All,3]],style->Opacity[0.5,Blue]],infoEvolution[parsHistory[[All,3]],start->800,offset->800,style->Opacity[0.995,Red]],FrameLabel->{"time step","(Subscript[N, B])"}]

Definizione delle funzioni

PriorRatio[prior_,thetaprimo_,theta_]:=((prior[thetaprimo])/(prior[theta]))
productPorbRatio[listPR_]:=Tr[listPR,Times]  (*Produttoria dei rapporti LLR *)

Options[generateNextPoint]={debug->False}
generateNextPoint[par_,step_,opt:OptionsPattern[]]:=Module[{u,dpar,res},
u=RandomReal[{-1,1},3];
dpar=u*step;
If[OptionValue[debug]==True,Print[dpar]];
res=par+dpar;
If[OptionValue[debug]==True,Print[res]];
res
]

Options[EvaluteLLR]={weights->None,debug->False};
EvaluteLLR[myExponentialPlusConts_,opt:OptionsPattern[]]:=Module[{chi2old,chi2new,$weights,LLR},
If[Not[OptionValue[weights]==None],$weights=OptionValue[weights],$weights=Sqrt[counts]];
chi2old=ChiSquareTestStatistics[bincenters,counts,$weights,myExponentialPlusConts,par];
chi2new=ChiSquareTestStatistics[bincenters,counts,$weights,myExponentialPlusConts,proposedParameters];
If[OptionValue[debug],Print[{chi2old,chi2new}]];
LLR=Exp[-0.5*(chi2new-chi2old)];
{LLR,chi2new}
];

AcceptReject[prob_,par_,proposedParameters_]:=Module[{r,res=par},
r=RandomReal[{0,1},1][[1]];
If[r<prob,
 res=proposedParameters;
(*chi2old=chi2new;*)
];
 res
]

Ciclo MCMC

NMC=100;
History={};
(*chiHistory={};*)
step={0.2,50,1};
parametersInitialValues={7.2,722,7.12};
par=parametersInitialValues;
chi2old=ChiSquareTestStatistics[bincenters,counts,Sqrt[counts],myExponentialPlusConts,par];
$parsHistory=<|"parameters"->par,"proposedChi2"->chi2old|>;
(*AppendTo[History,$parsHistory];*)

For[i=1,i<=NMC,i++,
(*  adatto step in base alla storia della catena  *)

(*  generazione del prossimo punto *)
proposedParameters=generateNextPoint[par,step,debug->False];
$parsHistory["proposedParameter"]=proposedParameters;
(* valutazione della LLR *)
{LLR,proposedChi2}=EvaluteLLR[myExponentialPlusConts,debug->False];

(* valutazione dei prior *)
priorRatio=PriorRatio[prior,proposedParameters,par];

(* valutazione del fattore di Metropolis *)
prob=productPorbRatio[{LLR,priorRatio}];

(*  AcceptReject par or proposedParameters *)
par=AcceptReject[prob,par,proposedParameters];

$parsHistory=<|"parameters"->par,"proposedChi2"->proposedChi2|>;
(* aumento la catena *)
AppendTo[History,$parsHistory];
Unset[$parsHistory];

]
Length[History]

Map[#["parameters"]&,History]
Map[#["proposedChi2"]&,History]