#include <GaussianSumUtilities1D.h>
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 SingleGaussianState1D & | mode () 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 MultiGaussianState1D & | theState |
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.
enum GaussianSumUtilities1D::ModeStatus [private] |
Definition at line 18 of file GaussianSumUtilities1D.h.
{ Valid, NotValid, NotComputed };
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.
{}
double GaussianSumUtilities1D::cdf | ( | const double & | x | ) | const |
value of the c.d.f.
Definition at line 280 of file GaussianSumUtilities1D.cc.
References gaussInt(), i, mean(), query::result, asciidump::s, size(), standardDeviation(), and weight().
double GaussianSumUtilities1D::combinedMean | ( | ) | const [private] |
Mean value of combined state.
Definition at line 455 of file GaussianSumUtilities1D.cc.
References components(), i, asciidump::s, size(), and weight().
const std::vector<SingleGaussianState1D>& GaussianSumUtilities1D::components | ( | ) | const [inline] |
components
Definition at line 30 of file GaussianSumUtilities1D.h.
References MultiGaussianState1D::components(), and theState.
Referenced by combinedMean(), computeMode(), mean(), pdfComponents(), size(), standardDeviation(), variance(), and weight().
{ return theState.components(); }
void GaussianSumUtilities1D::computeMode | ( | ) | const [private] |
calculation of mode
Definition at line 102 of file GaussianSumUtilities1D.cc.
References components(), ExpressReco_HICollisions_FallBack::e, findMode(), first, i, localVariance(), M_PI, mean(), mode(), NotValid, pdf(), asciidump::s, edm::second(), size(), mathSSE::sqrt(), standardDeviation(), theMode, theModeStatus, Valid, TrackValidation_HighPurity_cff::valid, variance(), weight(), ExpressReco_HICollisions_FallBack::x, and ExpressReco_HICollisions_FallBack::y.
Referenced by mode(), and modeIsValid().
{ // TimerStack tstack; // tstack.benchmark("GSU1D::benchmark",100000); // FastTimerStackPush(tstack,"GaussianSumUtilities1D::computeMode"); theModeStatus = NotValid; // // 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 314 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 414 of file GaussianSumUtilities1D.cc.
References d1Pdf(), f, min, pdf(), and query::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 368 of file GaussianSumUtilities1D.cc.
References i, mean(), query::result, asciidump::s, size(), and standardDeviation().
double GaussianSumUtilities1D::d1Pdf | ( | const double & | x | ) | const |
first derivative of the p.d.f.
Definition at line 290 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 320 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 425 of file GaussianSumUtilities1D.cc.
References d1LnPdf(), d2Pdf(), f, min, pdf(), and query::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 380 of file GaussianSumUtilities1D.cc.
References i, mean(), query::result, asciidump::s, size(), and standardDeviation().
double GaussianSumUtilities1D::d2Pdf | ( | const double & | x | ) | const |
second derivative of the p.d.f.
Definition at line 296 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 392 of file GaussianSumUtilities1D.cc.
References i, mean(), query::result, asciidump::s, size(), and standardDeviation().
double GaussianSumUtilities1D::d3Pdf | ( | const double & | x | ) | const |
third derivative of the p.d.f.
Definition at line 302 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 207 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 436 of file GaussianSumUtilities1D.cc.
Referenced by pdf(), and pdfComponents().
double GaussianSumUtilities1D::gaussInt | ( | double | x, |
double | mean, | ||
double | sigma | ||
) | [static, private] |
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 405 of file GaussianSumUtilities1D.cc.
References f, funct::log(), max(), min, pdf(), and query::result.
double GaussianSumUtilities1D::lnPdf | ( | const double & | x | ) | const |
ln(pdf)
Definition at line 308 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 469 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] |
mean value of a component
Definition at line 36 of file GaussianSumUtilities1D.h.
References components(), and i.
Referenced by MultiTrajectoryStateMode::chargeFromMode(), GsfTrackProducerBase::fillMode(), MultiTrajectoryStateMode::momentumFromModeLocal(), MultiTrajectoryStateMode::momentumFromModeP(), MultiTrajectoryStateMode::momentumFromModePPhiEta(), MultiTrajectoryStateMode::momentumFromModeQP(), and MultiTrajectoryStateMode::positionFromModeLocal().
{return components()[i].mean();}
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().
const SingleGaussianState1D & GaussianSumUtilities1D::mode | ( | void | ) | const |
Mode "state": mean = mode, variance = local variance at mode, weight chosen to have pdf(mode) equal to the one of the mixture
Definition at line 94 of file GaussianSumUtilities1D.cc.
References computeMode(), NotComputed, theMode, and theModeStatus.
Referenced by MultiTrajectoryStateMode::chargeFromMode(), computeMode(), PFGsfHelper::computeQpMode(), ElectronMomentumCorrector::correct(), GsfTrackProducerBase::fillMode(), GsfTrackProducerBase::localParametersFromQpMode(), MultiTrajectoryStateMode::momentumFromModeCartesian(), MultiTrajectoryStateMode::momentumFromModeLocal(), MultiTrajectoryStateMode::momentumFromModeP(), MultiTrajectoryStateMode::momentumFromModePPhiEta(), MultiTrajectoryStateMode::momentumFromModeQP(), PFGsfHelper::PFGsfHelper(), MultiTrajectoryStateMode::positionFromModeCartesian(), and MultiTrajectoryStateMode::positionFromModeLocal().
{ if ( theModeStatus==NotComputed ) computeMode(); // std::cout << "Mode calculation failed!!" << std::endl; return theMode; }
bool GaussianSumUtilities1D::modeIsValid | ( | ) | const |
mode status
Definition at line 87 of file GaussianSumUtilities1D.cc.
References computeMode(), NotComputed, theModeStatus, and Valid.
Referenced by MultiTrajectoryStateMode::chargeFromMode(), PFGsfHelper::computeQpMode(), GsfTrackProducerBase::fillMode(), GsfTrackProducerBase::localParametersFromQpMode(), MultiTrajectoryStateMode::momentumFromModeLocal(), MultiTrajectoryStateMode::momentumFromModeP(), MultiTrajectoryStateMode::momentumFromModePPhiEta(), MultiTrajectoryStateMode::momentumFromModeQP(), and MultiTrajectoryStateMode::positionFromModeLocal().
{ if ( theModeStatus==NotComputed ) computeMode(); return theModeStatus==Valid; }
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().
double GaussianSumUtilities1D::pdf | ( | double | x | ) | const |
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 362 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 344 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 326 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().
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] |
Definition at line 257 of file GaussianSumUtilities1D.cc.
References d1Pdf(), d2Pdf(), min, pdf(), pdfComponents(), GaussianSumUtilities1D::FinderState::pdfs, ExpressReco_HICollisions_FallBack::x, GaussianSumUtilities1D::FinderState::x, GaussianSumUtilities1D::FinderState::y, GaussianSumUtilities1D::FinderState::yd, and GaussianSumUtilities1D::FinderState::yd2.
Referenced by findMode().
{ state.x = x; pdfComponents(state.x, state.pdfs); state.y = pdf(state.x, state.pdfs); state.yd = 0; if (state.y>std::numeric_limits<double>::min()) state.yd= d1Pdf(state.x,state.pdfs)/state.y; state.yd2 = -state.yd*state.yd; if (state.y > std::numeric_limits<double>::min()) state.yd2 += d2Pdf(state.x,state.pdfs)/state.y; }
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().
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().
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();}
SingleGaussianState1D GaussianSumUtilities1D::theMode [mutable, private] |
Definition at line 139 of file GaussianSumUtilities1D.h.
Referenced by computeMode(), and mode().
ModeStatus GaussianSumUtilities1D::theModeStatus [mutable, private] |
Definition at line 138 of file GaussianSumUtilities1D.h.
Referenced by computeMode(), mode(), and modeIsValid().
const MultiGaussianState1D& GaussianSumUtilities1D::theState [private] |
Definition at line 135 of file GaussianSumUtilities1D.h.
Referenced by components(), mean(), variance(), and weight().