52 std::vector<Double_t> fWeights;
55 void ComputeAdaptiveWeights();
63 enum EIntegralResult{
kNorm, kMu, kSigma2, kUnitIntegration};
68 EIntegralResult fIntegralResult;
73 fKernelFunction(nullptr),
78 fApproximateBias(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)
102 fData = std::vector<Double_t>(events, 0.0);
103 fEvents = std::vector<Double_t>(events, 0.0);
114 fNBins = events < 10000 ? 100 : events / 10;
142 std::string options = opt.
Data();
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) {
151 *it = options.substr(pos + 1);
152 options = options.substr(0, pos);
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));
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) {
175 *it = options.substr(pos + 1);
176 options = options.substr(0, pos);
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);
188 foundPlotOPt =
kTRUE;
189 if (optionInstance.
Contains(
"estimate") || optionInstance.
Contains(
"errors") || optionInstance.
Contains(
"confidenceinterval"))
190 plotOpt = optionInstance;
192 this->
Warning(
"SetDrawOptions",
"Unknown plotting option: setting to KDE estimate plot.");
193 }
else if (optionType.
Contains(
"drawoptions")) {
194 foundDrawOPt =
kTRUE;
195 drawOpt = optionInstance;
199 this->
Warning(
"SetDrawOptions",
"No plotting option: setting to KDE estimate plot.");
200 plotOpt =
"estimate";
203 this->
Warning(
"SetDrawOptions",
"No drawing options: setting to default ones.");
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) {
223 this->
Warning(
"GetOptions",
"Unknown kernel type option: setting to Gaussian");
226 }
else if (optionType.compare(
"iteration") == 0) {
228 if (option.compare(
"adaptive") == 0) {
230 }
else if (option.compare(
"fixed") == 0) {
233 this->
Warning(
"GetOptions",
"Unknown iteration option: setting to Adaptive");
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) {
257 this->
Warning(
"GetOptions",
"Unknown mirror option: setting to NoMirror");
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) {
269 this->
Warning(
"GetOptions",
"Unknown binning option: setting to RelaxedBinning");
294 Error(
"CheckOptions",
"Illegal user kernel type input! Use template constructor for user defined kernel.");
297 Warning(
"CheckOptions",
"Illegal user iteration type input - use default value !");
301 Warning(
"CheckOptions",
"Illegal user mirroring type input - use default value !");
305 Warning(
"CheckOptions",
"Illegal user binning type input - use default value !");
309 Warning(
"CheckOptions",
"Tuning factor rho cannot be non-positive - use default value !");
355 Error(
"SetNBins",
"Number of bins must be greater than zero.");
365 Warning(
"SetNBins",
"Bin type using SetBinning must be set for using a binned evaluation");
367 Warning(
"SetNBins",
"Bin type using SetBinning or with SetUseBinsNEvents must be set for using a binned evaluation");
390 Error(
"SetRange",
"Minimum range cannot be bigger or equal than the maximum range! Present range values remain the same.");
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");
482 if (
fKernelFunction)
Error(
"ReInit",
"Kernel function pointer should be a nullptr when re-initializing after reading from a file");
485 Error(
"ReInit",
"TKDE does not contain any data !");
498 Error(
"InitFromNewData",
"Re-felling is not supported with binning");
522 std::vector<Double_t> originalEvents =
fEvents;
527 std::bind(std::minus<Double_t>(), 2 *
fXMin, std::placeholders::_1));
532 std::bind(std::minus<Double_t>(), 2 *
fXMax, std::placeholders::_1));
574 fKernel->ComputeAdaptiveWeights();
613 Error(
"SetKernelFunction",
"User kernel function is not defined !");
667 this->
Warning(
"Fill",
"Cannot fill data with data binned option. Data input ignored.");
670 fData.push_back(data);
678 this->
Warning(
"Fill",
"Cannot fill data with data binned option. Data input ignored.");
681 fData.push_back(data);
695 (
const_cast<TKDE*
>(
this))->ReInit();
704 if (
fNewData) (
const_cast<TKDE*
>(
this))->InitFromNewData();
710 if (
fNewData) (
const_cast<TKDE*
>(
this))->InitFromNewData();
723fNWeights(kde->fData.size()),
724fWeights(fNWeights, weight)
727void TKDE::TKernel::ComputeAdaptiveWeights() {
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);
735 for (
unsigned int i = 0; i <
n; ++i) {
737 if (useDataWeights && fKDE->fBinCount[i] <= 0)
continue;
738 f = (*fKDE->fKernel)(fKDE->fData[i]);
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);
746 Double_t kAPPROX_GEO_MEAN = 0.241970724519143365;
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));
755 return fWeights[fKDE->Index(
x)];
834 else if (plotOpt.
Contains(
"confidenceinterval") ||
840 const char *
s = strstr(plotOpt.
Data(),
"interval@");
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);
871 for (
UInt_t i = 0; i <=
n; ++i) {
873 y[i] = (*this)(
x[i]);
878 ge->
SetName(
"kde_graph_error");
892 upper->
Draw((
"same" + drawOpt).Data());
895 lower->
Draw((
"same" + drawOpt).Data());
906 this->
Warning(
"GetFixedWeight()",
"Fixed iteration option not enabled. Returning %f.", result);
908 result =
fKernel->GetFixedWeight();
916 this->
Warning(
"GetFixedWeight()",
"Adaptive iteration option not enabled. Returning a NULL pointer<");
919 if (
fNewData) (
const_cast<TKDE*
>(
this))->InitFromNewData();
920 return &(
fKernel->GetAdaptiveWeights()).front();
923Double_t TKDE::TKernel::GetFixedWeight()
const {
928const std::vector<Double_t> & TKDE::TKernel::GetAdaptiveWeights()
const {
938 Bool_t useBins = (fKDE->fBinCount.size() ==
n);
939 Double_t nSum = (useBins) ? fKDE->fSumOfCounts : fKDE->fNEvents;
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]);
948 if (fKDE->fAsymRight) {
949 result -= binCount / fWeights[i] * (*fKDE->fKernelFunction)((
x - (2. * fKDE->fXMax - fKDE->fData[i])) / fWeights[i]);
967 fKDE->Warning(
"operator()",
"Result is NaN for x %f \n",
x);
970 return result / nSum;
982 }
else if (bin <= 0) {
1029 valid = valid && unity == 1.;
1031 Error(
"CheckKernelValidity",
"Kernel's integral is %f",unity);
1034 valid = valid && mu == 0.;
1036 Error(
"CheckKernelValidity",
"Kernel's mu is %f" ,mu);
1039 valid = valid && sigma2 > 0 && sigma2 != std::numeric_limits<Double_t>::infinity();
1041 Error(
"CheckKernelValidity",
"Kernel's sigma2 is %f",sigma2);
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.");
1079 KernelIntegrand kernel(
this, TKDE::KernelIntegrand::kUnitIntegration);
1097 Double_t quantiles[2] = {0.0, 0.0};
1100 Double_t midspread = quantiles[1] - quantiles[0];
1117 Double_t quantiles[2] = {0.0, 0.0};
1120 Double_t lowerquartile = quantiles[0];
1121 Double_t upperquartile = quantiles[1];
1122 return upperquartile - lowerquartile;
1135TKDE::KernelIntegrand::KernelIntegrand(
const TKDE* kde, EIntegralResult intRes) : fKDE(kde), fIntegralResult(intRes) {
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);
1158 if (xMin >= xMax) { xMin =
fXMin; xMax =
fXMax; }
1161 TF1 * pdf =
new TF1(
name.Data(),
this, xMin, xMax, 0);
1163 if (npx > 0) pdf->
SetNpx(npx);
1171 name.Form(
"KDE_UpperCL%f5.3_%s",confidenceLevel,
GetName());
1172 if (xMin >= xMax) { xMin =
fXMin; xMax =
fXMax; }
1175 if (npx > 0) upperPDF->
SetNpx(npx);
1184 name.Form(
"KDE_LowerCL%f5.3_%s",confidenceLevel,
GetName());
1185 if (xMin >= xMax) { xMin =
fXMin; xMax =
fXMax; }
1188 if (npx > 0) lowerPDF->
SetNpx(npx);
1197 if (xMin >= xMax) { xMin =
fXMin; xMax =
fXMax; }
1199 if (npx > 0) approximateBias->
SetNpx(npx);
1201 delete approximateBias;
static const double x2[5]
static const double x1[5]
double pow(double, double)
TRObject operator()(const T1 &t1) const
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
User Class for performing numerical integration of a function in one dimension.
void SetFunction(Function &f)
method to set the a generic integration function
double Integral(Function &f, double a, double b)
evaluate the Integral of a function f over the defined interval (a,b)
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.
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...
virtual void SetNpx(Int_t npx=100)
Set the number of points used to draw the function.
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
static Bool_t DefaultAddToGlobalList(Bool_t on=kTRUE)
Static method to add/avoid to add automatically functions to the global list (gROOT->GetListOfFunctio...
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
virtual void SetParameter(Int_t param, Double_t value)
A TGraphErrors is a TGraph with error bars.
virtual void SetName(const char *name="")
Set graph name.
virtual void SetTitle(const char *title="")
Change (i.e.
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
1-D histogram with a double per channel (see TH1 documentation)}
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...
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.
Double_t GetRMS(Int_t axis=1) const
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.
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Kernel Density Estimation class.
TF1 * GetPDFUpperConfidenceInterval(Double_t confidenceLevel=0.95, UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
TF1 * GetKDEApproximateBias(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
void SetData(const Double_t *data, const Double_t *weights)
TF1 * fLowerPDF
Output Kernel Density Estimation upper confidence interval PDF function.
std::vector< Double_t > fKernelSigmas2
Double_t ComputeKernelL2Norm() const
TF1 * GetPDFLowerConfidenceInterval(Double_t confidenceLevel=0.95, UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
void SetKernelType(EKernelType kern)
std::vector< Double_t > fCanonicalBandwidths
void SetKernelFunction(KernelFunction_Ptr kernfunc=0)
Double_t UpperConfidenceInterval(const Double_t *x, const Double_t *p) const
Double_t ApproximateBias(const Double_t *x, const Double_t *) const
Double_t ComputeMidspread()
void DrawConfidenceInterval(TString &drawOpt, double cl=0.95)
void SetUserCanonicalBandwidth()
void CheckKernelValidity()
const Double_t * GetAdaptiveWeights() const
Double_t fAdaptiveBandwidthFactor
Double_t LowerConfidenceInterval(const Double_t *x, const Double_t *p) const
std::vector< Double_t > fBinCount
Double_t GetRAMISE() const
void SetIteration(EIteration iter)
Double_t ComputeKernelIntegral() const
Double_t CosineArchKernel(Double_t x) const
Double_t operator()(Double_t x) const
void SetUserKernelSigma2()
Double_t GetBias(Double_t x) const
std::vector< Double_t > fData
internal kernel class. Transient because it is recreated after reading from a file
TGraphErrors * GetGraphWithErrors(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
void SetUseBinsNEvents(UInt_t nEvents)
std::vector< Double_t > fEvents
Double_t GetError(Double_t x) const
TF1 * GetKDEFunction(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
void SetBinning(EBinning)
void GetOptions(std::string optionType, std::string option)
std::vector< Bool_t > fSettedOptions
Double_t GaussianKernel(Double_t x) const
void SetRange(Double_t xMin, Double_t xMax)
virtual void Draw(const Option_t *option="")
Double_t ComputeKernelSigma2() const
void SetOptions(const Option_t *option, Double_t rho)
Double_t GetFixedWeight() const
TF1 * GetFunction(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
TF1 * GetUpperFunction(Double_t confidenceLevel=0.95, UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
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)
void SetDrawOptions(const Option_t *option, TString &plotOpt, TString &drawOpt)
TGraphErrors * fGraph
Output Kernel Density Estimation approximate bias.
void SetCanonicalBandwidths()
TKDE()
default constructor used by I/O
void SetBinCentreData(Double_t xmin, Double_t xmax)
void SetTuneFactor(Double_t rho)
UInt_t Index(Double_t x) const
TF1 * fUpperPDF
Output Kernel Density Estimation PDF function.
Double_t ComputeKernelMu() const
void DrawErrors(TString &drawOpt)
void SetNBins(UInt_t nbins)
void CheckOptions(Bool_t isUserDefinedKernel=kFALSE)
Double_t EpanechnikovKernel(Double_t x) const
Double_t BiweightKernel(Double_t x) const
Bool_t fUseMinMaxFromData
Double_t GetSigma() const
EKernelType fKernelType
Graph with the errors.
TF1 * fApproximateBias
Output Kernel Density Estimation lower confidence interval PDF function.
void SetSigma(Double_t R)
std::vector< Double_t > fEventWeights
TF1 * GetApproximateBias(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
TF1 * GetLowerFunction(Double_t confidenceLevel=0.95, UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
KernelFunction_Ptr fKernelFunction
friend struct KernelIntegrand
virtual const char * GetTitle() const
Returns title of object.
virtual const char * GetName() const
Returns name of object.
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
void ToLower()
Change string to lower-case.
const char * Data() const
TString & ReplaceAll(const TString &s1, const TString &s2)
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
double normal_quantile(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the lower tail of the normal (Gaussian) distri...
double inner_product(const LAVector &, const LAVector &)
static constexpr double s
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754
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.