GSUtilities Class Reference

#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::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, theWeights, and hltDeepSecondaryVertexTagInfosPFPuppi_cfi::weights.

      : 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::~GSUtilities ( )

inline

Definition at line 30 of file GSUtilities.h.

References theErrors, theParameters, and theWeights.

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

Member Function Documentation

cdf()

double GSUtilities::cdf

(

const double &

x

)

const

value of integral(pdf)

Definition at line 145 of file GSUtilities.cc.

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

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

combinedMean()

double GSUtilities::combinedMean

(

)

const

mean value of combined state

Definition at line 185 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;
}

dpdf1()

double GSUtilities::dpdf1

(

const double &

x

)

const

first derivative of pdf

Definition at line 152 of file GSUtilities.cc.

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

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

dpdf2()

double GSUtilities::dpdf2

(

const double &

x

)

const

second derivative of pdf

Definition at line 161 of file GSUtilities.cc.

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

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

errorCombinedMean()

double GSUtilities::errorCombinedMean

(

)

const

Definition at line 195 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));
}

errorHighestWeight()

float GSUtilities::errorHighestWeight

(

)

const

Definition at line 69 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];
}

errorMode()

float GSUtilities::errorMode

(

)

Definition at line 217 of file GSUtilities.cc.

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

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

findMode()

double GSUtilities::findMode ( const double  xStart ) const

private

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

Definition at line 111 of file GSUtilities.cc.

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

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

gauss()

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

private

value of gaussian distribution

Definition at line 170 of file GSUtilities.cc.

References ztail::d, JetChargeProducer_cfi::exp, SiStripPI::mean, Pi, mps_fire::result, mathSSE::sqrt(), and x.

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

gaussInt()

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

private

integrated value of gaussian distribution

Definition at line 181 of file GSUtilities.cc.

References SiStripPI::mean, and x.

Referenced by cdf().

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

getMax()

float GSUtilities::getMax

(

float

x

)

Definition at line 260 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;
}

getMin()

float GSUtilities::getMin

(

float

x

)

Definition at line 246 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;
}

maxWeight()

float GSUtilities::maxWeight

(

)

const

Definition at line 203 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;
}

mode()

float GSUtilities::mode

(

void

)

const

mode

Definition at line 81 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;
}

pdf()

double GSUtilities::pdf

(

const double &

x

)

const

value of the pdf

Definition at line 138 of file GSUtilities.cc.

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

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

quantile()

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, MillePedeFileConverter_cfg::e, mps_fire::i, submitPVResolutionJobs::q, theErrors, theNComp, theParameters, theWeights, x, testProducerWithPsetDescEmpty_cfi::x1, testProducerWithPsetDescEmpty_cfi::x2, y, testProducerWithPsetDescEmpty_cfi::y1, and testProducerWithPsetDescEmpty_cfi::y2.

                                               {
  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);
  while (y1 * y2 > 0.) {
    x1 = x2;
    y1 = y2;
    x2 += dx;
    y2 = cdf(x2) - qq;
  }
  //   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

theErrors

float* GSUtilities::theErrors

private

Definition at line 77 of file GSUtilities.h.

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

theNComp

unsigned GSUtilities::theNComp

private

Definition at line 74 of file GSUtilities.h.

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

theParameters

float* GSUtilities::theParameters

private

Definition at line 76 of file GSUtilities.h.

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

theWeights

float* GSUtilities::theWeights

private

Definition at line 75 of file GSUtilities.h.

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