63 enum EIntegralResult{
kNorm, kMu, kSigma2, kUnitIntegration};
68 EIntegralResult fIntegralResult;
118 fData = std::vector<Double_t>(events, 0.0);
119 fEvents = std::vector<Double_t>(events, 0.0);
130 fNBins = events < 10000 ? 100 : events / 10;
157 std::string options = opt.
Data();
159 std::vector<std::string> voption(numOpt,
"");
160 for (std::vector<std::string>::iterator it = voption.begin(); it != voption.end() && !options.empty(); ++it) {
161 size_t pos = options.find_last_of(
';');
162 if (pos == std::string::npos) {
166 *it = options.substr(pos + 1);
167 options = options.substr(0, pos);
169 for (std::vector<std::string>::iterator it = voption.begin(); it != voption.end(); ++it) {
170 size_t pos = (*it).find(
':');
171 if (pos != std::string::npos) {
172 GetOptions((*it).substr(0, pos) , (*it).substr(pos + 1));
183 std::vector<std::string> voption(numOpt,
"");
184 for (std::vector<std::string>::iterator it = voption.begin(); it != voption.end() && !options.empty(); ++it) {
185 size_t pos = options.find_last_of(
';');
186 if (pos == std::string::npos) {
190 *it = options.substr(pos + 1);
191 options = options.substr(0, pos);
195 for (std::vector<std::string>::iterator it = voption.begin(); it != voption.end() && !options.empty(); ++it) {
196 size_t pos = (*it).find(
':');
197 if (pos == std::string::npos)
break;
198 TString optionType = (*it).substr(0, pos);
199 TString optionInstance = (*it).substr(pos + 1);
203 foundPlotOPt =
kTRUE;
204 if (optionInstance.
Contains(
"estimate") || optionInstance.
Contains(
"errors") || optionInstance.
Contains(
"confidenceinterval"))
205 plotOpt = optionInstance;
207 this->
Warning(
"SetDrawOptions",
"Unknown plotting option: setting to KDE estimate plot.");
208 }
else if (optionType.
Contains(
"drawoptions")) {
209 foundDrawOPt =
kTRUE;
210 drawOpt = optionInstance;
214 this->
Warning(
"SetDrawOptions",
"No plotting option: setting to KDE estimate plot.");
215 plotOpt =
"estimate";
218 this->
Warning(
"SetDrawOptions",
"No drawing options: setting to default ones.");
225 if (optionType.compare(
"kerneltype") == 0) {
227 if (option.compare(
"gaussian") == 0) {
229 }
else if (option.compare(
"epanechnikov") == 0) {
231 }
else if (option.compare(
"biweight") == 0) {
233 }
else if (option.compare(
"cosinearch") == 0) {
235 }
else if (option.compare(
"userdefined") == 0) {
238 this->
Warning(
"GetOptions",
"Unknown kernel type option: setting to Gaussian");
241 }
else if (optionType.compare(
"iteration") == 0) {
243 if (option.compare(
"adaptive") == 0) {
245 }
else if (option.compare(
"fixed") == 0) {
248 this->
Warning(
"GetOptions",
"Unknown iteration option: setting to Adaptive");
251 }
else if (optionType.compare(
"mirror") == 0) {
253 if (option.compare(
"nomirror") == 0) {
255 }
else if (option.compare(
"mirrorleft") == 0) {
257 }
else if (option.compare(
"mirrorright") == 0) {
259 }
else if (option.compare(
"mirrorboth") == 0) {
261 }
else if (option.compare(
"mirrorasymleft") == 0) {
263 }
else if (option.compare(
"mirrorasymleftright") == 0) {
265 }
else if (option.compare(
"mirrorasymright") == 0) {
267 }
else if (option.compare(
"mirrorleftasymright") == 0) {
269 }
else if (option.compare(
"mirrorasymboth") == 0) {
272 this->
Warning(
"GetOptions",
"Unknown mirror option: setting to NoMirror");
275 }
else if (optionType.compare(
"binning") == 0) {
277 if (option.compare(
"unbinned") == 0) {
279 }
else if (option.compare(
"relaxedbinning") == 0) {
281 }
else if (option.compare(
"forcedbinning") == 0) {
284 this->
Warning(
"GetOptions",
"Unknown binning option: setting to RelaxedBinning");
309 Error(
"CheckOptions",
"Illegal user kernel type input! Use template constructor for user defined kernel.");
312 Warning(
"CheckOptions",
"Illegal user iteration type input - use default value !");
316 Warning(
"CheckOptions",
"Illegal user mirroring type input - use default value !");
320 Warning(
"CheckOptions",
"Illegal user binning type input - use default value !");
324 Warning(
"CheckOptions",
"Tuning factor rho cannot be non-positive - use default value !");
371 Error(
"SetNBins",
"Number of bins must be greater than zero.");
380 Warning(
"SetNBins",
"Bin type using SetBinning must set for using a binned evaluation");
407 Error(
"SetRange",
"Minimum range cannot be bigger or equal than the maximum range! Present range values remain the same.");
462 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");
502 std::vector<Double_t> originalEvents =
fEvents;
552 fKernel->ComputeAdaptiveWeights();
590 Error(
"SetKernelFunction",
"User kernel function is not defined !");
644 this->
Warning(
"Fill",
"Cannot fill data with data binned option. Data input ignored.");
647 fData.push_back(data);
655 this->
Warning(
"Fill",
"Cannot fill data with data binned option. Data input ignored.");
658 fData.push_back(data);
671 if (
fNewData) (
const_cast<TKDE*
>(
this))->InitFromNewData();
677 if (
fNewData) (
const_cast<TKDE*
>(
this))->InitFromNewData();
683 if (
fNewData) (
const_cast<TKDE*
>(
this))->InitFromNewData();
702 std::vector<Double_t> weights =
fWeights;
703 Double_t minWeight = weights[0] * 0.05;
704 unsigned int n = fKDE->fData.size();
705 assert( n == weights.size() );
706 bool useDataWeights = (fKDE->fBinCount.size() ==
n);
708 for (
unsigned int i = 0; i <
n; ++i) {
710 if (useDataWeights && fKDE->fBinCount[i] <= 0)
continue;
711 f = (*fKDE->fKernel)(fKDE->fData[i]);
713 fKDE->Warning(
"ComputeAdativeWeights",
"function value is zero or negative for x = %f w = %f",
714 fKDE->fData[i],(useDataWeights) ? fKDE->fBinCount[i] : 1.);
716 fKDE->fAdaptiveBandwidthFactor +=
std::log(f);
719 Double_t kAPPROX_GEO_MEAN = 0.241970724519143365;
720 fKDE->fAdaptiveBandwidthFactor = fKDE->fUseMirroring ? kAPPROX_GEO_MEAN / fKDE->fSigmaRob :
std::sqrt(
std::exp(fKDE->fAdaptiveBandwidthFactor / fKDE->fData.size()));
721 transform(weights.begin(), weights.end(),
fWeights.begin(), std::bind2nd(std::multiplies<Double_t>(), fKDE->fAdaptiveBandwidthFactor));
735 fData[i] = xmin + (i + 0.5) * binWidth;
806 else if (plotOpt.
Contains(
"confidenceinterval") ||
812 const char * s = strstr(plotOpt.
Data(),
"interval@");
814 if (s != 0) sscanf(s,
"interval@%lf",&level);
815 if((level <= 0) || (level >= 1)) {
816 Warning(
"Draw",
"given confidence level %.3lf is invalid - use default 0.95",level);
836 if (xmin>= xmax) { xmin =
fXMin; xmax =
fXMax; }
843 for (
UInt_t i = 0; i <=
n; ++i) {
844 x[i] = xmin + i * (xmax -
xmin) / n;
845 y[i] = (*this)(x[i]);
850 ge->
SetName(
"kde_graph_error");
864 upper->
Draw((
"same" + drawOpt).
Data());
867 lower->
Draw((
"same" + drawOpt).
Data());
878 this->
Warning(
"GetFixedWeight()",
"Fixed iteration option not enabled. Returning %f.", result);
880 result =
fKernel->GetFixedWeight();
888 this->
Warning(
"GetFixedWeight()",
"Adaptive iteration option not enabled. Returning a NULL pointer<");
891 if (
fNewData) (
const_cast<TKDE*
>(
this))->InitFromNewData();
892 return &(
fKernel->GetAdaptiveWeights()).front();
908 UInt_t n = fKDE->fData.size();
910 Bool_t useBins = (fKDE->fBinCount.size() ==
n);
911 Double_t nSum = (useBins) ? fKDE->fSumOfCounts : fKDE->fNEvents;
915 Double_t binCount = (useBins) ? fKDE->fBinCount[i] : 1.0;
918 if (fKDE->fAsymLeft) {
919 result -= binCount /
fWeights[i] * (*fKDE->fKernelFunction)((x - (2. * fKDE->fXMin - fKDE->fData[i])) /
fWeights[i]);
921 if (fKDE->fAsymRight) {
922 result -= binCount /
fWeights[i] * (*fKDE->fKernelFunction)((x - (2. * fKDE->fXMax - fKDE->fData[i])) /
fWeights[i]);
938 fKDE->Warning(
"operator()",
"Result is NaN for x %f \n",x);
953 }
else if (bin <= 0) {
965 return f + z * sigma;
974 return f + z * sigma;
1000 valid = valid && unity == 1.;
1002 Error(
"CheckKernelValidity",
"Kernel's integral is %f",unity);
1005 valid = valid && mu == 0.;
1007 Error(
"CheckKernelValidity",
"Kernel's mu is %f" ,mu);
1012 Error(
"CheckKernelValidity",
"Kernel's sigma2 is %f",sigma2);
1015 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.");
1050 KernelIntegrand kernel(
this, TKDE::KernelIntegrand::kUnitIntegration);
1062 TH1D h1(
"temphist",
"", 500, x1, x2);
1068 Double_t quantiles[2] = {0.0, 0.0};
1071 Double_t midspread = quantiles[1] - quantiles[0];
1088 Double_t quantiles[2] = {0.0, 0.0};
1091 Double_t lowerquartile = quantiles[0];
1092 Double_t upperquartile = quantiles[1];
1093 return upperquartile - lowerquartile;
1106 TKDE::KernelIntegrand::KernelIntegrand(
const TKDE* kde, EIntegralResult intRes) : fKDE(kde), fIntegralResult(intRes) {
1112 if (fIntegralResult ==
kNorm) {
1113 return std::pow((*fKDE->fKernelFunction)(x), 2);
1114 }
else if (fIntegralResult == kMu) {
1115 return x * (*fKDE->fKernelFunction)(x);
1116 }
else if (fIntegralResult == kSigma2) {
1117 return std::pow(x, 2) * (*fKDE->fKernelFunction)(x);
1118 }
else if (fIntegralResult == kUnitIntegration) {
1119 return (*fKDE->fKernelFunction)(
x);
1129 if (xMin >= xMax) { xMin =
fXMin; xMax =
fXMax; }
1130 TF1 * pdf =
new TF1(name.
Data(),
this, xMin, xMax, 0);
1131 if (npx > 0) pdf->
SetNpx(npx);
1141 name.
Form(
"KDE_UpperCL%f5.3_%s",confidenceLevel,
GetName());
1142 if (xMin >= xMax) { xMin =
fXMin; xMax =
fXMax; }
1145 if (npx > 0) upperPDF->
SetNpx(npx);
1154 name.
Form(
"KDE_LowerCL%f5.3_%s",confidenceLevel,
GetName());
1155 if (xMin >= xMax) { xMin =
fXMin; xMax =
fXMax; }
1158 if (npx > 0) lowerPDF->
SetNpx(npx);
1167 if (xMin >= xMax) { xMin =
fXMin; xMax =
fXMax; }
1169 if (npx > 0) approximateBias->
SetNpx(npx);
1171 delete approximateBias;
double normal_quantile(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the lower tail of the normal (Gaussian) distri...
Double_t operator()(Double_t x) const
virtual const char * GetTitle() const
Returns title of object.
Double_t GaussianKernel(Double_t x) const
Bool_t fUseMinMaxFromData
static Vc_ALWAYS_INLINE int_v min(const int_v &x, const int_v &y)
void SetUserKernelSigma2()
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
virtual void SetNpx(Int_t npx=100)
Set the number of points used to draw the function.
TF1 * GetFunction(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
friend struct KernelIntegrand
Kernel Density Estimation class.
void SetRange(Double_t xMin, Double_t xMax)
TF1 * GetLowerFunction(Double_t confidenceLevel=0.95, UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
Double_t operator()(Double_t x) const
void GetOptions(std::string optionType, std::string option)
Double_t ComputeKernelSigma2() const
TString & ReplaceAll(const TString &s1, const TString &s2)
const std::vector< Double_t > & GetAdaptiveWeights() const
virtual void SetName(const char *name)
Change (i.e.
Double_t EpanechnikovKernel(Double_t x) const
Double_t ComputeKernelMu() const
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...
void SetIteration(EIteration iter)
double inner_product(const LAVector &, const LAVector &)
TKernel(Double_t weight, TKDE *kde)
void ToLower()
Change string to lower-case.
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
std::vector< Double_t > fKernelSigmas2
virtual void SetTitle(const char *title="")
Set graph title.
double Derivative2(double x)
Returns the second derivative of the function at point x, computed by Richardson's extrapolation meth...
Double_t GetFixedWeight() const
virtual void Draw(const Option_t *option="")
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
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
Template class to wrap any C++ callable object which takes one argument i.e.
const char * Data() const
TF1 * GetApproximateBias(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
void SetCanonicalBandwidths()
static const double x2[5]
Double_t GetFixedWeight() const
TGraphErrors * GetGraphWithErrors(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
std::vector< Double_t > fBinCount
void SetFunction(const IGenFunction &f)
Set function for derivative calculation (copy the function if option has been enabled in the construc...
TF1 * GetUpperFunction(Double_t confidenceLevel=0.95, UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
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 Parameters: x -the data sample n ...
std::map< std::string, std::string >::const_iterator iter
double pow(double, double)
void SetKernelFunction(KernelFunction_Ptr kernfunc=0)
Double_t ComputeKernelIntegral() const
std::vector< std::vector< double > > Data
Double_t BiweightKernel(Double_t x) const
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
void ComputeAdaptiveWeights()
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
std::vector< Double_t > fData
void SetUseBinsNEvents(UInt_t nEvents)
TF1 * GetKDEFunction(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
virtual void SetLineColor(Color_t lcolor)
double Integral(Function &f, double a, double b)
evaluate the Integral of a function f over the defined interval (a,b)
ClassImp(TKDE) class TKDE UInt_t fNWeights
void CheckKernelValidity()
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.
void SetOptions(const Option_t *option, Double_t rho)
void Form(const char *fmt,...)
Formats a string using a printf style format descriptor.
User Class for performing numerical integration of a function in one dimension.
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
std::vector< Double_t > fEventWeights
virtual const char * GetName() const
Returns name of object.
Double_t GetWeight(Double_t x) const
Double_t GetRMS(Int_t axis=1) const
Double_t CosineArchKernel(Double_t x) const
void SetSigma(Double_t R)
1-D histogram with a double per channel (see TH1 documentation)}
void SetNBins(UInt_t nbins)
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...
static const double x1[5]
Double_t ComputeKernelL2Norm() const
std::vector< Bool_t > fSettedOptions
void SetBinning(EBinning)
std::vector< Double_t > fEvents
TF1 * GetPDFUpperConfidenceInterval(Double_t confidenceLevel=0.95, UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
const Double_t * GetAdaptiveWeights() const
KernelFunction_Ptr fKernelFunction
void DrawErrors(TString &drawOpt)
Double_t GetBias(Double_t x) const
static Vc_ALWAYS_INLINE int_v max(const int_v &x, const int_v &y)
Double_t UpperConfidenceInterval(const Double_t *x, const Double_t *p) const
void CheckOptions(Bool_t isUserDefinedKernel=kFALSE)
void SetFunction(Function &f)
method to set the a generic integration function
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)
Double_t fAdaptiveBandwidthFactor
void SetUserCanonicalBandwidth()
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.
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
TF1 * GetKDEApproximateBias(UInt_t npx=100, Double_t xMin=1.0, Double_t xMax=0.0)
Double_t GetSigma() const
UInt_t Index(Double_t x) const
A TGraphErrors is a TGraph with error bars.
void SetTuneFactor(Double_t rho)
std::vector< Double_t > fWeights
Double_t ApproximateBias(const Double_t *x, const Double_t *) const
void SetData(const Double_t *data, const Double_t *weights)
Double_t LowerConfidenceInterval(const Double_t *x, const Double_t *p) const
virtual void SetParameter(Int_t param, Double_t value)
void DrawConfidenceInterval(TString &drawOpt, double cl=0.95)
void SetBinCentreData(Double_t xmin, Double_t xmax)
Double_t GetRAMISE() const
Template class to wrap any member function of a class taking a double and returning a double in a 1D ...
User class for calculating the derivatives of a function.
Double_t GetError(Double_t x) const
Double_t ComputeMidspread()
void SetDrawOptions(const Option_t *option, TString &plotOpt, TString &drawOpt)
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.