Logo ROOT   6.12/07
Reference Guide
SpecFuncCephesInv.cxx
Go to the documentation of this file.
1 // inverse of gamma and beta from Cephes library
2 // see: http://www.netlib.org/cephes
3 //
4 // Copyright 1985, 1987, 2000 by Stephen L. Moshier
5 
6 
7 
8 #include "Math/Error.h"
9 
10 #include "SpecFuncCephes.h"
11 
12 #include <cmath>
13 
14 #include <limits>
15 
16 namespace ROOT {
17 namespace Math {
18 
19 namespace Cephes {
20 
21 
22 
23 /*
24  *
25  * Inverse of Normal distribution function
26  *
27  *
28  *
29  * SYNOPSIS:
30  *
31  * double x, y, ndtri();
32  *
33  * x = ndtri( y );
34  *
35  *
36  *
37  * DESCRIPTION:
38  *
39  * Returns the argument, x, for which the area under the
40  * Gaussian probability density function (integrated from
41  * minus infinity to x) is equal to y.
42  *
43  *
44  * For small arguments 0 < y < exp(-2), the program computes
45  * z = sqrt( -2.0 * log(y) ); then the approximation is
46  * x = z - log(z)/z - (1/z) P(1/z) / Q(1/z).
47  * There are two rational functions P/Q, one for 0 < y < exp(-32)
48  * and the other for y up to exp(-2). For larger arguments,
49  * w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
50  *
51  *
52  * ACCURACY:
53  *
54  * Relative error:
55  * arithmetic domain # trials peak rms
56  * DEC 0.125, 1 5500 9.5e-17 2.1e-17
57  * DEC 6e-39, 0.135 3500 5.7e-17 1.3e-17
58  * IEEE 0.125, 1 20000 7.2e-16 1.3e-16
59  * IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17
60  *
61  *
62  * ERROR MESSAGES:
63  *
64  * message condition value returned
65  * ndtri domain x <= 0 -MAXNUM
66  * ndtri domain x >= 1 MAXNUM
67  *
68  */
69 
70 /*
71 Cephes Math Library Release 2.8: June, 2000
72 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
73 */
74 
75 
76 static double s2pi = 2.50662827463100050242E0;
77 
78 static double P0[5] = {
79 -5.99633501014107895267E1,
80  9.80010754185999661536E1,
81 -5.66762857469070293439E1,
82  1.39312609387279679503E1,
83 -1.23916583867381258016E0,
84 };
85 static double Q0[8] = {
86  1.95448858338141759834E0,
87  4.67627912898881538453E0,
88  8.63602421390890590575E1,
89 -2.25462687854119370527E2,
90  2.00260212380060660359E2,
91 -8.20372256168333339912E1,
92  1.59056225126211695515E1,
93 -1.18331621121330003142E0,
94 };
95 static double P1[9] = {
96  4.05544892305962419923E0,
97  3.15251094599893866154E1,
98  5.71628192246421288162E1,
99  4.40805073893200834700E1,
100  1.46849561928858024014E1,
101  2.18663306850790267539E0,
102 -1.40256079171354495875E-1,
103 -3.50424626827848203418E-2,
104 -8.57456785154685413611E-4,
105 };
106 static double Q1[8] = {
107  1.57799883256466749731E1,
108  4.53907635128879210584E1,
109  4.13172038254672030440E1,
110  1.50425385692907503408E1,
111  2.50464946208309415979E0,
112 -1.42182922854787788574E-1,
113 -3.80806407691578277194E-2,
114 -9.33259480895457427372E-4,
115 };
116 static double P2[9] = {
117  3.23774891776946035970E0,
118  6.91522889068984211695E0,
119  3.93881025292474443415E0,
120  1.33303460815807542389E0,
121  2.01485389549179081538E-1,
122  1.23716634817820021358E-2,
123  3.01581553508235416007E-4,
124  2.65806974686737550832E-6,
125  6.23974539184983293730E-9,
126 };
127 static double Q2[8] = {
128  6.02427039364742014255E0,
129  3.67983563856160859403E0,
130  1.37702099489081330271E0,
131  2.16236993594496635890E-1,
132  1.34204006088543189037E-2,
133  3.28014464682127739104E-4,
134  2.89247864745380683936E-6,
135  6.79019408009981274425E-9,
136 };
137 double ndtri( double y0 )
138 {
139  double x, y, z, y2, x0, x1;
140  int code;
141  if( y0 <= 0.0 )
142  return( - std::numeric_limits<double>::infinity() );
143  if( y0 >= 1.0 )
144  return( + std::numeric_limits<double>::infinity() );
145  code = 1;
146  y = y0;
147  if( y > (1.0 - 0.13533528323661269189) )
148  {
149  y = 1.0 - y;
150  code = 0;
151  }
152  if( y > 0.13533528323661269189 )
153  {
154  y = y - 0.5;
155  y2 = y * y;
156  x = y + y * (y2 * Polynomialeval( y2, P0, 4)/ Polynomial1eval( y2, Q0, 8 ));
157  x = x * s2pi;
158  return(x);
159  }
160  x = std::sqrt( -2.0 * std::log(y) );
161  x0 = x - std::log(x)/x;
162  z = 1.0/x;
163  if( x < 8.0 )
164  x1 = z * Polynomialeval( z, P1, 8 )/ Polynomial1eval ( z, Q1, 8 );
165  else
166  x1 = z * Polynomialeval( z, P2, 8 )/ Polynomial1eval( z, Q2, 8 );
167  x = x0 - x1;
168  if( code != 0 )
169  x = -x;
170  return( x );
171 }
172 
173 
174 
175 
176 /*
177  *
178  * Inverse of complemented imcomplete gamma integral
179  *
180  *
181  *
182  * SYNOPSIS:
183  *
184  * double a, x, p, igami();
185  *
186  * x = igami( a, p );
187  *
188  * DESCRIPTION:
189  *
190  * Given p, the function finds x such that
191  *
192  * igamc( a, x ) = p.
193  *
194  * Starting with the approximate value
195  *
196  * 3
197  * x = a t
198  *
199  * where
200  *
201  * t = 1 - d - ndtri(p) sqrt(d)
202  *
203  * and
204  *
205  * d = 1/9a,
206  *
207  * the routine performs up to 10 Newton iterations to find the
208  * root of igamc(a,x) - p = 0.
209  *
210  * ACCURACY:
211  *
212  * Tested at random a, p in the intervals indicated.
213  *
214  * a p Relative error:
215  * arithmetic domain domain # trials peak rms
216  * IEEE 0.5,100 0,0.5 100000 1.0e-14 1.7e-15
217  * IEEE 0.01,0.5 0,0.5 100000 9.0e-14 3.4e-15
218  * IEEE 0.5,10000 0,0.5 20000 2.3e-13 3.8e-14
219  */
220 /*
221 Cephes Math Library Release 2.8: June, 2000
222 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
223 */
224 
225 double igami( double a, double y0 )
226 {
227  double x0, x1, x, yl, yh, y, d, lgm, dithresh;
228  int i, dir;
229 
230  // check the domain
231  if (a<= 0) {
232  MATH_ERROR_MSG("Cephes::igami","Wrong domain for parameter a (must be > 0)");
233  return 0;
234  }
235  if (y0 <= 0) {
236  //if (y0<0) MATH_ERROR_MSG("Cephes::igami","Wrong domain for y (must be in [0,1])");
237  return std::numeric_limits<double>::infinity();
238  }
239  if (y0 >= 1) {
240  //if (y0>1) MATH_ERROR_MSG("Cephes::igami","Wrong domain for y (must be in [0,1])");
241  return 0;
242  }
243 
244 
245 /* bound the solution */
246  static double kMAXNUM = std::numeric_limits<double>::max();
247  x0 = kMAXNUM;
248  yl = 0;
249  x1 = 0;
250  yh = 1.0;
251  dithresh = 5.0 * kMACHEP;
252 
253 /* approximation to inverse function */
254  d = 1.0/(9.0*a);
255  y = ( 1.0 - d - ndtri(y0) * std::sqrt(d) );
256  x = a * y * y * y;
257 
258  lgm = lgam(a);
259 
260  for( i=0; i<10; i++ )
261  {
262  if( x > x0 || x < x1 )
263  goto ihalve;
264  y = igamc(a,x);
265  if( y < yl || y > yh )
266  goto ihalve;
267  if( y < y0 )
268  {
269  x0 = x;
270  yl = y;
271  }
272  else
273  {
274  x1 = x;
275  yh = y;
276  }
277 /* compute the derivative of the function at this point */
278  d = (a - 1.0) * std::log(x) - x - lgm;
279  if( d < -kMAXLOG )
280  goto ihalve;
281  d = -std::exp(d);
282 /* compute the step to the next approximation of x */
283  d = (y - y0)/d;
284  if( std::fabs(d/x) < kMACHEP )
285  goto done;
286  x = x - d;
287  }
288 
289 /* Resort to interval halving if Newton iteration did not converge. */
290 ihalve:
291 
292  d = 0.0625;
293  if( x0 == kMAXNUM )
294  {
295  if( x <= 0.0 )
296  x = 1.0;
297  while( x0 == kMAXNUM )
298  {
299  x = (1.0 + d) * x;
300  y = igamc( a, x );
301  if( y < y0 )
302  {
303  x0 = x;
304  yl = y;
305  break;
306  }
307  d = d + d;
308  }
309  }
310  d = 0.5;
311  dir = 0;
312 
313  for( i=0; i<400; i++ )
314  {
315  x = x1 + d * (x0 - x1);
316  y = igamc( a, x );
317  lgm = (x0 - x1)/(x1 + x0);
318  if( std::fabs(lgm) < dithresh )
319  break;
320  lgm = (y - y0)/y0;
321  if( std::fabs(lgm) < dithresh )
322  break;
323  if( x <= 0.0 )
324  break;
325  if( y >= y0 )
326  {
327  x1 = x;
328  yh = y;
329  if( dir < 0 )
330  {
331  dir = 0;
332  d = 0.5;
333  }
334  else if( dir > 1 )
335  d = 0.5 * d + 0.5;
336  else
337  d = (y0 - yl)/(yh - yl);
338  dir += 1;
339  }
340  else
341  {
342  x0 = x;
343  yl = y;
344  if( dir > 0 )
345  {
346  dir = 0;
347  d = 0.5;
348  }
349  else if( dir < -1 )
350  d = 0.5 * d;
351  else
352  d = (y0 - yl)/(yh - yl);
353  dir -= 1;
354  }
355  }
356 
357 // if( x == 0.0 )
358 // mtherr( "igami", UNDERFLOW );
359 
360 done:
361  return( x );
362 }
363 
364 
365 /*
366  *
367  * Inverse of imcomplete beta integral
368  *
369  *
370  *
371  * SYNOPSIS:
372  *
373  * double a, b, x, y, incbi();
374  *
375  * x = incbi( a, b, y );
376  *
377  *
378  *
379  * DESCRIPTION:
380  *
381  * Given y, the function finds x such that
382  *
383  * incbet( a, b, x ) = y .
384  *
385  * The routine performs interval halving or Newton iterations to find the
386  * root of incbet(a,b,x) - y = 0.
387  *
388  *
389  * ACCURACY:
390  *
391  * Relative error:
392  * x a,b
393  * arithmetic domain domain # trials peak rms
394  * IEEE 0,1 .5,10000 50000 5.8e-12 1.3e-13
395  * IEEE 0,1 .25,100 100000 1.8e-13 3.9e-15
396  * IEEE 0,1 0,5 50000 1.1e-12 5.5e-15
397  * VAX 0,1 .5,100 25000 3.5e-14 1.1e-15
398  * With a and b constrained to half-integer or integer values:
399  * IEEE 0,1 .5,10000 50000 5.8e-12 1.1e-13
400  * IEEE 0,1 .5,100 100000 1.7e-14 7.9e-16
401  * With a = .5, b constrained to half-integer or integer values:
402  * IEEE 0,1 .5,10000 10000 8.3e-11 1.0e-11
403  */
404 
405 /*
406 Cephes Math Library Release 2.8: June, 2000
407 Copyright 1984, 1996, 2000 by Stephen L. Moshier
408 */
409 
410 
411 double incbi( double aa, double bb, double yy0 )
412 {
413  double a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt;
414  int i, rflg, dir, nflg;
415 
416  // check the domain
417  if (aa<= 0) {
418  MATH_ERROR_MSG("Cephes::incbi","Wrong domain for parameter a (must be > 0)");
419  return 0;
420  }
421  if (bb<= 0) {
422  MATH_ERROR_MSG("Cephes::incbi","Wrong domain for parameter b (must be > 0)");
423  return 0;
424  }
425 
426 
427  i = 0;
428  if( yy0 <= 0 )
429  return(0.0);
430  if( yy0 >= 1.0 )
431  return(1.0);
432  x0 = 0.0;
433  yl = 0.0;
434  x1 = 1.0;
435  yh = 1.0;
436  nflg = 0;
437 
438  if( aa <= 1.0 || bb <= 1.0 )
439  {
440  dithresh = 1.0e-6;
441  rflg = 0;
442  a = aa;
443  b = bb;
444  y0 = yy0;
445  x = a/(a+b);
446  y = incbet( a, b, x );
447  goto ihalve;
448  }
449  else
450  {
451  dithresh = 1.0e-4;
452  }
453 /* approximation to inverse function */
454 
455  yp = -ndtri(yy0);
456 
457  if( yy0 > 0.5 )
458  {
459  rflg = 1;
460  a = bb;
461  b = aa;
462  y0 = 1.0 - yy0;
463  yp = -yp;
464  }
465  else
466  {
467  rflg = 0;
468  a = aa;
469  b = bb;
470  y0 = yy0;
471  }
472 
473  lgm = (yp * yp - 3.0)/6.0;
474  x = 2.0/( 1.0/(2.0*a-1.0) + 1.0/(2.0*b-1.0) );
475  d = yp * std::sqrt( x + lgm ) / x
476  - ( 1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0) )
477  * (lgm + 5.0/6.0 - 2.0/(3.0*x));
478  d = 2.0 * d;
479  if( d < kMINLOG )
480  {
481  x = 1.0;
482  goto under;
483  }
484  x = a/( a + b * std::exp(d) );
485  y = incbet( a, b, x );
486  yp = (y - y0)/y0;
487  if( std::fabs(yp) < 0.2 )
488  goto newt;
489 
490 /* Resort to interval halving if not close enough. */
491 ihalve:
492 
493  dir = 0;
494  di = 0.5;
495  for( i=0; i<100; i++ )
496  {
497  if( i != 0 )
498  {
499  x = x0 + di * (x1 - x0);
500  if( x == 1.0 )
501  x = 1.0 - kMACHEP;
502  if( x == 0.0 )
503  {
504  di = 0.5;
505  x = x0 + di * (x1 - x0);
506  if( x == 0.0 )
507  goto under;
508  }
509  y = incbet( a, b, x );
510  yp = (x1 - x0)/(x1 + x0);
511  if( std::fabs(yp) < dithresh )
512  goto newt;
513  yp = (y-y0)/y0;
514  if( std::fabs(yp) < dithresh )
515  goto newt;
516  }
517  if( y < y0 )
518  {
519  x0 = x;
520  yl = y;
521  if( dir < 0 )
522  {
523  dir = 0;
524  di = 0.5;
525  }
526  else if( dir > 3 )
527  di = 1.0 - (1.0 - di) * (1.0 - di);
528  else if( dir > 1 )
529  di = 0.5 * di + 0.5;
530  else
531  di = (y0 - y)/(yh - yl);
532  dir += 1;
533  if( x0 > 0.75 )
534  {
535  if( rflg == 1 )
536  {
537  rflg = 0;
538  a = aa;
539  b = bb;
540  y0 = yy0;
541  }
542  else
543  {
544  rflg = 1;
545  a = bb;
546  b = aa;
547  y0 = 1.0 - yy0;
548  }
549  x = 1.0 - x;
550  y = incbet( a, b, x );
551  x0 = 0.0;
552  yl = 0.0;
553  x1 = 1.0;
554  yh = 1.0;
555  goto ihalve;
556  }
557  }
558  else
559  {
560  x1 = x;
561  if( rflg == 1 && x1 < kMACHEP )
562  {
563  x = 0.0;
564  goto done;
565  }
566  yh = y;
567  if( dir > 0 )
568  {
569  dir = 0;
570  di = 0.5;
571  }
572  else if( dir < -3 )
573  di = di * di;
574  else if( dir < -1 )
575  di = 0.5 * di;
576  else
577  di = (y - y0)/(yh - yl);
578  dir -= 1;
579  }
580  }
581  //mtherr( "incbi", PLOSS );
582  if( x0 >= 1.0 )
583  {
584  x = 1.0 - kMACHEP;
585  goto done;
586  }
587  if( x <= 0.0 )
588  {
589  under:
590  //mtherr( "incbi", UNDERFLOW );
591  x = 0.0;
592  goto done;
593  }
594 
595 newt:
596 
597  if( nflg )
598  goto done;
599  nflg = 1;
600  lgm = lgam(a+b) - lgam(a) - lgam(b);
601 
602  for( i=0; i<8; i++ )
603  {
604  /* Compute the function at this point. */
605  if( i != 0 )
606  y = incbet(a,b,x);
607  if( y < yl )
608  {
609  x = x0;
610  y = yl;
611  }
612  else if( y > yh )
613  {
614  x = x1;
615  y = yh;
616  }
617  else if( y < y0 )
618  {
619  x0 = x;
620  yl = y;
621  }
622  else
623  {
624  x1 = x;
625  yh = y;
626  }
627  if( x == 1.0 || x == 0.0 )
628  break;
629  /* Compute the derivative of the function at this point. */
630  d = (a - 1.0) * std::log(x) + (b - 1.0) * std::log(1.0-x) + lgm;
631  if( d < kMINLOG )
632  goto done;
633  if( d > kMAXLOG )
634  break;
635  d = std::exp(d);
636  /* Compute the step to the next approximation of x. */
637  d = (y - y0)/d;
638  xt = x - d;
639  if( xt <= x0 )
640  {
641  y = (x - x0) / (x1 - x0);
642  xt = x0 + 0.5 * y * (x - x0);
643  if( xt <= 0.0 )
644  break;
645  }
646  if( xt >= x1 )
647  {
648  y = (x1 - x) / (x1 - x0);
649  xt = x1 - 0.5 * y * (x1 - x);
650  if( xt >= 1.0 )
651  break;
652  }
653  x = xt;
654  if( std::fabs(d/x) < 128.0 * kMACHEP )
655  goto done;
656  }
657 /* Did not converge. */
658  dithresh = 256.0 * kMACHEP;
659  goto ihalve;
660 
661 done:
662 
663  if( rflg )
664  {
665  if( x <= kMACHEP )
666  x = 1.0 - kMACHEP;
667  else
668  x = 1.0 - x;
669  }
670  return( x );
671 }
672 
673 } // end namespace Cephes
674 
675 } // end namespace Math
676 } // end namespace ROOT
static double P0[5]
#define kMAXLOG
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
static double P1[9]
double lgam(double x)
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
double incbi(double a, double b, double y)
double Polynomial1eval(double x, double *a, unsigned int N)
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
double Polynomialeval(double x, double *a, unsigned int N)
auto * a
Definition: textangle.C:12
double incbet(double aa, double bb, double xx)
DESCRIPTION:
double igamc(double a, double x)
incomplete complementary gamma function igamc(a, x) = 1 - igam(a, x)
double ndtri(double y)
static const double x1[5]
Double_t y[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
#define kMINLOG
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
static double Q0[8]
double igami(double a, double y)
double exp(double)
#define kMACHEP
double log(double)