#ifndef RooStats_RooStatsUtils
#include "RooStats/RooStatsUtils.h"
#endif
#ifndef ROOT_Rtypes
#include "Rtypes.h"
#endif
#include "RooRealVar.h"
#include "RooNLLVar.h"
#include "RooGlobalFunc.h"
#include "RooDataSet.h"
#include "RooArgSet.h"
#include "RooArgList.h"
#include "TRandom.h"
#include "TH1.h"
#include "TMath.h"
#include "TFile.h"
#include "RooStats/MCMCCalculator.h"
#include "RooStats/MCMCInterval.h"
ClassImp(RooStats::MCMCCalculator);
using namespace RooFit;
using namespace RooStats;
MCMCCalculator::MCMCCalculator()
{
fWS = NULL;
fPOI = NULL;
fNuisParams = NULL;
fOwnsWorkspace = false;
fPropFunc = NULL;
fPdfName = NULL;
fDataName = NULL;
fNumIters = 0;
fNumBurnInSteps = 0;
fNumBins = 0;
fAxes = NULL;
}
MCMCCalculator::MCMCCalculator(RooWorkspace& ws, RooAbsData& data,
RooAbsPdf& pdf, RooArgSet& paramsOfInterest,
ProposalFunction& proposalFunction, Int_t numIters,
RooArgList* axes, Double_t size)
{
fOwnsWorkspace = false;
SetWorkspace(ws);
SetData(data);
SetPdf(pdf);
SetParameters(paramsOfInterest);
SetTestSize(size);
SetProposalFunction(proposalFunction);
fNumIters = numIters;
fNumBurnInSteps = 0;
fNumBins = 0;
fAxes = axes;
}
MCMCCalculator::MCMCCalculator(RooAbsData& data, RooAbsPdf& pdf,
RooArgSet& paramsOfInterest, ProposalFunction& proposalFunction,
Int_t numIters, RooArgList* axes, Double_t size)
{
fWS = new RooWorkspace();
fOwnsWorkspace = true;
SetData(data);
SetPdf(pdf);
SetParameters(paramsOfInterest);
SetTestSize(size);
SetProposalFunction(proposalFunction);
fNumIters = numIters;
fNumBurnInSteps = 0;
fNumBins = 0;
fAxes = axes;
}
MCMCInterval* MCMCCalculator::GetInterval() const
{
RooAbsPdf* pdf = fWS->pdf(fPdfName);
RooAbsData* data = fWS->data(fDataName);
if (!data || !pdf || !fPOI) return 0;
RooArgSet x;
RooArgSet xPrime;
RooRealVar* w = new RooRealVar("w", "weight", 0);
RooArgSet* parameters = pdf->getParameters(data);
RemoveConstantParameters(parameters);
x.addClone(*parameters);
x.addOwned(*w);
xPrime.addClone(*parameters);
RooDataSet* points = new RooDataSet("points", "Markov Chain", x, WeightVar(*w));
TRandom gen;
RooArgSet* constrainedParams = pdf->getParameters(*data);
RooAbsReal* nll = pdf->createNLL(*data, Constrain(*constrainedParams));
delete constrainedParams;
RooArgSet* nllParams = nll->getParameters(*data);
Int_t weight = 0;
for (int i = 0; i < fNumIters; i++) {
if (i % 100 == 0){
fprintf(stdout, ".");
fflush(NULL);
}
fPropFunc->Propose(xPrime, x);
RooStats::SetParameters(&xPrime, nllParams);
Double_t xPrimeNLL = nll->getVal();
RooStats::SetParameters(&x, nllParams);
Double_t xNLL = nll->getVal();
Double_t diff = xPrimeNLL - xNLL;
if (!fPropFunc->IsSymmetric(xPrime, x))
diff += TMath::Log(fPropFunc->GetProposalDensity(xPrime, x)) -
TMath::Log(fPropFunc->GetProposalDensity(x, xPrime));
if (diff < 0.0) {
points->addFast(x, (Double_t)weight);
weight = 1;
RooStats::SetParameters(&xPrime, &x);
}
else {
Double_t rand = TMath::Log(gen.Uniform(1.0));
if (-1.0 * rand >= diff) {
points->addFast(x, (Double_t)weight);
weight = 1;
RooStats::SetParameters(&xPrime, &x);
} else {
weight++;
}
}
}
delete nllParams;
printf("\n");
points->addFast(x, (Double_t)weight);
MCMCInterval* interval = new MCMCInterval("mcmcinterval", "MCMCInterval",
*fPOI, *points);
if (fAxes != NULL)
interval->SetAxes(*fAxes);
if (fNumBins > 0)
interval->SetNumBins(fNumBins);
if (fNumBurnInSteps > 0)
interval->SetNumBurnInSteps(fNumBurnInSteps);
interval->SetConfidenceLevel(1.0 - fSize);
return interval;
}