BayesianMarginal.h
Go to the documentation of this file.
1 #pragma once
2 
3 
4 #include <map>
5 
6 #include "TMVA/MethodKNN.h"
7 #include "TMVA/DataLoader.h"
8 
9 #include "CAFAna/Core/Binning.h"
10 #include "CAFAna/Fit/MCMCSample.h"
11 class TDirectory;
12 
13 namespace ana
14 {
15  class IFitVar;
16  class ISyst;
17  class MCMCSamples;
18 
19  /// Type to specify quantile values, so that we don't run into floating-point matching issues
20  enum class Quantile
21  {
25  k0pc, // the very lowest LL point
26  k90pc,
27  k95pc,
28  k100pc // the very highest LL point
29  };
30 
32  {
33  public:
34  IFitVarOrISyst(const ISyst* s ) : isyst(s), ifitvar(0) {}
35  IFitVarOrISyst(const IFitVar* v) : isyst(0), ifitvar(v) {}
36  const ISyst* isyst;
37  const IFitVar* ifitvar;
38  };
39 
40  /// The result of marginalizing over a set of Bayesian MCMC samples down to a few dimensions
42  {
43  public:
44  /// How to construct marginal distributions
45  enum class MarginalMode
46  {
47  kHistogram,
48  kLLWgtdHistogram,
49  kKNN,
50  };
51 
52  // See eg the statistics section of the PDG for the values. Note these are two-tailed versions...
53  static const std::vector<std::pair<Quantile, double>> QuantileUpVals;
54  BayesianMarginal(const MCMCSamples &samples,
55  const std::vector<IFitVarOrISyst> & orderedNames,
56  MarginalMode mode = MarginalMode::kHistogram);
57 
58  /// Retrieve the value above which the remaining points comprise the given quantile.
59  ///
60  /// \param quantile The quantile value desired
61  /// \param pdf If fMode==MarginalMode::kHistogram, the returned value
62  /// depends on the histogram that results from the marginalization (see on return val, below).
63  /// Thus the relevant histogram must be supplied for that mode.
64  /// For the other modes it may be left as nullptr.
65  /// \return Interpretation of return value depends on the MarginalMode \a fMode:
66  /// * for \a MarginalMode::kHistogram, it's the associated pdf value;
67  /// * for \a MarginalMode::kLLWgtdHistogram or MarginalMode::kkNN, it's the associated LL.
68  /// In either case it should be suitable for use with the derived class methods's for this instance.
69  double QuantileThreshold(Quantile quantile, const TH1 * pdf=nullptr) const;
70 
71  /// See other variant of QuantileThreshold() for discussion.
72  /// This version computes the threshold for an arbitrary quantile,
73  /// rather than the prestored ones, and so requires the original MCMCSamples.
74  double QuantileThreshold(double quantile, const MCMCSamples & samples, const TH1 * pdf=nullptr) const;
75 
76  std::vector<const ISyst*> Systs() const { return fSysts; }
77  std::vector<const IFitVar*> Vars() const { return fVars; }
78 
79  protected:
81 
82  double EstimateLLFromKNN(const std::vector<float> &vals) const;
83 
84  /// Determine the threshold for the 'histogram' mode from the given histogram with a brute force method
85  static double ThresholdFromTH1(const TH1* pdf, double quantile);
86 
87  /// Note that the signature here is different from usual.
88  /// This version expects to be passed a derived class instance
89  /// created by a derived LoadFrom() as its argument...
90  static void LoadInto(BayesianMarginal * marg, TDirectory * dir);
91  void SaveTo(TDirectory * dir) const;
92 
93  std::unique_ptr<TH1> ToHistogram(const std::vector<Binning> & bins) const;
94 
95  /// What mode are we using?
97 
98  /// Ref to the MCMCSamples. Will <b>not</b> be persisted.
100 
101  /// the variables/systs to be marginalized over, in axis order (x-axis, y-axis, ...)
102  std::vector<std::string> fOrderedBrNames;
103 
104  private:
105  void SetupKNN(TMVA::DataLoader * loader = nullptr);
106 
107  /// The samples that are the boundaries for various commonly-used quantiles
108  std::map<Quantile, MCMCSample> fQuantileSampleMap;
109 
110  /// the kNN used for approximating the surface from the points when in kNN mode
111  mutable std::unique_ptr<TMVA::MethodKNN> fkNN;
112  /// Needed for kNN
113  mutable std::unique_ptr<TMVA::DataLoader> fdataLoader;
114 
115  /// Variables we want to marginalize to (see also fSysts)
116  std::vector<const IFitVar*> fVars;
117  /// Systs we want to marginalize to (see also fVars)
118  std::vector<const ISyst*> fSysts;
119 
120 
121  static const std::string sKNNName;
122 
123  };
124 
125 }
std::vector< const ISyst * > fSysts
Systs we want to marginalize to (see also fVars)
std::map< Quantile, MCMCSample > fQuantileSampleMap
The samples that are the boundaries for various commonly-used quantiles.
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
IFitVarOrISyst(const ISyst *s)
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
std::unique_ptr< TMVA::MethodKNN > fkNN
the kNN used for approximating the surface from the points when in kNN mode
std::vector< const IFitVar * > Vars() const
std::vector< std::string > fOrderedBrNames
the variables/systs to be marginalized over, in axis order (x-axis, y-axis, ...)
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
static const std::vector< std::pair< Quantile, double > > QuantileUpVals
const XML_Char * s
Definition: expat.h:262
MarginalMode
How to construct marginal distributions.
std::unique_ptr< TMVA::DataLoader > fdataLoader
Needed for kNN.
IFitVarOrISyst(const IFitVar *v)
Quantile
Type to specify quantile values, so that we don&#39;t run into floating-point matching issues...
loader
Definition: demo0.py:10
std::vector< const ISyst * > Systs() const
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
MarginalMode fMode
What mode are we using?
const Binning bins
Definition: NumuCC_CPiBin.h:8
const MCMCSamples * fMCMCSamples
Ref to the MCMCSamples. Will not be persisted.
static const std::string sKNNName
The result of marginalizing over a set of Bayesian MCMC samples down to a few dimensions.
TDirectory * dir
Definition: macro.C:5
Interface definition for fittable variables.
Definition: IFitVar.h:16
std::vector< const IFitVar * > fVars
Variables we want to marginalize to (see also fSysts)
const IFitVar * ifitvar
enum BeamMode string