MarkovChainMC Class Reference

#include <MarkovChainMC.h>

Inheritance diagram for MarkovChainMC:

Public Member Functions
virtual void	applyOptions (const boost::program_options::variables_map &vm)
	MarkovChainMC ()
virtual const std::string &	name () const
bool	run (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint)
Public Member Functions inherited from LimitAlgo
virtual void	applyDefaultOptions ()
	LimitAlgo ()
	LimitAlgo (const char *desc)
const boost::program_options::options_description &	options () const
virtual void	setNToys (const int)
virtual void	setToyNumber (const int)

Private Types
enum	ProposalType { FitP, UniformP, MultiGaussianP, TestP }

Private Member Functions
int	guessBurnInSteps (const RooStats::MarkovChain &chain) const
void	limitAndError (double &limit, double &limitErr, const std::vector< double > &limits) const
void	limitFromChain (double &limit, double &limitErr, const RooArgSet &poi, RooStats::MarkovChain &chain, int burnInSteps=-1)
RooStats::MarkovChain *	mergeChains (const RooArgSet &poi, const std::vector< double > &limits) const
void	readChains (const RooArgSet &poi, std::vector< double > &limits)
int	runOnce (RooWorkspace *w, RooStats::ModelConfig *mc_s, RooStats::ModelConfig *mc_b, RooAbsData &data, double &limit, double &limitErr, const double *hint) const
RooStats::MarkovChain *	slimChain (const RooArgSet &poi, const RooStats::MarkovChain &chain) const
int	stationarityTest (const RooStats::MarkovChain &chain, const RooArgSet &poi, int nchunks) const
	note: this is still being developed More...

Private Attributes
TList	chains_
float	mass_
	Mass of the Higgs boson (goes into the name of the saved chains) More...
int	modelNDF_
	Number of degrees of freedom of the problem, approximately. More...

Static Private Attributes
static bool	adaptiveBurnIn_ = false
	Adaptive burn-in (experimental!) More...
static bool	adaptiveTruncation_ = true
	do adaptive truncated mean More...
static bool	alwaysStepPoi_ = true
static float	burnInFraction_ = 0.25
	Discard these fraction of points. More...
static unsigned int	burnInSteps_ = 200
	Discard these points. More...
static float	cropNSigmas_ = 0
static int	debugProposal_ = false
static float	hintSafetyFactor_ = 5.
	Safety factor for hint (integrate up to this number of times the hinted limit) More...
static unsigned int	iterations_ = 10000
	Propose this number of points for the chain. More...
static bool	mergeChains_ = false
	Merge chains instead of averaging limits. More...
static bool	noReset_ = false
static bool	noSlimChain_ = false
	Leave all parameters in the markov chain, not just the POI. More...
static unsigned int	numberOfBins_
static unsigned int	proposalHelperCacheSize_
static float	proposalHelperUniformFraction_ = 0.0
static float	proposalHelperWidthRangeDivisor_ = 5.
static ProposalType	proposalType_ = TestP
static std::string	proposalTypeName_ = "ortho"
static bool	readChains_ = false
	Read chains from file instead of running them. More...
static bool	runMinos_ = false
static bool	saveChain_ = false
	Save Markov Chain in output file. More...
static unsigned int	tries_ = 10
	compute the limit N times More...
static float	truncatedMeanFraction_ = 0.0
	Ignore up to this fraction of results if they're too far from the median. More...
static bool	updateHint_ = false
static bool	updateProposalParams_ = false

Additional Inherited Members
Protected Attributes inherited from LimitAlgo
boost::program_options::options_description	options_

Detailed Description

Interface to RooStats MCMC, with handing of multiple chains

Author

Giovanni Petrucciani (UCSD)

Definition at line 16 of file MarkovChainMC.h.

Member Enumeration Documentation

enum MarkovChainMC::ProposalType

private

Enumerator
FitP
UniformP
MultiGaussianP
TestP

Definition at line 26 of file MarkovChainMC.h.

{ FitP, UniformP, MultiGaussianP, TestP };

Constructor & Destructor Documentation

MarkovChainMC::MarkovChainMC

(

)

Definition at line 56 of file MarkovChainMC.cc.

