CMS 3D CMS Logo

Classes | Public Member Functions | Private Types | Private Member Functions | Static Private Member Functions | Private Attributes

GaussianSumUtilities1D Class Reference

#include <GaussianSumUtilities1D.h>

List of all members.

Classes

struct  FinderState

Public Member Functions

double cdf (const double &) const
 value of the c.d.f.
const std::vector
< SingleGaussianState1D > & 
components () const
 components
double d1LnPdf (const double &) const
 first derivative of ln(pdf)
double d1Pdf (const double &) const
 first derivative of the p.d.f.
double d2LnPdf (const double &) const
 second derivative of ln(pdf)
double d2Pdf (const double &) const
 second derivative of the p.d.f.
double d3Pdf (const double &) const
 third derivative of the p.d.f.
 GaussianSumUtilities1D (const MultiGaussianState1D &state)
double lnPdf (const double &) const
 ln(pdf)
double mean () const
 combined mean
double mean (unsigned int i) const
 mean value of a component
const SingleGaussianState1Dmode () const
bool modeIsValid () const
 mode status
double pdf (unsigned int i, double x) const
 pdf of a single component at x
double pdf (double) const
 value of the p.d.f.
double quantile (const double) const
 Quantile (i.e. x for a given value of the c.d.f.)
unsigned int size () const
 number of components
double standardDeviation (unsigned int i) const
 standard deviation of a component
double variance (unsigned int i) const
 variance of a component
double variance () const
 combined covariance
double weight () const
 combined weight
double weight (unsigned int i) const
 weight of a component
 ~GaussianSumUtilities1D ()

Private Types

enum  ModeStatus { Valid, NotValid, NotComputed }

Private Member Functions

double combinedMean () const
 Mean value of combined state.
void computeMode () const
 calculation of mode
double d1LnPdf (double, const std::vector< double > &) const
 first derivative of ln(pdf) using the pdf components at the evaluation point
double d1Pdf (double, const std::vector< double > &) const
 first derivative of the p.d.f. using the pdf components at the evaluation point
double d2LnPdf (double, const std::vector< double > &) const
 second derivative of ln(pdf) using the pdf components at the evaluation point
double d2Pdf (double, const std::vector< double > &) const
 second derivative of the p.d.f. using the pdf components at the evaluation point
double d3Pdf (double, const std::vector< double > &) const
 third derivative of the p.d.f. using the pdf components at the evaluation point
bool findMode (double &mode, double &pdfAtMode, const double &xStart, const double &scale) const
double localVariance (double x) const
void pdfComponents (double, std::vector< double > &) const
 pdf components
std::vector< double > pdfComponents (const double &) const
 pdf components
void update (FinderState &state, double x) const

Static Private Member Functions

static double gauss (double, double, double)
 Value of gaussian distribution.
static double gaussInt (double, double, double)
 Integrated value of gaussian distribution.
static double lnPdf (double, const std::vector< double > &)
 ln(pdf) using the pdf components at the evaluation point
static double pdf (double, const std::vector< double > &)
 value of the p.d.f. using the pdf components at the evaluation point

Private Attributes

SingleGaussianState1D theMode
ModeStatus theModeStatus
const MultiGaussianState1DtheState

Detailed Description

Utility class for the analysis of one-dimensional Gaussian mixtures. The input state is assumed to exist for the lifetime of this object.

Definition at line 16 of file GaussianSumUtilities1D.h.


Member Enumeration Documentation

Enumerator:
Valid 
NotValid 
NotComputed 

Definition at line 18 of file GaussianSumUtilities1D.h.


Constructor & Destructor Documentation

GaussianSumUtilities1D::GaussianSumUtilities1D ( const MultiGaussianState1D state) [inline]

Definition at line 21 of file GaussianSumUtilities1D.h.

                                                             :
    theState(state),
//     theStates(state.components()),
    theModeStatus(NotComputed) {} 
GaussianSumUtilities1D::~GaussianSumUtilities1D ( ) [inline]

