/**************************************************************************** * chi2.c * * This module contains the function for the chi square distribution. * * All functions in this module are taken from the Cephes math library. * * Cephes Math Library Release 2.0: April, 1987 * Copyright 1984, 1987 by Stephen L. Moshier * Direct inquiries to 30 Frost Street, Cambridge, MA 02140 * * from Persistence of Vision(tm) Ray Tracer * Copyright 1996,1999 Persistence of Vision Team *--------------------------------------------------------------------------- * NOTICE: This source code file is provided so that users may experiment * with enhancements to POV-Ray and to port the software to platforms other * than those supported by the POV-Ray Team. There are strict rules under * which you are permitted to use this file. The rules are in the file * named POVLEGAL.DOC which should be distributed with this file. * If POVLEGAL.DOC is not available or for more info please contact the POV-Ray * Team Coordinator by email to team-coord@povray.org or visit us on the web at * http://www.povray.org. The latest version of POV-Ray may be found at this site. * * This program is based on the popular DKB raytracer version 2.12. * DKBTrace was originally written by David K. Buck. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins. * *****************************************************************************/ #include "frame.h" #include "povproto.h" #include "chi2.h" /* Cephes Math Library Release 2.0: April, 1987 Copyright 1984, 1987 by Stephen L. Moshier Direct inquiries to 30 Frost Street, Cambridge, MA 02140 */ /***************************************************************************** * Local preprocessor defines ******************************************************************************/ #define MAXLGM 2.556348e305 #define BIG 1.44115188075855872E+17 #define MACHEP 1.38777878078144567553E-17 /* 2**-56 */ #define MAXLOG 8.8029691931113054295988E1 /* log(2**127) */ #define MAXNUM 1.701411834604692317316873e38 /* 2**127 */ #define LOGPI 1.14472988584940017414 /***************************************************************************** * Local typedefs ******************************************************************************/ /***************************************************************************** * Local variables ******************************************************************************/ static int sgngam = 0; /* * A[]: Stirling's formula expansion of log gamma * B[], C[]: log gamma function between 2 and 3 */ static DBL A[] = { 8.11614167470508450300E-4, -5.95061904284301438324E-4, 7.93650340457716943945E-4, -2.77777777730099687205E-3, 8.33333333333331927722E-2 }; static DBL B[] = { -1.37825152569120859100E3, -3.88016315134637840924E4, -3.31612992738871184744E5, -1.16237097492762307383E6, -1.72173700820839662146E6, -8.53555664245765465627E5 }; static DBL C[] = { 1.00000000000000000000E0, -3.51815701436523470549E2, -1.70642106651881159223E4, -2.20528590553854454839E5, -1.13933444367982507207E6, -2.53252307177582951285E6, -2.01889141433532773231E6 }; /* log(sqrt(2pi)) */ static DBL LS2PI = 0.91893853320467274178; /* sqrt(2pi) */ static DBL s2pi = 2.50662827463100050242E0; /* approximation for 0 <= |y - 0.5| <= 3/8 */ static DBL P0[5] = { -5.99633501014107895267E1, 9.80010754185999661536E1, -5.66762857469070293439E1, 1.39312609387279679503E1, -1.23916583867381258016E0, }; static DBL Q0[8] = { /* 1.00000000000000000000E0,*/ 1.95448858338141759834E0, 4.67627912898881538453E0, 8.63602421390890590575E1, -2.25462687854119370527E2, 2.00260212380060660359E2, -8.20372256168333339912E1, 1.59056225126211695515E1, -1.18331621121330003142E0, }; /* * Approximation for interval z = sqrt(-2 log y ) between 2 and 8 * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14. */ static DBL P1[9] = { 4.05544892305962419923E0, 3.15251094599893866154E1, 5.71628192246421288162E1, 4.40805073893200834700E1, 1.46849561928858024014E1, 2.18663306850790267539E0, -1.40256079171354495875E-1, -3.50424626827848203418E-2, -8.57456785154685413611E-4, }; static DBL Q1[8] = { /* 1.00000000000000000000E0,*/ 1.57799883256466749731E1, 4.53907635128879210584E1, 4.13172038254672030440E1, 1.50425385692907503408E1, 2.50464946208309415979E0, -1.42182922854787788574E-1, -3.80806407691578277194E-2, -9.33259480895457427372E-4, }; /* * Approximation for interval z = sqrt(-2 log y ) between 8 and 64 * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890. */ static DBL P2[9] = { 3.23774891776946035970E0, 6.91522889068984211695E0, 3.93881025292474443415E0, 1.33303460815807542389E0, 2.01485389549179081538E-1, 1.23716634817820021358E-2, 3.01581553508235416007E-4, 2.65806974686737550832E-6, 6.23974539184983293730E-9, }; static DBL Q2[8] = { /* 1.00000000000000000000E0,*/ 6.02427039364742014255E0, 3.67983563856160859403E0, 1.37702099489081330271E0, 2.16236993594496635890E-1, 1.34204006088543189037E-2, 3.28014464682127739104E-4, 2.89247864745380683936E-6, 6.79019408009981274425E-9, }; /***************************************************************************** * Static functions ******************************************************************************/ static DBL igami (DBL a, DBL y0); static DBL lgam (DBL x); static DBL polevl (DBL x, DBL * coef, int N); static DBL p1evl (DBL x, DBL * coef, int N); static DBL igamc (DBL a, DBL x); static DBL igam (DBL a, DBL x); static DBL ndtri (DBL y0); /* chdtri() * * Inverse of complemented Chi-square distribution * * * * SYNOPSIS: * * DBL df, x, y, chdtri(); * * x = chdtri( df, y ); * * * * * DESCRIPTION: * * Finds the Chi-square argument x such that the integral * from x to infinity of the Chi-square density is equal * to the given cumulative probability y. * * This is accomplished using the inverse gamma integral * function and the relation * * x/2 = igami( df/2, y ); * * * * * ACCURACY: * * See igami.c. * * ERROR MESSAGES: * * message condition value returned * chdtri domain y < 0 or y > 1 0.0 * v < 1 * */ DBL chdtri(DBL df, DBL y) { DBL x; if ((y < 0.0) || (y > 1.0) || (df < 1.0)) { Error("Illegal values fd=%f and y=%f in chdtri().\n", df, y); return (0.0); } x = igami(0.5 * df, y); return (2.0 * x); } /* lgam() * * Natural logarithm of gamma function * * * * SYNOPSIS: * * DBL x, y, lgam(); * extern int sgngam; * * y = lgam( x ); * * * * DESCRIPTION: * * Returns the base e (2.718...) logarithm of the absolute * value of the gamma function of the argument. * The sign (+1 or -1) of the gamma function is returned in a * global (extern) variable named sgngam. * * For arguments greater than 13, the logarithm of the gamma * function is approximated by the logarithmic version of * Stirling's formula using a polynomial approximation of * degree 4. Arguments between -33 and +33 are reduced by * recurrence to the interval [2,3] of a rational approximation. * The cosecant reflection formula is employed for arguments * less than -33. * * Arguments greater than MAXLGM return MAXNUM and an error * message. MAXLGM = 2.035093e36 for DEC * arithmetic or 2.556348e305 for IEEE arithmetic. * * * * ACCURACY: * * * arithmetic domain # trials peak rms * DEC 0, 3 7000 5.2e-17 1.3e-17 * DEC 2.718, 2.035e36 5000 3.9e-17 9.9e-18 * IEEE 0, 3 28000 5.4e-16 1.1e-16 * IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17 * The error criterion was relative when the function magnitude * was greater than one but absolute when it was less than one. * * The following test used the relative error criterion, though * at certain points the relative error could be much higher than * indicated. * IEEE -200, -4 10000 4.8e-16 1.3e-16 * */ static DBL lgam(DBL x) { DBL p, q, w, z; int i; sgngam = 1; if (x < -34.0) { q = -x; w = lgam(q); /* note this modifies sgngam! */ p = floor(q); if (p == q) { goto loverf; } i = p; if ((i & 1) == 0) { sgngam = -1; } else { sgngam = 1; } z = q - p; if (z > 0.5) { p += 1.0; z = p - q; } z = q * sin(M_PI * z); if (z == 0.0) { goto loverf; } /* z = log(M_PI) - log( z ) - w;*/ z = LOGPI - log(z) - w; return (z); } if (x < 13.0) { z = 1.0; while (x >= 3.0) { x -= 1.0; z *= x; } while (x < 2.0) { if (x == 0.0) { goto loverf; } z /= x; x += 1.0; } if (z < 0.0) { sgngam = -1; z = -z; } else { sgngam = 1; } if (x == 2.0) { return (log(z)); } x -= 2.0; p = x * polevl(x, B, 5) / p1evl(x, C, 6); return (log(z) + p); } if (x > MAXLGM) { loverf: /* mtherr("lgam", OVERFLOW); */ return (sgngam * MAXNUM); } q = (x - 0.5) * log(x) - x + LS2PI; if (x > 1.0e8) { return (q); } p = 1.0 / (x * x); if (x >= 1000.0) { q += ((7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3) * p + 0.0833333333333333333333) / x; } else { q += polevl(p, A, 4) / x; } return (q); } /* igamc() * * Complemented incomplete gamma integral * * * * SYNOPSIS: * * DBL a, x, y, igamc(); * * y = igamc( a, x ); * * * * DESCRIPTION: * * The function is defined by * * * igamc(a,x) = 1 - igam(a,x) * * inf. * - * 1 | | -t a-1 * = ----- | e t dt. * - | | * | (a) - * x * * * In this implementation both arguments must be positive. * The integral is evaluated by either a power series or * continued fraction expansion, depending on the relative * values of a and x. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * DEC 0,30 2000 2.7e-15 4.0e-16 * IEEE 0,30 60000 1.4e-12 6.3e-15 * */ static DBL igamc(DBL a, DBL x) { DBL ans, c, yc, ax, y, z; DBL pk, pkm1, pkm2, qk, qkm1, qkm2; DBL r, t; if ((x <= 0) || (a <= 0)) { return (1.0); } if ((x < 1.0) || (x < a)) { return (1.0 - igam(a, x)); } ax = a * log(x) - x - lgam(a); if (ax < -MAXLOG) { /* mtherr("igamc", UNDERFLOW); */ return (0.0); } ax = exp(ax); /* continued fraction */ y = 1.0 - a; z = x + y + 1.0; c = 0.0; pkm2 = 1.0; qkm2 = x; pkm1 = x + 1.0; qkm1 = z * x; ans = pkm1 / qkm1; do { c += 1.0; y += 1.0; z += 2.0; yc = y * c; pk = pkm1 * z - pkm2 * yc; qk = qkm1 * z - qkm2 * yc; if (qk != 0) { r = pk / qk; t = fabs((ans - r) / r); ans = r; } else { t = 1.0; } pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; if (fabs(pk) > BIG) { pkm2 /= BIG; pkm1 /= BIG; qkm2 /= BIG; qkm1 /= BIG; } } while (t > MACHEP); return (ans * ax); } /* igam.c * * Incomplete gamma integral * * * * SYNOPSIS: * * DBL a, x, y, igam(); * * y = igam( a, x ); * * * * DESCRIPTION: * * The function is defined by * * x * - * 1 | | -t a-1 * igam(a,x) = ----- | e t dt. * - | | * | (a) - * 0 * * * In this implementation both arguments must be positive. * The integral is evaluated by either a power series or * continued fraction expansion, depending on the relative * values of a and x. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * DEC 0,30 4000 4.4e-15 6.3e-16 * IEEE 0,30 10000 3.6e-14 5.1e-15 * */ /* left tail of incomplete gamma function: * * inf. k * a -x - x * x e > ---------- * - - * k=0 | (a+k+1) * */ static DBL igam(DBL a, DBL x) { DBL ans, ax, c, r; if ((x <= 0) || (a <= 0)) { return (0.0); } if ((x > 1.0) && (x > a)) { return (1.0 - igamc(a, x)); } /* Compute x**a * exp(-x) / gamma(a) */ ax = a * log(x) - x - lgam(a); if (ax < -MAXLOG) { /* mtherr("igam", UNDERFLOW); */ return (0.0); } ax = exp(ax); /* power series */ r = a; c = 1.0; ans = 1.0; do { r += 1.0; c *= x / r; ans += c; } while (c / ans > MACHEP); return (ans * ax / a); } /* igami() * * Inverse of complemented imcomplete gamma integral * * * * SYNOPSIS: * * DBL a, x, y, igami(); * * x = igami( a, y ); * * * * DESCRIPTION: * * Given y, the function finds x such that * * igamc( a, x ) = y. * * Starting with the approximate value * * 3 * x = a t * * where * * t = 1 - d - ndtri(y) sqrt(d) * * and * * d = 1/9a, * * the routine performs up to 10 Newton iterations to find the * root of igamc(a,x) - y = 0. * * * ACCURACY: * * Tested for a ranging from 0.5 to 30 and x from 0 to 0.5. * * Relative error: * arithmetic domain # trials peak rms * DEC 0,0.5 3400 8.8e-16 1.3e-16 * IEEE 0,0.5 10000 1.1e-14 1.0e-15 * */ static DBL igami(DBL a, DBL y0) { DBL d, y, x0, lgm; int i; /* approximation to inverse function */ d = 1.0 / (9.0 * a); y = (1.0 - d - ndtri(y0) * sqrt(d)); x0 = a * y * y * y; lgm = lgam(a); for (i = 0; i < 10; i++) { if (x0 <= 0.0) { /* mtherr("igami", UNDERFLOW); */ return (0.0); } y = igamc(a, x0); /* compute the derivative of the function at this point */ d = (a - 1.0) * log(x0) - x0 - lgm; if (d < -MAXLOG) { /* mtherr("igami", UNDERFLOW); */ goto done; } d = -exp(d); /* compute the step to the next approximation of x */ if (d == 0.0) { goto done; } d = (y - y0) / d; x0 = x0 - d; if (i < 3) { continue; } if (fabs(d / x0) < 2.0 * MACHEP) { goto done; } } done: return (x0); } /* ndtri.c * * Inverse of Normal distribution function * * * * SYNOPSIS: * * DBL x, y, ndtri(); * * x = ndtri( y ); * * * * DESCRIPTION: * * Returns the argument, x, for which the area under the * Gaussian probability density function (integrated from * minus infinity to x) is equal to y. * * * For small arguments 0 < y < exp(-2), the program computes * z = sqrt( -2.0 * log(y) ); then the approximation is * x = z - log(z)/z - (1/z) P(1/z) / Q(1/z). * There are two rational functions P/Q, one for 0 < y < exp(-32) * and the other for y up to exp(-2). For larger arguments, * w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)). * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * DEC 0.125, 1 5500 9.5e-17 2.1e-17 * DEC 6e-39, 0.135 3500 5.7e-17 1.3e-17 * IEEE 0.125, 1 20000 7.2e-16 1.3e-16 * IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17 * * * ERROR MESSAGES: * * message condition value returned * ndtri domain x <= 0 -MAXNUM * ndtri domain x >= 1 MAXNUM * */ static DBL ndtri(DBL y0) { DBL x, y, z, y2, x0, x1; int code; if (y0 <= 0.0) { /* mtherr("ndtri", DOMAIN); */ return (-MAXNUM); } if (y0 >= 1.0) { /* mtherr("ndtri", DOMAIN); */ return (MAXNUM); } code = 1; y = y0; if (y > (1.0 - 0.13533528323661269189)) /* 0.135... = exp(-2) */ { y = 1.0 - y; code = 0; } if (y > 0.13533528323661269189) { y = y - 0.5; y2 = y * y; x = y + y * (y2 * polevl(y2, P0, 4) / p1evl(y2, Q0, 8)); x = x * s2pi; return (x); } x = sqrt(-2.0 * log(y)); x0 = x - log(x) / x; z = 1.0 / x; if (x < 8.0) /* y > exp(-32) = 1.2664165549e-14 */ { x1 = z * polevl(z, P1, 8) / p1evl(z, Q1, 8); } else { x1 = z * polevl(z, P2, 8) / p1evl(z, Q2, 8); } x = x0 - x1; if (code != 0) { x = -x; } return (x); } /* polevl.c * p1evl.c * * Evaluate polynomial * * * * SYNOPSIS: * * int N; * DBL x, y, coef[N+1], polevl[]; * * y = polevl( x, coef, N ); * * * * DESCRIPTION: * * Evaluates polynomial of degree N: * * 2 N * y = C + C x + C x +...+ C x * 0 1 2 N * * Coefficients are stored in reverse order: * * coef[0] = C , ..., coef[N] = C . * N 0 * * The function p1evl() assumes that coef[N] = 1.0 and is * omitted from the array. Its calling arguments are * otherwise the same as polevl(). * * * SPEED: * * In the interest of speed, there are no checks for out * of bounds arithmetic. This routine is used by most of * the functions in the library. Depending on available * equipment features, the user may wish to rewrite the * program in microcode or assembly language. * */ static DBL polevl(DBL x, DBL coef[], int N) { DBL ans; int i; DBL *p; p = coef; ans = *p++; i = N; do { ans = ans * x + *p++; } while (--i); return (ans); } /* p1evl() */ /* N * Evaluate polynomial when coefficient of x is 1.0. * Otherwise same as polevl. */ static DBL p1evl(DBL x, DBL coef[], int N) { DBL ans; DBL *p; int i; p = coef; ans = x + *p++; i = N - 1; do { ans = ans * x + *p++; } while (--i); return (ans); }