References adaptiveBurnIn_, adaptiveTruncation_, alwaysStepPoi_, burnInFraction_, burnInSteps_, cropNSigmas_, debugProposal_, hintSafetyFactor_, iterations_, LimitAlgo::options_, proposalHelperUniformFraction_, proposalHelperWidthRangeDivisor_, proposalTypeName_, tries_, truncatedMeanFraction_, and updateProposalParams_.

                             : 
    LimitAlgo("Markov Chain MC specific options") 
{
    options_.add_options()
        ("iteration,i", boost::program_options::value<unsigned int>(&iterations_)->default_value(iterations_), "Number of iterations")
        ("tries", boost::program_options::value<unsigned int>(&tries_)->default_value(tries_), "Number of times to run the MCMC on the same data")
        ("burnInSteps,b", boost::program_options::value<unsigned int>(&burnInSteps_)->default_value(burnInSteps_), "Burn in steps (absolute number)")
        ("burnInFraction", boost::program_options::value<float>(&burnInFraction_)->default_value(burnInFraction_), "Burn in steps (fraction of total accepted steps)")
        ("adaptiveBurnIn", boost::program_options::value<bool>(&adaptiveBurnIn_)->default_value(adaptiveBurnIn_), "Adaptively determine burn in steps (experimental!).")
        ("proposal", boost::program_options::value<std::string>(&proposalTypeName_)->default_value(proposalTypeName_), 
                              "Proposal function to use: 'fit', 'uniform', 'gaus', 'ortho' (also known as 'test')")
        ("runMinos",          "Run MINOS when fitting the data")
        ("noReset",           "Don't reset variable state after fit")
        ("updateHint",        "Update hint with the results")
        ("updateProposalParams", 
                boost::program_options::value<bool>(&updateProposalParams_)->default_value(updateProposalParams_), 
                "Control ProposalHelper::SetUpdateProposalParameters")
        ("propHelperWidthRangeDivisor", 
                boost::program_options::value<float>(&proposalHelperWidthRangeDivisor_)->default_value(proposalHelperWidthRangeDivisor_), 
                "Sets the fractional size of the gaussians in the proposal")
        ("alwaysStepPOI", boost::program_options::value<bool>(&alwaysStepPoi_)->default_value(alwaysStepPoi_),
                            "When using 'ortho' proposal, always step also the parameter of interest. On by default, as it improves convergence, but you can turn it off (e.g. if you turn off --optimizeSimPdf)")
        ("propHelperUniformFraction", 
                boost::program_options::value<float>(&proposalHelperUniformFraction_)->default_value(proposalHelperUniformFraction_), 
                "Add a fraction of uniform proposals to the algorithm")
        ("debugProposal", boost::program_options::value<int>(&debugProposal_)->default_value(debugProposal_), "Printout the first N proposals")
        ("cropNSigmas", 
                boost::program_options::value<float>(&cropNSigmas_)->default_value(cropNSigmas_),
                "crop range of all parameters to N times their uncertainty") 
        ("truncatedMeanFraction", 
                boost::program_options::value<float>(&truncatedMeanFraction_)->default_value(truncatedMeanFraction_), 
                "Discard this fraction of the results before computing the mean and rms")
        ("adaptiveTruncation", boost::program_options::value<bool>(&adaptiveTruncation_)->default_value(adaptiveTruncation_),
                            "When averaging multiple runs, ignore results that are more far away from the median than the inter-quartile range")
        ("hintSafetyFactor",
                boost::program_options::value<float>(&hintSafetyFactor_)->default_value(hintSafetyFactor_),
                "set range of integration equal to this number of times the hinted limit")
        ("saveChain", "Save MarkovChain to output file")
        ("noSlimChain", "Include also nuisance parameters in the chain that is saved to file")
        ("mergeChains", "Merge MarkovChains instead of averaging limits")
        ("readChains", "Just read MarkovChains from toysFile instead of running MCMC directly")
    ;
}

Member Function Documentation

void MarkovChainMC::applyOptions ( const boost::program_options::variables_map &  vm )

virtual

Reimplemented from LimitAlgo.

Definition at line 100 of file MarkovChainMC.cc.

References dtNoiseDBValidation_cfg::cerr, chains_, FitP, mass_, mergeChains_, MultiGaussianP, noReset_, noSlimChain_, LimitAlgo::options_, proposalType_, proposalTypeName_, readChains_, runMinos_, saveChain_, TestP, UniformP, and updateHint_.

                                                                            {
    if      (proposalTypeName_ == "fit")     proposalType_ = FitP;
    else if (proposalTypeName_ == "uniform") proposalType_ = UniformP;
    else if (proposalTypeName_ == "gaus")    proposalType_ = MultiGaussianP;
    else if (proposalTypeName_ == "ortho")   proposalType_ = TestP;
    else if (proposalTypeName_ == "test")    proposalType_ = TestP;
    else {
        std::cerr << "MarkovChainMC: proposal type " << proposalTypeName_ << " not known." << "\n" << options_ << std::endl;
        throw std::invalid_argument("MarkovChainMC: unsupported proposal");
    }
        
    runMinos_ = vm.count("runMinos");
    noReset_  = vm.count("noReset");
    updateHint_  = vm.count("updateHint");

    mass_ = vm["mass"].as<float>();
    saveChain_   = vm.count("saveChain");
    noSlimChain_   = vm.count("noSlimChain");
    mergeChains_ = vm.count("mergeChains");
    readChains_  = vm.count("readChains");

    if (mergeChains_ && !saveChain_ && !readChains_) chains_.SetOwner(true);
}

int MarkovChainMC::guessBurnInSteps ( const RooStats::MarkovChain &  chain ) const

private

Definition at line 457 of file MarkovChainMC.cc.

References i, modelNDF_, n, mathSSE::sqrt(), and errorMatrix2Lands_multiChannel::start.

Referenced by limitFromChain(), mergeChains(), and runOnce().

