jepler/povray
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity
povray/source/backend/math/chi2.cpp
chris20 86f4639d37 Initial import
2013-11-06 13:07:19 -05:00

1009 lines
18 KiB
C++
Raw Blame History

/*******************************************************************************
* chi2.cpp
*
* This module contains the function for the chi square distribution.
*
* This software is derived from the Cephes Math Library and is
* incorporated herein by permission of the author, Stephen L. Moshier.
*
* Cephes Math Library Release 2.0: April, 1987
* Copyright 1984, 1987 by Stephen L. Moshier
*
* The author reserves the right to distribute this material elsewhere under
* different terms. This file is provided within POV-Ray under the following
* terms:
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
* 3. Neither the names of the copyright holders nor the names of contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
* ---------------------------------------------------------------------------
* $File: //depot/public/povray/3.x/source/backend/math/chi2.cpp $
* $Revision: #1 $
* $Change: 6069 $
* $DateTime: 2013/11/06 11:59:40 $
* $Author: chrisc $
*******************************************************************************/
// frame.h must always be the first POV file included (pulls in platform config)
#include "backend/frame.h"
#include "backend/math/chi2.h"
#include "base/pov_err.h"
// this must be the last file included
#include "base/povdebug.h"
namespace pov
{
/*
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
******************************************************************************/
const DBL MAXLGM = 2.556348e305;
const DBL BIG = 1.44115188075855872E+17;
const DBL MACHEP = 1.38777878078144567553E-17; /* 2**-56 */
const DBL MAXLOG = 8.8029691931113054295988E1; /* log(2**127) */
const DBL MAXNUM = 1.701411834604692317316873e38; /* 2**127 */
const DBL LOGPI = 1.14472988584940017414;
/*
* A[]: Stirling's formula expansion of log gamma
* B[], C[]: log gamma function between 2 and 3
*/
const DBL A[] =
{
8.11614167470508450300E-4,
-5.95061904284301438324E-4,
7.93650340457716943945E-4,
-2.77777777730099687205E-3,
8.33333333333331927722E-2
};
const DBL B[] =
{
-1.37825152569120859100E3,
-3.88016315134637840924E4,
-3.31612992738871184744E5,
-1.16237097492762307383E6,
-1.72173700820839662146E6,
-8.53555664245765465627E5
};
const DBL C[] =
{
1.00000000000000000000E0,
-3.51815701436523470549E2,
-1.70642106651881159223E4,
-2.20528590553854454839E5,
-1.13933444367982507207E6,
-2.53252307177582951285E6,
-2.01889141433532773231E6
};
/* log(sqrt(2pi)) */
const DBL LS2PI = 0.91893853320467274178;
/* sqrt(2pi) */
const DBL s2pi = 2.50662827463100050242E0;
/* approximation for 0 <= |y - 0.5| <= 3/8 */
const DBL P0[5] =
{
-5.99633501014107895267E1,
9.80010754185999661536E1,
-5.66762857469070293439E1,
1.39312609387279679503E1,
-1.23916583867381258016E0,
};
const 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.
*/
const DBL P1[9] =
{
4.05544892305962419923E0,
3.15251094599893866154E1,
5.71628192246421288162E1,
4.40805073893200834700E1,
1.46849561928858024014E1,
2.18663306850790267539E0,
-1.40256079171354495875E-1,
-3.50424626827848203418E-2,
-8.57456785154685413611E-4,
};
const 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.
*/
const 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,
};
const 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, int *sgngam);
static DBL polevl (DBL x, const DBL * coef, int N);
static DBL p1evl (DBL x, const 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))
{
throw POV_EXCEPTION_STRING("Illegal values in chdtri().");
return (0.0);
}
x = igami(0.5 * df, y);
return (2.0 * x);
}
/* lgam()
*
* Natural logarithm of gamma function
*
*
*
* SYNOPSIS:
*
* DBL x, y, lgam();
*
* 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
* 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, int *sgngam)
{
DBL p, q, w, z;
int i;
*sgngam = 1;
if (x < -34.0)
{
q = -x;
w = lgam(q, sgngam); /* 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;
int sgngam = 0;
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, &sgngam);
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;
int sgngam = 0;
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, &sgngam);
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;
int sgngam = 0;
/* 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, &sgngam);
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, const DBL coef[], int N)
{
DBL ans;
int i;
DBL const *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, const DBL coef[], int N)
{
DBL ans;
DBL const *p;
int i;
p = coef;
ans = x + *p++;
i = N - 1;
do
{
ans = ans * x + *p++;
}
while (--i);
return (ans);
}
}
Reference in a new issue View git blame Copy permalink