CMSSW Reference Manual

enum GaussianSumUtilities1D::ModeStatus [private]

Enumerator:

Valid
NotValid
NotComputed

Definition at line 18 of file GaussianSumUtilities1D.h.

00018 { Valid, NotValid, NotComputed };

GaussianSumUtilities1D::GaussianSumUtilities1D ( const MultiGaussianState1D & state ) [inline]

Definition at line 21 of file GaussianSumUtilities1D.h.

00021                                                              :
00022     theState(state),
00023 //     theStates(state.components()),
00024     theModeStatus(NotComputed) {} 
  ~GaussianSumUtilities1D () {}

GaussianSumUtilities1D::~GaussianSumUtilities1D ( ) [inline]

Definition at line 25 of file GaussianSumUtilities1D.h.

00025 {}

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

value of the c.d.f.

Definition at line 397 of file GaussianSumUtilities1D.cc.

References gaussInt(), i, mean(), HLT_VtxMuL3::result, size(), standardDeviation(), and weight().

Referenced by quantile().

00398 {
00399   double result(0.);
00400   for ( unsigned int i=0; i<size(); i++ )
00401     result += weight(i)*gaussInt(x,mean(i),standardDeviation(i));
00402   return result;
00403 }

double GaussianSumUtilities1D::combinedMean ( ) const [private]

Mean value of combined state.

Definition at line 546 of file GaussianSumUtilities1D.cc.

References i, mean(), s1, size(), and weight().

00547 {
00548   double s0(0.);
00549   double s1(0.);
00550   for ( unsigned int i=0; i<size(); i++ ) {
00551     s0 += weight(i);
00552     s1 += weight(i)*mean(i);
00553   }
00554   return s1/s0;
00555 }

const std::vector<SingleGaussianState1D>& GaussianSumUtilities1D::components ( ) const [inline]

components

Definition at line 30 of file GaussianSumUtilities1D.h.

References MultiGaussianState1D::components(), and theState.

Referenced by computeMode(), mean(), size(), standardDeviation(), variance(), and weight().

00030                                                                      {
00031     return theState.components();
00032   }

void GaussianSumUtilities1D::computeMode ( ) const [private]

calculation of mode

Definition at line 92 of file GaussianSumUtilities1D.cc.

References components(), GenMuonPlsPt100GeV_cfg::cout, e, lat::endl(), findMode(), first, g, i, localVariance(), mean(), SingleGaussianState1D::mean(), mode(), NotValid, np, pdf(), Pi, edm::second(), size(), funct::sqrt(), standardDeviation(), theMode, theModeStatus, Valid, TrackValidation_HighPurity_cff::valid, variance(), SingleGaussianState1D::variance(), w, weight(), x, and y.

Referenced by mode(), and modeIsValid().

