#include <GSUtilities.h>

Public Member Functions
double	cdf (const double &) const
	value of integral(pdf) More...

double	combinedMean () const
	mean value of combined state More...

double	dpdf1 (const double &) const
	first derivative of pdf More...

double	dpdf2 (const double &) const
	second derivative of pdf More...

double	errorCombinedMean () const

float	errorHighestWeight () const

float	errorMode ()

float	getMax (float)

float	getMin (float)

	GSUtilities (const unsigned nComp, const float weights, const float parameters, const float *errors)
	constructor from arrays of weights, parameters and standard deviations More...

float	maxWeight () const

float	mode () const
	mode More...

double	pdf (const double &) const
	value of the pdf More...

float	quantile (const float) const

	~GSUtilities ()

Private Member Functions
double	findMode (const double) const

double	gauss (const double &, const double &, const double &) const
	value of gaussian distribution More...

double	gaussInt (const double &, const double &, const double &) const
	integrated value of gaussian distribution More...

Private Attributes
float *	theErrors

unsigned	theNComp

float *	theParameters

float *	theWeights

Detailed Description

Some utilities for analysing 1D Gaussian mixtures. Copied from ORCA's EgammaGSUtilities.

Definition at line 8 of file GSUtilities.h.

Constructor & Destructor Documentation

GSUtilities::GSUtilities	(	const unsigned	nComp,
		const float *	weights,
		const float *	parameters,
		const float *	errors
	)

inline

constructor from arrays of weights, parameters and standard deviations

Definition at line 11 of file GSUtilities.h.

References mps_fire::i, theErrors, theNComp, theParameters, and theWeights.

                                                          :
     theNComp(nComp),
     theWeights(nullptr),
     theParameters(nullptr),
     theErrors(nullptr)
   {
     if ( theNComp ) {
       theWeights = new float[theNComp];
       theParameters = new float[theNComp];
       theErrors = new float[theNComp];
     }
     const float* wPtr1(weights);
     const float* pPtr1(parameters);
     const float* ePtr1(errors);
     float* wPtr2(theWeights);
     float* pPtr2(theParameters);
     float* ePtr2(theErrors);
     for ( unsigned i=0; i<theNComp; i++ ) {
       *(wPtr2++) = weights ? *(wPtr1++) : 1.;
       *(pPtr2++) = *(pPtr1++);
       *(ePtr2++) = *(ePtr1++);
     }
   } 

GSUtilities::~GSUtilities ( )

inline

Definition at line 35 of file GSUtilities.h.

References cdf(), combinedMean(), dpdf1(), dpdf2(), errorCombinedMean(), errorHighestWeight(), errorMode(), findMode(), gauss(), gaussInt(), getMax(), getMin(), maxWeight(), mode(), pdf(), quantile(), theErrors, theParameters, and theWeights.

   {
     delete [] theWeights;
     delete [] theParameters;
     delete [] theErrors;
   }

Member Function Documentation

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

value of integral(pdf)

Definition at line 157 of file GSUtilities.cc.

References gaussInt(), mps_fire::i, mps_fire::result, theErrors, theNComp, theParameters, and theWeights.

Referenced by errorMode(), getMax(), getMin(), quantile(), and ~GSUtilities().

 {
   double result(0.);
   for ( unsigned i=0; i<theNComp; i++ )
     result += theWeights[i]*gaussInt(x,theParameters[i],theErrors[i]);
   return result;
 }

double GSUtilities::combinedMean ( ) const

mean value of combined state

Definition at line 210 of file GSUtilities.cc.

References mps_fire::i, theNComp, theParameters, and theWeights.

Referenced by ~GSUtilities().

 {
   double s0(0.);
   double s1(0.);
   for ( unsigned i=0; i<theNComp; i++ ) {
     s0 += theWeights[i];
     s1 += theWeights[i]*theParameters[i];
   }
   return s1/s0;
 }

double GSUtilities::dpdf1 ( const double & x ) const

first derivative of pdf

Definition at line 166 of file GSUtilities.cc.

References PVValHelper::dx, gauss(), mps_fire::i, mps_fire::result, theErrors, theNComp, theParameters, and theWeights.

Referenced by findMode(), and ~GSUtilities().

 {
   double result(0.);
   for ( unsigned i=0; i<theNComp; i++ ) {
     double dx = (x-theParameters[i])/theErrors[i];
     result += -theWeights[i]*dx/theErrors[i]*
       gauss(x,theParameters[i],theErrors[i]);
   }
   return result;
 }

double GSUtilities::dpdf2 ( const double & x ) const

second derivative of pdf

Definition at line 178 of file GSUtilities.cc.

References PVValHelper::dx, gauss(), mps_fire::i, mps_fire::result, theErrors, theNComp, theParameters, and theWeights.