{
    int n = chain.Size();
    std::vector<double> nll(n);
    for (int i = 0; i < n; ++i) {
        nll[i] = chain.NLL(i);
    }
    // get minimum of nll
    double minnll = nll[0];
    for (int i = 0; i < n; ++i) { if (nll[i] < minnll) minnll = nll[i]; }
    // subtract it from all
    for (int i = 0; i < n; ++i) nll[i] -= minnll;
    // the NLL looks like 0.5 * a chi2 with nparam degrees of freedom, plus an arbitrary constant
    // so it should have a sigma of 0.5*sqrt(2*nparams)
    // anything sensible should be between minimum and minimum + 10*sigma
    double maxcut = 5*sqrt(2.0*modelNDF_); 
    // now loop backwards until you find something that's outside
    int start = 0;
    for (start = n-1; start >= 0; --start) {
       if (nll[start] > maxcut) break;
    }
    return start;
}

void MarkovChainMC::limitAndError ( double &  limit, double &  limitErr, const std::vector< double > &  limits  ) const

private

Definition at line 298 of file MarkovChainMC.cc.

References abs, adaptiveTruncation_, gather_cfg::cout, end, QuantileCalculator::Jacknife, gen::k, MessageLogger_cff::limit, max(), min, QuantileCalculator::quantileAndError(), QuantileCalculator::Sectioning, QuantileCalculator::Simple, python.multivaluedict::sort(), mathSSE::sqrt(), errorMatrix2Lands_multiChannel::start, std::swap(), and truncatedMeanFraction_.

Referenced by run().

                                                                                                          {
  std::vector<double> limits(limitsIn);
  int num = limits.size();
  // possibly remove outliers before computing mean
  if (adaptiveTruncation_ && num >= 10) {
      std::sort(limits.begin(), limits.end());
      // determine location and size of the sample
      double median = (num % 2 ? limits[num/2] : 0.5*(limits[num/2] + limits[num/2+1]));
      double iqr = limits[3*num/4] - limits[num/4];
      // determine range of plausible values
      double min = median - iqr, max = median + iqr; 
      int start = 0, end = num-1; 
      while (start < end   && limits[start] < min) ++start;
      while (end   > start && limits[end]   > max) --end;
      num = end-start+1;
      // compute mean and rms of accepted part
      limit = 0; limitErr = 0;
      for (int k = start; k <= end; ++k) limit += limits[k];
      limit /= num;
      for (int k = start; k <= end; k++) limitErr += (limits[k]-limit)*(limits[k]-limit);
      limitErr = (num > 1 ? sqrt(limitErr/(num*(num-1))) : 0);
      std::cout << "Result from truncated mean: " << limit << " +/- " << limitErr << std::endl;
#if 0
      QuantileCalculator qc(limits);
      std::pair<double,double> qn = qc.quantileAndError(0.5, QuantileCalculator::Simple);
      std::pair<double,double> qs = qc.quantileAndError(0.5, QuantileCalculator::Sectioning);
      std::pair<double,double> qj = qc.quantileAndError(0.5, QuantileCalculator::Jacknife);
      std::cout << "Median of limits (simple):     " << qn.first << " +/- " << qn.second << std::endl;
      std::cout << "Median of limits (sectioning): " << qs.first << " +/- " << qs.second << std::endl;
      std::cout << "Median of limits (jacknife):   " << qj.first << " +/- " << qj.second << std::endl;
#endif
  } else {
      int noutl = floor(truncatedMeanFraction_ * num);
      if (noutl >= 1) { 
          std::sort(limits.begin(), limits.end());
          double median = (num % 2 ? limits[num/2] : 0.5*(limits[num/2] + limits[num/2+1]));
          for (int k = 0; k < noutl; ++k) {
              // make sure outliers are all at the end
              if (std::abs(limits[0]-median) > std::abs(limits[num-k-1]-median)) {
                  std::swap(limits[0], limits[num-k-1]);
              }
          }
          num -= noutl;
      }
      // compute mean and rms
      limit = 0; limitErr = 0;
      for (int k = 0; k < num; k++) limit += limits[k];
      limit /= num;
      for (int k = 0; k < num; k++) limitErr += (limits[k]-limit)*(limits[k]-limit);
      limitErr = (num > 1 ? sqrt(limitErr/(num*(num-1))) : 0);
  }
}

void MarkovChainMC::limitFromChain ( double &  limit, double &  limitErr, const RooArgSet &  poi, RooStats::MarkovChain &  chain, int  burnInSteps = -1  )

private

Definition at line 415 of file MarkovChainMC.cc.

References adaptiveBurnIn_, burnInFraction_, burnInSteps_, cl, guessBurnInSteps(), MergeJob_cfg::interval, mergeChains_, QuantileCalculator::quantileAndError(), and QuantileCalculator::Simple.

Referenced by readChains(), and run().