Definition at line 25 of file GaussianSumUtilities1D.h.

{}

Member Function Documentation

double GaussianSumUtilities1D::cdf ( const double &  x) const

value of the c.d.f.

Definition at line 288 of file GaussianSumUtilities1D.cc.

References gaussInt(), i, mean(), query::result, asciidump::s, size(), standardDeviation(), and weight().

{
  double result(0.);
  size_t s=size();
  for ( unsigned int i=0; i<s; i++ )
    result += weight(i)*gaussInt(x,mean(i),standardDeviation(i));
  return result;
}
double GaussianSumUtilities1D::combinedMean ( ) const [private]

Mean value of combined state.

Definition at line 463 of file GaussianSumUtilities1D.cc.

References components(), i, asciidump::s, size(), and weight().

{
  double s0(0.);
  double s1(0.);
  int s = size();
  SingleGaussianState1D const * __restrict__ sgv = &components().front();
  for (int i=0; i<s; i++ )
    s0 += sgv[i].weight();
  for (int i=0; i<s; i++ )
    s1 += sgv[i].weight()*sgv[i].mean();
  return s1/s0;
}
const std::vector<SingleGaussianState1D>& GaussianSumUtilities1D::components ( ) const [inline]
void GaussianSumUtilities1D::computeMode ( ) const [private]

calculation of mode

Definition at line 101 of file GaussianSumUtilities1D.cc.

References components(), findMode(), first, i, localVariance(), M_PI, MultiGaussianState1D::mean(), mean(), mode(), NotValid, pdf(), asciidump::s, edm::second(), size(), mathSSE::sqrt(), standardDeviation(), theMode, theModeStatus, theState, Valid, TrackValidation_HighPurity_cff::valid, MultiGaussianState1D::variance(), variance(), weight(), MultiGaussianState1D::weight(), x, and detailsBasic3DVector::y.

Referenced by mode(), and modeIsValid().

{
//   TimerStack tstack;
//   tstack.benchmark("GSU1D::benchmark",100000);
//   FastTimerStackPush(tstack,"GaussianSumUtilities1D::computeMode");
  theModeStatus = NotValid;
  //
  // check for "degenerate" states (identical values with vanishing error)
  //
  if ( theState.variance()<1.e-14 ) {
    theMode = SingleGaussianState1D(theState.mean(),theState.variance(),
                                    theState.weight());
    theModeStatus = Valid;
    return;
  }
  //
  // Use means of individual components as start values.
  // Sort by value of pdf.
  //
  typedef std::multimap<double, int, std::greater<double> > StartMap;
  StartMap xStart;
  for ( unsigned int i=0; i<size(); i++ ) {
    xStart.insert(std::make_pair(pdf(mean(i)),i));
  }
  //
  // Now try with each start value
  //
  int iRes(-1);     // index of start component for current estimate
  double xRes(mean((*xStart.begin()).second)); // current estimate of mode
  double yRes(-1.); // pdf at current estimate of mode
//   std::pair<double,double> result(-1.,mean((*xStart.begin()).second));
  for ( StartMap::const_iterator i=xStart.begin(); i!=xStart.end(); i++ ) {
    //
    // Convergence radius for a single Gaussian = 1 sigma: don't try
    // start values within 1 sigma of the current solution
    //
    if ( theModeStatus==Valid &&
         fabs(mean((*i).second)-mean(iRes))/standardDeviation(iRes)<1. )  continue;
    //
    // If a solution exists: drop as soon as the pdf at
    // start value drops to < 75% of maximum (try to avoid
    // unnecessary searches for the mode)
    //
    if ( theModeStatus==Valid && 
         (*i).first/(*xStart.begin()).first<0.75 )  break;
    //
    // Try to find mode
    //
    double x;
    double y;
    bool valid = findMode(x,y,mean((*i).second),standardDeviation((*i).second));
    //
    // consider only successful searches
    //
    if ( valid ) { //...
      //
      // update result only for significant changes in pdf(solution)
      //
      if ( yRes<0. || (y-yRes)/(y+yRes)>1.e-10 ) {
      iRes = (*i).second;               // store index
      theModeStatus = Valid;            // update status
      xRes = x;                         // current solution
      yRes = y;                         // and its pdf value
//       result = std::make_pair(y,x);     // store solution and pdf(solution)
      }
    } //...
  } 
  //
  // check (existance of) solution
  //
  if ( theModeStatus== Valid ) {
    //
    // Construct single Gaussian state with 
    //  mean = mode
    //  variance = local variance at mode
    //  weight such that the pdf's of the mixture and the
    //    single state are equal at the mode
    //
    double mode = xRes;
    double varMode = localVariance(mode);
    double wgtMode = pdf(mode)*sqrt(2*M_PI*varMode);
    theMode = SingleGaussianState1D(mode,varMode,wgtMode);
  }
  else {
    //
    // mode finding failed: set solution to highest component
    //  (alternative would be global mean / variance ..?)
    //
    //two many log warnings to actually be useful - comment out
    //edm::LogWarning("GaussianSumUtilities") << "1D mode calculation failed";
//     double x = mean();
//     double y = pdf(x);
//     result = std::make_pair(y,x);
//     theMode = SingleGaussianState1D(mean(),variance(),weight());
    //
    // look for component with highest value at mean
    //
    int icMax(0);
    double ySqMax(0.);
    int s = size();
    for (int ic=0; ic<s; ++ic ) {
      double w = weight(ic);
      double ySq = w*w/variance(ic);
      if ( ySq>ySqMax ) {
        icMax = ic;
        ySqMax = ySq;
      }
    }
    theMode = SingleGaussianState1D(components()[icMax]);
  }
  
}
double GaussianSumUtilities1D::d1LnPdf ( const double &  x) const

