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