ROOT   6.10/09 Reference Guide
testStatFunc.cxx
Go to the documentation of this file.
1 /// test of the statistical functions cdf and quantiles
2
3 //#define NO_MATHCORE
4
6 #ifndef NO_MATHCORE
7 #include "Math/GSLIntegrator.h"
8 #include "Math/WrappedFunction.h"
10 #endif
11
12
13 #include <iostream>
14 #include <limits>
15
16 using namespace ROOT::Math;
17
18 int compare( std::string name, double v1, double v2, double scale = 2.0) {
19  // ntest = ntest + 1;
20
21  //std::cout << std::setw(50) << std::left << name << ":\t";
22
23  // numerical double limit for epsilon
24  double eps = scale* std::numeric_limits<double>::epsilon();
25  int iret = 0;
26  double delta = v2 - v1;
27  double d = 0;
28  if (delta < 0 ) delta = - delta;
29  if (v1 == 0 || v2 == 0) {
30  if (delta > eps ) {
31  iret = 1;
32  }
33  }
34  // skip case v1 or v2 is infinity
35  else {
36  d = v1;
37
38  if ( v1 < 0) d = -d;
39  // add also case when delta is small by default
40  if ( delta/d > eps && delta > eps )
41  iret = 1;
42  }
43
44  if (iret == 0)
45  std::cout <<".";
46  else {
47  int pr = std::cout.precision (18);
48  std::cout << "\nDiscrepancy in " << name.c_str() << "() :\n " << v1 << " != " << v2 << " discr = " << int(delta/d/eps)
49  << " (Allowed discrepancy is " << eps << ")\n\n";
50  std::cout.precision (pr);
51  //nfail = nfail + 1;
52  }
53  return iret;
54 }
55
56
57
58 // for functions with one parameters
59
60 struct TestFunc1 {
61
62  typedef double ( * Pdf) ( double, double, double);
63  typedef double ( * Cdf) ( double, double, double);
64  typedef double ( * Quant) ( double, double);
65
66  // wrappers to pdf to be integrated
67  struct PdfFunction {
68  PdfFunction( Pdf pdf, double p1 ) : fPdf(pdf), fP1(p1) {}
69  double operator() ( double x) const { return fPdf(x,fP1,0.0); }
70
71  Pdf fPdf;
72  double fP1;
73  };
74
75
76  TestFunc1( double p1, double s = 2) :
77  scale(s)
78  {
79  p[0] = p1;
80  }
81
82 #ifndef NO_MATHCORE // skip cdf test when Mathcore is missing
83
84  int testCdf( Pdf f_pdf, Cdf f_cdf, bool c = false) {
85  int iret = 0;
86  double val = 1.0; // value used for testing
87  double q1 = f_cdf(val, p[0],0.0);
88  // calculate integral of pdf
89  PdfFunction f(f_pdf,p[0]);
91  GSLIntegrator ig(1.E-12,1.E-12,100000);
92  ig.SetFunction(wf);
93  if (!c) {
94  // lower intergal (cdf)
95  double q2 = ig.IntegralLow(val);
96  // use a larger scale (integral error is 10-9)
97  iret |= compare("test _cdf", q1, q2, 1.0E6);
98  }
99  else {
100  // upper integral (cdf_c)
101  double q2 = ig.IntegralUp(val);
102  iret |= compare("test _cdf_c", q1, q2, 1.0E6);
103  }
104  return iret;
105  }
106
107 #endif
108
109  int testQuantile (Cdf f_cdf, Quant f_quantile, bool c= false) {
110  int iret = 0;
111  double z1,z2,q;
112  z1 = 1.0E-6;
113  for (int i = 0; i < 10; ++i) {
114  q = f_quantile(z1,p[0]);
115  z2 = f_cdf(q, p[0],0.);
116  if (!c)
117  iret |= compare("test quantile", z1, z2, scale);
118  else
119  iret |= compare("test quantile_c", z1, z2, scale);
120  z1 += 0.1;
121  }
122  return iret;
123  }
124
125  double p[1]; // parameters
126  double scale;
127 };
128
129 // for functions with two parameters
130
131
132 struct TestFunc2 {
133
134  typedef double ( * Pdf) ( double, double, double, double);
135  typedef double ( * Cdf) ( double, double, double, double);
136  typedef double ( * Quant) ( double, double, double);
137
138
139  struct PdfFunction {
140  PdfFunction( Pdf pdf, double p1, double p2 ) : fPdf(pdf), fP1(p1), fP2(p2) {}
141  double operator() ( double x) const { return fPdf(x,fP1,fP2,0.0); }
142
143  Pdf fPdf;
144  double fP1;
145  double fP2;
146  };
147
148  TestFunc2( double p1, double p2, double s = 2) :
149  scale(s)
150  {
151  p[0] = p1,
152  p[1] = p2;
153  }
154
155 #ifndef NO_MATHCORE // skip cdf test when Mathcore is missing
156
157  int testCdf( Pdf f_pdf, Cdf f_cdf, bool c = false) {
158  int iret = 0;
159  double val = 1.0; // value used for testing
160  double q1 = f_cdf(val, p[0],p[1],0.0);
161  // calculate integral of pdf
162  PdfFunction f(f_pdf,p[0],p[1]);
164  GSLIntegrator ig(1.E-12,1.E-12,100000);
165  ig.SetFunction(wf);
166  if (!c) {
167  // lower intergal (cdf)
168  double q2 = ig.IntegralLow(val);
169  // use a larger scale (integral error is 10-9)
170  iret |= compare("test _cdf", q1, q2, 1.0E6);
171  }
172  else {
173  // upper integral (cdf_c)
174  double q2 = ig.IntegralUp(val);
175  iret |= compare("test _cdf_c", q1, q2, 1.0E6);
176  }
177  return iret;
178  }
179
180 #endif
181
182  int testQuantile (Cdf f_cdf, Quant f_quantile, bool c=false) {
183  int iret = 0;
184  double z1,z2,q;
185  z1 = 1.0E-6;
186  for (int i = 0; i < 10; ++i) {
187  q = f_quantile(z1,p[0],p[1]);
188  z2 = f_cdf(q, p[0],p[1],0.0);
189  if (!c)
190  iret |= compare("test quantile", z1, z2, scale);
191  else
192  iret |= compare("test quantile_c", z1, z2, scale);
193  z1 += 0.1;
194  }
195  return iret;
196  }
197
198  double p[2]; // parameters
199  double scale;
200 };
201
202 void printStatus(int iret) {
203  if (iret == 0)
204  std::cout <<"\t\t\t\t OK" << std::endl;
205  else
206  std::cout <<"\t\t\t\t FAILED " << std::endl;
207 }
208
209 #ifndef NO_MATHCORE
210
211 #define TESTDIST1(name, p1, s) {\
212  int ir = 0; \
213  TestFunc1 t(p1,s);\
214  ir |= t.testCdf( name ## _pdf , name ## _cdf );\
215  ir |= t.testCdf( name ##_pdf , name ##_cdf_c,true);\
216  ir |= t.testQuantile( name ## _cdf, name ##_quantile);\
217  ir |= t.testQuantile( name ##_cdf_c, name ##_quantile_c,true);\
218  printStatus(ir);\
219  iret |= ir; }
220
221 #define TESTDIST2(name, p1, p2, s) {\
222  int ir = 0; \
223  TestFunc2 t(p1,p2,s);\
224  ir |= t.testCdf( name ## _pdf , name ## _cdf );\
225  ir |= t.testCdf( name ##_pdf , name ##_cdf_c,true);\
226  ir |= t.testQuantile( name ## _cdf, name ##_quantile);\
227  ir |= t.testQuantile( name ##_cdf_c, name ##_quantile_c,true);\
228  printStatus(ir);\
229  iret |= ir; }
230
231
232
233
234 #else
235 // without mathcore pdf are missing so skip cdf test
236 #define TESTDIST1(name, p1, s) {\
237  int ir = 0; \
238  TestFunc1 t(p1,s);\
239  ir |= t.testQuantile( name ## _cdf, name ##_quantile);\
240  ir |= t.testQuantile( name ##_cdf_c, name ##_quantile_c,true);\
241  printStatus(ir);\
242  iret |= ir; }
243
244 #define TESTDIST2(name, p1, p2, s) {\
245  int ir = 0; \
246  TestFunc2 t(p1,p2,s);\
247  ir |= t.testQuantile( name ## _cdf, name ##_quantile);\
248  ir |= t.testQuantile( name ##_cdf_c, name ##_quantile_c,true);\
249  printStatus(ir);\
250  iret |= ir; }
251
252 #endif
253
254 // wrapper for the beta
255 double mbeta_pdf(double x, double a, double b, double){
256  return beta_pdf(x,a,b);
257 }
258 double mbeta_cdf(double x, double a, double b, double){
259  return beta_cdf(x,a,b);
260 }
261 double mbeta_cdf_c(double x, double a, double b, double){
262  return beta_cdf_c(x,a,b);
263 }
264 double mbeta_quantile(double x, double a, double b){
265  return beta_quantile(x,a,b);
266 }
267 double mbeta_quantile_c(double x, double a, double b){
268  return beta_quantile_c(x,a,b);
269 }
270
271
272 #ifndef NO_MATHCORE
273
274 int testPoissonCdf(double mu,double tol) {
275  int iret = 0;
276  for (int i = 0; i < 12; ++i) {
277  double q1 = poisson_cdf(i,mu);
278  double q1c = 1.- poisson_cdf_c(i,mu);
279  double q2 = 0;
280  for (int j = 0; j <= i; ++j) {
281  q2 += poisson_pdf(j,mu);
282  }
283  iret |= compare("test cdf",q1,q2,tol);
284  iret |= compare("test cdf_c",q1c,q2,tol);
285  }
286  printStatus(iret);
287  return iret;
288 }
289
290 int testBinomialCdf(double p, int n, double tol) {
291  int iret = 0;
292  for (int i = 0; i < 12; ++i) {
293  double q1 = binomial_cdf(i,p,n);
294  double q1c = 1.- binomial_cdf_c(i,p,n);
295  double q2 = 0;
296  for (int j = 0; j <= i; ++j) {
297  q2 += binomial_pdf(j,p,n);
298  }
299  iret |= compare("test cdf",q1,q2,tol);
300  iret |= compare("test cdf_c",q1c,q2,tol);
301  }
302  printStatus(iret);
303  return iret;
304 }
305
306 #endif
307
309
310  int iret = 0;
311  double tol = 2;
312
313
314 #ifndef NO_MATHCORE
315
316  tol = 8;
317  std::cout << "Poisson distrib. \t: ";
318  double mu = 5;
319  iret |= testPoissonCdf(mu,tol);
320
321  tol = 32;
322  std::cout << "Binomial distrib. \t: ";
323  double p = 0.5; int nt = 9;
324  iret |= testBinomialCdf(p,nt,tol);
325
326  tol = 2;
327  std::cout << "BreitWigner distrib\t: ";
328  TESTDIST1(breitwigner,1.0,tol);
329
330  std::cout << "Cauchy distribution\t: ";
331  TESTDIST1(cauchy,2.0,tol);
332
333  std::cout << "Exponential distrib\t: ";
334  TESTDIST1(exponential,1.0,tol);
335
336  std::cout << "Gaussian distribution\t: ";
337  TESTDIST1(gaussian,1.0,tol);
338
339  std::cout << "Log-normal distribution\t: ";
340  TESTDIST2(lognormal,1.0,1.0,tol);
341
342  std::cout << "Normal distribution\t: ";
343  TESTDIST1(normal,1.0,tol);
344
345  std::cout << "Uniform distribution\t: ";
346  TESTDIST2(uniform,0.0,10.0,tol);
347
348
349
350 #endif
351
352  std::cout << "Chisquare distribution\t: ";
353  tol = 8;
354  TESTDIST1(chisquared,9.,tol);
355
356  tol = 2;
357
358  std::cout << "F distribution\t\t: ";
359  double n = 5; double m = 6;
360  TESTDIST2(fdistribution,n,m,tol);
361
362  std::cout << "Gamma distribution\t: ";
363  double a = 1; double b = 2;
364  TESTDIST2(gamma,a,b,tol);
365
366  std::cout << "t distribution\t\t: ";
367  double nu = 10;
368  TESTDIST1(tdistribution,nu,tol);
369
370  std::cout << "Beta distribution\t: ";
371  a = 2; b = 1;
372  TESTDIST2(mbeta,a,b,tol);
373
374
375
376  return iret;
377 }
378 int main() {
379
380  return testStatFunc();
381 }
double binomial_cdf(unsigned int k, double p, unsigned int n)
Cumulative distribution function of the Binomial distribution Lower tail of the integral of the binom...
double gaussian(const double *x, const double *p)
double IntegralUp(const IGenFunction &f, double a)
evaluate the Integral of a function f over the semi-infinite interval (a,+inf)
double beta_pdf(double x, double a, double b)
Probability density function of the beta distribution.
double mbeta_quantile_c(double x, double a, double b)
double poisson_cdf(unsigned int n, double mu)
Cumulative distribution function of the Poisson distribution Lower tail of the integral of the poisso...
double poisson_cdf_c(unsigned int n, double mu)
Complement of the cumulative distribution function of the Poisson distribution.
double binomial_cdf_c(unsigned int k, double p, unsigned int n)
Complement of the cumulative distribution function of the Binomial distribution.
double mbeta_cdf_c(double x, double a, double b, double)
TArc * a
Definition: textangle.C:12
int main()
double mbeta_pdf(double x, double a, double b, double)
Float_t delta
TRObject operator()(const T1 &t1) const
Template class to wrap any C++ callable object which takes one argument i.e.
Double_t x[n]
Definition: legend1.C:17
static double p2(double t, double a, double b, double c)
int testStatFunc()
Class for performing numerical integration of a function in one dimension.
Definition: GSLIntegrator.h:90
double beta_cdf(double x, double a, double b)
Cumulative distribution function of the beta distribution Upper tail of the integral of the beta_pdf...
int compare(std::string name, double v1, double v2, double scale=2.0)
double IntegralLow(const IGenFunction &f, double b)
evaluate the Integral of a function f over the over the semi-infinite interval (-inf,b)
double gamma(double x)
const double tol
int testBinomialCdf(double p, int n, double tol)
double binomial_pdf(unsigned int k, double p, unsigned int n)
Probability density function of the binomial distribution.
void SetFunction(const IGenFunction &f)
method to set the a generic integration function
TMarker * m
Definition: textangle.C:8
double mbeta_quantile(double x, double a, double b)
static double p1(double t, double a, double b)
REAL epsilon
Definition: triangle.c:617
constexpr Double_t E()
Definition: TMath.h:74
double beta_quantile_c(double x, double a, double b)
Inverse ( ) of the cumulative distribution function of the lower tail of the beta distribution (beta_...
double f(double x)
int testPoissonCdf(double mu, double tol)
double mbeta_cdf(double x, double a, double b, double)
void printStatus(int iret)
double poisson_pdf(unsigned int n, double mu)
Probability density function of the Poisson distribution.
double beta_cdf_c(double x, double a, double b)
Complement of the cumulative distribution function of the beta distribution.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
def normal(shape, name=None)
#define TESTDIST2(name, p1, p2, s)
float * q
Definition: THbookFile.cxx:87
const Int_t n
Definition: legend1.C:16
double beta_quantile(double x, double a, double b)
Inverse ( ) of the cumulative distribution function of the upper tail of the beta distribution (beta_...
#define TESTDIST1(name, p1, s)