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 <iomanip>
42#include <cmath>
43#include <limits>
44
45
46namespace ROOT {
47namespace Math {
48
49VavilovAccurate *VavilovAccurate::fgInstance = 0;
50
51
52VavilovAccurate::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
63void VavilovAccurate::SetKappaBeta2 (double kappa, double beta2) {
64 Set (kappa, beta2);
65}
66
67void 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)
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
218double 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
252double VavilovAccurate::Pdf (double x, double kappa, double beta2) {
253 if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
254 return Pdf (x);
255}
256
257double 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
292double VavilovAccurate::Cdf (double x, double kappa, double beta2) {
293 if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
294 return Cdf (x);
295}
296
297double 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
332double 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
337double VavilovAccurate::Quantile (double z) const {
338 if (z < 0 || z > 1) return std::numeric_limits<double>::signaling_NaN();
339
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
384double VavilovAccurate::Quantile (double z, double kappa, double beta2) {
385 if (kappa != fKappa || beta2 != fBeta2) Set (kappa, beta2);
386 return Quantile (z);
387}
388
389double VavilovAccurate::Quantile_c (double z) const {
390 if (z < 0 || z > 1) return std::numeric_limits<double>::signaling_NaN();
391
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
446double 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
456VavilovAccurate *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
462double vavilov_accurate_pdf (double x, double kappa, double beta2) {
463 VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
464 return vavilov->Pdf (x);
465}
466
467double 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
472double vavilov_accurate_cdf (double x, double kappa, double beta2) {
473 VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
474 return vavilov->Cdf (x);
475}
476
477double vavilov_accurate_quantile (double z, double kappa, double beta2) {
478 VavilovAccurate *vavilov = VavilovAccurate::GetInstance (kappa, beta2);
479 return vavilov->Quantile (z);
480}
481
482double 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
487double 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
496double 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
505int 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)
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
676double 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
695double 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
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:87
int * lq
Definition: THbookFile.cxx:86
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)
VSD Structures.
Definition: StringConv.hxx:21
static constexpr double s
static constexpr double pi
static constexpr double pi2
constexpr Double_t E()
Base of natural log:
Definition: TMath.h:97
auto * l
Definition: textangle.C:4
auto * a
Definition: textangle.C:12
REAL epsilon
Definition: triangle.c:617