{
    if (burnInSteps < 0) {
        if (adaptiveBurnIn_) burnInSteps = guessBurnInSteps(chain);
        else burnInSteps = max<int>(burnInSteps_, chain.Size() * burnInFraction_);
    }
#if 1   // This is much faster and gives the same result
    QuantileCalculator qc(*chain.GetAsConstDataSet(), poi.first()->GetName(), burnInSteps);
    limit = qc.quantileAndError(cl, QuantileCalculator::Simple).first;
#else
    MCMCInterval interval("",poi,chain);
    RooArgList axes(poi);
    interval.SetConfidenceLevel(cl);
    interval.SetIntervalType(MCMCInterval::kTailFraction);
    interval.SetLeftSideTailFraction(0);
    interval.SetNumBurnInSteps(burnInSteps);
    interval.SetAxes(axes);
    limit = interval.UpperLimit((RooRealVar&)*poi.first());
    if (mergeChains_) {
        // must avoid that MCMCInterval deletes the chain
        interval.SetChain(*(RooStats::MarkovChain *)0);
    }
#endif
}

RooStats::MarkovChain * MarkovChainMC::mergeChains ( const RooArgSet &  poi, const std::vector< double > &  limits  ) const

private

Definition at line 351 of file MarkovChainMC.cc.

References adaptiveBurnIn_, burnInFraction_, burnInSteps_, chains_, gather_cfg::cout, guessBurnInSteps(), i, getHLTprescales::index, n, point, python.multivaluedict::sort(), and CommonMethods::weight().

Referenced by run().

{
    std::vector<double> limitsSorted(limits); std::sort(limitsSorted.begin(), limitsSorted.end());
    double lmin = limitsSorted.front(), lmax = limitsSorted.back();
    if (limitsSorted.size() > 5) {
        int n = limitsSorted.size();
        double lmedian = limitsSorted[n/2];
        lmin = lmedian - 2*(+lmedian - limitsSorted[1*n/4]);
        lmax = lmedian + 2*(-lmedian + limitsSorted[3*n/4]);
    }
    if (chains_.GetSize() == 0) throw std::runtime_error("No chains to merge");
    if (verbose > 1) std::cout << "Will merge " << chains_.GetSize() << " chains." << std::endl;
    RooArgSet pars(poi);
    RooStats::MarkovChain *merged = new RooStats::MarkovChain("Merged","",pars);
    TIter iter(&chains_);
    int index = 0;
    for (RooStats::MarkovChain *other = (RooStats::MarkovChain *) iter.Next();
         other != 0;
         other = (RooStats::MarkovChain *) iter.Next(), ++index) {
        if (limits[index] < lmin || limits[index] > lmax) continue;
        int burninSteps = adaptiveBurnIn_ ? guessBurnInSteps(*other) : max<int>(burnInSteps_, other->Size() * burnInFraction_);
        if (verbose > 1) std::cout << "Adding chain of " << other->Size() << " entries, skipping the first " <<  burninSteps << "; individual limit " << limits[index] << std::endl;
        for (int i = burninSteps, n = other->Size(); i < n; ++i) {
            RooArgSet point(*other->Get(i));
            double nllval = other->NLL();
            double weight = other->Weight();
            merged->Add(point,nllval,weight);
            if (verbose > 2 && (i % 500 == 0)) std::cout << "   added " << i << "/" << other->Size() << " entries." << std::endl;
        }
    }
    return merged;
}

virtual const std::string& MarkovChainMC::name ( void  ) const

inlinevirtual

Implements LimitAlgo.

Definition at line 21 of file MarkovChainMC.h.

Referenced by BeautifulSoup.Tag::_invert().

                                         {
    static const std::string name("MarkovChainMC");
    return name;
  }

void MarkovChainMC::readChains ( const RooArgSet &  poi, std::vector< double > &  limits  )

private

Definition at line 384 of file MarkovChainMC.cc.

References chains_, gather_cfg::cout, runtimedef::get(), gen::k, limitFromChain(), mass_, GetRecoTauVFromDQM_MC_cff::next, prof2calltree::prefix, readToysFromHere, and stationarityTest().

Referenced by run().

{
    double mylim, myerr;
    chains_.Clear();
    chains_.SetOwner(false);
    if (!readToysFromHere) throw std::logic_error("Cannot use readChains: option toysFile not specified, or input file empty");
    TDirectory *toyDir = readToysFromHere->GetDirectory("toys");
    if (!toyDir) throw std::logic_error("Cannot use readChains: empty toy dir in input file empty");
    TString prefix = TString::Format("MarkovChain_mh%g_",mass_);
    TIter next(toyDir->GetListOfKeys()); TKey *k;
    while ((k = (TKey *) next()) != 0) {
        if (TString(k->GetName()).Index(prefix) != 0) continue;
        RooStats::MarkovChain *toy = dynamic_cast<RooStats::MarkovChain *>(toyDir->Get(k->GetName()));
        if (toy == 0) continue;
        limitFromChain(mylim, myerr, poi, *toy);
        if (verbose > 1) std::cout << " limit " << mylim << " +/- " << myerr << std::endl;
        // vvvvv ---- begin convergence test, still being developed, not recommended yet.
        if (runtimedef::get("MCMC_STATIONARITY")) {
            if (!stationarityTest(*toy, poi, 30)) {
                if (verbose > 1) std::cout << " ---> rejecting chain!" << std::endl; 
                continue;
            }
        }
        // ^^^^^ ---- end of convergence test
        chains_.Add(toy);
        limits.push_back(mylim);
    }
    if (verbose) { std::cout << "Read " << chains_.GetSize() << " Markov Chains from input file." << std::endl; }
}