Referenced by findMode(), and ~GSUtilities().

 {
   double result(0.);
   for ( unsigned i=0; i<theNComp; i++ ) {
     double dx = (x-theParameters[i])/theErrors[i];
     result += theWeights[i]/theErrors[i]/theErrors[i]*
       (dx*dx-1)*gauss(x,theParameters[i],theErrors[i]);
   }
   return result;
 }

double GSUtilities::errorCombinedMean ( ) const

Definition at line 222 of file GSUtilities.cc.

References mps_fire::i, mathSSE::sqrt(), theNComp, and theWeights.

Referenced by ~GSUtilities().

 {
   double s0(0.);
   for ( unsigned i=0; i<theNComp; i++ ) {
     s0 += theWeights[i];
   }
   return 1./(sqrt(s0));
 }

float GSUtilities::errorHighestWeight ( ) const

Definition at line 73 of file GSUtilities.cc.

References mps_fire::i, theErrors, theNComp, and theWeights.

Referenced by ~GSUtilities().

 {
     int iwMax(-1);
     float wMax(0.);
     for ( unsigned i=0; i<theNComp; i++ ) {
         if ( theWeights[i]>wMax ) {
             iwMax = i;
             wMax = theWeights[i];
         }
     }
     return theErrors[iwMax];
 }

float GSUtilities::errorMode ( )

Definition at line 249 of file GSUtilities.cc.

References cdf(), PVValHelper::dx, getMax(), getMin(), Exhume::I, SiStripPI::max, min(), mod(), mode(), pdf(), globals_cff::x1, and globals_cff::x2.

Referenced by ~GSUtilities().

 {
     float mod = mode();
     float min = getMin(mod);
     float max = getMax(mod);
     int nBins = 1000;
     float dx = (max - min )/(float)nBins;
 
     float x1=mod, x2=mod, x1f=mod, x2f=mod;
     float I=0;
     int cnt=0;
     while (I<.68) {
         x1 -= dx;
         x2 += dx;
         if (pdf(x1)>=pdf(x2)) {
             x1f = x1;
             x2 -= dx;
         } else {
             x2f = x2;
             x1 += dx;
         }
         I = cdf(x2f) - cdf(x1f);
         cnt++;
         // for crazy pdf's return a crazy value
         if (cnt>2500) return 100000.;
     }
     return (x2f - x1f)/2;
 }

double GSUtilities::findMode ( const double xStart ) const

private

mean value of combined state double combinedMean() const; mode from starting value

Definition at line 122 of file GSUtilities.cc.

References dpdf1(), dpdf2(), MillePedeFileConverter_cfg::e, pdf(), and x.

Referenced by mode(), and ~GSUtilities().

 {
   //
   // try with Newton
   //
   double y1(0.);
   double x(xStart);
   double y2(pdf(xStart));
   double yd(dpdf1(xStart));
   int nLoop(0);
   if ( (y1+y2)<10*DBL_MIN )  return xStart;
   while ( nLoop++<20 && fabs(y2-y1)/(y2+y1)>1.e-6 ) {
 //     std::cout << "dy = " << y2-y1 << std::endl;
     double yd2 = dpdf2(x);
     if ( fabs(yd2)<10*DBL_MIN )  return xStart;
     x -= yd/dpdf2(x);
     yd = dpdf1(x);
     y1 = y2;
     y2 = pdf(x);
 //     std::cout << "New x / yd = " << x << " / " << yd << std::endl;
   }
   if ( nLoop >= 20 )  return xStart;
   return x;
 }

double GSUtilities::gauss	(	const double &	x,
		const double &	mean,
		const double &	sigma
	)		const

private

value of gaussian distribution

Definition at line 190 of file GSUtilities.cc.

References edmIntegrityCheck::d, JetChargeProducer_cfi::exp, Pi, mps_fire::result, and mathSSE::sqrt().

Referenced by dpdf1(), dpdf2(), pdf(), and ~GSUtilities().

 {
   const double fNorm(1./sqrt(2*TMath::Pi()));
   double result(0.);
 
   double d((x-mean)/sigma);
   if ( fabs(d)<20. )  result = exp(-d*d/2.);
   result *= fNorm/sigma;
   return result;
 }

double GSUtilities::gaussInt	(	const double &	x,
		const double &	mean,
		const double &	sigma
	)		const

private

integrated value of gaussian distribution

Definition at line 203 of file GSUtilities.cc.

Referenced by cdf(), and ~GSUtilities().

 {
   return TMath::Freq((x-mean)/sigma);
 }

float GSUtilities::getMax ( float x )

Definition at line 291 of file GSUtilities.cc.

References cdf(), PVValHelper::dx, and x.

Referenced by errorMode(), and ~GSUtilities().

 {
     int cnt=0;
     float dx;
     if (fabs(x)<2) dx = .5;
     else dx = fabs(x)/10;
     while( cdf(x)<.9 && cnt<1000) {
         x += dx;
         cnt++;
     }
     return x;
 }