00093 {
00094 //   TimerStack tstack;
00095 //   tstack.benchmark("GSU1D::benchmark",100000);
00096 //   FastTimerStackPush(tstack,"GaussianSumUtilities1D::computeMode");
00097 #ifdef DRAW_GS1D
00098   {
00099   gPad->Clear();
00100   std::cout << "gPad = " << gPad << std::endl;
00101   std::multimap<double,double> drawMap;
00102   const int npDraw(1024);
00103   double xpDraw[npDraw];
00104   double ypDraw[npDraw];
00105   for ( unsigned int i=0; i<size(); i++ ) {
00106     double ave = mean(i);
00107     double sig = standardDeviation(i);
00108     for ( double xi=-3.; xi<3.; ) {
00109       double x = ave + xi*sig;
00110       drawMap.insert(std::make_pair(x,pdf(x)));
00111       xi += 0.2;
00112     }
00113   }
00114   int np(0);
00115   double xMin(FLT_MAX);
00116   double xMax(-FLT_MAX);
00117   double yMax(-FLT_MAX);
00118   for ( std::multimap<double,double>::const_iterator im=drawMap.begin();
00119         im!=drawMap.end(); ++im,++np ) {
00120     if ( np>=1024 )  break;
00121     xpDraw[np] = (*im).first;
00122     ypDraw[np] = (*im).second;
00123     if ( xMin>(*im).first )  xMin = (*im).first;
00124     if ( xMax<(*im).first )  xMax = (*im).first;
00125     if ( yMax<(*im).second )  yMax = (*im).second;
00126   }
00127 //   for ( int i=0; i<np; ++i )
00128 //     std::cout << i << " " << xpDraw[i] << " " << ypDraw[i] << std::endl;
00129   gPad->DrawFrame(xMin,0.,xMax,1.05*yMax);  
00130   TGraph* g = new TGraph(np,xpDraw,ypDraw);
00131   g->SetLineWidth(2);
00132   g->Draw("C");
00133   TGraph* gc = new TGraph();
00134   gc->SetLineStyle(2);
00135   int np2(0);
00136   double xpDraw2[1024];
00137   double ypDraw2[1024];
00138   for ( unsigned int i=0; i<size(); i++ ) {
00139     double ave = mean(i);
00140     double sig = standardDeviation(i);
00141     SingleGaussianState1D sgs(ave,variance(i),weight(i));
00142     std::vector<SingleGaussianState1D> sgsv(1,sgs);
00143     MultiGaussianState1D mgs(sgsv);
00144     GaussianSumUtilities1D gsu(mgs);
00145     np2 = 0;
00146     for ( double xi=-3.; xi<3.; ) {
00147       double x = ave + xi*sig;
00148       xpDraw2[np2] = x;
00149       ypDraw2[np2] = gsu.pdf(x);
00150       ++np2;
00151       xi += 0.2;
00152     }
00153 //    for ( int i=0; i<np2; ++i )
00154 //      std::cout << i << " " << xpDraw2[i] << " " << ypDraw2[i] << std::endl;
00155 //     std::cout << "np2 = " << np2 
00156 //       << " ave = " << ave 
00157 //       << " sig = " << sig 
00158 //       << " weight = " << weight(i)
00159 //       << " pdf = " << gsu.pdf(ave) << std::endl;
00160     if ( np2>0 )  gc->DrawGraph(np2,xpDraw2,ypDraw2);
00161 //     break;
00162   }
00163   gPad->Modified();
00164   gPad->Update();
00165   }
00166 #endif
00167   theModeStatus = NotValid;
00168   //
00169   // Use means of individual components as start values.
00170   // Sort by value of pdf.
00171   //
00172   typedef std::multimap<double, int, std::greater<double> > StartMap;
00173   StartMap xStart;
00174   for ( unsigned int i=0; i<size(); i++ ) {
00175     xStart.insert(std::make_pair(pdf(mean(i)),i));
00176   }
00177   //
00178   // Now try with each start value
00179   //
00180 #ifdef DRAW_GS1D
00181   int ind(0);
00182 #endif
00183   int iRes(-1);     // index of start component for current estimate
00184   double xRes(mean((*xStart.begin()).second)); // current estimate of mode
00185   double yRes(-1.); // pdf at current estimate of mode
00186 //   std::pair<double,double> result(-1.,mean((*xStart.begin()).second));
00187   for ( StartMap::const_iterator i=xStart.begin(); i!=xStart.end(); i++ ) {
00188 #ifdef DRAW_GS1D
00189     double xp[2];
00190     double yp[2];
00191     TPolyMarker* g = new TPolyMarker();
00192     g->SetMarkerStyle(20+ind/4);
00193     g->SetMarkerColor(1+ind%4);
00194     ++ind;
00195 #endif
00196     //
00197     // Convergence radius for a single Gaussian = 1 sigma: don't try
00198     // start values within 1 sigma of the current solution
00199     //
00200     if ( theModeStatus==Valid &&
00201          fabs(mean((*i).second)-mean(iRes))/standardDeviation(iRes)<1. )  continue;
00202     //
00203     // If a solution exists: drop as soon as the pdf at
00204     // start value drops to < 75% of maximum (try to avoid
00205     // unnecessary searches for the mode)
00206     //
00207     if ( theModeStatus==Valid && 
00208          (*i).first/(*xStart.begin()).first<0.75 )  break;
00209     //
00210     // Try to find mode
00211     //
00212     double x;
00213     double y;
00214     bool valid = findMode(x,y,mean((*i).second),standardDeviation((*i).second));
00215     //
00216     // consider only successful searches
00217     //
00218     if ( valid ) { //...
00219       //
00220       // update result only for significant changes in pdf(solution)
00221       //
00222       if ( yRes<0. || (y-yRes)/(y+yRes)>1.e-10 ) {
00223 #ifdef DBG_GS1D
00224       if ( iRes>=0 ) 
00225         std::cout << "dxStart = " << (xRes-mean((*i).second))/standardDeviation(iRes) << std::endl;
00226 #endif
00227       iRes = (*i).second;               // store index
00228       theModeStatus = Valid;            // update status
00229       xRes = x;                         // current solution
00230       yRes = y;                         // and its pdf value
00231 //       result = std::make_pair(y,x);     // store solution and pdf(solution)
00232       }
00233     } //...
00234 #ifdef DRAW_GS1D
00235     xp[0] = xp[1] = mean((*i).second);
00236     yp[0] = yp[1] = (*i).first;
00237     if ( valid ) {
00238       xp[1] = x;
00239       yp[1] = y;
00240     }
00241     g->DrawPolyMarker(2,xp,yp);
00242     ++ind; //...
00243 #endif
00244   } 
00245   //
00246   // check (existance of) solution
00247   //
00248   if ( theModeStatus== Valid ) {
00249 #ifdef DBG_GS1D
00250     std::cout << "Ratio of pdfs at  first / real maximum = " << iRes << " " << mean(iRes) << " "
00251          << pdf(mean(iRes))/(*xStart.begin()).first << std::endl;
00252 #endif
00253     //
00254     // Construct single Gaussian state with 
00255     //  mean = mode
00256     //  variance = local variance at mode
00257     //  weight such that the pdf's of the mixture and the
00258     //    single state are equal at the mode
00259     //
00260     double mode = xRes;
00261     double varMode = localVariance(mode);
00262     double wgtMode = pdf(mode)*sqrt(2*TMath::Pi()*varMode);
00263     theMode = SingleGaussianState1D(mode,varMode,wgtMode);
00264   }
00265   else {
00266     //
00267     // mode finding failed: set solution to highest component
00268     //  (alternative would be global mean / variance ..?)
00269     //
00270     edm::LogWarning("GaussianSumUtilities") << "1D mode calculation failed";
00271 //     double x = mean();
00272 //     double y = pdf(x);
00273 //     result = std::make_pair(y,x);
00274 //     theMode = SingleGaussianState1D(mean(),variance(),weight());
00275     //
00276     // look for component with highest value at mean
00277     //
00278     unsigned int icMax(0);
00279     double ySqMax(0.);
00280     for ( unsigned int ic=0; ic<size(); ++ic ) {
00281       double w = weight(ic);
00282       double ySq = w*w/variance(ic);
00283       if ( ic==0 || ySqMax ) {
00284         icMax = ic;
00285         ySqMax = ySq;
00286       }
00287     }
00288     theMode = SingleGaussianState1D(components()[icMax]);
00289   }
00290 #ifdef DRAW_GS1D
00291   gPad->Modified();
00292   gPad->Update();
00293 #endif
00294 #ifdef DRAW_GS1D
00295   {
00296     TGraph* gm = new TGraph();
00297     gm->SetLineColor(2);
00298     int np(0);
00299     double xp[1024];
00300     double yp[1024];
00301     double mode = theMode.mean();
00302     double var = theMode.variance();
00303     double sig = sqrt(var);
00304     SingleGaussianState1D sgs(mode,var,pdf(mode)*sqrt(2*TMath::Pi())*sig);
00305     std::vector<SingleGaussianState1D> sgsv(1,sgs);
00306     MultiGaussianState1D mgs(sgsv);
00307     GaussianSumUtilities1D gsu(mgs);
00308     for ( double xi=-3.; xi<3.; ) {
00309       double x = mode + xi*sig;
00310       xp[np] = x;
00311       yp[np] = gsu.pdf(x);
00312       ++np;
00313       xi += 0.2;
00314     }
00315     gm->DrawGraph(np,xp,yp);
00316     gPad->Modified();
00317     gPad->Update();
00318   }
00319 #endif
00320   
00321 }

