#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)
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
Private Attributes
TList	chains_
float	mass_
	Mass of the Higgs boson (goes into the name of the saved chains)
int	modelNDF_
	Number of degrees of freedom of the problem, approximately.
static bool	noReset_ = false
static float	proposalHelperUniformFraction_ = 0.0
static bool	updateHint_ = false
static bool	updateProposalParams_ = false
Static Private Attributes
static bool	adaptiveBurnIn_ = false
	Adaptive burn-in (experimental!)
static bool	adaptiveTruncation_ = true
	do adaptive truncated mean
static bool	alwaysStepPoi_ = true
static float	burnInFraction_ = 0.25
	Discard these fraction of points.
static unsigned int	burnInSteps_ = 200
	Discard these points.
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)
static unsigned int	iterations_ = 10000
	Propose this number of points for the chain.
static bool	mergeChains_ = false
	Merge chains instead of averaging limits.
static bool	noSlimChain_ = false
	Leave all parameters in the markov chain, not just the POI.
static unsigned int	numberOfBins_
static unsigned int	proposalHelperCacheSize_
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.
static bool	runMinos_ = false
static bool	saveChain_ = false
	Save Markov Chain in output file.
static unsigned int	tries_ = 10
	compute the limit N times
static float	truncatedMeanFraction_ = 0.0
	Ignore up to this fraction of results if they're too far from the median.

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, 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 [inline, virtual]

Implements LimitAlgo.

Definition at line 21 of file MarkovChainMC.h.

                                         {
    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_, 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;
}