float GSUtilities::getMin ( float x )

Definition at line 278 of file GSUtilities.cc.

References cdf(), PVValHelper::dx, and x.

Referenced by errorMode(), and ~GSUtilities().

 {
     int cnt=0;
     float dx;
     if (fabs(x)<2) dx = .5;
     else dx = fabs(x)/10;
     while( cdf(x)>.1 && cnt<1000) {
         x -= dx;
         cnt++;
     }
     return x;
 }

float GSUtilities::maxWeight ( ) const

Definition at line 232 of file GSUtilities.cc.

References mps_fire::i, theNComp, and theWeights.

Referenced by ~GSUtilities().

 {
   // Look for the highest weight component
   //
   //int iwMax(-1);
   float wMax(0.);
   for ( unsigned i=0; i<theNComp; i++ ) {
     if ( theWeights[i]>wMax ) {
       //iwMax = i;
       wMax = theWeights[i];
     }
   }
   return wMax;
 }

float GSUtilities::mode ( void ) const

mode

Definition at line 88 of file GSUtilities.cc.

References findMode(), mps_fire::i, pdf(), theNComp, theParameters, and x.

Referenced by errorMode(), and ~GSUtilities().

 {
   //
   // start values = means of components
   //
   typedef std::multimap<double, double, std::greater<double> > StartMap;
   StartMap xStart;
   for ( unsigned i=0; i<theNComp; i++ ) {
     xStart.insert(std::pair<double,float>(pdf(theParameters[i]),
                      theParameters[i]));
   }
   //
   // now try with each start value
   //
   typedef std::multimap<double, double, std::greater<double> > ResultMap;
   ResultMap xFound;
   for ( StartMap::const_iterator i=xStart.begin();
     i!=xStart.end(); i++ ) {
     double x = findMode((*i).second);
     xFound.insert(std::pair<double,double>(pdf(x),x));
 //     std::cout << "Started at " << (*i).second 
 //   << " , found " << x << std::endl;
   }  
   //
   // results
   //
 //   for ( ResultMap::const_iterator i=xFound.begin();
 //  i!=xFound.end(); i++ ) {
 //     std::cout << "pdf at " << (*i).second << " = " << (*i).first << std::endl;
 //   }
   return xFound.begin()->second;
 }

double GSUtilities::pdf ( const double & x ) const

value of the pdf

Definition at line 148 of file GSUtilities.cc.

References gauss(), mps_fire::i, mps_fire::result, theErrors, theNComp, theParameters, and theWeights.

Referenced by errorMode(), findMode(), mode(), and ~GSUtilities().

 {
   double result(0.);
   for ( unsigned i=0; i<theNComp; i++ )
     result += theWeights[i]*gauss(x,theParameters[i],theErrors[i]);
   return result;
 }

float GSUtilities::quantile ( const float q ) const

normalised integral from -inf to x (taking into account under- & overflows)

Definition at line 12 of file GSUtilities.cc.

References cdf(), PVValHelper::dx, MillePedeFileConverter_cfg::e, mps_fire::i, lumiQueryAPI::q, theErrors, theNComp, theParameters, theWeights, x, globals_cff::x1, globals_cff::x2, and y.

Referenced by ~GSUtilities().

 {
   float qq = q;
   if (q>1) qq = 1;
   else if (q<0) qq = 0;
   //
   // mean and sigma of highest weight component
   //
   int iwMax(-1);
   float wMax(0.);
   for ( unsigned i=0; i<theNComp; i++ ) {
     if ( theWeights[i]>wMax ) {
       iwMax = i;
       wMax = theWeights[i];
     }
   }
   //
   // Find start values: begin with mean and
   // go towards f(x)=0 until two points on
   // either side are found (assumes monotonously
   // increasing function)
   //
   double x1(theParameters[iwMax]);
   double y1(cdf(x1)-qq);
   double dx = y1>0. ? -theErrors[iwMax] : theErrors[iwMax];
   double x2(x1+dx);
   double y2(cdf(x2)-qq);
   int cnt =0;
   while ( y1*y2>0. ) {
     x1 = x2;
     y1 = y2;
     x2 += dx;
     y2 = cdf(x2) - qq;
     cnt++;
   }
 //   std::cout << "(" << x1 << "," << y1 << ") / (" 
 //        << x2 << "," << y2 << ")" << std::endl;
   //
   // now use bisection to find final value
   //
   double x(0.);
   while ( true ) {
     // use linear interpolation
     x = -(x2*y1-x1*y2)/(y2-y1);
     double y = cdf(x) - qq;
 //     std::cout << x << " " << y << std::endl;
     if ( fabs(y)<1.e-6 )  break;
     if ( y*y1>0. ) {
       y1 = y;
       x1 = x;
     }
     else {
       y2 = y;
       x2 = x;
     }
   }
   return x;
 }