ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
VavilovAccurate.cxx
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Authors: B. List 29.4.2010
3 
4 
5  /**********************************************************************
6  * *
7  * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
8  * *
9  * This library is free software; you can redistribute it and/or *
10  * modify it under the terms of the GNU General Public License *
11  * as published by the Free Software Foundation; either version 2 *
12  * of the License, or (at your option) any later version. *
13  * *
14  * This library is distributed in the hope that it will be useful, *
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
17  * General Public License for more details. *
18  * *
19  * You should have received a copy of the GNU General Public License *
20  * along with this library (see file COPYING); if not, write *
21  * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
22  * 330, Boston, MA 02111-1307 USA, or contact the author. *
23  * *
24  **********************************************************************/
25 
26 // Implementation file for class VavilovAccurate
27 //
28 // Created by: blist at Thu Apr 29 11:19:00 2010
29 //
30 // Last update: Thu Apr 29 11:19:00 2010
31 //
32 
33 
34 #include "Math/VavilovAccurate.h"
35 #include "Math/SpecFuncMathCore.h"
36 #include "Math/SpecFuncMathMore.h"
37 #include "Math/QuantFuncMathCore.h"
38 
39 #include <cassert>
40 #include <iostream>
41 #include <iomanip>
42 #include <cmath>
43 #include <limits>
44 
45 
46 namespace ROOT {
47 namespace Math {
48 
49 VavilovAccurate *VavilovAccurate::fgInstance = 0;
50 
51 
52 VavilovAccurate::VavilovAccurate(double kappa, double beta2, double epsilonPM, double epsilon)
53 {
54  Set (kappa, beta2, epsilonPM, epsilon);
55 }
56 
57 
59 {
60  // desctructor (clean up resources)
61 }
62 
63 void VavilovAccurate::SetKappaBeta2 (double kappa, double beta2) {
64  Set (kappa, beta2);
65 }
66 
67 void VavilovAccurate::Set(double kappa, double beta2, double epsilonPM, double epsilon) {
68  // Method described in
69  // B. Schorr, Programs for the Landau and the Vavilov distributions and the corresponding random numbers,
70  // <A HREF="http://dx.doi.org/10.1016/0010-4655(74)90091-5">Computer Phys. Comm. 7 (1974) 215-224</A>.
71  fQuantileInit = false;
72 
73  fKappa = kappa;
74  fBeta2 = beta2;
75  fEpsilonPM = epsilonPM; // epsilon_+ = epsilon_-: determines support (T0, T1)
76  fEpsilon = epsilon;
77 
78  static const double eu = 0.577215664901532860606; // Euler's constant
79  static const double pi2 = 6.28318530717958647693, // 2pi
80  rpi = 0.318309886183790671538, // 1/pi
81  pih = 1.57079632679489661923; // pi/2
82  double h1 = -std::log(fEpsilon)-1.59631259113885503887; // -ln(fEpsilon) + ln(2/pi**2)
83  double deltaEpsilon = 0.001;
84  static const double logdeltaEpsilon = -std::log(deltaEpsilon); // 3 ln 10 = -ln(.001);
85  double logEpsilonPM = std::log(fEpsilonPM);
86  static const double eps = 1e-5; // accuracy of root finding for x0
87 
88  double xp[9] = {0,
89  9.29, 2.47, 0.89, 0.36, 0.15, 0.07, 0.03, 0.02};
90  double xq[7] = {0,
91  0.012, 0.03, 0.08, 0.26, 0.87, 3.83};
92 
93  if (kappa < 0.001) {
94  std::cerr << "VavilovAccurate::Set: kappa = " << kappa << " - out of range" << std::endl;
95  if (kappa < 0.001) kappa = 0.001;
96  }
97  if (beta2 < 0 || beta2 > 1) {
98  std::cerr << "VavilovAccurate::Set: beta2 = " << beta2 << " - out of range" << std::endl;
99  if (beta2 < 0) beta2 = -beta2;
100  if (beta2 > 1) beta2 = 1;
101  }
102 
103  // approximation of x_-
104  fH[5] = 1-beta2*(1-eu)-logEpsilonPM/kappa; // eq. 3.9
105  fH[6] = beta2;
106  fH[7] = 1-beta2;
107  double h4 = logEpsilonPM/kappa-(1+beta2*eu);
108  double logKappa = std::log(kappa);
109  double kappaInv = 1/kappa;
110  // Calculate T0 from Eq. (3.6), using x_- = fH[5]
111 // double e1h5 = (fH[5] > 40 ) ? 0 : -ROOT::Math::expint (-fH[5]);
112 // fT0 = (h4-fH[5]*logKappa-(fH[5]+beta2)*(std::log(fH[5])+e1h5)+std::exp(-fH[5]))/fH[5];
113  fT0 = (h4-fH[5]*logKappa-(fH[5]+beta2)*E1plLog(fH[5])+std::exp(-fH[5]))/fH[5];
114  int lp = 1;
115  while (lp < 9 && kappa < xp[lp]) ++lp;
116  int lq = 1;
117  while (lq < 7 && kappa >= xq[lq]) ++lq;
118  // Solve eq. 3.7 to get fH[0] = x_+
119  double delta = 0;
120  int ifail = 0;
121  do {
122  ifail = Rzero(-lp-0.5-delta,lq-7.5+delta,fH[0],eps,1000,&ROOT::Math::VavilovAccurate::G116f2);
123  delta += 0.5;
124  } while (ifail == 2);
125 
126  double q = 1/fH[0];
127  // Calculate T1 from Eq. (3.6)
128 // double e1h0 = (fH[0] > 40 ) ? 0 : -ROOT::Math::expint (-fH[0]);
129 // fT1 = h4*q-logKappa-(1+beta2*q)*(std::log(std::fabs(fH[0]))+e1h0)+std::exp(-fH[0])*q;
130  fT1 = h4*q-logKappa-(1+beta2*q)*E1plLog(fH[0])+std::exp(-fH[0])*q;
131 
132  fT = fT1-fT0; // Eq. (2.5)
133  fOmega = pi2/fT; // Eq. (2.5)
134  fH[1] = kappa*(2+beta2*eu)+h1;
135  if(kappa >= 0.07) fH[1] += logdeltaEpsilon; // reduce fEpsilon by a factor .001 for large kappa
136  fH[2] = beta2*kappa;
137  fH[3] = kappaInv*fOmega;
138  fH[4] = pih*fOmega;
139 
140  // Solve log(eq. (4.10)) to get fX0 = N
142 // if (ifail) {
143 // std::cerr << "Rzero failed for x0: ifail=" << ifail << ", kappa=" << kappa << ", beta2=" << beta2 << std::endl;
144 // std::cerr << "G116f1(" << 5. << ")=" << G116f1(5.) << ", G116f1(" << MAXTERMS << ")=" << G116f1(MAXTERMS) << std::endl;
145 // std::cerr << "fH[0]=" << fH[0] << ", fH[1]=" << fH[1] << ", fH[2]=" << fH[2] << ", fH[3]=" << fH[3] << ", fH[4]=" << fH[4] << std::endl;
146 // std::cerr << "fH[5]=" << fH[5] << ", fH[6]=" << fH[6] << ", fH[7]=" << fH[7] << std::endl;
147 // std::cerr << "fT0=" << fT0 << ", fT1=" << fT1 << std::endl;
148 // std::cerr << "G116f2(" << fH[0] << ")=" << G116f2(fH[0]) << std::endl;
149 // }
150  if (ifail == 2) {
151  fX0 = (G116f1(5) > G116f1(MAXTERMS)) ? MAXTERMS : 5;
152  }
153  if (fX0 < 5) fX0 = 5;
154  else if (fX0 > MAXTERMS) fX0 = MAXTERMS;
155  int n = int(fX0+1);
156  // logKappa=log(kappa)
157  double d = rpi*std::exp(kappa*(1+beta2*(eu-logKappa)));
158  fA_pdf[n] = rpi*fOmega;
159  fA_cdf[n] = 0;
160  q = -1;
161  double q2 = 2;
162  for (int k = 1; k < n; ++k) {
163  int l = n-k;
164  double x = fOmega*k;
165  double x1 = kappaInv*x;
166  double c1 = std::log(x)-ROOT::Math::cosint(x1);
167  double c2 = ROOT::Math::sinint(x1);
168  double c3 = std::sin(x1);
169  double c4 = std::cos(x1);
170  double xf1 = kappa*(beta2*c1-c4)-x*c2;
171  double xf2 = x*(c1 + fT0) + kappa*(c3+beta2*c2);
172  double d1 = q*d*fOmega*std::exp(xf1);
173  double s = std::sin(xf2);
174  double c = std::cos(xf2);
175  fA_pdf[l] = d1*c;
176  fB_pdf[l] = -d1*s;
177  d1 = q*d*std::exp(xf1)/k;
178  fA_cdf[l] = d1*s;
179  fB_cdf[l] = d1*c;
180  fA_cdf[n] += q2*fA_cdf[l];
181  q = -q;
182  q2 = -q2;
183  }
184 }
185 
187  fQuantileInit = true;
188 
189  fNQuant = 16;
190  // for kappa<0.02: use landau_quantile as first approximation
191  if (fKappa < 0.02) return;
192  else if (fKappa < 0.05) fNQuant = 32;
193 
194  // crude approximation for the median:
195 
196  double estmedian = -4.22784335098467134e-01-std::log(fKappa)-fBeta2;
197  if (estmedian>1.3) estmedian = 1.3;
198 
199  // distribute test values evenly below and above the median
200  for (int i = 1; i < fNQuant/2; ++i) {
201  double x = fT0 + i*(estmedian-fT0)/(fNQuant/2);
202  fQuant[i] = Cdf(x);
203  fLambda[i] = x;
204  }
205  for (int i = fNQuant/2; i < fNQuant-1; ++i) {
206  double x = estmedian + (i-fNQuant/2)*(fT1-estmedian)/(fNQuant/2-1);
207  fQuant[i] = Cdf(x);
208  fLambda[i] = x;
209  }
210 
211  fQuant[0] = 0;
212  fLambda[0] = fT0;
213  fQuant[fNQuant-1] = 1;
214  fLambda[fNQuant-1] = fT1;
215 
216 }
217 
218 double VavilovAccurate::Pdf (double x) const {
219 
220  static const double pi = 3.14159265358979323846; // pi
221 
222  int n = int(fX0);
223  double f;
224  if (x < fT0) {
225  f = 0;
226  } else if (x <= fT1) {
227  double y = x-fT0;
228  double u = fOmega*y-pi;
229  double cof = 2*cos(u);
230  double a1 = 0;
231  double a0 = fA_pdf[1];
232  double a2 = 0;
233  for (int k = 2; k <= n+1; ++k) {
234  a2 = a1;
235  a1 = a0;
236  a0 = fA_pdf[k]+cof*a1-a2;
237  }
238  double b1 = 0;
239  double b0 = fB_pdf[1];
240  for (int k = 2; k <= n; ++k) {
241  double b2 = b1;
242  b1 = b0;
243  b0 = fB_pdf[k]+cof*b1-b2;
244  }
245  f = 0.5*(a0-a2)+b0*sin(u);
246  } else {
247  f = 0;
248  }
249  return f;
250 }
251 
252 double VavilovAccurate::Pdf (double x, double kappa, double beta2) {
253  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
254  return Pdf (x);
255 }
256 
257 double VavilovAccurate::Cdf (double x) const {
258 
259  static const double pi = 3.14159265358979323846; // pi
260 
261  int n = int(fX0);
262  double f;
263  if (x < fT0) {
264  f = 0;
265  } else if (x <= fT1) {
266  double y = x-fT0;
267  double u = fOmega*y-pi;
268  double cof = 2*cos(u);
269  double a1 = 0;
270  double a0 = fA_cdf[1];
271  double a2 = 0;
272  for (int k = 2; k <= n+1; ++k) {
273  a2 = a1;
274  a1 = a0;
275  a0 = fA_cdf[k]+cof*a1-a2;
276  }
277  double b1 = 0;
278  double b0 = fB_cdf[1];
279  for (int k = 2; k <= n; ++k) {
280  double b2 = b1;
281  b1 = b0;
282  b0 = fB_cdf[k]+cof*b1-b2;
283  }
284  f = 0.5*(a0-a2)+b0*sin(u);
285  f += y/fT;
286  } else {
287  f = 1;
288  }
289  return f;
290 }
291 
292 double VavilovAccurate::Cdf (double x, double kappa, double beta2) {
293  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
294  return Cdf (x);
295 }
296 
297 double VavilovAccurate::Cdf_c (double x) const {
298 
299  static const double pi = 3.14159265358979323846; // pi
300 
301  int n = int(fX0);
302  double f;
303  if (x < fT0) {
304  f = 1;
305  } else if (x <= fT1) {
306  double y = fT1-x;
307  double u = fOmega*y-pi;
308  double cof = 2*cos(u);
309  double a1 = 0;
310  double a0 = fA_cdf[1];
311  double a2 = 0;
312  for (int k = 2; k <= n+1; ++k) {
313  a2 = a1;
314  a1 = a0;
315  a0 = fA_cdf[k]+cof*a1-a2;
316  }
317  double b1 = 0;
318  double b0 = fB_cdf[1];
319  for (int k = 2; k <= n; ++k) {
320  double b2 = b1;
321  b1 = b0;
322  b0 = fB_cdf[k]+cof*b1-b2;
323  }
324  f = -0.5*(a0-a2)+b0*sin(u);
325  f += y/fT;
326  } else {
327  f = 0;
328  }
329  return f;
330 }
331 
332 double VavilovAccurate::Cdf_c (double x, double kappa, double beta2) {
333  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
334  return Cdf_c (x);
335 }
336 
337 double VavilovAccurate::Quantile (double z) const {
338  if (z < 0 || z > 1) return std::numeric_limits<double>::signaling_NaN();
339 
340  if (!fQuantileInit) InitQuantile();
341 
342  double x;
343  if (fKappa < 0.02) {
345  if (x < fT0+5*fEpsilon) x = fT0+5*fEpsilon;
346  else if (x > fT1-10*fEpsilon) x = fT1-10*fEpsilon;
347  }
348  else {
349  // yes, I know what a binary search is, but linear search is faster for small n!
350  int i = 1;
351  while (z > fQuant[i]) ++i;
352  assert (i < fNQuant);
353 
354  assert (i >= 1);
355  assert (i < fNQuant);
356 
357  // initial solution
358  double f = (z-fQuant[i-1])/(fQuant[i]-fQuant[i-1]);
359  assert (f >= 0);
360  assert (f <= 1);
361  assert (fQuant[i] > fQuant[i-1]);
362 
363  x = f*fLambda[i] + (1-f)*fLambda[i-1];
364  }
365  if (fabs(x-fT0) < fEpsilon || fabs(x-fT1) < fEpsilon) return x;
366 
367  assert (x > fT0 && x < fT1);
368  double dx;
369  int n = 0;
370  do {
371  ++n;
372  double y = Cdf(x)-z;
373  double y1 = Pdf (x);
374  dx = - y/y1;
375  x = x + dx;
376  // protect against shooting beyond the support
377  if (x < fT0) x = 0.5*(fT0+x-dx);
378  else if (x > fT1) x = 0.5*(fT1+x-dx);
379  assert (x > fT0 && x < fT1);
380  } while (fabs(dx) > fEpsilon && n < 100);
381  return x;
382 }
383 
384 double VavilovAccurate::Quantile (double z, double kappa, double beta2) {
385  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
386  return Quantile (z);
387 }
388 
389 double VavilovAccurate::Quantile_c (double z) const {
390  if (z < 0 || z > 1) return std::numeric_limits<double>::signaling_NaN();
391 
392  if (!fQuantileInit) InitQuantile();
393 
394  double z1 = 1-z;
395 
396  double x;
397  if (fKappa < 0.02) {
399  if (x < fT0+5*fEpsilon) x = fT0+5*fEpsilon;
400  else if (x > fT1-10*fEpsilon) x = fT1-10*fEpsilon;
401  }
402  else {
403  // yes, I know what a binary search is, but linear search is faster for small n!
404  int i = 1;
405  while (z1 > fQuant[i]) ++i;
406  assert (i < fNQuant);
407 
408 // int i0=0, i1=fNQuant, i;
409 // for (int it = 0; it < LOG2fNQuant; ++it) {
410 // i = (i0+i1)/2;
411 // if (z > fQuant[i]) i0 = i;
412 // else i1 = i;
413 // }
414 // assert (i1-i0 == 1);
415 
416  assert (i >= 1);
417  assert (i < fNQuant);
418 
419  // initial solution
420  double f = (z1-fQuant[i-1])/(fQuant[i]-fQuant[i-1]);
421  assert (f >= 0);
422  assert (f <= 1);
423  assert (fQuant[i] > fQuant[i-1]);
424 
425  x = f*fLambda[i] + (1-f)*fLambda[i-1];
426  }
427  if (fabs(x-fT0) < fEpsilon || fabs(x-fT1) < fEpsilon) return x;
428 
429  assert (x > fT0 && x < fT1);
430  double dx;
431  int n = 0;
432  do {
433  ++n;
434  double y = Cdf_c(x)-z;
435  double y1 = -Pdf (x);
436  dx = - y/y1;
437  x = x + dx;
438  // protect against shooting beyond the support
439  if (x < fT0) x = 0.5*(fT0+x-dx);
440  else if (x > fT1) x = 0.5*(fT1+x-dx);
441  assert (x > fT0 && x < fT1);
442  } while (fabs(dx) > fEpsilon && n < 100);
443  return x;
444 }
445 
446 double VavilovAccurate::Quantile_c (double z, double kappa, double beta2) {
447  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
448  return Quantile_c (z);
449 }
450 
452  if (!fgInstance) fgInstance = new VavilovAccurate (1, 1);
453  return fgInstance;
454 }
455 
456 VavilovAccurate *VavilovAccurate::GetInstance(double kappa, double beta2) {
457  if (!fgInstance) fgInstance = new VavilovAccurate (kappa, beta2);
458  else if (kappa != fgInstance->fKappa || beta2 != fgInstance->fBeta2) fgInstance->Set (kappa, beta2);
459  return fgInstance;
460 }
461 
462 double vavilov_accurate_pdf (double x, double kappa, double beta2) {
463  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
464  return vavilov->Pdf (x);
465 }
466 
467 double vavilov_accurate_cdf_c (double x, double kappa, double beta2) {
468  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
469  return vavilov->Cdf_c (x);
470 }
471 
472 double vavilov_accurate_cdf (double x, double kappa, double beta2) {
473  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
474  return vavilov->Cdf (x);
475 }
476 
477 double vavilov_accurate_quantile (double z, double kappa, double beta2) {
478  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
479  return vavilov->Quantile (z);
480 }
481 
482 double vavilov_accurate_quantile_c (double z, double kappa, double beta2) {
483  VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
484  return vavilov->Quantile_c (z);
485 }
486 
487 double VavilovAccurate::G116f1 (double x) const {
488  // fH[1] = kappa*(2+beta2*eu) -ln(fEpsilon) + ln(2/pi**2)
489  // fH[2] = beta2*kappa
490  // fH[3] = omwga/kappa
491  // fH[4] = pi/2 *fOmega
492  // log of Eq. (4.10)
493  return fH[1]+fH[2]*std::log(fH[3]*x)-fH[4]*x;
494 }
495 
496 double VavilovAccurate::G116f2 (double x) const {
497  // fH[5] = 1-beta2*(1-eu)-logEpsilonPM/kappa; // eq. 3.9
498  // fH[6] = beta2;
499  // fH[7] = 1-beta2;
500  // Eq. 3.7 of Schorr
501 // return fH[5]-x+fH[6]*(std::log(std::fabs(x))-ROOT::Math::expint (-x))-fH[7]*std::exp(-x);
502  return fH[5]-x+fH[6]*E1plLog(x)-fH[7]*std::exp(-x);
503 }
504 
505 int VavilovAccurate::Rzero (double a, double b, double& x0,
506  double eps, int mxf, double (VavilovAccurate::*f)(double)const) const {
507 
508  double xa, xb, fa, fb, r;
509 
510  if (a <= b) {
511  xa = a;
512  xb = b;
513  } else {
514  xa = b;
515  xb = a;
516  }
517  fa = (this->*f)(xa);
518  fb = (this->*f)(xb);
519 
520  if(fa*fb > 0) {
521  r = -2*(xb-xa);
522  x0 = 0;
523 // std::cerr << "VavilovAccurate::Rzero: fail=2, f(" << a << ")=" << (this->*f) (a)
524 // << ", f(" << b << ")=" << (this->*f) (b) << std::endl;
525  return 2;
526  }
527  int mc = 0;
528 
529  bool recalcF12 = true;
530  bool recalcFab = true;
531  bool fail = false;
532 
533 
534  double x1=0, x2=0, f1=0, f2=0, fx=0, ee=0;
535  do {
536  if (recalcF12) {
537  x0 = 0.5*(xa+xb);
538  r = x0-xa;
539  ee = eps*(std::abs(x0)+1);
540  if(r <= ee) break;
541  f1 = fa;
542  x1 = xa;
543  f2 = fb;
544  x2 = xb;
545  }
546  if (recalcFab) {
547  fx = (this->*f)(x0);
548  ++mc;
549  if(mc > mxf) {
550  fail = true;
551  break;
552  }
553  if(fx*fa > 0) {
554  xa = x0;
555  fa = fx;
556  } else {
557  xb = x0;
558  fb = fx;
559  }
560  }
561  recalcF12 = true;
562  recalcFab = true;
563 
564  double u1 = f1-f2;
565  double u2 = x1-x2;
566  double u3 = f2-fx;
567  double u4 = x2-x0;
568  if(u2 == 0 || u4 == 0) continue;
569  double f3 = fx;
570  double x3 = x0;
571  u1 = u1/u2;
572  u2 = u3/u4;
573  double ca = u1-u2;
574  double cb = (x1+x2)*u2-(x2+x0)*u1;
575  double cc = (x1-x0)*f1-x1*(ca*x1+cb);
576  if(ca == 0) {
577  if(cb == 0) continue;
578  x0 = -cc/cb;
579  } else {
580  u3 = cb/(2*ca);
581  u4 = u3*u3-cc/ca;
582  if(u4 < 0) continue;
583  x0 = -u3 + (x0+u3 >= 0 ? +1 : -1)*std::sqrt(u4);
584  }
585  if(x0 < xa || x0 > xb) continue;
586 
587  recalcF12 = false;
588 
589  r = std::abs(x0-x3) < std::abs(x0-x2) ? std::abs(x0-x3) : std::abs(x0-x2);
590  ee = eps*(std::abs(x0)+1);
591  if (r > ee) {
592  f1 = f2;
593  x1 = x2;
594  f2 = f3;
595  x2 = x3;
596  continue;
597  }
598 
599  recalcFab = false;
600 
601  fx = (this->*f) (x0);
602  if (fx == 0) break;
603  double xx, ff;
604  if (fx*fa < 0) {
605  xx = x0-ee;
606  if (xx <= xa) break;
607  ff = (this->*f)(xx);
608  fb = ff;
609  xb = xx;
610  } else {
611  xx = x0+ee;
612  if(xx >= xb) break;
613  ff = (this->*f)(xx);
614  fa = ff;
615  xa = xx;
616  }
617  if (fx*ff <= 0) break;
618  mc += 2;
619  if (mc > mxf) {
620  fail = true;
621  break;
622  }
623  f1 = f3;
624  x1 = x3;
625  f2 = fx;
626  x2 = x0;
627  x0 = xx;
628  fx = ff;
629  }
630  while (true);
631 
632  if (fail) {
633  r = -0.5*std::abs(xb-xa);
634  x0 = 0;
635  std::cerr << "VavilovAccurate::Rzero: fail=" << fail << ", f(" << a << ")=" << (this->*f) (a)
636  << ", f(" << b << ")=" << (this->*f) (b) << std::endl;
637  return 1;
638  }
639 
640  r = ee;
641  return 0;
642 }
643 
644 // Calculates log(|x|)+E_1(x)
645 double VavilovAccurate::E1plLog (double x) {
646  static const double eu = 0.577215664901532860606; // Euler's constant
647  double absx = std::fabs(x);
648  if (absx < 1E-4) {
649  return (x-0.25*x)*x-eu;
650  }
651  else if (x > 35) {
652  return log (x);
653  }
654  else if (x < -50) {
655  return -ROOT::Math::expint (-x);
656  }
657  return log(absx) -ROOT::Math::expint (-x);
658 }
659 
661  return fT0;
662 }
663 
665  return fT1;
666 }
667 
669  return fKappa;
670 }
671 
673  return fBeta2;
674 }
675 
676 double VavilovAccurate::Mode() const {
677  double x = -4.22784335098467134e-01-std::log(fKappa)-fBeta2;
678  if (x>-0.223172) x = -0.223172;
679  double eps = 0.01;
680  double dx;
681 
682  do {
683  double p0 = Pdf (x - eps);
684  double p1 = Pdf (x);
685  double p2 = Pdf (x + eps);
686  double y1 = 0.5*(p2-p0)/eps;
687  double y2 = (p2-2*p1+p0)/(eps*eps);
688  dx = - y1/y2;
689  x += dx;
690  if (fabs(dx) < eps) eps = 0.1*fabs(dx);
691  } while (fabs(dx) > 1E-5);
692  return x;
693 }
694 
695 double VavilovAccurate::Mode(double kappa, double beta2) {
696  if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
697  return Mode();
698 }
699 
701  return fEpsilonPM;
702 }
703 
705  return fEpsilon;
706 }
707 
709  return fX0;
710 }
711 
712 
713 
714 } // namespace Math
715 } // namespace ROOT
double cosint(double x)
Calculates the real part of the cosine integral (Ci).
double G116f1(double x) const
const double pi
return c
TCanvas * c1
Definition: legend1.C:2
#define assert(cond)
Definition: unittest.h:542
double sinint(double x)
Calculates the sine integral.
double Cdf_c(double x) const
Evaluate the Vavilov complementary cummulative probability density function.
tuple f2
Definition: surfaces.py:24
int * lq
Definition: THbookFile.cxx:86
void f3()
Definition: na49.C:50
int Rzero(double a, double b, double &x0, double eps, int mxf, double(VavilovAccurate::*f)(double) const) const
TArc * a
Definition: textangle.C:12
double cos(double)
TFile * f
double vavilov_accurate_cdf_c(double x, double kappa, double beta2)
The Vavilov complementary cummulative probability density function.
double sqrt(double)
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
virtual void SetKappaBeta2(double kappa, double beta2)
Change and and recalculate coefficients if necessary.
double G116f2(double x) const
int d
Definition: tornado.py:11
static double p2(double t, double a, double b, double c)
static Vc_ALWAYS_INLINE Vector< T > abs(const Vector< T > &x)
Definition: vector.h:450
double GetEpsilonPM() const
Return the current value of .
TH1F * h1
Definition: legend1.C:5
double sin(double)
Float_t z[5]
Definition: Ifit.C:16
double Pdf(double x) const
Evaluate the Vavilov probability density function.
double Cdf(double x) const
Evaluate the Vavilov cummulative probability density function.
virtual double GetLambdaMax() const
Return the maximum value of for which is nonzero in the current approximation.
double expint(double x)
Calculates the exponential integral.
void Set(double kappa, double beta2, double epsilonPM=5E-4, double epsilon=1E-5)
(Re)Initialize the object
virtual double GetLambdaMin() const
Return the minimum value of for which is nonzero in the current approximation.
static VavilovAccurate * GetInstance()
Returns a static instance of class VavilovFast.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
ROOT::R::TRInterface & r
Definition: Object.C:4
static VavilovAccurate * fgInstance
Double_t E()
Definition: TMath.h:54
TLine * l
Definition: textangle.C:4
double Quantile_c(double z) const
Evaluate the inverse of the complementary Vavilov cummulative probability density function...
static double p1(double t, double a, double b)
VavilovAccurate(double kappa=1, double beta2=1, double epsilonPM=5E-4, double epsilon=1E-5)
Initialize an object to calculate the Vavilov distribution.
REAL epsilon
Definition: triangle.c:617
return c2
Definition: legend2.C:14
static const double x1[5]
double vavilov_accurate_pdf(double x, double kappa, double beta2)
The Vavilov probability density function.
Class describing a Vavilov distribution.
Double_t y[n]
Definition: legend1.C:17
double Quantile(double z) const
Evaluate the inverse of the Vavilov cummulative probability density function.
virtual double GetBeta2() const
Return the current value of .
static const double eu
Definition: Vavilov.cxx:52
virtual double Mode() const
Return the value of where the pdf is maximal.
static double E1plLog(double x)
double landau_quantile(double z, double xi=1)
Inverse ( ) of the cumulative distribution function of the lower tail of the Landau distribution (lan...
double GetNTerms() const
Return the number of terms used in the series expansion.
double fLambda[kNquantMax]
virtual ~VavilovAccurate()
Destructor.
TF1 * f1
Definition: legend1.C:11
virtual double GetKappa() const
Return the current value of .
double GetEpsilon() const
Return the current value of .
double vavilov_accurate_quantile_c(double z, double kappa, double beta2)
The inverse of the complementary Vavilov cummulative probability density function.
double exp(double)
double vavilov_accurate_quantile(double z, double kappa, double beta2)
The inverse of the Vavilov cummulative probability density function.
float * q
Definition: THbookFile.cxx:87
double vavilov_accurate_cdf(double x, double kappa, double beta2)
The Vavilov cummulative probability density function.
return c3
Definition: legend3.C:15
const Int_t n
Definition: legend1.C:16
double log(double)
static const double x3[11]