bool MarkovChainMC::run ( RooWorkspace *  w, RooStats::ModelConfig *  mc_s, RooStats::ModelConfig *  mc_b, RooAbsData &  data, double &  limit, double &  limitErr, const double *  hint  )

virtual

Implements LimitAlgo.

Definition at line 124 of file MarkovChainMC.cc.

References dtNoiseDBValidation_cfg::cerr, cl, CloseCoutSentry::clear(), gather_cfg::cout, i, dtDQMClient_cfg::INFO, iterations_, MessageLogger_cff::limit, limitAndError(), limitFromChain(), mergeChains(), mergeChains_, modelNDF_, MultiGaussianP, proposalType_, alignCSCRings::r, readChains(), readChains_, runOnce(), tries_, UniformP, updateHint_, dqm::qstatus::WARNING, and withSystematics.

                                                                                                                                                                    {
  if (proposalType_ == MultiGaussianP && !withSystematics) { 
      std::cerr << "Sorry, the multi-gaussian proposal does not work without systematics.\n" << 
                   "Uniform proposal will be used instead.\n" << std::endl;
      proposalType_ = UniformP;
  }

  RooFitGlobalKillSentry silence(verbose > 0 ? RooFit::INFO : RooFit::WARNING);

  CloseCoutSentry coutSentry(verbose <= 0); // close standard output and error, so that we don't flood them with minuit messages

  // Get degrees of freedom
  modelNDF_ = mc_s->GetParametersOfInterest()->getSize(); 
  if (withSystematics) modelNDF_ += mc_s->GetNuisanceParameters()->getSize();

  double suma = 0; int num = 0;
  double savhint = (hint ? *hint : -1); const double *thehint = hint;
  std::vector<double> limits;
  if (readChains_)  {
      readChains(*mc_s->GetParametersOfInterest(), limits);
  } else {
      for (unsigned int i = 0; i < tries_; ++i) {
          if (int nacc = runOnce(w,mc_s,mc_b,data,limit,limitErr,thehint)) {
              suma += nacc;
              if (verbose > 1) std::cout << "Limit from this run: " << limit << std::endl;
              limits.push_back(limit);
              if (updateHint_ && tries_ > 1 && limit > savhint) { 
                if (verbose > 0) std::cout << "Updating hint from " << savhint << " to " << limit << std::endl;
                savhint = limit; thehint = &savhint; 
              }
          }
      }
  } 
  num = limits.size();
  if (num == 0) return false;
  // average acceptance
  suma  = suma / (num * double(iterations_));
  limitAndError(limit, limitErr, limits);
  if (mergeChains_) {
    std::cout << "Limit from averaging:    " << limit << " +/- " << limitErr << std::endl;
    // copy constructors don't work, so we just have to leak memory :-(
    RooStats::MarkovChain *merged = mergeChains(*mc_s->GetParametersOfInterest(), limits);
    // set burn-in to zero, since steps have already been discarded when merging
    limitFromChain(limit, limitErr, *mc_s->GetParametersOfInterest(), *merged, 0);
    std::cout << "Limit from merged chain: " << limit << " +/- " << limitErr << std::endl;
  }
  coutSentry.clear();

  if (verbose >= 0) {
      std::cout << "\n -- MarkovChainMC -- " << "\n";
      RooRealVar *r = dynamic_cast<RooRealVar *>(mc_s->GetParametersOfInterest()->first());
      if (num > 1) {
          std::cout << "Limit: " << r->GetName() <<" < " << limit << " +/- " << limitErr << " @ " << cl * 100 << "% CL (" << num << " tries)" << std::endl;
          if (verbose > 0 && !readChains_) std::cout << "Average chain acceptance: " << suma << std::endl;
      } else {
          std::cout << "Limit: " << r->GetName() <<" < " << limit << " @ " << cl * 100 << "% CL" << std::endl;
      }
  }
  return true;
}

int MarkovChainMC::runOnce ( RooWorkspace *  w, RooStats::ModelConfig *  mc_s, RooStats::ModelConfig *  mc_b, RooAbsData &  data, double &  limit, double &  limitErr, const double *  hint  ) const

private

Definition at line 184 of file MarkovChainMC.cc.

References adaptiveBurnIn_, alwaysStepPoi_, burnInFraction_, burnInSteps_, dtNoiseDBValidation_cfg::cerr, chains_, cl, CloseCoutSentry::clear(), gather_cfg::cout, cropNSigmas_, data, debugProposal_, FitP, guessBurnInSteps(), hintSafetyFactor_, i, iterations_, mass_, max(), mergeChains_, min, MultiGaussianP, n, noReset_, proposalHelperUniformFraction_, proposalHelperWidthRangeDivisor_, proposalType_, alignCSCRings::r, runMinos_, saveChain_, slimChain(), TestP, UniformP, updateProposalParams_, v, withSystematics, and writeToysHere.

