GsfTools/interface/GaussianSumUtilities1D.h

0001 #ifndef GaussianSumUtilities1D_h_
0002 #define GaussianSumUtilities1D_h_
0003
0004 // #include "TROOT.h"
0005
0006 #include "TrackingTools/GsfTools/interface/SingleGaussianState1D.h"
0007 #include "TrackingTools/GsfTools/interface/MultiGaussianState1D.h"
0008
0009 #include <vector>
0010
0011 /** Utility class for the analysis of one-dimensional Gaussian
0012  *  mixtures. The input state is assumed to exist for
0013  *  the lifetime of this object.
0014  */
0015
0016 class GaussianSumUtilities1D {
0017 private:
0018   enum ModeStatus { Valid, NotValid, NotComputed };
0019
0020 public:
0021   GaussianSumUtilities1D(const MultiGaussianState1D& state)
0022       : theState(state),
0023         //     theStates(state.components()),
0024         theModeStatus(NotComputed) {}
0025   ~GaussianSumUtilities1D() {}
0026
0027   /// number of components
0028   inline unsigned int size() const { return components().size(); }
0029   /// components
0030   inline const std::vector<SingleGaussianState1D>& components() const { return theState.components(); }
0031   /// weight of a component
0032   inline double weight(unsigned int i) const { return components()[i].weight(); }
0033   /// mean value of a component
0034   inline double mean(unsigned int i) const { return components()[i].mean(); }
0035   /// standard deviation of a component
0036   inline double standardDeviation(unsigned int i) const {
0037     //     return sqrt(components()[i].variance());
0038     return components()[i].standardDeviation();
0039   }
0040   /// variance of a component
0041   inline double variance(unsigned int i) const { return components()[i].variance(); }
0042   /// pdf of a single component at x
0043   double pdf(unsigned int i, double x) const;
0044   /// Quantile (i.e. x for a given value of the c.d.f.)
0045   double quantile(const double) const;
0046   /// mode status
0047   bool modeIsValid() const;
0048   /** Mode "state": mean = mode, variance = local variance at mode,
0049    *  weight chosen to have pdf(mode) equal to the one of the mixture */
0050   const SingleGaussianState1D& mode() const;
0051   /// value of the p.d.f.
0052   double pdf(double) const;
0053   /// value of the c.d.f.
0054   double cdf(const double&) const;
0055   /// first derivative of the p.d.f.
0056   double d1Pdf(const double&) const;
0057   /// second derivative of the p.d.f.
0058   double d2Pdf(const double&) const;
0059   /// third derivative of the p.d.f.
0060   double d3Pdf(const double&) const;
0061   /// ln(pdf)
0062   double lnPdf(const double&) const;
0063   /// first derivative of ln(pdf)
0064   double d1LnPdf(const double&) const;
0065   /// second derivative of ln(pdf)
0066   double d2LnPdf(const double&) const;
0067
0068   /// combined weight
0069   double weight() const { return theState.weight(); }
0070   /// combined mean
0071   double mean() const { return theState.mean(); }
0072   /// combined covariance
0073   double variance() const { return theState.variance(); }
0074
0075 private:
0076   /** Finds mode. Input: start value and typical scale.
0077    *  Output: mode and pdf(mode). Return value is true on success.
0078    */
0079   bool findMode(double& mode, double& pdfAtMode, const double& xStart, const double& scale) const;
0080   /// Value of gaussian distribution
0081   static double gauss(double, double, double);
0082   /// Integrated value of gaussian distribution
0083   static double gaussInt(double, double, double);
0084   /// Mean value of combined state
0085   double combinedMean() const;
0086   /// calculation of mode
0087   void computeMode() const;
0088   /** Local variance from Hessian matrix.
0089    *  Only valid if x corresponds to a (local) maximum! */
0090   double localVariance(double x) const;
0091
0092   // the state of the mode finder
0093   struct FinderState {
0094     FinderState() {}
0095     FinderState(size_t n) : pdfs(n) {}
0096     double x;
0097     double y;
0098     double yd;   // d1LnPdf
0099     double yd2;  // d2LnPdf
0100     std::vector<double> pdfs;
0101   };
0102
0103   // update tre state at x
0104   void update(FinderState& state, double x) const;
0105
0106   /// pdf components
0107   std::vector<double> pdfComponents(const double&) const;
0108   /// pdf components
0109   void pdfComponents(double, std::vector<double>&) const;
0110   /// value of the p.d.f. using the pdf components at the evaluation point
0111   static double pdf(double, const std::vector<double>&);
0112   /// first derivative of the p.d.f. using the pdf components at the evaluation point
0113   double d1Pdf(double, const std::vector<double>&) const;
0114   /// second derivative of the p.d.f. using the pdf components at the evaluation point
0115   double d2Pdf(double, const std::vector<double>&) const;
0116   /// third derivative of the p.d.f. using the pdf components at the evaluation point
0117   double d3Pdf(double, const std::vector<double>&) const;
0118   /// ln(pdf) using the pdf components at the evaluation point
0119   static double lnPdf(double, const std::vector<double>&);
0120   /// first derivative of ln(pdf) using the pdf components at the evaluation point
0121   double d1LnPdf(double, const std::vector<double>&) const;
0122   /// second derivative of ln(pdf) using the pdf components at the evaluation point
0123   double d2LnPdf(double, const std::vector<double>&) const;
0124
0125 private:
0126   const MultiGaussianState1D& theState;
0127   //   std::vector<SingleGaussianState1D> theStates;
0128
0129   mutable ModeStatus theModeStatus;
0130   mutable SingleGaussianState1D theMode;
0131   //   mutable double theMode;
0132 };
0133 #endif