first derivative of ln(pdf)

Definition at line 322 of file GaussianSumUtilities1D.cc.

References pdfComponents().

Referenced by d2LnPdf().

{
  return d1LnPdf(x,pdfComponents(x));
}
double GaussianSumUtilities1D::d1LnPdf ( double  x,
const std::vector< double > &  pdfs 
) const [private]

first derivative of ln(pdf) using the pdf components at the evaluation point

Definition at line 422 of file GaussianSumUtilities1D.cc.

References d1Pdf(), f, min, pdf(), and query::result.

{

  double f = pdf(x,pdfs);
  double result(d1Pdf(x,pdfs));
  if ( f>std::numeric_limits<double>::min() )  result /= f;
  else  result = 0.;
  return result;
}
double GaussianSumUtilities1D::d1Pdf ( double  x,
const std::vector< double > &  pdfs 
) const [private]

first derivative of the p.d.f. using the pdf components at the evaluation point

Definition at line 376 of file GaussianSumUtilities1D.cc.

References i, mean(), query::result, asciidump::s, size(), and standardDeviation().

{
  double result(0.);
  size_t s=size();
  for ( unsigned int i=0; i<s; i++ ) {
    double dx = (x-mean(i))/standardDeviation(i);
    result += -pdfs[i]*dx/standardDeviation(i);
  }
  return result;
}
double GaussianSumUtilities1D::d1Pdf ( const double &  x) const

first derivative of the p.d.f.

Definition at line 298 of file GaussianSumUtilities1D.cc.

References pdfComponents().

Referenced by d1LnPdf(), and update().

{
  return d1Pdf(x,pdfComponents(x));
}
double GaussianSumUtilities1D::d2LnPdf ( const double &  x) const

second derivative of ln(pdf)

Definition at line 328 of file GaussianSumUtilities1D.cc.

References pdfComponents().

{
  return d2LnPdf(x,pdfComponents(x));
}
double GaussianSumUtilities1D::d2LnPdf ( double  x,
const std::vector< double > &  pdfs 
) const [private]

second derivative of ln(pdf) using the pdf components at the evaluation point

Definition at line 433 of file GaussianSumUtilities1D.cc.