Referenced by run().

                                                                                                                                                                             {
  RooArgList poi(*mc_s->GetParametersOfInterest());
  RooRealVar *r = dynamic_cast<RooRealVar *>(poi.first());

  if ((hint != 0) && (*hint > r->getMin())) {
    r->setMax(hintSafetyFactor_*(*hint));
  }

  if (withSystematics && (mc_s->GetNuisanceParameters() == 0)) {
    throw std::logic_error("MarkovChainMC: running with systematics enabled, but nuisances not defined.");
  }
  
  w->loadSnapshot("clean");
  std::auto_ptr<RooFitResult> fit(0);
  if (proposalType_ == FitP || (cropNSigmas_ > 0)) {
      CloseCoutSentry coutSentry(verbose <= 1); // close standard output and error, so that we don't flood them with minuit messages
      fit.reset(mc_s->GetPdf()->fitTo(data, RooFit::Save(), RooFit::Minos(runMinos_)));
      coutSentry.clear();
      if (fit.get() == 0) { std::cerr << "Fit failed." << std::endl; return false; }
      if (verbose > 1) fit->Print("V");
      if (!noReset_) w->loadSnapshot("clean");
  }

  if (cropNSigmas_ > 0) {
      const RooArgList &fpf = fit->floatParsFinal();
      for (int i = 0, n = fpf.getSize(); i < n; ++i) {
          RooRealVar *fv = dynamic_cast<RooRealVar *>(fpf.at(i));
          if (std::string(r->GetName()) == fv->GetName()) continue;
          RooRealVar *v  = w->var(fv->GetName());
          double min = v->getMin(), max = v->getMax();
          if (fv->hasAsymError(false)) {
              min = (std::max(v->getMin(), fv->getVal() + cropNSigmas_ * fv->getAsymErrorLo()));
              max = (std::min(v->getMax(), fv->getVal() + cropNSigmas_ * fv->getAsymErrorHi()));
          } else if (fv->hasError(false)) {
              min = (std::max(v->getMin(), fv->getVal() - cropNSigmas_ * fv->getError()));
              max = (std::min(v->getMax(), fv->getVal() + cropNSigmas_ * fv->getError()));
          }
          if (verbose > 1) {
              std::cout << "  " << fv->GetName() << "[" << v->getMin() << ", " << v->getMax() << "] -> [" << min << ", " << max << "]" << std::endl;
          }
          v->setMin(min); v->setMax(max);
      }
  }

  std::auto_ptr<ProposalFunction> ownedPdfProp; 
  ProposalFunction* pdfProp = 0;
  ProposalHelper ph;
  switch (proposalType_) {
    case UniformP:  
        if (verbose) std::cout << "Using uniform proposal" << std::endl;
        ownedPdfProp.reset(new UniformProposal());
        pdfProp = ownedPdfProp.get();
        break;
    case FitP:
        if (verbose) std::cout << "Using fit proposal" << std::endl;
        ph.SetVariables(fit->floatParsFinal());
        ph.SetCovMatrix(fit->covarianceMatrix());
        pdfProp = ph.GetProposalFunction();
        break;
    case MultiGaussianP:
        if (verbose) std::cout << "Using multi-gaussian proposal" << std::endl;
        ph.SetVariables(*mc_s->GetNuisanceParameters());
        ph.SetWidthRangeDivisor(proposalHelperWidthRangeDivisor_);
        pdfProp = ph.GetProposalFunction();
        break;
    case TestP:
        ownedPdfProp.reset(new TestProposal(proposalHelperWidthRangeDivisor_, alwaysStepPoi_ ? poi : RooArgList()));
        pdfProp = ownedPdfProp.get();
        break;
  }
  if (proposalType_ != UniformP) {
      ph.SetUpdateProposalParameters(updateProposalParams_);
      if (proposalHelperUniformFraction_ > 0) ph.SetUniformFraction(proposalHelperUniformFraction_);
  }

  std::auto_ptr<DebugProposal> pdfDebugProp(debugProposal_ > 0 ? new DebugProposal(pdfProp, mc_s->GetPdf(), &data, debugProposal_) : 0);
  
  MCMCCalculator mc(data, *mc_s);
  mc.SetNumIters(iterations_); 
  mc.SetConfidenceLevel(cl);
  mc.SetNumBurnInSteps(burnInSteps_); 
  mc.SetProposalFunction(debugProposal_ > 0 ? *pdfDebugProp : *pdfProp);
  mc.SetLeftSideTailFraction(0);

  if (typeid(*mc_s->GetPriorPdf()) == typeid(RooUniform)) {
    mc.SetPriorPdf(*((RooAbsPdf *)0));
  }

  std::auto_ptr<MCMCInterval> mcInt;
  try {  
      mcInt.reset((MCMCInterval*)mc.GetInterval()); 
  } catch (std::length_error &ex) {
      mcInt.reset(0);
  }
  if (mcInt.get() == 0) return false;
  if (adaptiveBurnIn_) {
    mcInt->SetNumBurnInSteps(guessBurnInSteps(*mcInt->GetChain()));
  } else if (mcInt->GetChain()->Size() * burnInFraction_ > burnInSteps_) {
    mcInt->SetNumBurnInSteps(mcInt->GetChain()->Size() * burnInFraction_);
  }
  limit = mcInt->UpperLimit(*r);

  if (saveChain_ || mergeChains_) {
      // Copy-constructors don't work properly, so we just have to leak memory.
      //RooStats::MarkovChain *chain = new RooStats::MarkovChain(*mcInt->GetChain());
      RooStats::MarkovChain *chain = slimChain(*mc_s->GetParametersOfInterest(), *mcInt->GetChain());
      if (mergeChains_) chains_.Add(chain);
      if (saveChain_)  writeToysHere->WriteTObject(chain,  TString::Format("MarkovChain_mh%g_%u",mass_, RooRandom::integer(std::numeric_limits<UInt_t>::max() - 1)));
      return chain->Size();
  } else {
      return mcInt->GetChain()->Size();
  }
}

