npstat is hosted by Hepforge, IPPP Durham
NPStat  5.10.0
LocalQuantileRegression.hh
Go to the documentation of this file.
1 #ifndef NPSTAT_LOCALQUANTILEREGRESSION_HH_
2 #define NPSTAT_LOCALQUANTILEREGRESSION_HH_
3 
4 /*!
5 // \file LocalQuantileRegression.hh
6 //
7 // \brief Facilities for performing local linear and quadratic quantile
8 // regression.
9 //
10 // Gradient calculation is not supported.
11 //
12 // Classes declared in this header can be used with the
13 // top-level driver functions minuitQuantileRegression
14 // and minuitQuantileRegressionIncrBW declared in the
15 // minuitQuantileRegression.hh header (interfaces section).
16 //
17 // Author: I. Volobouev
18 //
19 // October 2011
20 */
21 
22 #include "npstat/nm/KDTree.hh"
23 
24 #include "npstat/stat/HistoND.hh"
27 
29 
30 namespace npstat {
31  /** Base class for quantile regression */
32  template <class Numeric>
34  {
35  public:
37  double quantile);
38  virtual ~QuantileRegressionBase() {}
39 
40  //@{
41  /** Inspect object properties */
42  inline unsigned dim() const {return mydim_;}
43  inline double quantile() const {return quantile_;}
44  inline const QuadraticOrthoPolyND& getPoly() const {return poly_;}
45  inline const std::vector<double>& lastCoeffs() const {return coeffs_;}
46  //@}
47 
48  /**
49  // Only the data inside the regression box will be used
50  // in the fit (and some points may be subsequently weighted
51  // down to zero by the weigh function).
52  //
53  // Note that the "setRegressionBox" method automatically
54  // sets up the standard mapping: the edges of the
55  // regression box are mapped into the boundaries which
56  // were used to create the orthogonal polynomials. If the
57  // standard mapping is not what you want (for example,
58  // this might be the case near the boundaries of the
59  // fitted region), you must call "setLinearMapping" after
60  // calling "setRegressionBox" in order to override the
61  // standard mapping.
62  //
63  // This method must be called at least once after the object
64  // is constructed and before any call to "linearLoss".
65  */
66  void setRegressionBox(const BoxND<Numeric>& box);
67 
68  /**
69  // This method sets up a mapping from the fitted region
70  // into the region which was used to create the orthogonal
71  // polynomials. The mapping is linear in every coordinate
72  // separately: ortho[i] = fitted[i]*scale[i] + location[i].
73  // Note that the "setRegressionBox" method automatically sets
74  // up the default mapping which will be adequate for most
75  // applications.
76  */
77  void setLinearMapping(const double* location, const double* scale,
78  unsigned locationAndScaleArraySize);
79 
80  /**
81  // This method will provide a guess for the 0th order
82  // coefficient of the regression polynomial. It will also
83  // calculate the effective number of points inside the current
84  // regression box.
85  */
86  void empiricalC0(double* c0, double* c0Uncertainty, double* npoints);
87 
88  /**
89  // This method calculates the function that has to be
90  // minimized; the minimization itself is up to some external
91  // optimization engine.
92  */
93  double linearLoss(const double* coeffs, unsigned nCoeffs);
94 
95  protected:
96  void resetAccumulators();
97 
98  const QuadraticOrthoPolyND& poly_;
99 
100  // Various vectors with size equal to the space dimensionality
101  BoxND<Numeric> regressionBox_;
102  std::vector<double> location_;
103  std::vector<double> scale_;
104  std::vector<double> coords_;
105 
106  // Various vectors with size equal to the number of orthogonal
107  // polynomials
108  std::vector<double> coeffs_;
109 
110  // Main accumulator for the loss function
111  long double loss_;
112 
113  // Stuff which does not change
114  const double quantile_;
115  const double onemq_;
116  const unsigned mydim_;
117 
118  private:
122 
123  virtual void makeStandardMapping();
124  virtual double empiricalQuantile(double* err, double* npoints) = 0;
125  virtual void actualLossCalculation() = 0;
126 
127  bool mappingIsDone_;
128  };
129 
130 
131  /** Quantile regression on data samples arranged into k-d trees */
132  template<class Point, class Numeric>
134  public AbsVisitor<Point, double>
135  {
136  public:
137  /**
138  // Constructor arguments are as follows:
139  //
140  // dataTree -- the tree of data points.
141  //
142  // regressedValue -- a functor that provides the observed
143  // value for the given input point. Typically,
144  // this value will be just one of the Point
145  // coordinates not used in k-d tree construction.
146  //
147  // poly -- the set of orthogonal polynomials used to
148  // construct the local regression surface.
149  //
150  // quantile -- the target quantile (between 0.0 and 1.0).
151  //
152  // This object will not own "dataTree", "regressedValue",
153  // or "poly" objects. These objects must still exist when the
154  // QuantileRegressionOnKDTree object is in use.
155  */
157  const KDTree<Point,Numeric>& dataTree,
158  const Functor1<double,Point>& regressedValue,
159  const QuadraticOrthoPolyND& poly,
160  double quantile);
161  virtual ~QuantileRegressionOnKDTree() {}
162 
163  /** Examine the data */
164  inline const KDTree<Point,Numeric>& getDataTree() const
165  {return dataTree_;}
166 
167  //@{
168  /** Method from AbsVisitor that we have to implement */
169  inline virtual void clear() {this->resetAccumulators();}
170  inline virtual double result()
171  {return static_cast<double>(this->loss_);}
172  virtual void process(const Point& value);
173  //@}
174 
175  protected:
176  const KDTree<Point,Numeric>& dataTree_;
177  const Functor1<double,Point>& regressedValue_;
178 
179  private:
180  virtual double empiricalQuantile(double* err, double* npoints);
181 
182  inline virtual void actualLossCalculation()
183  {dataTree_.visitInBox(*this, this->regressionBox_);}
184  };
185 
186 
187  /** Quantile regression on histograms */
188  template <typename Numeric>
190  public AbsArrayProjector<Numeric, double>
191  {
192  public:
193  /**
194  // Constructor arguments are as follows:
195  //
196  // histo -- Histogrammed data. It is assumed that the histogram
197  // is filled with actual event counts and that
198  // the quantity which is being regressed is used for
199  // the last histogram axis. The dimensionality of
200  // the histogram therefore must be larger by 1 than
201  // the dimensionality of the QuadraticOrthoPolyND object.
202  // Histogram overflow bins will be ignored.
203  //
204  // poly -- the set of orthogonal polynomial used to construct
205  // the local regression surface.
206  //
207  // quantile -- The quantile to determine (between 0.0 and 1.0).
208  //
209  // This object will not own "histo" or "poly" objects. These objects
210  // must still exist when this QuantileRegressionOnHisto object is
211  // in use.
212  */
214  const QuadraticOrthoPolyND& poly,
215  double quantile);
216  virtual ~QuantileRegressionOnHisto() {}
217 
218  /** Examine the data */
219  inline const HistoND<Numeric>& getHisto() const
220  {return histo_;}
221 
222  //@{
223  /** Method from AbsArrayProjector we have to implement */
224  inline virtual void clear() {this->resetAccumulators();}
225  inline virtual double result()
226  {return static_cast<double>(this->loss_);}
227  virtual void process(const unsigned *index, unsigned indexLen,
228  unsigned long linearIndex, const Numeric& value);
229  //@}
230 
231  protected:
232  const HistoAxis& lastAxis_;
233  const double halfbin_;
234  const Numeric zero_;
235  const HistoND<Numeric>& histo_;
236  BoxND<int> binBox_;
237 
238  private:
239  virtual void makeStandardMapping();
240  virtual double empiricalQuantile(double* err, double* npoints);
241 
242  inline virtual void actualLossCalculation()
243  {histo_.binContents().processSubrange(*this, binBox_);}
244  };
245 }
246 
247 #include "npstat/stat/LocalQuantileRegression.icc"
248 
249 #endif // NPSTAT_LOCALQUANTILEREGRESSION_HH_
Arbitrary-dimensional histogram template.
k-d tree template
Multivariate quadratic orthogonal polynomials with arbitrary weight functions.
Interface definitions and concrete simple functors for a variety of functor-based calculations.
Accumulates a sample of weighted points and calculates various descriptive statistics.
void processSubrange(AbsArrayProjector< Numeric, Num2 > &f, const BoxND< Integer > &subrange) const
Definition: HistoAxis.hh:31
Definition: HistoND.hh:46
const ArrayND< Numeric > & binContents() const
Definition: HistoND.hh:197
void visitInBox(AbsVisitor< Point, Result > &visitor, const BoxND< Numeric > &box) const
Definition: QuadraticOrthoPolyND.hh:24
Definition: LocalQuantileRegression.hh:34
unsigned dim() const
Definition: LocalQuantileRegression.hh:42
double linearLoss(const double *coeffs, unsigned nCoeffs)
void empiricalC0(double *c0, double *c0Uncertainty, double *npoints)
void setLinearMapping(const double *location, const double *scale, unsigned locationAndScaleArraySize)
void setRegressionBox(const BoxND< Numeric > &box)
Definition: LocalQuantileRegression.hh:191
virtual void process(const unsigned *index, unsigned indexLen, unsigned long linearIndex, const Numeric &value)
virtual double result()
Definition: LocalQuantileRegression.hh:225
const HistoND< Numeric > & getHisto() const
Definition: LocalQuantileRegression.hh:219
QuantileRegressionOnHisto(const HistoND< Numeric > &histo, const QuadraticOrthoPolyND &poly, double quantile)
virtual void clear()
Definition: LocalQuantileRegression.hh:224
Definition: LocalQuantileRegression.hh:135
QuantileRegressionOnKDTree(const KDTree< Point, Numeric > &dataTree, const Functor1< double, Point > &regressedValue, const QuadraticOrthoPolyND &poly, double quantile)
const KDTree< Point, Numeric > & getDataTree() const
Definition: LocalQuantileRegression.hh:164
virtual double result()
Definition: LocalQuantileRegression.hh:170
virtual void clear()
Definition: LocalQuantileRegression.hh:169
virtual void process(const Point &value)
Definition: AbsArrayProjector.hh:14
Definition: AbsArrayProjector.hh:21
Definition: AbsVisitor.hh:20
Definition: BoxND.hh:25