References d1LnPdf(), d2Pdf(), f, min, pdf(), and query::result.

{

  double f = pdf(x,pdfs);
  double df = d1LnPdf(x,pdfs);
  double result(-df*df);
  if ( f>std::numeric_limits<double>::min() )  result += d2Pdf(x,pdfs)/f;
  return result;
}
double GaussianSumUtilities1D::d2Pdf ( double  x,
const std::vector< double > &  pdfs 
) const [private]

second derivative of the p.d.f. using the pdf components at the evaluation point

Definition at line 388 of file GaussianSumUtilities1D.cc.

References i, mean(), query::result, asciidump::s, size(), and standardDeviation().

{
  double result(0.);
  size_t s=size();
  for ( unsigned int i=0; i<s; i++ ) {
    double dx = (x-mean(i))/standardDeviation(i);
    result += pdfs[i]/standardDeviation(i)/standardDeviation(i)*(dx*dx-1);
  }
  return result;
}
double GaussianSumUtilities1D::d2Pdf ( const double &  x) const

second derivative of the p.d.f.

Definition at line 304 of file GaussianSumUtilities1D.cc.

References pdfComponents().

Referenced by d2LnPdf(), localVariance(), and update().

{
  return d2Pdf(x,pdfComponents(x));
}
double GaussianSumUtilities1D::d3Pdf ( double  x,
const std::vector< double > &  pdfs 
) const [private]

third derivative of the p.d.f. using the pdf components at the evaluation point

Definition at line 400 of file GaussianSumUtilities1D.cc.

References i, mean(), query::result, asciidump::s, size(), and standardDeviation().

{
  double result(0.);
  size_t s=size();
  for ( unsigned int i=0; i<s; i++ ) {
    double dx = (x-mean(i))/standardDeviation(i);
    result += pdfs[i]/standardDeviation(i)/standardDeviation(i)/standardDeviation(i)*
      (-dx*dx+3)*dx;
  }
  return result;
}
double GaussianSumUtilities1D::d3Pdf ( const double &  x) const

third derivative of the p.d.f.

Definition at line 310 of file GaussianSumUtilities1D.cc.

References pdfComponents().

{
  return d3Pdf(x,pdfComponents(x));
}
bool GaussianSumUtilities1D::findMode ( double &  mode,
double &  pdfAtMode,
const double &  xStart,
const double &  scale 
) const [private]

Finds mode. Input: start value and typical scale. Output: mode and pdf(mode). Return value is true on success.

Definition at line 215 of file GaussianSumUtilities1D.cc.

References min, query::result, size(), evf::utils::state, update(), GaussianSumUtilities1D::FinderState::x, GaussianSumUtilities1D::FinderState::y, GaussianSumUtilities1D::FinderState::yd, and GaussianSumUtilities1D::FinderState::yd2.

Referenced by computeMode().

{
  //
  // try with Newton on (lnPdf)'(x)
  //
  double x1(0.);
  double y1(0.);
  FinderState  state(size());
  update(state,xStart);

  double xmin(xStart-1.*scale);
  double xmax(xStart+1.*scale);
  //
  // preset result
  //
  bool result(false);
  xMode = state.x;
  yMode = state.y;
  //
  // Iterate
  //
  int nLoop(0);
  while ( nLoop++<20 ) {
    if ( nLoop>1 && state.yd2<0. &&  
         ( fabs(state.yd*scale)<1.e-10 || fabs(state.y-y1)/(state.y+y1)<1.e-14 ) ) {
      result = true;
      break;
    }
    if ( fabs(state.yd2)<std::numeric_limits<float>::min() )  
      state.yd2 = state.yd2>0. ? std::numeric_limits<float>::min() : -std::numeric_limits<float>::min();
    double dx = -state.yd/state.yd2;
    x1 = state.x;
    y1 = state.y;
    double x2 = x1 + dx;
    if ( state.yd2>0. && (x2<xmin||x2>xmax) )  return false;
    update(state,x2);
  }
  //
  // result
  //
  if ( result ) {
    xMode = state.x;
    yMode = state.y;
  }
  return result;
}
double GaussianSumUtilities1D::gauss ( double  x,
double  mean,
double  sigma 
) [static, private]