RooStats::MarkovChain * MarkovChainMC::slimChain ( const RooArgSet &  poi, const RooStats::MarkovChain &  chain  ) const

private

Definition at line 441 of file MarkovChainMC.cc.

References i, n, noSlimChain_, run_regression::ret, and CommonMethods::weight().

Referenced by runOnce().

{
    RooArgSet poilvalue(poi); // wtf they don't take a const & ??
    if (noSlimChain_) poilvalue.add(*chain.Get());
    RooStats::MarkovChain * ret = new RooStats::MarkovChain("","",poilvalue);
    for (int i = 0, n = chain.Size(); i < n; ++i) {
        RooArgSet entry(*chain.Get(i));
        Double_t nll = chain.NLL();
        Double_t weight = chain.Weight();
        if (i) ret->AddFast(entry, nll, weight);
        else   ret->Add(entry, nll, weight);
    }    
    return ret;
}

int MarkovChainMC::stationarityTest ( const RooStats::MarkovChain &  chain, const RooArgSet &  poi, int  nchunks  ) const

private

note: this is still being developed

Definition at line 482 of file MarkovChainMC.cc.

References trackerHits::c, data, python.tagInventory::entries, i, timingPdfMaker::mean, n, alignCSCRings::r, and python.multivaluedict::sort().

Referenced by readChains().

{
    std::vector<int>    entries(nchunks, 0);
    std::vector<double> mean(nchunks, .0);
    const RooDataSet *data = chain.GetAsConstDataSet();
    const RooRealVar *r = dynamic_cast<const RooRealVar *>(data->get()->find(poi.first()->GetName()));
    int  n = data->numEntries(), chunksize = ceil(n/double(nchunks));
    for (int i = 0, chunk = 0; i < n; i++) {
        data->get(i);
        if (i > 0 && i % chunksize == 0) chunk++;
        entries[chunk]++;
        mean[chunk] += r->getVal();
    }
    for (int c = 0; c < nchunks; ++c) { mean[c] /= entries[c]; }

    std::vector<double> dists, dists25;
    for (int c = 0; c < nchunks; ++c) {
        for (int c2 = 0; c2 < nchunks; ++c2) {
            if (c2 != c) dists.push_back(fabs(mean[c]-mean[c2]));
        }
        std::sort(dists.begin(), dists.end());
        dists25.push_back(dists[ceil(0.25*nchunks)]/mean[c]);
        dists.clear();
        //printf("chunk %3d: mean  %9.5f  dist25 %9.5f abs, %9.5f real\n", c, mean[c], mean[c]*dists25.back(), dists25.back());
    }
    std::sort(dists25.begin(), dists25.end());
    double tolerance = 10*dists25[ceil(0.25*nchunks)];
    //printf("Q(25) is %9.5f\n", dists25[ceil(0.25*nchunks)]);
    //printf("Q(50) is %9.5f\n", dists25[ceil(0.50*nchunks)]);
    //printf("Tolerance set to %9.5f\n", tolerance);
    bool converged = true;
    std::vector<int> friends(nchunks, 0), foes(nchunks, 0);
    for (int c = 0; c < nchunks; ++c) {
        for (int c2 = c+1; c2 < nchunks; ++c2) {
            if (c2 == c) continue;
            if (fabs(mean[c] - mean[c2]) < tolerance*mean[c]) {
                friends[c]++;
            } else {
                foes[c]++;
            }
        }   
        //printf("chunk %3d: mean  %9.5f  friends %3d  foes %3d \n", c, mean[c], friends[c], foes[c]);
        if (friends[c] >= 2 && foes[c] > 1) {
            converged = false;
        }
        //fflush(stdout);
    }
    return converged;
}

Member Data Documentation

bool MarkovChainMC::adaptiveBurnIn_ = false

staticprivate

Adaptive burn-in (experimental!)

Definition at line 37 of file MarkovChainMC.h.

Referenced by limitFromChain(), MarkovChainMC(), mergeChains(), and runOnce().

bool MarkovChainMC::adaptiveTruncation_ = true

staticprivate

do adaptive truncated mean

Definition at line 43 of file MarkovChainMC.h.

Referenced by limitAndError(), and MarkovChainMC().

bool MarkovChainMC::alwaysStepPoi_ = true

staticprivate

