Logo ROOT  
Reference Guide
TKDE.cxx
Go to the documentation of this file.
1// // @(#)root/hist:$Id$
2// Authors: Bartolomeu Rabacal 07/2010
3/**********************************************************************
4 * *
5 * Copyright (c) 2006 , ROOT MathLib Team *
6 * *
7 * For the licensing terms see $ROOTSYS/LICENSE. *
8 * For the list of contributors see $ROOTSYS/README/CREDITS. *
9 * *
10 **********************************************************************/
11//////////////////////////////////////////////////////////////////////////////
12/** \class TKDE
13 \ingroup Hist
14 Kernel Density Estimation class.
15 The three main references are:
16 1. "Scott DW, Multivariate Density Estimation. Theory, Practice and Visualization. New York: Wiley",
17 2. "Jann Ben - ETH Zurich, Switzerland -, Univariate kernel density estimation document for KDENS:
18 Stata module for univariate kernel density estimation."
19 3. "Hardle W, Muller M, Sperlich S, Werwatz A, Nonparametric and Semiparametric Models. Springer."
20 4. "Cranmer KS, Kernel Estimation in High-Energy
21 Physics. Computer Physics Communications 136:198-207,2001" - e-Print Archive: hep ex/0011057.
22
23 The algorithm is briefly described in (4). A binned version is also implemented to address the
24 performance issue due to its data size dependance.
25 */
26
27
28#include <functional>
29#include <algorithm>
30#include <numeric>
31#include <limits>
32#include <cassert>
33
34#include "Math/Error.h"
35#include "TMath.h"
36#include "Math/Functor.h"
37#include "Math/Integrator.h"
40#include "TGraphErrors.h"
41#include "TF1.h"
42#include "TH1.h"
43#include "TVirtualPad.h"
44#include "TKDE.h"
45
46
48
49class TKDE::TKernel {
50 TKDE* fKDE;
51 UInt_t fNWeights; // Number of kernel weights (bandwidth as vectorized for binning)
52 std::vector<Double_t> fWeights; // Kernel weights (bandwidth)
53public:
54 TKernel(Double_t weight, TKDE* kde);
55 void ComputeAdaptiveWeights();
57 Double_t GetWeight(Double_t x) const;
59 const std::vector<Double_t> & GetAdaptiveWeights() const;
60};
61
63 enum EIntegralResult{kNorm, kMu, kSigma2, kUnitIntegration};
64 KernelIntegrand(const TKDE* kde, EIntegralResult intRes);
66private:
67 const TKDE* fKDE;
68 EIntegralResult fIntegralResult;
69};
70
71/// default constructor used by I/O
73 fKernelFunction(nullptr),
74 fKernel(nullptr),
75 fPDF(nullptr),
76 fUpperPDF(nullptr),
77 fLowerPDF(nullptr),
78 fApproximateBias(nullptr),
79 fGraph(nullptr),
80 fUseMirroring(false), fMirrorLeft(false), fMirrorRight(false), fAsymLeft(false), fAsymRight(false),
81 fUseBins(false), fNewData(false), fUseMinMaxFromData(false),
82 fNBins(0), fNEvents(0), fSumOfCounts(0), fUseBinsNEvents(0),
83 fMean(0.),fSigma(0.), fSigmaRob(0.), fXMin(0.), fXMax(0.),
84 fRho(0.), fAdaptiveBandwidthFactor(0.), fWeightSize(0)
85{
86}
87
89 //Class destructor
90 if (fPDF) delete fPDF;
91 if (fUpperPDF) delete fUpperPDF;
92 if (fLowerPDF) delete fLowerPDF;
93 if (fGraph) delete fGraph;
96 if (fKernel) delete fKernel;
97}
98
99void TKDE::Instantiate(KernelFunction_Ptr kernfunc, UInt_t events, const Double_t* data, const Double_t* dataWeights, Double_t xMin, Double_t xMax, const Option_t* option, Double_t rho) {
100 // Template's constructor surrogate
101
102 fData = std::vector<Double_t>(events, 0.0);
103 fEvents = std::vector<Double_t>(events, 0.0);
104 fPDF = 0;
105 fKernel = 0;
106 fKernelFunction = 0;
107 fUpperPDF = 0;
108 fLowerPDF = 0;
110 fGraph = 0;
111 fNewData = false;
112 fUseMirroring = false; fMirrorLeft = false; fMirrorRight = false;
113 fAsymLeft = false; fAsymRight = false;
114 fNBins = events < 10000 ? 100 : events / 10;
115 fNEvents = events;
116 fUseBinsNEvents = 10000;
117 fMean = 0.0;
118 fSigma = 0.0;
119 fXMin = xMin;
120 fXMax = xMax;
122 fSumOfCounts = 0;
124 fRho = rho;
125 fWeightSize = 0;
126 fCanonicalBandwidths = std::vector<Double_t>(kTotalKernels, 0.0);
127 fKernelSigmas2 = std::vector<Double_t>(kTotalKernels, -1.0);
128 fSettedOptions = std::vector<Bool_t>(4, kFALSE);
129 SetOptions(option, rho);
131 SetMirror();
132 SetUseBins();
133 SetData(data, dataWeights);
134 SetKernelFunction(kernfunc);
135
136}
137
138void TKDE::SetOptions(const Option_t* option, Double_t rho) {
139 //Sets User defined construction options
140 TString opt = option;
141 opt.ToLower();
142 std::string options = opt.Data();
143 size_t numOpt = 4;
144 std::vector<std::string> voption(numOpt, "");
145 for (std::vector<std::string>::iterator it = voption.begin(); it != voption.end() && !options.empty(); ++it) {
146 size_t pos = options.find_last_of(';');
147 if (pos == std::string::npos) {
148 *it = options;
149 break;
150 }
151 *it = options.substr(pos + 1);
152 options = options.substr(0, pos);
153 }
154 for (std::vector<std::string>::iterator it = voption.begin(); it != voption.end(); ++it) {
155 size_t pos = (*it).find(':');
156 if (pos != std::string::npos) {
157 GetOptions((*it).substr(0, pos) , (*it).substr(pos + 1));
158 }
159 }
161 fRho = rho;
162}
163
164void TKDE::SetDrawOptions(const Option_t* option, TString& plotOpt, TString& drawOpt) {
165 // Sets User defined drawing options
166 size_t numOpt = 2;
167 std::string options = TString(option).Data();
168 std::vector<std::string> voption(numOpt, "");
169 for (std::vector<std::string>::iterator it = voption.begin(); it != voption.end() && !options.empty(); ++it) {
170 size_t pos = options.find_last_of(';');
171 if (pos == std::string::npos) {
172 *it = options;
173 break;
174 }
175 *it = options.substr(pos + 1);
176 options = options.substr(0, pos);
177 }
178 Bool_t foundPlotOPt = kFALSE;
179 Bool_t foundDrawOPt = kFALSE;
180 for (std::vector<std::string>::iterator it = voption.begin(); it != voption.end() && !options.empty(); ++it) {
181 size_t pos = (*it).find(':');
182 if (pos == std::string::npos) break;
183 TString optionType = (*it).substr(0, pos);
184 TString optionInstance = (*it).substr(pos + 1);
185 optionType.ToLower();
186 optionInstance.ToLower();
187 if (optionType.Contains("plot")) {
188 foundPlotOPt = kTRUE;
189 if (optionInstance.Contains("estimate") || optionInstance.Contains("errors") || optionInstance.Contains("confidenceinterval"))
190 plotOpt = optionInstance;
191 else
192 this->Warning("SetDrawOptions", "Unknown plotting option: setting to KDE estimate plot.");
193 } else if (optionType.Contains("drawoptions")) {
194 foundDrawOPt = kTRUE;
195 drawOpt = optionInstance;
196 }
197 }
198 if (!foundPlotOPt) {
199 this->Warning("SetDrawOptions", "No plotting option: setting to KDE estimate plot.");
200 plotOpt = "estimate";
201 }
202 if (!foundDrawOPt) {
203 this->Warning("SetDrawOptions", "No drawing options: setting to default ones.");
204 drawOpt = "apl4";
205 }
206}
207
208void TKDE::GetOptions(std::string optionType, std::string option) {
209 // Gets User defined KDE construction options
210 if (optionType.compare("kerneltype") == 0) {
212 if (option.compare("gaussian") == 0) {
214 } else if (option.compare("epanechnikov") == 0) {
216 } else if (option.compare("biweight") == 0) {
218 } else if (option.compare("cosinearch") == 0) {
220 } else if (option.compare("userdefined") == 0) {
222 } else {
223 this->Warning("GetOptions", "Unknown kernel type option: setting to Gaussian");
225 }
226 } else if (optionType.compare("iteration") == 0) {
228 if (option.compare("adaptive") == 0) {
230 } else if (option.compare("fixed") == 0) {
232 } else {
233 this->Warning("GetOptions", "Unknown iteration option: setting to Adaptive");
235 }
236 } else if (optionType.compare("mirror") == 0) {
238 if (option.compare("nomirror") == 0) {
240 } else if (option.compare("mirrorleft") == 0) {
242 } else if (option.compare("mirrorright") == 0) {
244 } else if (option.compare("mirrorboth") == 0) {
246 } else if (option.compare("mirrorasymleft") == 0) {
248 } else if (option.compare("mirrorasymleftright") == 0) {
250 } else if (option.compare("mirrorasymright") == 0) {
252 } else if (option.compare("mirrorleftasymright") == 0) {
254 } else if (option.compare("mirrorasymboth") == 0) {
256 } else {
257 this->Warning("GetOptions", "Unknown mirror option: setting to NoMirror");
259 }
260 } else if (optionType.compare("binning") == 0) {
262 if (option.compare("unbinned") == 0) {
264 } else if (option.compare("relaxedbinning") == 0) {
266 } else if (option.compare("forcedbinning") == 0) {
268 } else {
269 this->Warning("GetOptions", "Unknown binning option: setting to RelaxedBinning");
271 }
272 }
273}
274
276 // Sets missing construction options to default ones
277 if (!fSettedOptions[0]) {
279 }
280 if (!fSettedOptions[1]) {
282 }
283 if (!fSettedOptions[2]) {
285 }
286 if (!fSettedOptions[3]) {
288 }
289}
290
291void TKDE::CheckOptions(Bool_t isUserDefinedKernel) {
292 // Sets User global options
293 if (!(isUserDefinedKernel) && !(fKernelType >= kGaussian && fKernelType < kUserDefined)) {
294 Error("CheckOptions", "Illegal user kernel type input! Use template constructor for user defined kernel.");
295 }
296 if (fIteration != kAdaptive && fIteration != kFixed) {
297 Warning("CheckOptions", "Illegal user iteration type input - use default value !");
299 }
300 if (!(fMirror >= kNoMirror && fMirror <= kMirrorAsymBoth)) {
301 Warning("CheckOptions", "Illegal user mirroring type input - use default value !");
303 }
304 if (!(fBinning >= kUnbinned && fBinning <= kForcedBinning)) {
305 Warning("CheckOptions", "Illegal user binning type input - use default value !");
307 }
308 if (fRho <= 0.0) {
309 Warning("CheckOptions", "Tuning factor rho cannot be non-positive - use default value !");
310 fRho = 1.0;
311 }
312}
313
315 // Sets User option for the choice of kernel estimator
317 delete fKernelFunction;
318 fKernelFunction = 0;
319 }
320 fKernelType = kern;
321 CheckOptions();
323}
324
326 // Sets User option for fixed or adaptive iteration
327 fIteration = iter;
328 CheckOptions();
329 SetKernel();
330}
331
332
334 // Sets User option for mirroring the data
335 fMirror = mir;
336 CheckOptions();
337 SetMirror();
338 if (fUseMirroring) {
340 }
341 SetKernel();
342}
343
344
346 // Sets User option for binning the weights
347 fBinning = bin;
348 CheckOptions();
349 SetUseBins();
350}
351
353 // Sets User option for number of bins
354 if (!nbins) {
355 Error("SetNBins", "Number of bins must be greater than zero.");
356 return;
357 }
358
359 fNBins = nbins;
361
362 SetUseBins();
363 if (!fUseBins) {
364 if (fBinning == kUnbinned)
365 Warning("SetNBins", "Bin type using SetBinning must be set for using a binned evaluation");
366 else
367 Warning("SetNBins", "Bin type using SetBinning or with SetUseBinsNEvents must be set for using a binned evaluation");
368 }
369}
370
372 // Sets User option for the minimum number of events for allowing automatic binning
373 fUseBinsNEvents = nEvents;
374 SetUseBins();
375}
376
378 // Factor which can be used to tune the smoothing.
379 // It is used as multiplicative factor for the fixed and adaptive bandwidth.
380 // A value < 1 will reproduce better the tails but oversmooth the peak
381 // while a factor > 1 will overestimate the tail
382 fRho = rho;
383 CheckOptions();
384 SetKernel();
385}
386
388 // Sets minimum range value and maximum range value
389 if (xMin >= xMax) {
390 Error("SetRange", "Minimum range cannot be bigger or equal than the maximum range! Present range values remain the same.");
391 return;
392 }
393 fXMin = xMin;
394 fXMax = xMax;
395 fUseMinMaxFromData = false;
396 SetKernel();
397}
398
399// private methods
400
402 // Sets User option for using binned weights
403 switch (fBinning) {
404 default:
405 case kRelaxedBinning:
406 if (fNEvents >= fUseBinsNEvents) {
407 fUseBins = kTRUE;
408 } else {
410 }
411 break;
412 case kForcedBinning:
413 fUseBins = kTRUE;
414 break;
415 case kUnbinned:
417 }
418
419 // set the binning
420 if (fUseBins && !fEvents.empty() ) {
423 SetKernel();
424 }
425}
426
428 // Sets the mirroring
434}
435
436void TKDE::SetData(const Double_t* data, const Double_t* wgts) {
437 // Sets the data events input sample or bin centres for binned option and computes basic estimators
438 if (!data) {
439 if (fNEvents) fData.reserve(fNEvents);
440 return;
441 }
442 fEvents.assign(data, data + fNEvents);
443 if (wgts) fEventWeights.assign(wgts, wgts + fNEvents);
444
445 if (fUseMinMaxFromData) {
446 fXMin = *std::min_element(fEvents.begin(), fEvents.end());
447 fXMax = *std::max_element(fEvents.begin(), fEvents.end());
448 }
449
450 if (fUseBins) {
451 if (fNBins >= fNEvents) {
452 this->Warning("SetData", "Default number of bins is greater or equal to number of events. Use SetNBins(UInt_t) to set the appropriate number of bins");
453 }
456 } else {
458 fData = fEvents;
459 }
460 // to set fBinCount and fSumOfCounts
462
463
465 if (fUseMirroring) {
467 }
468}
469
471 // re-initialize the KDE in case of new data or in case of reading from a file
472
473 // in case of new data re-compute the statistics (works only for unbinned case)
474 if (fNewData) {
476 return;
477 }
478 // case of reading from a file.
479 // we need to recreate Kernel class (with adaptive weights if needed) and
480 // recreate kernel function pointer
481
482 if (fKernelFunction) Error("ReInit","Kernel function pointer should be a nullptr when re-initializing after reading from a file");
483
484 if (fEvents.size() == 0) {
485 Error("ReInit","TKDE does not contain any data !");
486 return;
487 }
488
490
491 SetKernel();
492}
493
495 // re-initialize when new data have been filled in TKDE
496 // re-compute kernel quantities and statistics
497 if (fUseBins) {
498 Error("InitFromNewData","Re-felling is not supported with binning");
499 return;
500 }
501
502 fNewData = false;
503 // in case of unbinned data
504 if (!fUseBins) fEvents = fData;
505 if (fUseMinMaxFromData) {
506 fXMin = *std::min_element(fEvents.begin(), fEvents.end());
507 fXMax = *std::max_element(fEvents.begin(), fEvents.end());
508 }
510 // if (fUseBins) {
511 // } // bin usage is not supported in this case
512 //
514 if (fUseMirroring) {
516 }
517 // in case of I/O reset the kernel
518}
519
521 // Mirrors the data
522 std::vector<Double_t> originalEvents = fEvents;
523 std::vector<Double_t> originalWeights = fEventWeights;
524 if (fMirrorLeft) {
525 fEvents.resize(2 * fNEvents, 0.0);
526 transform(fEvents.begin(), fEvents.begin() + fNEvents, fEvents.begin() + fNEvents,
527 std::bind(std::minus<Double_t>(), 2 * fXMin, std::placeholders::_1));
528 }
529 if (fMirrorRight) {
530 fEvents.resize((fMirrorLeft + 2) * fNEvents, 0.0);
531 transform(fEvents.begin(), fEvents.begin() + fNEvents, fEvents.begin() + (fMirrorLeft + 1) * fNEvents,
532 std::bind(std::minus<Double_t>(), 2 * fXMax, std::placeholders::_1));
533 }
534 if (!fEventWeights.empty() && (fMirrorLeft || fMirrorRight)) {
535 // copy weights too
536 fEventWeights.insert(fEventWeights.end(), fEventWeights.begin(), fEventWeights.end() );
537 }
538
539 if(fUseBins) {
544 } else {
545 fData = fEvents;
546 }
548
549 fEvents = originalEvents;
550 fEventWeights = originalWeights;
551}
552
554 // Computes input data's mean
555 fMean = std::accumulate(fEvents.begin(), fEvents.end(), 0.0) / fEvents.size();
556}
557
559 // Computes input data's sigma
560 fSigma = std::sqrt(1. / (fEvents.size() - 1.) * (std::inner_product(fEvents.begin(), fEvents.end(), fEvents.begin(), 0.0) - fEvents.size() * std::pow(fMean, 2.)));
561 fSigmaRob = std::min(fSigma, R / 1.349); // Sigma's robust estimator
562}
563
565 // Sets the kernel density estimator
566 UInt_t n = fData.size();
567 if (n == 0) return;
568 // Optimal bandwidth (Silverman's rule of thumb with assumed Gaussian density)
569 Double_t weight = fCanonicalBandwidths[kGaussian] * fSigmaRob * std::pow(3. / (8. * std::sqrt(M_PI)) * n, -0.2);
571 if (fKernel) delete fKernel;
572 fKernel = new TKernel(weight, this);
573 if (fIteration == kAdaptive) {
574 fKernel->ComputeAdaptiveWeights();
575 }
576 //std::cout << "setting the kernel - n = " << n << " weight is " << weight << " " << fRho << " " << fSigmaRob << " " << fSigma << " " << fMean << " " << fCanonicalBandwidths[kGaussian] << std::endl;
577}
578
580
581 assert(fKernelFunction == 0); // to avoid memory leaks
582 switch (fKernelType) {
583 case kGaussian :
585 break;
586 case kEpanechnikov :
588 break;
589 case kBiweight :
591 break;
592 case kCosineArch :
594 break;
595 case kUserDefined :
596 fKernelFunction = kernfunc;
598 break;
599 case kTotalKernels :
600 default:
601 /// for user defined kernels
602 fKernelFunction = kernfunc;
604 }
605
606 if (fKernelType == kUserDefined) {
607 if (fKernelFunction) {
611 }
612 else {
613 Error("SetKernelFunction", "User kernel function is not defined !");
614 return;
615 }
616 }
617 assert(fKernelFunction);
620 SetKernel();
621}
622
624 // Sets the canonical bandwidths according to the kernel type
625 fCanonicalBandwidths[kGaussian] = 0.7764; // Checked in Mathematica
626 fCanonicalBandwidths[kEpanechnikov] = 1.7188; // Checked in Mathematica
627 fCanonicalBandwidths[kBiweight] = 2.03617; // Checked in Mathematica
628 fCanonicalBandwidths[kCosineArch] = 1.7663; // Checked in Mathematica
629 fCanonicalBandwidths[kUserDefined] = 1.0; // To be Checked
630}
631
633 // Sets the kernel sigmas2 according to the kernel type
635 fKernelSigmas2[kEpanechnikov] = 1.0 / 5.0;
636 fKernelSigmas2[kBiweight] = 1.0 / 7.0;
637 fKernelSigmas2[kCosineArch] = 1.0 - 8.0 / std::pow(M_PI, 2);
638}
639
641 // Returns the PDF estimate as a function sampled in npx between xMin and xMax
642 // the KDE is not re-normalized to the xMin/xMax range.
643 // The user manages the returned function
644 // For getting a non-sampled TF1, one can create a TF1 directly from the TKDE by doing
645 // TF1 * f1 = new TF1("f1",kde,xMin,xMax,0);
646 return GetKDEFunction(npx,xMin,xMax);
647}
648
649TF1* TKDE::GetUpperFunction(Double_t confidenceLevel, UInt_t npx, Double_t xMin, Double_t xMax) {
650 // Returns the PDF upper estimate (upper confidence interval limit)
651 return GetPDFUpperConfidenceInterval(confidenceLevel,npx,xMin,xMax);
652}
653
654TF1* TKDE::GetLowerFunction(Double_t confidenceLevel, UInt_t npx, Double_t xMin, Double_t xMax) {
655 // Returns the PDF lower estimate (lower confidence interval limit)
656 return GetPDFLowerConfidenceInterval(confidenceLevel,npx,xMin,xMax);
657}
658
660 // Returns the PDF estimate bias
661 return GetKDEApproximateBias(npx,xMin,xMax);
662}
663
665 // Fills data member with User input data event for the unbinned option
666 if (fUseBins) {
667 this->Warning("Fill", "Cannot fill data with data binned option. Data input ignored.");
668 return;
669 }
670 fData.push_back(data);
671 fNEvents++;
672 fNewData = kTRUE;
673}
674
675void TKDE::Fill(Double_t data, Double_t weight) {
676 // Fills data member with User input data event for the unbinned option
677 if (fUseBins) {
678 this->Warning("Fill", "Cannot fill data with data binned option. Data input ignored.");
679 return;
680 }
681 fData.push_back(data); // should not be here fEvent ??
682 fEventWeights.push_back(weight);
683 fNEvents++;
684 fNewData = kTRUE;
685}
686
688 // The class's unary function: returns the kernel density estimate
689 return (*this)(*x);
690}
691
693 // The class's unary function: returns the kernel density estimate
694 if (!fKernel) {
695 (const_cast<TKDE*>(this))->ReInit();
696 // in case of failed re-initialization
697 if (!fKernel) return TMath::QuietNaN();
698 }
699 return (*fKernel)(x);
700}
701
703 // return the mean of the data
704 if (fNewData) (const_cast<TKDE*>(this))->InitFromNewData();
705 return fMean;
706}
707
709 // return the standard deviation of the data
710 if (fNewData) (const_cast<TKDE*>(this))->InitFromNewData();
711 return fSigma;
712}
713
715 // Returns the Root Asymptotic Mean Integrated Squared Error according to Silverman's rule of thumb with assumed Gaussian density
717 return std::sqrt(result);
718}
719
720TKDE::TKernel::TKernel(Double_t weight, TKDE* kde) :
721// Internal class constructor
722fKDE(kde),
723fNWeights(kde->fData.size()),
724fWeights(fNWeights, weight)
725{}
726
727void TKDE::TKernel::ComputeAdaptiveWeights() {
728 // Gets the adaptive weights (bandwidths) for TKernel internal computation
729 std::vector<Double_t> weights = fWeights;
730 Double_t minWeight = weights[0] * 0.05;
731 unsigned int n = fKDE->fData.size();
732 assert( n == weights.size() );
733 bool useDataWeights = (fKDE->fBinCount.size() == n);
734 Double_t f = 0.0;
735 for (unsigned int i = 0; i < n; ++i) {
736// for (; weight != weights.end(); ++weight, ++data, ++dataW) {
737 if (useDataWeights && fKDE->fBinCount[i] <= 0) continue; // skip negative or null weights
738 f = (*fKDE->fKernel)(fKDE->fData[i]);
739 if (f <= 0)
740 fKDE->Warning("ComputeAdativeWeights","function value is zero or negative for x = %f w = %f",
741 fKDE->fData[i],(useDataWeights) ? fKDE->fBinCount[i] : 1.);
742 weights[i] = std::max(weights[i] /= std::sqrt(f), minWeight);
743 fKDE->fAdaptiveBandwidthFactor += std::log(f);
744 //printf("(f = %f w = %f af = %f ),",f,*weight,fKDE->fAdaptiveBandwidthFactor);
745 }
746 Double_t kAPPROX_GEO_MEAN = 0.241970724519143365; // 1 / TMath::Power(2 * TMath::Pi(), .5) * TMath::Exp(-.5). Approximated geometric mean over pointwise data (the KDE function is substituted by the "real Gaussian" pdf) and proportional to sigma. Used directly when the mirroring is enabled, otherwise computed from the data
747 fKDE->fAdaptiveBandwidthFactor = fKDE->fUseMirroring ? kAPPROX_GEO_MEAN / fKDE->fSigmaRob : std::sqrt(std::exp(fKDE->fAdaptiveBandwidthFactor / fKDE->fData.size()));
748 transform(weights.begin(), weights.end(), fWeights.begin(),
749 std::bind(std::multiplies<Double_t>(), std::placeholders::_1, fKDE->fAdaptiveBandwidthFactor));
750 //printf("adaptive bandwidth factor % f weight 0 %f , %f \n",fKDE->fAdaptiveBandwidthFactor, weights[0],fWeights[0] );
751}
752
753Double_t TKDE::TKernel::GetWeight(Double_t x) const {
754 // Returns the bandwidth
755 return fWeights[fKDE->Index(x)];
756}
757
759 // Returns the bins' centres from the data for using with the binned option
760 fData.assign(fNBins, 0.0);
761 Double_t binWidth((xmax - xmin) / fNBins);
762 for (UInt_t i = 0; i < fNBins; ++i) {
763 fData[i] = xmin + (i + 0.5) * binWidth;
764 }
765}
766
768 // Returns the bins' count from the data for using with the binned option
769 // or set the bin count to the weights in case of weighted data
770 if (fUseBins) {
771 fBinCount.resize(fNBins);
772 fSumOfCounts = 0;
773 // case of weighted events
774 if (!fEventWeights.empty() ) {
775 for (UInt_t i = 0; i < fNEvents; ++i) {
776 if (fEvents[i] >= fXMin && fEvents[i] < fXMax) {
779 //printf("sum of counts %f - bin count %d - %f \n",fSumOfCounts, Index(fEvents[i]), fBinCount[Index(fEvents[i])] );
780 }
781 }
782 }
783 // case of unweighted data
784 else {
785 for (UInt_t i = 0; i < fNEvents; ++i) {
786 if (fEvents[i] >= fXMin && fEvents[i] < fXMax) {
787 fBinCount[Index(fEvents[i])] += 1;
788 fSumOfCounts += 1;
789 }
790 }
791 }
792 }
793 else if (!fEventWeights.empty() ) {
795 fSumOfCounts = 0;
796 for (UInt_t i = 0; i < fNEvents; ++i)
798 }
799 else {
801 fBinCount.clear();
802 }
803}
804
805void TKDE::Draw(const Option_t* opt) {
806 // Draws either the KDE functions or its errors
807 // Possible options:
808 // "" (default) - draw just the kde
809 // "same" draw on top of existing pad
810 // "Errors" draw a TGraphErrors with the point and errors
811 // "confidenceinterval" draw KDE + conf interval functions (default is 95%)
812 // "confidenceinterval@0.90" draw KDE + conf interval functions at 90%
813 // Extra options can be passed in opt for drawing the TF1 or the TGraph
814 //
815 //NOTE: The functions GetDrawnFunction(), GetDrawnUpperFunction(), GetDrawnLowerFunction()
816 // and GetGraphWithErrors() return the corresponding drawn objects (which are maneged by the TKDE)
817 // They can be used to changes style, color, etc...
818
819 // TString plotOpt = "";
820 // TString drawOpt = "";
821 // LM : this is too complicates - skip it - not needed for just
822 // three options
823 // SetDrawOptions(opt, plotOpt, drawOpt);
824 TString plotOpt = opt;
825 plotOpt.ToLower();
826 TString drawOpt = plotOpt;
827 if(gPad && !plotOpt.Contains("same")) {
828 gPad->Clear();
829 }
830 if (plotOpt.Contains("errors")) {
831 drawOpt.ReplaceAll("errors","");
832 DrawErrors(drawOpt);
833 }
834 else if (plotOpt.Contains("confidenceinterval") ||
835 plotOpt.Contains("confinterval")) {
836 // parse level option
837 drawOpt.ReplaceAll("confidenceinterval","");
838 drawOpt.ReplaceAll("confinterval","");
839 Double_t level = 0.95;
840 const char * s = strstr(plotOpt.Data(),"interval@");
841 // coverity [secure_coding : FALSE]
842 if (s != 0) sscanf(s,"interval@%lf",&level);
843 if((level <= 0) || (level >= 1)) {
844 Warning("Draw","given confidence level %.3lf is invalid - use default 0.95",level);
845 level = 0.95;
846 }
847 DrawConfidenceInterval(drawOpt,level);
848 }
849 else {
850 if (fPDF) delete fPDF;
852 fPDF->Draw(drawOpt);
853 }
854}
855
856void TKDE::DrawErrors(TString& drawOpt) {
857 // Draws a TGraphErrors for the KDE errors
858 if (fGraph) delete fGraph;
860 fGraph->Draw(drawOpt.Data());
861}
862
864 if (xmin>= xmax) { xmin = fXMin; xmax = fXMax; }
865 // return a TGraphErrors for the KDE errors
866 UInt_t n = npx;
867 Double_t* x = new Double_t[n + 1];
868 Double_t* ex = new Double_t[n + 1];
869 Double_t* y = new Double_t[n + 1];
870 Double_t* ey = new Double_t[n + 1];
871 for (UInt_t i = 0; i <= n; ++i) {
872 x[i] = xmin + i * (xmax - xmin) / n;
873 y[i] = (*this)(x[i]);
874 ex[i] = 0;
875 ey[i] = this->GetError(x[i]);
876 }
877 TGraphErrors* ge = new TGraphErrors(n, &x[0], &y[0], &ex[0], &ey[0]);
878 ge->SetName("kde_graph_error");
879 ge->SetTitle("Errors");
880 delete [] x;
881 delete [] ex;
882 delete [] y;
883 delete [] ey;
884 return ge;
885}
886
887void TKDE::DrawConfidenceInterval(TString& drawOpt,double cl) {
888 // Draws the KDE and its confidence interval
889 GetKDEFunction()->Draw(drawOpt.Data());
891 upper->SetLineColor(kBlue);
892 upper->Draw(("same" + drawOpt).Data());
894 lower->SetLineColor(kRed);
895 lower->Draw(("same" + drawOpt).Data());
896 if (fUpperPDF) delete fUpperPDF;
897 if (fLowerPDF) delete fLowerPDF;
898 fUpperPDF = upper;
899 fLowerPDF = lower;
900}
901
903 // Returns the bandwidth for the non adaptive KDE
904 Double_t result = -1.;
906 this->Warning("GetFixedWeight()", "Fixed iteration option not enabled. Returning %f.", result);
907 } else {
908 result = fKernel->GetFixedWeight();
909 }
910 return result;
911}
912
914 // Returns the bandwidths for the adaptive KDE
916 this->Warning("GetFixedWeight()", "Adaptive iteration option not enabled. Returning a NULL pointer<");
917 return 0;
918 }
919 if (fNewData) (const_cast<TKDE*>(this))->InitFromNewData();
920 return &(fKernel->GetAdaptiveWeights()).front();
921}
922
923Double_t TKDE::TKernel::GetFixedWeight() const {
924 // Returns the bandwidth for the non adaptive KDE
925 return fWeights[0];
926}
927
928const std::vector<Double_t> & TKDE::TKernel::GetAdaptiveWeights() const {
929 // Returns the bandwidth for the non adaptive KDE
930 return fWeights;
931}
932
934 // The internal class's unary function: returns the kernel density estimate
935 Double_t result(0.0);
936 UInt_t n = fKDE->fData.size();
937 // case of bins or weighted data
938 Bool_t useBins = (fKDE->fBinCount.size() == n);
939 Double_t nSum = (useBins) ? fKDE->fSumOfCounts : fKDE->fNEvents;
940 // double dmin = 1.E10;
941 // double xmin,bmin,wmin;
942 for (UInt_t i = 0; i < n; ++i) {
943 Double_t binCount = (useBins) ? fKDE->fBinCount[i] : 1.0;
944 result += binCount / fWeights[i] * (*fKDE->fKernelFunction)((x - fKDE->fData[i]) / fWeights[i]);
945 if (fKDE->fAsymLeft) {
946 result -= binCount / fWeights[i] * (*fKDE->fKernelFunction)((x - (2. * fKDE->fXMin - fKDE->fData[i])) / fWeights[i]);
947 }
948 if (fKDE->fAsymRight) {
949 result -= binCount / fWeights[i] * (*fKDE->fKernelFunction)((x - (2. * fKDE->fXMax - fKDE->fData[i])) / fWeights[i]);
950 }
951 // if ( TMath::IsNaN(result) ) {
952 // printf("event %i count %f weight %f data % f x %f \n",i,binCount,fWeights[i],fKDE->fData[i],x );
953 // }
954 // if ( result <= 0 ) {
955 // printf("event %i count %f weight %f data % f x %f \n",i,binCount,fWeights[i],fKDE->fData[i],x );
956 // }
957 // if (std::abs(x - fKDE->fData[i]) < dmin ) {
958 // xmin = x;
959 // bmin = binCount;
960 // wmin = fWeights[i];
961 // dmin = std::abs(x - fKDE->fData[i]);
962 // }
963 // if (i < fKDE->fEvents.size() )
964 // printf("data point %i %f %f count %f weight % f result % f\n",i,fKDE->fData[i],fKDE->fEvents[i],binCount,fWeights[i], result);
965 }
966 if ( TMath::IsNaN(result) ) {
967 fKDE->Warning("operator()","Result is NaN for x %f \n",x);
968 //xmin % f , %f, %f \n",result,x,xmin,bmin,wmin );
969 }
970 return result / nSum;
971}
972
974 // Returns the indices (bins) for the binned weights
975 Int_t bin = Int_t((x - fXMin) * fWeightSize);
976 if (bin == (Int_t)fData.size()) return --bin;
978 bin += fData.size() / (fMirrorLeft + fMirrorRight + 1);
979 }
980 if (bin > (Int_t)fData.size()) {
981 return (Int_t)(fData.size()) - 1;
982 } else if (bin <= 0) {
983 return 0;
984 }
985 return bin;
986}
987
989 // Returns the pointwise upper estimated density
990 Double_t f = (*this)(x);
992 Double_t prob = 1. - (1.-*p)/2; // this is 1.-alpha/2
994 return f + z * sigma;
995}
996
998 // Returns the pointwise lower estimated density
999 Double_t f = (*this)(x);
1001 Double_t prob = (1.-*p)/2; // this is alpha/2
1003 return f + z * sigma;
1004}
1005
1006
1008 // Returns the pointwise approximate estimated density bias
1011 rd.SetFunction(kern);
1012 Double_t df2 = rd.Derivative2(x);
1013 Double_t weight = fKernel->GetWeight(x); // Bandwidth
1014 return 0.5 * fKernelSigmas2[fKernelType] * std::pow(weight, 2) * df2;
1015}
1017 // Returns the pointwise sigma of estimated density
1018 Double_t kernelL2Norm = ComputeKernelL2Norm();
1019 Double_t f = (*this)(x);
1020 Double_t weight = fKernel->GetWeight(x); // Bandwidth
1021 Double_t result = f * kernelL2Norm / (fNEvents * weight);
1022 return std::sqrt(result);
1023}
1024
1026 // Checks if kernel has unit integral, mu = 0 and positive finite sigma conditions
1027 Double_t valid = kTRUE;
1029 valid = valid && unity == 1.;
1030 if (!valid) {
1031 Error("CheckKernelValidity", "Kernel's integral is %f",unity);
1032 }
1034 valid = valid && mu == 0.;
1035 if (!valid) {
1036 Error("CheckKernelValidity", "Kernel's mu is %f" ,mu);
1037 }
1038 Double_t sigma2 = ComputeKernelSigma2();
1039 valid = valid && sigma2 > 0 && sigma2 != std::numeric_limits<Double_t>::infinity();
1040 if (!valid) {
1041 Error("CheckKernelValidity", "Kernel's sigma2 is %f",sigma2);
1042 }
1043 if (!valid) {
1044 Error("CheckKernelValidity", "Validation conditions: the kernel's integral must be 1, the kernel's mu must be zero and the kernel's sigma2 must be finite positive to be a suitable kernel.");
1045 //exit(EXIT_FAILURE);
1046 }
1047}
1048
1050 // Computes the kernel's L2 norm
1053 ig.SetFunction(kernel);
1054 Double_t result = ig.Integral();
1055 return result;
1056}
1057
1059 // Computes the kernel's sigma squared
1061 KernelIntegrand kernel(this, TKDE::KernelIntegrand::kSigma2);
1062 ig.SetFunction(kernel);
1063 Double_t result = ig.Integral();
1064 return result;
1065}
1066
1068 // Computes the kernel's mu
1070 KernelIntegrand kernel(this, TKDE::KernelIntegrand::kMu);
1071 ig.SetFunction(kernel);
1072 Double_t result = ig.Integral();
1073 return result;
1074}
1075
1077 // Computes the kernel's integral which ought to be unity
1079 KernelIntegrand kernel(this, TKDE::KernelIntegrand::kUnitIntegration);
1080 ig.SetFunction(kernel);
1081 Double_t result = ig.Integral();
1082 return result;
1083}
1084
1086 /// in case of weights use
1087 if (!fEventWeights.empty() ) {
1088 // weighted data
1089 double x1 = fXMin - 0.001*(fXMax-fXMin);
1090 double x2 = fXMax + 0.001*(fXMax-fXMin);
1091 TH1D h1("temphist","", 500, x1, x2);
1092 h1.FillN(fEvents.size(), fEvents.data(), fEventWeights.data() );
1093 assert (h1.GetSumOfWeights() > 0) ;
1094 fMean = h1.GetMean();
1095 fSigma = h1.GetRMS();
1096 // compute robust sigma using midspread
1097 Double_t quantiles[2] = {0.0, 0.0};
1098 Double_t prob[2] = {0.25, 0.75};
1099 h1.GetQuantiles(2, quantiles, prob);
1100 Double_t midspread = quantiles[1] - quantiles[0];
1101 fSigmaRob = std::min(fSigma, midspread / 1.349); // Sigma's robust estimator
1102 //printf("weight case - stat: m = %f, s= %f, sr = %f \n",fMean, fSigma, midspread);
1103 return;
1104 }
1105 else {
1106 // compute statistics using the data
1107 SetMean();
1108 Double_t midspread = ComputeMidspread();
1109 SetSigma(midspread);
1110 //printf("un-weight case - stat: m = %f, s= %f, sr = %f \n",fMean, fSigma, fSigmaRob);
1111 }
1112}
1113
1115 // Computes the inter-quartile range from the data
1116 std::sort(fEvents.begin(), fEvents.end());
1117 Double_t quantiles[2] = {0.0, 0.0};
1118 Double_t prob[2] = {0.25, 0.75};
1119 TMath::Quantiles(fEvents.size(), 2, &fEvents[0], quantiles, prob);
1120 Double_t lowerquartile = quantiles[0];
1121 Double_t upperquartile = quantiles[1];
1122 return upperquartile - lowerquartile;
1123}
1124
1126 // Computes the user's input kernel function canonical bandwidth
1128}
1129
1131 // Computes the user's input kernel function sigma2
1133}
1134
1135TKDE::KernelIntegrand::KernelIntegrand(const TKDE* kde, EIntegralResult intRes) : fKDE(kde), fIntegralResult(intRes) {
1136 // Internal class constructor
1137}
1138
1140 // Internal class unary function
1141 if (fIntegralResult == kNorm) {
1142 return std::pow((*fKDE->fKernelFunction)(x), 2);
1143 } else if (fIntegralResult == kMu) {
1144 return x * (*fKDE->fKernelFunction)(x);
1145 } else if (fIntegralResult == kSigma2) {
1146 return std::pow(x, 2) * (*fKDE->fKernelFunction)(x);
1147 } else if (fIntegralResult == kUnitIntegration) {
1148 return (*fKDE->fKernelFunction)(x);
1149 } else {
1150 return -1;
1151 }
1152}
1153
1155 //Returns the estimated density
1156 TString name = "KDEFunc_"; name+= GetName();
1157 TString title = "KDE "; title+= GetTitle();
1158 if (xMin >= xMax) { xMin = fXMin; xMax = fXMax; }
1159 //Do not register the TF1 to global list
1160 bool previous = TF1::DefaultAddToGlobalList(kFALSE);
1161 TF1 * pdf = new TF1(name.Data(), this, xMin, xMax, 0);
1163 if (npx > 0) pdf->SetNpx(npx);
1164 pdf->SetTitle(title);
1165 return pdf;
1166}
1167
1169 // Returns the upper estimated density
1170 TString name;
1171 name.Form("KDE_UpperCL%f5.3_%s",confidenceLevel,GetName());
1172 if (xMin >= xMax) { xMin = fXMin; xMax = fXMax; }
1173 TF1 * upperPDF = new TF1(name, this, &TKDE::UpperConfidenceInterval, xMin, xMax, 1);
1174 upperPDF->SetParameter(0, confidenceLevel);
1175 if (npx > 0) upperPDF->SetNpx(npx);
1176 TF1 * f = (TF1*)upperPDF->Clone();
1177 delete upperPDF;
1178 return f;
1179}
1180
1182 // Returns the upper estimated density
1183 TString name;
1184 name.Form("KDE_LowerCL%f5.3_%s",confidenceLevel,GetName());
1185 if (xMin >= xMax) { xMin = fXMin; xMax = fXMax; }
1186 TF1 * lowerPDF = new TF1(name, this, &TKDE::LowerConfidenceInterval, xMin, xMax, 1);
1187 lowerPDF->SetParameter(0, confidenceLevel);
1188 if (npx > 0) lowerPDF->SetNpx(npx);
1189 TF1 * f = (TF1*)lowerPDF->Clone();
1190 delete lowerPDF;
1191 return f;
1192}
1193
1195 // Returns the approximate bias
1196 TString name = "KDE_Bias_"; name += GetName();
1197 if (xMin >= xMax) { xMin = fXMin; xMax = fXMax; }
1198 TF1 * approximateBias = new TF1(name, this, &TKDE::ApproximateBias, xMin, xMax, 0);
1199 if (npx > 0) approximateBias->SetNpx(npx);
1200 TF1 * f = (TF1*)approximateBias->Clone();
1201 delete approximateBias;
1202 return f;
1203}
#define f(i)
Definition: RSha256.hxx:104
static const double x2[5]
static const double x1[5]
#define M_PI
Definition: Rotated.cxx:105
int Int_t
Definition: RtypesCore.h:43
unsigned int UInt_t
Definition: RtypesCore.h:44
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
const char Option_t
Definition: RtypesCore.h:64
#define ClassImp(name)
Definition: Rtypes.h:361
@ kRed
Definition: Rtypes.h:64
@ kBlue
Definition: Rtypes.h:64
char name[80]
Definition: TGX11.cxx:109
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
double pow(double, double)
double sqrt(double)
double exp(double)
double log(double)
TRObject operator()(const T1 &t1) const
#define gPad
Definition: TVirtualPad.h:287
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
User Class for performing numerical integration of a function in one dimension.
Definition: Integrator.h:94
void SetFunction(Function &f)
method to set the a generic integration function
Definition: Integrator.h:489
double Integral(Function &f, double a, double b)
evaluate the Integral of a function f over the defined interval (a,b)
Definition: Integrator.h:496
User class for calculating the derivatives of a function.
double Derivative2(double x)
Returns the second derivative of the function at point x, computed by Richardson's extrapolation meth...
void SetFunction(const IGenFunction &f)
Set function for derivative calculation (copy the function if option has been enabled in the construc...
Template class to wrap any C++ callable object which takes one argument i.e.
Template class to wrap any member function of a class taking a double and returning a double in a 1D ...
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
1-Dim function class
Definition: TF1.h:210
virtual void SetTitle(const char *title="")
Set function title if title has the form "fffffff;xxxx;yyyy", it is assumed that the function title i...
Definition: TF1.cxx:3551
virtual void SetNpx(Int_t npx=100)
Set the number of points used to draw the function.
Definition: TF1.cxx:3435
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF1.cxx:1320
static Bool_t DefaultAddToGlobalList(Bool_t on=kTRUE)
Static method to add/avoid to add automatically functions to the global list (gROOT->GetListOfFunctio...
Definition: TF1.cxx:824
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TF1.cxx:1051
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:628
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:26
virtual void SetName(const char *name="")
Set graph name.
Definition: TGraph.cxx:2308
virtual void SetTitle(const char *title="")
Change (i.e.
Definition: TGraph.cxx:2324
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:760
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
virtual Int_t GetQuantiles(Int_t nprobSum, Double_t *q, const Double_t *probSum=0)
Compute Quantiles for this histogram Quantile x_q of a probability distribution Function F is defined...
Definition: TH1.cxx:4459
virtual Double_t GetMean(Int_t axis=1) const
For axis = 1,2 or 3 returns the mean value of the histogram along X,Y or Z axis.
Definition: TH1.cxx:7086
Double_t GetRMS(Int_t axis=1) const
Definition: TH1.h:311
virtual void FillN(Int_t ntimes, const Double_t *x, const Double_t *w, Int_t stride=1)
Fill this histogram with an array x and weights w.
Definition: TH1.cxx:3378
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition: TH1.cxx:7426
Kernel Density Estimation class.
Definition: TKDE.h:31
TF1 * GetPDFUpperConfidenceInterval(Double_t confidenceLevel=0.95, UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
Definition: TKDE.cxx:1168
TF1 * GetKDEApproximateBias(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
Definition: TKDE.cxx:1194
void SetData(const Double_t *data, const Double_t *weights)
Definition: TKDE.cxx:436
TF1 * fLowerPDF
Output Kernel Density Estimation upper confidence interval PDF function.
Definition: TKDE.h:151
virtual ~TKDE()
Definition: TKDE.cxx:88
std::vector< Double_t > fKernelSigmas2
Definition: TKDE.h:182
Double_t ComputeKernelL2Norm() const
Definition: TKDE.cxx:1049
TF1 * fPDF
Definition: TKDE.h:149
TF1 * GetPDFLowerConfidenceInterval(Double_t confidenceLevel=0.95, UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
Definition: TKDE.cxx:1181
void SetKernelType(EKernelType kern)
Definition: TKDE.cxx:314
std::vector< Double_t > fCanonicalBandwidths
Definition: TKDE.h:181
void SetKernelFunction(KernelFunction_Ptr kernfunc=0)
Definition: TKDE.cxx:579
UInt_t fNEvents
Definition: TKDE.h:167
void ComputeDataStats()
Definition: TKDE.cxx:1085
Double_t fXMax
Definition: TKDE.h:175
Double_t UpperConfidenceInterval(const Double_t *x, const Double_t *p) const
Definition: TKDE.cxx:988
friend class TKernel
Definition: TKDE.h:140
void ReInit()
Definition: TKDE.cxx:470
Double_t ApproximateBias(const Double_t *x, const Double_t *) const
Definition: TKDE.h:212
Double_t ComputeMidspread()
Definition: TKDE.cxx:1114
Bool_t fNewData
Definition: TKDE.h:163
void SetMirror()
Definition: TKDE.cxx:427
Bool_t fUseMirroring
Definition: TKDE.h:161
void DrawConfidenceInterval(TString &drawOpt, double cl=0.95)
Definition: TKDE.cxx:887
void SetMirroredEvents()
Definition: TKDE.cxx:520
EMirror fMirror
Definition: TKDE.h:157
void SetUserCanonicalBandwidth()
Definition: TKDE.cxx:1125
Bool_t fMirrorLeft
Definition: TKDE.h:161
EIteration fIteration
Definition: TKDE.h:156
void CheckKernelValidity()
Definition: TKDE.cxx:1025
const Double_t * GetAdaptiveWeights() const
Definition: TKDE.cxx:913
Double_t fAdaptiveBandwidthFactor
Definition: TKDE.h:177
Double_t LowerConfidenceInterval(const Double_t *x, const Double_t *p) const
Definition: TKDE.cxx:997
Double_t fSigmaRob
Definition: TKDE.h:173
std::vector< Double_t > fBinCount
Definition: TKDE.h:184
EBinning fBinning
Definition: TKDE.h:158
EIteration
Definition: TKDE.h:43
@ kAdaptive
Definition: TKDE.h:44
@ kFixed
Definition: TKDE.h:45
Double_t GetRAMISE() const
Definition: TKDE.cxx:714
void SetIteration(EIteration iter)
Definition: TKDE.cxx:325
Double_t ComputeKernelIntegral() const
Definition: TKDE.cxx:1076
Double_t CosineArchKernel(Double_t x) const
Definition: TKDE.h:206
Double_t fXMin
Definition: TKDE.h:174
Double_t operator()(Double_t x) const
Definition: TKDE.cxx:692
void SetUserKernelSigma2()
Definition: TKDE.cxx:1130
Double_t GetBias(Double_t x) const
Definition: TKDE.cxx:1007
std::vector< Double_t > fData
internal kernel class. Transient because it is recreated after reading from a file
Definition: TKDE.h:145
Double_t fSumOfCounts
Definition: TKDE.h:168
UInt_t fUseBinsNEvents
Definition: TKDE.h:169
void SetKernel()
Definition: TKDE.cxx:564
Double_t fSigma
Definition: TKDE.h:172
TGraphErrors * GetGraphWithErrors(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
Definition: TKDE.cxx:863
Double_t GetMean() const
Definition: TKDE.cxx:702
Double_t fRho
Definition: TKDE.h:176
Bool_t fAsymRight
Definition: TKDE.h:161
Double_t fWeightSize
Definition: TKDE.h:179
TKernel * fKernel
Definition: TKDE.h:143
void SetUseBinsNEvents(UInt_t nEvents)
Definition: TKDE.cxx:371
std::vector< Double_t > fEvents
Definition: TKDE.h:146
Double_t GetError(Double_t x) const
Definition: TKDE.cxx:1016
TF1 * GetKDEFunction(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
Definition: TKDE.cxx:1154
void SetBinning(EBinning)
Definition: TKDE.cxx:345
void GetOptions(std::string optionType, std::string option)
Definition: TKDE.cxx:208
Bool_t fAsymLeft
Definition: TKDE.h:161
std::vector< Bool_t > fSettedOptions
Definition: TKDE.h:186
Double_t GaussianKernel(Double_t x) const
Definition: TKDE.h:194
void SetRange(Double_t xMin, Double_t xMax)
Definition: TKDE.cxx:387
void SetMean()
Definition: TKDE.cxx:553
virtual void Draw(const Option_t *option="")
Definition: TKDE.cxx:805
Double_t ComputeKernelSigma2() const
Definition: TKDE.cxx:1058
void SetOptions(const Option_t *option, Double_t rho)
Definition: TKDE.cxx:138
void SetUseBins()
Definition: TKDE.cxx:401
Double_t GetFixedWeight() const
Definition: TKDE.cxx:902
void AssureOptions()
Definition: TKDE.cxx:275
TF1 * GetFunction(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
Definition: TKDE.cxx:640
void SetBinCountData()
Definition: TKDE.cxx:767
TF1 * GetUpperFunction(Double_t confidenceLevel=0.95, UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
Definition: TKDE.cxx:649
void InitFromNewData()
Definition: TKDE.cxx:494
void Instantiate(KernelFunction_Ptr kernfunc, UInt_t events, const Double_t *data, const Double_t *weight, Double_t xMin, Double_t xMax, const Option_t *option, Double_t rho)
Definition: TKDE.cxx:99
void SetDrawOptions(const Option_t *option, TString &plotOpt, TString &drawOpt)
Definition: TKDE.cxx:164
EMirror
Definition: TKDE.h:48
@ kMirrorLeft
Definition: TKDE.h:50
@ kMirrorAsymLeftRight
Definition: TKDE.h:54
@ kMirrorAsymRight
Definition: TKDE.h:55
@ kMirrorAsymBoth
Definition: TKDE.h:57
@ kNoMirror
Definition: TKDE.h:49
@ kMirrorRight
Definition: TKDE.h:51
@ kMirrorLeftAsymRight
Definition: TKDE.h:56
@ kMirrorBoth
Definition: TKDE.h:52
@ kMirrorAsymLeft
Definition: TKDE.h:53
TGraphErrors * fGraph
Output Kernel Density Estimation approximate bias.
Definition: TKDE.h:153
void SetCanonicalBandwidths()
Definition: TKDE.cxx:623
TKDE()
default constructor used by I/O
Definition: TKDE.cxx:72
void SetBinCentreData(Double_t xmin, Double_t xmax)
Definition: TKDE.cxx:758
void SetTuneFactor(Double_t rho)
Definition: TKDE.cxx:377
UInt_t Index(Double_t x) const
Definition: TKDE.cxx:973
void Fill(Double_t data)
Definition: TKDE.cxx:664
TF1 * fUpperPDF
Output Kernel Density Estimation PDF function.
Definition: TKDE.h:150
Bool_t fMirrorRight
Definition: TKDE.h:161
UInt_t fNBins
Definition: TKDE.h:166
Double_t ComputeKernelMu() const
Definition: TKDE.cxx:1067
EKernelType
Definition: TKDE.h:34
@ kGaussian
Definition: TKDE.h:35
@ kEpanechnikov
Definition: TKDE.h:36
@ kCosineArch
Definition: TKDE.h:38
@ kBiweight
Definition: TKDE.h:37
@ kTotalKernels
Definition: TKDE.h:40
@ kUserDefined
Definition: TKDE.h:39
Double_t fMean
Definition: TKDE.h:171
void DrawErrors(TString &drawOpt)
Definition: TKDE.cxx:856
void SetNBins(UInt_t nbins)
Definition: TKDE.cxx:352
void CheckOptions(Bool_t isUserDefinedKernel=kFALSE)
Definition: TKDE.cxx:291
Double_t EpanechnikovKernel(Double_t x) const
Definition: TKDE.h:199
Double_t BiweightKernel(Double_t x) const
Definition: TKDE.h:202
void SetKernelSigmas2()
Definition: TKDE.cxx:632
Bool_t fUseMinMaxFromData
Definition: TKDE.h:164
Double_t GetSigma() const
Definition: TKDE.cxx:708
EKernelType fKernelType
Graph with the errors.
Definition: TKDE.h:155
TF1 * fApproximateBias
Output Kernel Density Estimation lower confidence interval PDF function.
Definition: TKDE.h:152
void SetSigma(Double_t R)
Definition: TKDE.cxx:558
std::vector< Double_t > fEventWeights
Definition: TKDE.h:147
TF1 * GetApproximateBias(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
Definition: TKDE.cxx:659
TF1 * GetLowerFunction(Double_t confidenceLevel=0.95, UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
Definition: TKDE.cxx:654
Bool_t fUseBins
Definition: TKDE.h:162
KernelFunction_Ptr fKernelFunction
Definition: TKDE.h:138
friend struct KernelIntegrand
Definition: TKDE.h:188
EBinning
Definition: TKDE.h:60
@ kUnbinned
Definition: TKDE.h:61
@ kRelaxedBinning
Definition: TKDE.h:62
@ kForcedBinning
Definition: TKDE.h:63
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:877
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:891
Basic string class.
Definition: TString.h:131
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
const char * Data() const
Definition: TString.h:364
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:687
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
double normal_quantile(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the lower tail of the normal (Gaussian) distri...
const Double_t sigma
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Double_t ey[n]
Definition: legend1.C:17
for(Int_t i=0;i< n;i++)
Definition: legend1.C:18
Double_t ex[n]
Definition: legend1.C:17
TH1F * h1
Definition: legend1.C:5
double inner_product(const LAVector &, const LAVector &)
static constexpr double s
Bool_t IsNaN(Double_t x)
Definition: TMath.h:882
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754
Definition: TMath.h:891
void Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob, Bool_t isSorted=kTRUE, Int_t *index=0, Int_t type=7)
Computes sample quantiles, corresponding to the given probabilities.
Definition: TMath.cxx:1183