Value of gaussian distribution.

Definition at line 444 of file GaussianSumUtilities1D.cc.

Referenced by pdf(), and pdfComponents().

{
//   const double fNorm(1./sqrt(2*M_PI));
//   double result(0.);

//   double d((x-mean)/sigma);
//   if ( fabs(d)<20. )  result = exp(-d*d/2.);
//   result *= fNorm/sigma;
//   return result;
  return ROOT::Math::gaussian_pdf(x,sigma,mean);
}
double GaussianSumUtilities1D::gaussInt ( double  x,
double  mean,
double  sigma 
) [static, private]

Integrated value of gaussian distribution.

Definition at line 457 of file GaussianSumUtilities1D.cc.

Referenced by cdf().

{
  return ROOT::Math::normal_cdf(x,sigma,mean);
}
double GaussianSumUtilities1D::lnPdf ( double  x,
const std::vector< double > &  pdfs 
) [static, private]

ln(pdf) using the pdf components at the evaluation point

Definition at line 413 of file GaussianSumUtilities1D.cc.

References f, funct::log(), max(), min, pdf(), and query::result.

{
  double f(pdf(x,pdfs));
  double result(-std::numeric_limits<float>::max());
  if ( f>std::numeric_limits<double>::min() )  result = log(f);
  return result;
}
double GaussianSumUtilities1D::lnPdf ( const double &  x) const

ln(pdf)

Definition at line 316 of file GaussianSumUtilities1D.cc.

References pdfComponents().

{
  return lnPdf(x,pdfComponents(x));
}
double GaussianSumUtilities1D::localVariance ( double  x) const [private]

Local variance from Hessian matrix. Only valid if x corresponds to a (local) maximum!

Definition at line 477 of file GaussianSumUtilities1D.cc.

References d2Pdf(), pdf(), pdfComponents(), and query::result.

Referenced by computeMode().

{
  std::vector<double> pdfs;
  pdfComponents(x,pdfs);
  double result = -pdf(x,pdfs)/d2Pdf(x,pdfs);
  // FIXME: wrong curvature seems to be non-existant but should add a proper recovery
  if ( result<0. )
    edm::LogWarning("GaussianSumUtilities") << "1D variance at mode < 0";    
  return result;
}
double GaussianSumUtilities1D::mean ( unsigned int  i) const [inline]
double GaussianSumUtilities1D::mean ( ) const [inline]

combined mean

Definition at line 75 of file GaussianSumUtilities1D.h.

References MultiGaussianState1D::mean(), and theState.

Referenced by cdf(), computeMode(), d1Pdf(), d2Pdf(), d3Pdf(), pdf(), and pdfComponents().

                       {
    return theState.mean();
  }
const SingleGaussianState1D & GaussianSumUtilities1D::mode ( void  ) const
bool GaussianSumUtilities1D::modeIsValid ( ) const
double GaussianSumUtilities1D::pdf ( unsigned int  i,
double  x 
) const

pdf of a single component at x

Definition at line 22 of file GaussianSumUtilities1D.cc.

References gauss(), mean(), standardDeviation(), and weight().

Referenced by computeMode(), d1LnPdf(), d2LnPdf(), lnPdf(), localVariance(), pdf(), and update().

                                                                  {
  return weight(i)*gauss(x,mean(i),standardDeviation(i));
}
double GaussianSumUtilities1D::pdf ( double  x) const

value of the p.d.f.

Definition at line 278 of file GaussianSumUtilities1D.cc.

References i, pdf(), query::result, asciidump::s, and size().

{
  double result(0.);
  size_t s=size();
  for ( unsigned int i=0; i<s; i++ )
    result += pdf(i,x);
  return result;
}
double GaussianSumUtilities1D::pdf ( double  ,
const std::vector< double > &  pdfs 
) [static, private]

