#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 30 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 147 of file GSUtilities.cc.

References gaussInt(), mps_fire::i, mps_fire::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 187 of file GSUtilities.cc.

References mps_fire::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 154 of file GSUtilities.cc.

References PVValHelper::dx, gauss(), mps_fire::i, mps_fire::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 163 of file GSUtilities.cc.

References PVValHelper::dx, gauss(), mps_fire::i, mps_fire::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 197 of file GSUtilities.cc.

References mps_fire::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 71 of file GSUtilities.cc.

References mps_fire::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 219 of file GSUtilities.cc.

References cdf(), PVValHelper::dx, getMax(), getMin(), Exhume::I, SiStripPI::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 113 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 172 of file GSUtilities.cc.

References ztail::d, funct::exp(), Pi, mps_fire::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 183 of file GSUtilities.cc.

Referenced by cdf().

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

float GSUtilities::getMax ( float x )

Definition at line 262 of file GSUtilities.cc.

References cdf(), PVValHelper::dx, 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 248 of file GSUtilities.cc.

References cdf(), PVValHelper::dx, 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 205 of file GSUtilities.cc.

References mps_fire::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 83 of file GSUtilities.cc.

References findMode(), mps_fire::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 140 of file GSUtilities.cc.

References gauss(), mps_fire::i, mps_fire::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 11 of file GSUtilities.cc.

References cdf(), PVValHelper::dx, alignCSCRings::e, mps_fire::i, submitPVResolutionJobs::q, theErrors, theNComp, theParameters, theWeights, x, and 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;
 }