double GaussianSumUtilities1D::d1LnPdf	(	const 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 504 of file GaussianSumUtilities1D.cc.

References d1Pdf(), f, pdf(), and HLT_VtxMuL3::result.

00505 {
00506 
00507   double f = pdf(x,pdfs);
00508   double result(d1Pdf(x,pdfs));
00509   if ( f>DBL_MIN )  result /= f;
00510   else  result = 0.;
00511   return result;
00512 }

double GaussianSumUtilities1D::d1LnPdf ( const double & x ) const

first derivative of ln(pdf)

Definition at line 430 of file GaussianSumUtilities1D.cc.

References pdfComponents().

Referenced by d2LnPdf(), and findMode().

00431 {
00432   return d1LnPdf(x,pdfComponents(x));
00433 }

double GaussianSumUtilities1D::d1Pdf	(	const 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 461 of file GaussianSumUtilities1D.cc.

References i, mean(), HLT_VtxMuL3::result, size(), and standardDeviation().

00462 {
00463   double result(0.);
00464   for ( unsigned int i=0; i<size(); i++ ) {
00465     double dx = (x-mean(i))/standardDeviation(i);
00466     result += -pdfs[i]*dx/standardDeviation(i);
00467   }
00468   return result;
00469 }

double GaussianSumUtilities1D::d1Pdf ( const double & x ) const

first derivative of the p.d.f.

Definition at line 406 of file GaussianSumUtilities1D.cc.

References pdfComponents().

Referenced by d1LnPdf().

00407 {
00408   return d1Pdf(x,pdfComponents(x));
00409 }

double GaussianSumUtilities1D::d2LnPdf	(	const 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 515 of file GaussianSumUtilities1D.cc.

References d1LnPdf(), d2Pdf(), f, pdf(), and HLT_VtxMuL3::result.

00516 {
00517 
00518   double f = pdf(x,pdfs);
00519   double df = d1LnPdf(x,pdfs);
00520   double result(-df*df);
00521   if ( f>DBL_MIN )  result += d2Pdf(x,pdfs)/f;
00522   return result;
00523 }

double GaussianSumUtilities1D::d2LnPdf ( const double & x ) const

second derivative of ln(pdf)

Definition at line 436 of file GaussianSumUtilities1D.cc.

References pdfComponents().

Referenced by findMode().

00437 {
00438   return d2LnPdf(x,pdfComponents(x));
00439 }

double GaussianSumUtilities1D::d2Pdf	(	const 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 472 of file GaussianSumUtilities1D.cc.

References i, mean(), HLT_VtxMuL3::result, size(), and standardDeviation().

00473 {
00474   double result(0.);
00475   for ( unsigned int i=0; i<size(); i++ ) {
00476     double dx = (x-mean(i))/standardDeviation(i);
00477     result += pdfs[i]/standardDeviation(i)/standardDeviation(i)*(dx*dx-1);
00478   }
00479   return result;
00480 }

double GaussianSumUtilities1D::d2Pdf ( const double & x ) const

second derivative of the p.d.f.

Definition at line 412 of file GaussianSumUtilities1D.cc.

References pdfComponents().

Referenced by d2LnPdf(), and localVariance().

00413 {
00414   return d2Pdf(x,pdfComponents(x));
00415 }

double GaussianSumUtilities1D::d3Pdf	(	const 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 483 of file GaussianSumUtilities1D.cc.

References i, mean(), HLT_VtxMuL3::result, size(), and standardDeviation().

00484 {
00485   double result(0.);
00486   for ( unsigned int i=0; i<size(); i++ ) {
00487     double dx = (x-mean(i))/standardDeviation(i);
00488     result += pdfs[i]/standardDeviation(i)/standardDeviation(i)/standardDeviation(i)*
00489       (-dx*dx+3)*dx;
00490   }
00491   return result;
00492 }

double GaussianSumUtilities1D::d3Pdf ( const double & x ) const

third derivative of the p.d.f.

Definition at line 418 of file GaussianSumUtilities1D.cc.

References pdfComponents().

00419 {
00420   return d3Pdf(x,pdfComponents(x));
00421 }

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 324 of file GaussianSumUtilities1D.cc.

References GenMuonPlsPt100GeV_cfg::cout, d1LnPdf(), d2LnPdf(), e, lat::endl(), pdf(), pdfComponents(), pdfs, and HLT_VtxMuL3::result.

Referenced by computeMode().

00327 {
00328   //
00329   // try with Newton on (lnPdf)'(x)
00330   //
00331   double x1(0.);
00332   double y1(0.);
00333   std::vector<double> pdfs(pdfComponents(xStart));
00334   double x2(xStart);
00335   double y2(pdf(xStart,pdfs));
00336   double yd(d1LnPdf(xStart,pdfs));
00337   double yd2(d2LnPdf(xStart,pdfs));
00338   double xmin(xStart-1.*scale);
00339   double xmax(xStart+1.*scale);
00340   //
00341   // preset result
00342   //
00343   bool result(false);
00344   xMode = x2;
00345   yMode = y2;
00346   //
00347   // Iterate
00348   //
00349   int nLoop(0);
00350   while ( nLoop++<20 ) {
00351 #ifdef DBG_GS1D
00352     std::cout << "Loop " << nLoop << " " << yd2 << " "
00353               << (yd*scale) << " " << (y2-y1)/(y2+y1) << std::endl;
00354 #endif
00355     if ( nLoop>1 && yd2<0. &&  
00356          ( fabs(yd*scale)<1.e-10 || fabs(y2-y1)/(y2+y1)<1.e-14 ) ) {
00357       result = true;
00358       break;
00359     }
00360     if ( fabs(yd2)<FLT_MIN )  yd2 = yd2>0. ? FLT_MIN : -FLT_MIN;
00361     double dx = -yd/yd2;
00362     x1 = x2;
00363     y1 = y2;
00364     x2 += dx;
00365 #ifdef DBG_GS1D
00366     std::cout << yd << " " << yd2 << " " << x2 << " " << xmin << " " << xmax << std::endl;
00367 #endif
00368     if ( yd2>0. && (x2<xmin||x2>xmax) )  return false;
00369 
00370     pdfs = pdfComponents(x2);
00371     y2 = pdf(x2,pdfs);
00372     yd = d1LnPdf(x2,pdfs);
00373     yd2 = d2LnPdf(x2,pdfs);
00374   }
00375   //
00376   // result
00377   //
00378   if ( result ) {
00379     xMode = x2;
00380     yMode = y2;
00381   }
00382 #ifdef DBG_GS1D
00383   std::cout << "Started from " << xStart << " " << pdf(xStart)
00384             << " ; ended at " << xMode << " " << yMode << " after " 
00385             << nLoop << " iterations " << result << std::endl;
00386 #endif
00387   return result;
00388 }

double GaussianSumUtilities1D::gauss	(	const double &	x,
		const double &	mean,
		const double &	sigma
	)			const `[private]`

Value of gaussian distribution.

Definition at line 526 of file GaussianSumUtilities1D.cc.

References d, funct::exp(), Pi, HLT_VtxMuL3::result, and funct::sqrt().

Referenced by pdfComponents().

00528 {
00529   const double fNorm(1./sqrt(2*TMath::Pi()));
00530   double result(0.);
00531 
00532   double d((x-mean)/sigma);
00533   if ( fabs(d)<20. )  result = exp(-d*d/2.);
00534   result *= fNorm/sigma;
00535   return result;
00536 }

double GaussianSumUtilities1D::gaussInt	(	const double &	x,
		const double &	mean,
		const double &	sigma
	)			const `[private]`

Integrated value of gaussian distribution.

Definition at line 539 of file GaussianSumUtilities1D.cc.

Referenced by cdf().

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

double GaussianSumUtilities1D::lnPdf	(	const double &	x,
		const std::vector< double > &	pdfs
	)			const `[private]`

ln(pdf) using the pdf components at the evaluation point

Definition at line 495 of file GaussianSumUtilities1D.cc.

References f, funct::log(), pdf(), and HLT_VtxMuL3::result.

00496 {
00497   double f(pdf(x,pdfs));
00498   double result(-FLT_MAX);
00499   if ( result>DBL_MIN )  result = log(f);
00500   return result;
00501 }

double GaussianSumUtilities1D::lnPdf ( const double & x ) const

ln(pdf)

Definition at line 424 of file GaussianSumUtilities1D.cc.

References pdfComponents().

00425 {
00426   return lnPdf(x,pdfComponents(x));
00427 }

double GaussianSumUtilities1D::localVariance ( const double & x ) const [private]

Local variance from Hessian matrix.

Only valid if x corresponds to a (local) maximum!

Definition at line 558 of file GaussianSumUtilities1D.cc.

References d2Pdf(), pdf(), and HLT_VtxMuL3::result.

Referenced by computeMode().

00559 {
00560   double result = -pdf(x)/d2Pdf(x);
00561   // FIXME: wrong curvature seems to be non-existant but should add a proper recovery
00562   if ( result<0. )
00563     edm::LogWarning("GaussianSumUtilities") << "1D variance at mode < 0";    
00564   return result;
00565 }

double GaussianSumUtilities1D::mean ( ) const [inline]

combined mean

Definition at line 72 of file GaussianSumUtilities1D.h.

References MultiGaussianState1D::mean(), and theState.

Referenced by cdf(), combinedMean(), computeMode(), d1Pdf(), d2Pdf(), d3Pdf(), pdfComponents(), and quantile().

00072                        {
00073     return theState.mean();
00074   }

double GaussianSumUtilities1D::mean ( unsigned int i ) const [inline]

mean value of a component

Definition at line 36 of file GaussianSumUtilities1D.h.

References components().

Referenced by GsfTrackProducerBase::fillMode().

00036 {return components()[i].mean();}

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 84 of file GaussianSumUtilities1D.cc.

References computeMode(), NotComputed, theMode, and theModeStatus.

Referenced by GlobalGsfElectronAlgo::computeMode(), GsfElectronAlgo::computeMode(), computeMode(), PFGsfHelper::computeQpMode(), ElectronMomentumCorrector::correct(), GsfTrackProducerBase::fillMode(), and PFGsfHelper::PFGsfHelper().

00085 {
00086   if ( theModeStatus==NotComputed )  computeMode();
00087 //     std::cout << "Mode calculation failed!!" << std::endl;
00088   return theMode;
00089 }

bool GaussianSumUtilities1D::modeIsValid ( ) const

mode status

Definition at line 77 of file GaussianSumUtilities1D.cc.

References computeMode(), NotComputed, theModeStatus, and Valid.

Referenced by PFGsfHelper::computeQpMode(), and GsfTrackProducerBase::fillMode().

00078 {
00079   if ( theModeStatus==NotComputed )  computeMode();
00080   return theModeStatus==Valid;
00081 }

double GaussianSumUtilities1D::pdf	(	const double &	x,
		const std::vector< double > &	pdfs
	)			const `[private]`

value of the p.d.f. using the pdf components at the evaluation point

Definition at line 452 of file GaussianSumUtilities1D.cc.

References i, HLT_VtxMuL3::result, and size().

00453 {
00454   double result(0.);
00455   for ( unsigned int i=0; i<size(); i++ )
00456     result += pdfs[i];
00457   return result;
00458 }

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

value of the p.d.f.

Definition at line 391 of file GaussianSumUtilities1D.cc.

References pdfComponents().

Referenced by computeMode(), d1LnPdf(), d2LnPdf(), findMode(), lnPdf(), and localVariance().

00392 {
00393   return pdf(x,pdfComponents(x));
00394 }

std::vector< double > GaussianSumUtilities1D::pdfComponents ( const double & x ) const [private]

pdf components

Definition at line 442 of file GaussianSumUtilities1D.cc.

References gauss(), i, mean(), HLT_VtxMuL3::result, size(), standardDeviation(), and weight().

Referenced by d1LnPdf(), d1Pdf(), d2LnPdf(), d2Pdf(), d3Pdf(), findMode(), lnPdf(), and pdf().

00443 {
00444   std::vector<double> result;
00445   result.reserve(size());
00446   for ( unsigned int i=0; i<size(); i++ )
00447     result.push_back(weight(i)*gauss(x,mean(i),standardDeviation(i)));
00448   return result;
00449 }

double GaussianSumUtilities1D::quantile ( const double ) const

Quantile (i.e. x for a given value of the c.d.f.).

Definition at line 25 of file GaussianSumUtilities1D.cc.

References cdf(), e, i, mean(), size(), standardDeviation(), weight(), x, and y.

00026 {
00027   //
00028   // mean and sigma of highest weight component
00029   //
00030   int iwMax(-1);
00031   float wMax(0.);
00032   for ( unsigned int i=0; i<size(); i++ ) {
00033     if ( weight(i)>wMax ) {
00034       iwMax = i;
00035       wMax = weight(i);
00036     }
00037   }
00038   //
00039   // Find start values: begin with mean and
00040   // go towards f(x)=0 until two points on
00041   // either side are found (assumes monotonously
00042   // increasing function)
00043   //
00044   double x1(mean(iwMax));
00045   double y1(cdf(x1)-q);
00046   double dx = y1>0. ? -standardDeviation(iwMax) : standardDeviation(iwMax);
00047   double x2(x1+dx);
00048   double y2(cdf(x2)-q);
00049   while ( y1*y2>0. ) {
00050     x1 = x2;
00051     y1 = y2;
00052     x2 += dx;
00053     y2 = cdf(x2) - q;
00054   }
00055   //
00056   // now use bisection to find final value
00057   //
00058   double x(0.);
00059   while ( true ) {
00060     // use linear interpolation
00061     x = -(x2*y1-x1*y2)/(y2-y1);
00062     double y = cdf(x) - q;
00063     if ( fabs(y)<1.e-6 )  break;
00064     if ( y*y1>0. ) {
00065       y1 = y;
00066       x1 = x;
00067     }
00068     else {
00069       y2 = y;
00070       x2 = x;
00071     }
00072   }
00073   return x;
00074 }

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(), pdf(), pdfComponents(), and quantile().

00028 {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(), funct::sqrt(), and variance().

Referenced by cdf(), computeMode(), d1Pdf(), d2Pdf(), d3Pdf(), pdfComponents(), and quantile().

00038                                                          {
00039     return sqrt(components()[i].variance());
00040   }

double GaussianSumUtilities1D::variance ( ) const [inline]

combined covariance

Definition at line 76 of file GaussianSumUtilities1D.h.

References theState, and MultiGaussianState1D::variance().

Referenced by computeMode(), and standardDeviation().

00076                            {
00077     return theState.variance();
00078   }

double GaussianSumUtilities1D::variance ( unsigned int i ) const [inline]

variance of a component

Definition at line 42 of file GaussianSumUtilities1D.h.

References components().

Referenced by GsfTrackProducerBase::fillMode().

00042 {return components()[i].variance();}

double GaussianSumUtilities1D::weight ( ) const [inline]

combined weight

Definition at line 68 of file GaussianSumUtilities1D.h.

References theState, and MultiGaussianState1D::weight().

Referenced by cdf(), combinedMean(), computeMode(), pdfComponents(), and quantile().

00068                          {
00069     return theState.weight();
00070   }

double GaussianSumUtilities1D::weight ( unsigned int i ) const [inline]

weight of a component

Definition at line 34 of file GaussianSumUtilities1D.h.

References components().

00034 {return components()[i].weight();}

SingleGaussianState1D GaussianSumUtilities1D::theMode [mutable, private]

Definition at line 120 of file GaussianSumUtilities1D.h.

Referenced by computeMode(), and mode().

ModeStatus GaussianSumUtilities1D::theModeStatus [mutable, private]

Definition at line 119 of file GaussianSumUtilities1D.h.

Referenced by computeMode(), mode(), and modeIsValid().

const MultiGaussianState1D& GaussianSumUtilities1D::theState [private]

Definition at line 116 of file GaussianSumUtilities1D.h.

Referenced by components(), mean(), variance(), and weight().

GaussianSumUtilities1D Class Reference

Public Member Functions

Private Types

Private Member Functions

Private Attributes

Detailed Description

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


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
	Mode "state": mean = mode, variance = local variance at mode, weight chosen to have pdf(mode) equal to the one of the mixture.
bool	modeIsValid () const
	mode status
double	pdf (const 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 () const
	combined covariance
double	variance (unsigned int i) const
	variance of a component
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 (const double &, const std::vector< double > &) const
	first derivative of ln(pdf) using the pdf components at the evaluation point
double	d1Pdf (const double &, const std::vector< double > &) const
	first derivative of the p.d.f. using the pdf components at the evaluation point
double	d2LnPdf (const double &, const std::vector< double > &) const
	second derivative of ln(pdf) using the pdf components at the evaluation point
double	d2Pdf (const double &, const std::vector< double > &) const
	second derivative of the p.d.f. using the pdf components at the evaluation point
double	d3Pdf (const 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
	Finds mode.
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.
double	lnPdf (const double &, const std::vector< double > &) const
	ln(pdf) using the pdf components at the evaluation point
double	localVariance (const double &x) const
	Local variance from Hessian matrix.
double	pdf (const double &, const std::vector< double > &) const
	value of the p.d.f. using the pdf components at the evaluation point
std::vector< double >	pdfComponents (const double &) const
	pdf components
Private Attributes
SingleGaussianState1D	theMode
ModeStatus	theModeStatus
const MultiGaussianState1D &	theState