GSUtilities Class Reference

#include <GSUtilities.h>

List of all members.

Public Member Functions
double	cdf (const double &) const
	value of integral(pdf)
double	combinedMean () const
	mean value of combined state
double	dpdf1 (const double &) const
	first derivative of pdf
double	dpdf2 (const double &) const
	second derivative of pdf
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
float	maxWeight () const
float	mode () const
	mode
double	pdf (const double &) const
	value of the pdf
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
double	gaussInt (const double &, const double &, const double &) const
	integrated value of gaussian distribution
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 i, theErrors, theNComp, theParameters, and theWeights.

                                                             :
    theNComp(nComp),
    theWeights(0),
    theParameters(0),
    theErrors(0)
  {
    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 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(), i, query::result, theErrors, theNComp, theParameters, and theWeights.

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

{
  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 i, theNComp, theParameters, and theWeights.

{
  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 gauss(), i, query::result, theErrors, theNComp, theParameters, and theWeights.

Referenced by findMode().

{
  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 gauss(), i, query::result, theErrors, theNComp, theParameters, and theWeights.

Referenced by findMode().

{
  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 i, mathSSE::sqrt(), theNComp, and theWeights.

{
  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 i, theErrors, theNComp, and theWeights.

{
        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(), getMax(), getMin(), Exhume::I, max(), min, mod(), mode(), and pdf().

{
        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(), alignCSCRings::e, pdf(), and x.

Referenced by mode().

{
  //
  // 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 funct::exp(), Pi, query::result, and mathSSE::sqrt().

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

{
  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().

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

float GSUtilities::getMax

(

float

x

)

Definition at line 291 of file GSUtilities.cc.

References cdf(), and x.

Referenced by errorMode().

{
        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(), and x.

Referenced by errorMode().

{
        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 i, theNComp, and theWeights.

{
  // 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(), i, pdf(), theNComp, theParameters, and x.

Referenced by errorMode().

{
  //
  // 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(), i, query::result, theErrors, theNComp, theParameters, and theWeights.

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

{
  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(), alignCSCRings::e, i, lumiQueryAPI::q, theErrors, theNComp, theParameters, theWeights, x, and detailsBasic3DVector::y.

{
  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;
}

---

Member Data Documentation

float* GSUtilities::theErrors [private]

Definition at line 83 of file GSUtilities.h.

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

unsigned GSUtilities::theNComp [private]

Definition at line 80 of file GSUtilities.h.

Referenced by cdf(), combinedMean(), dpdf1(), dpdf2(), errorCombinedMean(), errorHighestWeight(), GSUtilities(), maxWeight(), mode(), pdf(), and quantile().

float* GSUtilities::theParameters [private]

Definition at line 82 of file GSUtilities.h.

Referenced by cdf(), combinedMean(), dpdf1(), dpdf2(), GSUtilities(), mode(), pdf(), quantile(), and ~GSUtilities().

float* GSUtilities::theWeights [private]

Definition at line 81 of file GSUtilities.h.

Referenced by cdf(), combinedMean(), dpdf1(), dpdf2(), errorCombinedMean(), errorHighestWeight(), GSUtilities(), maxWeight(), pdf(), quantile(), and ~GSUtilities().