Definition at line 61 of file MarkovChainMC.h.

Referenced by MarkovChainMC(), and runOnce().

float MarkovChainMC::burnInFraction_ = 0.25

staticprivate

Discard these fraction of points.

Definition at line 35 of file MarkovChainMC.h.

Referenced by limitFromChain(), MarkovChainMC(), mergeChains(), and runOnce().

unsigned int MarkovChainMC::burnInSteps_ = 200

staticprivate

Discard these points.

Definition at line 33 of file MarkovChainMC.h.

Referenced by limitFromChain(), MarkovChainMC(), mergeChains(), and runOnce().

TList MarkovChainMC::chains_

mutableprivate

Definition at line 66 of file MarkovChainMC.h.

Referenced by applyOptions(), mergeChains(), readChains(), and runOnce().

float MarkovChainMC::cropNSigmas_ = 0

staticprivate

Definition at line 63 of file MarkovChainMC.h.

Referenced by MarkovChainMC(), and runOnce().

int MarkovChainMC::debugProposal_ = false

staticprivate

Definition at line 64 of file MarkovChainMC.h.

Referenced by MarkovChainMC(), and runOnce().

float MarkovChainMC::hintSafetyFactor_ = 5.

staticprivate

Safety factor for hint (integrate up to this number of times the hinted limit)

Definition at line 45 of file MarkovChainMC.h.

Referenced by MarkovChainMC(), and runOnce().

unsigned int MarkovChainMC::iterations_ = 10000

staticprivate

Propose this number of points for the chain.

Definition at line 31 of file MarkovChainMC.h.

Referenced by MarkovChainMC(), run(), and runOnce().

float MarkovChainMC::mass_

private

Mass of the Higgs boson (goes into the name of the saved chains)

Definition at line 55 of file MarkovChainMC.h.

Referenced by applyOptions(), readChains(), and runOnce().

bool MarkovChainMC::mergeChains_ = false

staticprivate

Merge chains instead of averaging limits.

Definition at line 51 of file MarkovChainMC.h.

Referenced by applyOptions(), limitFromChain(), run(), and runOnce().

int MarkovChainMC::modelNDF_

private

Number of degrees of freedom of the problem, approximately.

Definition at line 57 of file MarkovChainMC.h.

Referenced by guessBurnInSteps(), and run().

bool MarkovChainMC::noReset_ = false

staticprivate

Definition at line 29 of file MarkovChainMC.h.

Referenced by applyOptions(), and runOnce().

bool MarkovChainMC::noSlimChain_ = false

staticprivate

Leave all parameters in the markov chain, not just the POI.

Definition at line 49 of file MarkovChainMC.h.

Referenced by applyOptions(), and slimChain().

unsigned int MarkovChainMC::numberOfBins_

staticprivate

Definition at line 59 of file MarkovChainMC.h.

unsigned int MarkovChainMC::proposalHelperCacheSize_

staticprivate

Definition at line 60 of file MarkovChainMC.h.

float MarkovChainMC::proposalHelperUniformFraction_ = 0.0

staticprivate

Definition at line 62 of file MarkovChainMC.h.

Referenced by MarkovChainMC(), and runOnce().

float MarkovChainMC::proposalHelperWidthRangeDivisor_ = 5.

staticprivate

Definition at line 62 of file MarkovChainMC.h.

Referenced by MarkovChainMC(), and runOnce().

MarkovChainMC::ProposalType MarkovChainMC::proposalType_ = TestP

staticprivate

Definition at line 28 of file MarkovChainMC.h.

Referenced by applyOptions(), run(), and runOnce().

std::string MarkovChainMC::proposalTypeName_ = "ortho"

staticprivate

Definition at line 27 of file MarkovChainMC.h.

Referenced by applyOptions(), and MarkovChainMC().

bool MarkovChainMC::readChains_ = false

staticprivate

Read chains from file instead of running them.

Definition at line 53 of file MarkovChainMC.h.

Referenced by applyOptions(), and run().

bool MarkovChainMC::runMinos_ = false

staticprivate

Definition at line 29 of file MarkovChainMC.h.

Referenced by applyOptions(), and runOnce().

bool MarkovChainMC::saveChain_ = false

staticprivate

Save Markov Chain in output file.

Definition at line 47 of file MarkovChainMC.h.

Referenced by applyOptions(), and runOnce().

unsigned int MarkovChainMC::tries_ = 10

staticprivate

compute the limit N times

Definition at line 39 of file MarkovChainMC.h.

Referenced by MarkovChainMC(), and run().

float MarkovChainMC::truncatedMeanFraction_ = 0.0

staticprivate

Ignore up to this fraction of results if they're too far from the median.

Definition at line 41 of file MarkovChainMC.h.

Referenced by limitAndError(), and MarkovChainMC().

bool MarkovChainMC::updateHint_ = false

staticprivate

Definition at line 29 of file MarkovChainMC.h.

Referenced by applyOptions(), and run().

bool MarkovChainMC::updateProposalParams_ = false

staticprivate

Definition at line 29 of file MarkovChainMC.h.

Referenced by MarkovChainMC(), and runOnce().