value of the p.d.f. using the pdf components at the evaluation point

Definition at line 370 of file GaussianSumUtilities1D.cc.

{
  return std::accumulate(pdfs.begin(),pdfs.end(),0.);
}
void GaussianSumUtilities1D::pdfComponents ( double  x,
std::vector< double > &  result 
) const [private]

pdf components

Definition at line 352 of file GaussianSumUtilities1D.cc.

References begin, components(), end, asciidump::s, and size().

                                                                                      {
  size_t s = size();
  if (s!=result.size()) result.resize(s);  
  std::transform(components().begin(),components().end(),result.begin(),PDF(x));
}
std::vector< double > GaussianSumUtilities1D::pdfComponents ( const double &  x) const [private]

pdf components

Definition at line 334 of file GaussianSumUtilities1D.cc.

References gauss(), i, mean(), query::result, size(), standardDeviation(), and weight().

Referenced by d1LnPdf(), d1Pdf(), d2LnPdf(), d2Pdf(), d3Pdf(), lnPdf(), localVariance(), and update().

{
  std::vector<double> result;
  result.reserve(size());
  for ( unsigned int i=0; i<size(); i++ )
    result.push_back(weight(i)*gauss(x,mean(i),standardDeviation(i)));
  return result;
}
double GaussianSumUtilities1D::quantile ( const double  q) const

Quantile (i.e. x for a given value of the c.d.f.)

Definition at line 28 of file GaussianSumUtilities1D.cc.

{
  return ROOT::Math::gaussian_quantile(q,1.);
}
unsigned int GaussianSumUtilities1D::size ( void  ) const [inline]

number of components

Definition at line 28 of file GaussianSumUtilities1D.h.

References components().

Referenced by cdf(), combinedMean(), computeMode(), d1Pdf(), d2Pdf(), d3Pdf(), findMode(), pdf(), and pdfComponents().

{return components().size();}
double GaussianSumUtilities1D::standardDeviation ( unsigned int  i) const [inline]

standard deviation of a component

Definition at line 38 of file GaussianSumUtilities1D.h.

References components(), and i.

Referenced by cdf(), computeMode(), d1Pdf(), d2Pdf(), d3Pdf(), pdf(), and pdfComponents().

                                                         {
//     return sqrt(components()[i].variance());
    return components()[i].standardDeviation();
  }
void GaussianSumUtilities1D::update ( FinderState state,
double  x 
) const [private]
double GaussianSumUtilities1D::variance ( unsigned int  i) const [inline]

variance of a component

Definition at line 43 of file GaussianSumUtilities1D.h.

References components(), and i.

Referenced by GsfTrackProducerBase::fillMode().

{return components()[i].variance();}
double GaussianSumUtilities1D::variance ( ) const [inline]

combined covariance

Definition at line 79 of file GaussianSumUtilities1D.h.

References theState, and MultiGaussianState1D::variance().

Referenced by computeMode().

                           {
    return theState.variance();
  }
double GaussianSumUtilities1D::weight ( ) const [inline]

combined weight

Definition at line 71 of file GaussianSumUtilities1D.h.

References theState, and MultiGaussianState1D::weight().

Referenced by cdf(), combinedMean(), computeMode(), pdf(), and pdfComponents().

                         {
    return theState.weight();
  }
double GaussianSumUtilities1D::weight ( unsigned int  i) const [inline]

weight of a component

Definition at line 34 of file GaussianSumUtilities1D.h.

References components(), and i.

{return components()[i].weight();}

Member Data Documentation

Definition at line 139 of file GaussianSumUtilities1D.h.

Referenced by computeMode(), and mode().

Definition at line 138 of file GaussianSumUtilities1D.h.

Referenced by computeMode(), mode(), and modeIsValid().

Definition at line 135 of file GaussianSumUtilities1D.h.

Referenced by components(), computeMode(), mean(), variance(), and weight().