Index: trunk/npstat/nm/SpecialFunctions.hh =================================================================== --- trunk/npstat/nm/SpecialFunctions.hh (revision 741) +++ trunk/npstat/nm/SpecialFunctions.hh (revision 742) @@ -1,56 +1,59 @@ #ifndef NPSTAT_SPECIALFUNCTIONS_HH_ #define NPSTAT_SPECIALFUNCTIONS_HH_ /*! // \file SpecialFunctions.hh // // \brief Mathematical special functions // // A number of mathematical special functions needed by this package. // The implementations are not optimized in any systematic way, and // may be slow. // // Author: I. Volobouev // // November 2009 */ namespace npstat { /** Inverse cumulative distribition function for 1-d Gaussian */ double inverseGaussCdf(double cdf); /** Regularized incomplete beta function */ double incompleteBeta(double a, double b, double x); /** Inverse regularized incomplete beta function */ double inverseIncompleteBeta(double a, double b, double x); /** The gamma function for positive real arguments */ double Gamma(double x); /** Incomplete gamma ratio */ double incompleteGamma(double a, double x); /** Incomplete gamma ratio complement */ double incompleteGammaC(double a, double x); /** Inverse incomplete gamma ratio */ double inverseIncompleteGamma(double a, double x); /** Inverse incomplete gamma ratio complement */ double inverseIncompleteGammaC(double a, double x); /** Dawson's integral exp(-x^2) Int_0^x exp(t^2) dt */ long double dawsonIntegral(long double x); /** Inverse of the integral Int_0^x exp(t^2) dt */ long double inverseExpsqIntegral(long double x); /** Order n derivative of Gaussian density with mean 0 and sigma 1 */ long double normalDensityDerivative(unsigned n, long double x); /** Bivariate normal cumulative probability */ double bivariateNormalIntegral(double rho, double x, double y); + + /** Modified Bessel function of the second kind */ + double besselK(double nu, double z); } #endif // NPSTAT_SPECIALFUNCTIONS_HH_ Index: trunk/npstat/nm/SpecialFunctions.cc =================================================================== --- trunk/npstat/nm/SpecialFunctions.cc (revision 741) +++ trunk/npstat/nm/SpecialFunctions.cc (revision 742) @@ -1,1977 +1,2131 @@ // Most of the code below is cherry-picked from the Cephes library // written by Stephen L. Moshier #include #include #include #include #include #include #include "npstat/nm/SpecialFunctions.hh" #include "npstat/nm/MathUtils.hh" #include "npstat/nm/Matrix.hh" #include "npstat/rng/permutation.hh" #define SQRTTWOPIL 2.506628274631000502415765L #define SQTPI 2.50662827463100050242 #define LOGPI 1.14472988584940017414 #define MAXGAM 171.624376956302725 #define MAXSTIR 143.01608 #define MAXLOG 709.782712893383996732224 #define MINLOG (-708.396418532264106224) #define MACHEP 1.2e-16 #define MAXLGM 2.556348e305 #define MAXNUM 1.79769313486231570815e308 #define LS2PI 0.91893853320467274178 #define GAUSS_MAX_SIGMA 37.615 static double igam(double a, double x); static double invgauss(double x); static double polevl(const double x, const double *p, int i) { double ans = *p++; do { ans = ans * x + *p++; } while (--i); return ans; } static double p1evl(double x, const double *p, const int N) { double ans = x + *p++; int i = N - 1; do { ans = ans * x + *p++; } while (--i); return ans; } /* Gamma function computed by Stirling's formula. * The polynomial STIR is valid for 33 <= x <= 172. */ static double stirf(const double x) { static const double STIR[5] = { 7.87311395793093628397E-4, -2.29549961613378126380E-4, -2.68132617805781232825E-3, 3.47222221605458667310E-3, 8.33333333333482257126E-2, }; double y, w; w = 1.0/x; w = 1.0 + w * polevl( w, STIR, 4 ); y = exp(x); if( x > MAXSTIR ) { // Avoid overflow in pow() const double v = pow( x, 0.5 * x - 0.25 ); y = v * (v / y); } else { y = pow( x, x - 0.5 ) / y; } return SQTPI * y * w; } // Logarithm of the gamma function static double lgam(double x, int* sgngam_out=0) { int sgngam = 1; if (sgngam_out) *sgngam_out = sgngam; double p, q, u, w, z; int i; static const double A[] = { 8.11614167470508450300E-4, -5.95061904284301438324E-4, 7.93650340457716943945E-4, -2.77777777730099687205E-3, 8.33333333333331927722E-2 }; static const double B[] = { -1.37825152569120859100E3, -3.88016315134637840924E4, -3.31612992738871184744E5, -1.16237097492762307383E6, -1.72173700820839662146E6, -8.53555664245765465627E5 }; static const double C[] = { -3.51815701436523470549E2, -1.70642106651881159223E4, -2.20528590553854454839E5, -1.13933444367982507207E6, -2.53252307177582951285E6, -2.01889141433532773231E6 }; if( x < -34.0 ) { q = -x; w = lgam(q, &sgngam); if (sgngam_out) *sgngam_out = sgngam; p = floor(q); if( p == q ) goto loverf; i = static_cast(p); if( (i & 1) == 0 ) sgngam = -1; else sgngam = 1; if (sgngam_out) *sgngam_out = sgngam; 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(PI) - log( z ) - w; z = LOGPI - log( z ) - w; return( z ); } if( x < 13.0 ) { z = 1.0; p = 0.0; u = x; while( u >= 3.0 ) { p -= 1.0; u = x + p; z *= u; } while( u < 2.0 ) { if( u == 0.0 ) goto loverf; z /= u; p += 1.0; u = x + p; } if( z < 0.0 ) { sgngam = -1; z = -z; } else sgngam = 1; if (sgngam_out) *sgngam_out = sgngam; if( u == 2.0 ) return( log(z) ); p -= 2.0; x = x + p; p = x * polevl( x, B, 5 ) / p1evl( x, C, 6); return( log(z) + p ); } if( x > MAXLGM ) { loverf: assert(!"Overflow in lgam"); 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 ); } // Complementary incomplete gamma static double igamc(double a, double x) { const double big = 4.503599627370496e15; const double biginv = 2.22044604925031308085e-16; double ans, ax, c, yc, r, t, y, z; double pk, pkm1, pkm2, qk, qkm1, qkm2; if( (x <= 0) || ( a <= 0) ) return( 1.0 ); if( (x < 1.0) || (x < a) ) return( 1.0 - igam(a,x) ); // Compute x**a * exp(-x) / gamma(a) ax = a * log(x) - x - lgam(a); if( ax < -MAXLOG ) { // mtherr( "igamc", UNDERFLOW ); // assert(!"Underflow in igamc"); 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 *= biginv; pkm1 *= biginv; qkm2 *= biginv; qkm1 *= biginv; } } while( t > MACHEP ); return( ans * ax ); } // Incomplete gamma static double igam(double a, double x) { double 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 ); assert(!"Underflow in igam"); 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 ); } // Inverse of complemented incomplete gamma integral. // Work only for y0 <= 0.5. static double igami(double a, double y0) { double x0, x1, x, yl, yh, y, d, lgm, dithresh; int i, dir; assert(y0 <= 0.5); /* bound the solution */ x0 = MAXNUM; yl = 0; x1 = 0; yh = 1.0; dithresh = 5.0 * MACHEP; /* approximation to inverse function */ d = 1.0/(9.0*a); y = ( 1.0 - d - invgauss(y0) * sqrt(d) ); x = a * y * y * y; lgm = lgam(a); for( i=0; i<10; i++ ) { if( x > x0 || x < x1 ) goto ihalve; y = igamc(a,x); if( y < yl || y > yh ) goto ihalve; if( y < y0 ) { x0 = x; yl = y; } else { x1 = x; yh = y; } /* compute the derivative of the function at this point */ d = (a - 1.0) * log(x) - x - lgm; if( d < -MAXLOG ) goto ihalve; d = -exp(d); /* compute the step to the next approximation of x */ d = (y - y0)/d; if( fabs(d/x) < MACHEP ) goto done; x = x - d; } /* Resort to interval halving if Newton iteration did not converge. */ ihalve: d = 0.0625; if( x0 == MAXNUM ) { if( x <= 0.0 ) x = 1.0; while( x0 == MAXNUM ) { x = (1.0 + d) * x; y = igamc( a, x ); if( y < y0 ) { x0 = x; yl = y; break; } d = d + d; } } d = 0.5; dir = 0; for( i=0; i<400; i++ ) { x = x1 + d * (x0 - x1); y = igamc( a, x ); lgm = (x0 - x1)/(x1 + x0); if( fabs(lgm) < dithresh ) break; lgm = (y - y0)/y0; if( fabs(lgm) < dithresh ) break; if( x <= 0.0 ) break; if( y >= y0 ) { x1 = x; yh = y; if( dir < 0 ) { dir = 0; d = 0.5; } else if( dir > 1 ) d = 0.5 * d + 0.5; else d = (y0 - yl)/(yh - yl); dir += 1; } else { x0 = x; yl = y; if( dir > 0 ) { dir = 0; d = 0.5; } else if( dir < -1 ) d = 0.5 * d; else d = (y0 - yl)/(yh - yl); dir -= 1; } } if( x == 0.0 ) // mtherr( "igami", UNDERFLOW ); assert(!"Underflow in igami"); done: return( x ); } static double invgauss(const double P) { assert(P > 0.0); assert(P < 1.0); // Translated from PPND16 algorithm of Wichura (originally in Fortran). // This was not taken from Cephes. static const double ZERO = 0., ONE = 1., HALF = 0.5, SPLIT1 = 0.425, SPLIT2 = 5., CONST1 = 0.180625, CONST2 = 1.6; static const double A0 = 3.3871328727963666080, A1 = 1.3314166789178437745E+2, A2 = 1.9715909503065514427E+3, A3 = 1.3731693765509461125E+4, A4 = 4.5921953931549871457E+4, A5 = 6.7265770927008700853E+4, A6 = 3.3430575583588128105E+4, A7 = 2.5090809287301226727E+3, B1 = 4.2313330701600911252E+1, B2 = 6.8718700749205790830E+2, B3 = 5.3941960214247511077E+3, B4 = 2.1213794301586595867E+4, B5 = 3.9307895800092710610E+4, B6 = 2.8729085735721942674E+4, B7 = 5.2264952788528545610E+3; static const double C0 = 1.42343711074968357734, C1 = 4.63033784615654529590, C2 = 5.76949722146069140550, C3 = 3.64784832476320460504, C4 = 1.27045825245236838258, C5 = 2.41780725177450611770E-1, C6 = 2.27238449892691845833E-2, C7 = 7.74545014278341407640E-4, D1 = 2.05319162663775882187, D2 = 1.67638483018380384940, D3 = 6.89767334985100004550E-1, D4 = 1.48103976427480074590E-1, D5 = 1.51986665636164571966E-2, D6 = 5.47593808499534494600E-4, D7 = 1.05075007164441684324E-9; static const double E0 = 6.65790464350110377720, E1 = 5.46378491116411436990, E2 = 1.78482653991729133580, E3 = 2.96560571828504891230E-1, E4 = 2.65321895265761230930E-2, E5 = 1.24266094738807843860E-3, E6 = 2.71155556874348757815E-5, E7 = 2.01033439929228813265E-7, F1 = 5.99832206555887937690E-1, F2 = 1.36929880922735805310E-1, F3 = 1.48753612908506148525E-2, F4 = 7.86869131145613259100E-4, F5 = 1.84631831751005468180E-5, F6 = 1.42151175831644588870E-7, F7 = 2.04426310338993978564E-15; const double Q = P - HALF; double R, PPND16; if (fabs(Q) <= SPLIT1) { R = CONST1 - Q * Q; PPND16 = Q * (((((((A7 * R + A6) * R + A5) * R + A4) * R + A3) * R + A2) * R + A1) * R + A0) / (((((((B7 * R + B6) * R + B5) * R + B4) * R + B3) * R + B2) * R + B1) * R + ONE); } else { if (Q < ZERO) R = P; else R = ONE - P; R = sqrt(-log(R)); if (R <= SPLIT2) { R = R - CONST2; PPND16 = (((((((C7 * R + C6) * R + C5) * R + C4) * R + C3) * R + C2) * R + C1) * R + C0) / (((((((D7 * R + D6) * R + D5) * R + D4) * R + D3) * R + D2) * R + D1) * R + ONE); } else { R = R - SPLIT2 ; PPND16 = (((((((E7 * R + E6) * R + E5) * R + E4) * R + E3) * R + E2) * R + E1) * R + E0) / (((((((F7 * R + F6) * R + F5) * R + F4) * R + F3) * R + F2) * R + F1) * R + ONE); } if (Q < ZERO) PPND16 = -PPND16; } return PPND16; } /* Continued fraction expansion #1 * for incomplete beta integral */ static double incbcf( double a, double b, double x ) { static const double big = 4.503599627370496e15; static const double biginv = 2.22044604925031308085e-16; double xk, pk, pkm1, pkm2, qk, qkm1, qkm2; double k1, k2, k3, k4, k5, k6, k7, k8; double r, t, ans, thresh; int n; k1 = a; k2 = a + b; k3 = a; k4 = a + 1.0; k5 = 1.0; k6 = b - 1.0; k7 = k4; k8 = a + 2.0; pkm2 = 0.0; qkm2 = 1.0; pkm1 = 1.0; qkm1 = 1.0; ans = 1.0; r = 1.0; n = 0; thresh = 3.0 * MACHEP; do { xk = -( x * k1 * k2 )/( k3 * k4 ); pk = pkm1 + pkm2 * xk; qk = qkm1 + qkm2 * xk; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; xk = ( x * k5 * k6 )/( k7 * k8 ); pk = pkm1 + pkm2 * xk; qk = qkm1 + qkm2 * xk; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; if( qk != 0 ) r = pk/qk; if( r != 0 ) { t = fabs( (ans - r)/r ); ans = r; } else t = 1.0; if( t < thresh ) goto cdone; k1 += 1.0; k2 += 1.0; k3 += 2.0; k4 += 2.0; k5 += 1.0; k6 -= 1.0; k7 += 2.0; k8 += 2.0; if( (fabs(qk) + fabs(pk)) > big ) { pkm2 *= biginv; pkm1 *= biginv; qkm2 *= biginv; qkm1 *= biginv; } if( (fabs(qk) < biginv) || (fabs(pk) < biginv) ) { pkm2 *= big; pkm1 *= big; qkm2 *= big; qkm1 *= big; } } while( ++n < 300 ); cdone: return(ans); } /* Continued fraction expansion #2 * for incomplete beta integral */ static double incbd( double a, double b, double x ) { static const double big = 4.503599627370496e15; static const double biginv = 2.22044604925031308085e-16; double xk, pk, pkm1, pkm2, qk, qkm1, qkm2; double k1, k2, k3, k4, k5, k6, k7, k8; double r, t, ans, z, thresh; int n; k1 = a; k2 = b - 1.0; k3 = a; k4 = a + 1.0; k5 = 1.0; k6 = a + b; k7 = a + 1.0;; k8 = a + 2.0; pkm2 = 0.0; qkm2 = 1.0; pkm1 = 1.0; qkm1 = 1.0; z = x / (1.0-x); ans = 1.0; r = 1.0; n = 0; thresh = 3.0 * MACHEP; do { xk = -( z * k1 * k2 )/( k3 * k4 ); pk = pkm1 + pkm2 * xk; qk = qkm1 + qkm2 * xk; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; xk = ( z * k5 * k6 )/( k7 * k8 ); pk = pkm1 + pkm2 * xk; qk = qkm1 + qkm2 * xk; pkm2 = pkm1; pkm1 = pk; qkm2 = qkm1; qkm1 = qk; if( qk != 0 ) r = pk/qk; if( r != 0 ) { t = fabs( (ans - r)/r ); ans = r; } else t = 1.0; if( t < thresh ) goto cdone; k1 += 1.0; k2 -= 1.0; k3 += 2.0; k4 += 2.0; k5 += 1.0; k6 += 1.0; k7 += 2.0; k8 += 2.0; if( (fabs(qk) + fabs(pk)) > big ) { pkm2 *= biginv; pkm1 *= biginv; qkm2 *= biginv; qkm1 *= biginv; } if( (fabs(qk) < biginv) || (fabs(pk) < biginv) ) { pkm2 *= big; pkm1 *= big; qkm2 *= big; qkm1 *= big; } } while( ++n < 300 ); cdone: return(ans); } /* Power series for incomplete beta integral. Use when b*x is small and x not too close to 1. */ static double pseries( double a, double b, double x ) { double s, t, u, v, n, t1, z, ai; ai = 1.0 / a; u = (1.0 - b) * x; v = u / (a + 1.0); t1 = v; t = u; n = 2.0; s = 0.0; z = MACHEP * ai; while( fabs(v) > z ) { u = (n - b) * x / n; t *= u; v = t / (a + n); s += v; n += 1.0; } s += t1; s += ai; u = a * log(x); if( (a+b) < MAXGAM && fabs(u) < MAXLOG ) { t = npstat::Gamma(a+b)/(npstat::Gamma(a)*npstat::Gamma(b)); s = s * t * pow(x,a); } else { t = lgam(a+b) - lgam(a) - lgam(b) + u + log(s); if( t < MINLOG ) s = 0.0; else s = exp(t); } return(s); } /* Incomplete beta function */ static double incbet( double aa, double bb, double xx ) { double a, b, t, x, xc, w, y; int flag; if( aa <= 0.0 || bb <= 0.0 ) goto domerr; if( (xx <= 0.0) || ( xx >= 1.0) ) { if( xx == 0.0 ) return(0.0); if( xx == 1.0 ) return( 1.0 ); domerr: // mtherr( "incbet", DOMAIN ); assert(!"Domain error in incbet"); return( 0.0 ); } flag = 0; if( (bb * xx) <= 1.0 && xx <= 0.95) { t = pseries(aa, bb, xx); goto done; } w = 1.0 - xx; /* Reverse a and b if x is greater than the mean. */ if( xx > (aa/(aa+bb)) ) { flag = 1; a = bb; b = aa; xc = xx; x = w; } else { a = aa; b = bb; xc = w; x = xx; } if( flag == 1 && (b * x) <= 1.0 && x <= 0.95) { t = pseries(a, b, x); goto done; } /* Choose expansion for better convergence. */ y = x * (a+b-2.0) - (a-1.0); if( y < 0.0 ) w = incbcf( a, b, x ); else w = incbd( a, b, x ) / xc; /* Multiply w by the factor a b _ _ _ x (1-x) | (a+b) / ( a | (a) | (b) ) . */ y = a * log(x); t = b * log(xc); if( (a+b) < MAXGAM && fabs(y) < MAXLOG && fabs(t) < MAXLOG ) { t = pow(xc,b); t *= pow(x,a); t /= a; t *= w; t *= npstat::Gamma(a+b) / (npstat::Gamma(a) * npstat::Gamma(b)); goto done; } /* Resort to logarithms. */ y += t + lgam(a+b) - lgam(a) - lgam(b); y += log(w/a); if( y < MINLOG ) t = 0.0; else t = exp(y); done: if( flag == 1 ) { if( t <= MACHEP ) t = 1.0 - MACHEP; else t = 1.0 - t; } return( t ); } /* Inverse incomplete beta function */ static double incbi( double aa, double bb, double yy0 ) { double a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt; int i, rflg, dir, nflg; i = 0; if( yy0 <= 0.0 ) return(0.0); if( yy0 >= 1.0 ) return(1.0); x0 = 0.0; yl = 0.0; x1 = 1.0; yh = 1.0; nflg = 0; if( aa <= 1.0 || bb <= 1.0 ) { dithresh = 1.0e-6; rflg = 0; a = aa; b = bb; y0 = yy0; x = a/(a+b); y = incbet( a, b, x ); goto ihalve; } else { dithresh = 1.0e-4; } /* approximation to inverse function */ yp = -npstat::inverseGaussCdf(yy0); if( yy0 > 0.5 ) { rflg = 1; a = bb; b = aa; y0 = 1.0 - yy0; yp = -yp; } else { rflg = 0; a = aa; b = bb; y0 = yy0; } lgm = (yp * yp - 3.0)/6.0; x = 2.0/( 1.0/(2.0*a-1.0) + 1.0/(2.0*b-1.0) ); d = yp * sqrt( x + lgm ) / x - ( 1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0) ) * (lgm + 5.0/6.0 - 2.0/(3.0*x)); d = 2.0 * d; if( d < MINLOG ) { x = 1.0; goto under; } x = a/( a + b * exp(d) ); y = incbet( a, b, x ); yp = (y - y0)/y0; if( fabs(yp) < 0.2 ) goto newt; /* Resort to interval halving if not close enough. */ ihalve: dir = 0; di = 0.5; for( i=0; i<100; i++ ) { if( i != 0 ) { x = x0 + di * (x1 - x0); if( x == 1.0 ) x = 1.0 - MACHEP; if( x == 0.0 ) { di = 0.5; x = x0 + di * (x1 - x0); if( x == 0.0 ) goto under; } y = incbet( a, b, x ); yp = (x1 - x0)/(x1 + x0); if( fabs(yp) < dithresh ) goto newt; yp = (y-y0)/y0; if( fabs(yp) < dithresh ) goto newt; } if( y < y0 ) { x0 = x; yl = y; if( dir < 0 ) { dir = 0; di = 0.5; } else if( dir > 3 ) di = 1.0 - (1.0 - di) * (1.0 - di); else if( dir > 1 ) di = 0.5 * di + 0.5; else di = (y0 - y)/(yh - yl); dir += 1; if( x0 > 0.75 ) { if( rflg == 1 ) { rflg = 0; a = aa; b = bb; y0 = yy0; } else { rflg = 1; a = bb; b = aa; y0 = 1.0 - yy0; } x = 1.0 - x; y = incbet( a, b, x ); x0 = 0.0; yl = 0.0; x1 = 1.0; yh = 1.0; goto ihalve; } } else { x1 = x; if( rflg == 1 && x1 < MACHEP ) { x = 0.0; goto done; } yh = y; if( dir > 0 ) { dir = 0; di = 0.5; } else if( dir < -3 ) di = di * di; else if( dir < -1 ) di = 0.5 * di; else di = (y - y0)/(yh - yl); dir -= 1; } } // Partial loss of precision. Hm... Will ignore for now - igv // mtherr( "incbi", PLOSS ); if( x0 >= 1.0 ) { x = 1.0 - MACHEP; goto done; } if( x <= 0.0 ) { under: // mtherr( "incbi", UNDERFLOW ); assert(!"Underflow in incbi"); x = 0.0; goto done; } newt: if( nflg ) goto done; nflg = 1; lgm = lgam(a+b) - lgam(a) - lgam(b); for( i=0; i<8; i++ ) { /* Compute the function at this point. */ if( i != 0 ) y = incbet(a,b,x); if( y < yl ) { x = x0; y = yl; } else if( y > yh ) { x = x1; y = yh; } else if( y < y0 ) { x0 = x; yl = y; } else { x1 = x; yh = y; } if( x == 1.0 || x == 0.0 ) break; /* Compute the derivative of the function at this point. */ d = (a - 1.0) * log(x) + (b - 1.0) * log(1.0-x) + lgm; if( d < MINLOG ) goto done; if( d > MAXLOG ) break; d = exp(d); /* Compute the step to the next approximation of x. */ d = (y - y0)/d; xt = x - d; if( xt <= x0 ) { y = (x - x0) / (x1 - x0); xt = x0 + 0.5 * y * (x - x0); if( xt <= 0.0 ) break; } if( xt >= x1 ) { y = (x1 - x) / (x1 - x0); xt = x1 - 0.5 * y * (x1 - x); if( xt >= 1.0 ) break; } x = xt; if( fabs(d/x) < 128.0 * MACHEP ) goto done; } /* Did not converge. */ dithresh = 256.0 * MACHEP; goto ihalve; done: if( rflg ) { if( x <= MACHEP ) x = 1.0 - MACHEP; else x = 1.0 - x; } return( x ); } static void buildGaussDeriCoeffs(const unsigned order, std::vector* coeffs) { assert(coeffs); const unsigned dim = order + 1U; coeffs->resize(dim); if (order == 0U) { (*coeffs)[0] = 1LL; return; } npstat::Matrix deriOp(dim, dim, 0); for (unsigned row=0; row(row+1U); for (unsigned row=1; row& result = deriOp.pow(order); for (unsigned row=0; row= 0 // px = Phi( x ) (cumulative std. normal) // pxs = Phi( x * sqrt( ( 1 - rho ) / ( 1 + rho ) ) ) // Phi2diag() should not be used directly but Phi2() instead. static double Phi2diag( const double x, const double a, const double px, const double pxs ) { if( a <= 0.0 ) return px; // rho == 1 if( a >= 1.0 ) return px * px; // rho == 0 double b = 2.0 - a, sqrt_ab = sqrt( a * b ); // b == 1 + rho const double c1 = 6.36619772367581343e-1; // == 2 / PI const double c2 = 1.25331413731550025; // == sqrt( PI / 2 ) const double c3 = 1.57079632679489662; // == PI / 2 const double c4 = 1.591549430918953358e-1; // == 1 / ( 2 * PI ) // avoid asin() for values close to 1 (vertical tangent) double asr = ( a > 0.1 ? asin( 1.0 - a ) : acos( sqrt_ab ) ); // return upper bound if absolute error within double precision double comp = px * pxs; if( comp * ( 1.0 - a - c1 * asr ) < 5e-17 ) return b * comp; // initialize odd and even terms for recursion double tmp = c2 * x; double alpha = a * x * x / b; double a_even = -tmp * a; double a_odd = -sqrt_ab * alpha; double beta = x * x; double b_even = tmp * sqrt_ab; double b_odd = sqrt_ab * beta; double delta = 2.0 * x * x / b; double d_even = ( 1.0 - a ) * c3 - asr; double d_odd = tmp * ( sqrt_ab - a ); // recursion; update odd and even terms in each step double res = 0.0, res_new = d_even + d_odd; int k = 2; while( res != res_new ) { d_even = ( a_odd + b_odd + delta * d_even ) / k; a_even *= alpha / k; b_even *= beta / k; k++; a_odd *= alpha / k; b_odd *= beta / k; d_odd = ( a_even + b_even + delta * d_odd ) / k; k++; res = res_new; res_new += d_even + d_odd; } // transform; check against upper and lower bounds res *= exp( -x * x / b ) * c4; return std::max( ( 1.0 + c1 * asr ) * comp, b * comp - std::max( 0.0, res ) ); } static double squarethis(const double v) { return v * v; } // Auxiliary function Phi2help() will be called (twice, with x and y interchanged) by Phi2() // In particular, certain pathological cases will have been dealt with in Phi2() before // Axis is transformed into diagonal, and Phi2diag() is called static double Phi2help( const double x, const double y, const double rho ) { // note: case x == y == 0.0 is dealt with in Phi2() if( fabs(x) <= DBL_MIN ) return ( y >= 0.0 ? 0.0 : 0.5 ); double s = sqrt( ( 1.0 - rho ) * ( 1.0 + rho ) ); double a = 0.0, b1 = -fabs( x ), b2 = 0.0; // avoid cancellation by treating the cases rho -> +-1 // separately (cutoff at |0.99| might be optimized); if( rho > 0.99 ) { double tmp = sqrt( ( 1.0 - rho ) / ( 1.0 + rho ) ); b2 = -fabs( ( x - y ) / s - x * tmp ); a = squarethis( ( x - y ) / x / s - tmp ); } else if( rho < -0.99 ) { double tmp = sqrt( ( 1.0 + rho ) / ( 1.0 - rho ) ); b2 = -fabs( ( x + y ) / s - x * tmp ); a = squarethis( ( x + y ) / x / s - tmp ); } else { b2 = -fabs( rho * x - y ) / s; a = squarethis( b2 / x ); } // Phi(): cumulative std. normal (has to be defined elsewhere) double p1 = Phi( b1 ), p2 = Phi( b2 ); // reduction to the diagonal double q = 0.0; if( a <= 1.0 ) q = 0.5 * Phi2diag( b1, 2.0 * a / ( 1.0 + a ), p1, p2 ); else q = p1 * p2 - 0.5 * Phi2diag( b2, 2.0 / ( 1.0 + a ), p2, p1 ); // finally, transformation according to conditions on x, y, rho int c1 = ( y / x >= rho ), c2 = ( x < 0.0 ), c3 = c2 && ( y >= 0.0 ); return ( c1 && c3 ? q - 0.5 : c1 && c2 ? q : c1 ? 0.5 - p1 + q : c3 ? p1 - q - 0.5 : c2 ? p1 - q : 0.5 - q ); } // Evaluation of the cumulative bivariate normal distribution // for arbitrary (x,y) and correlation rho. This function is // described in the article "Recursive Numerical Evaluation of // the Cumulative Bivariate Normal Distribution" by Christian Meyer, // Journal of Statistical Software, Vol. 52, Issue 10, Mar 2013. // // It appears not to work very well for large rho (for example, // Phi2(-0, 0.5, -0.98865243049242746) is wrong), so for large // values of rho implementation by Genz is used (which has its // own problems). // static double Phi2( const double x, const double y, const double rho ) { // special case |rho| == 1 => lower or upper Frechet bound // Phi(): cumulative std. normal (has to be defined elsewhere) if( ( 1.0 - rho ) * ( 1.0 + rho ) <= 0.0 ) { if( rho > 0.0 ) return Phi( std::min( x, y ) ); else return std::max( 0.0, std::min( 1.0, Phi( x ) + Phi( y ) - 1.0 ) ); } // special case x == y == 0, avoids complications in Phi2help() if( x == 0.0 && y == 0.0 ) { if( rho > 0.0 ) return Phi2diag( 0.0, 1.0 - rho, 0.5, 0.5 ); else return 0.5 - Phi2diag( 0.0, 1.0 + rho, 0.5, 0.5 ); } // standard case by reduction formula const double result = std::max( 0.0, std::min( 1.0, Phi2help( x, y, rho ) + Phi2help( y, x, rho ) ) ); // std::cout.precision(17); // std::cout << "Phi2(" << x << ", " << y << ", " << rho << ") = " << result << std::endl; return result; } /* * A function for computing bivariate normal probabilities. * * Alan Genz * Department of Mathematics * Washington State University * Pullman, WA 99164-3113 * Email : alangenz@wsu.edu * * This function is based on the method described by * Drezner, Z and G.O. Wesolowsky, (1989), * On the computation of the bivariate normal integral, * Journal of Statist. Comput. Simul. 35, pp. 101-107, * with major modifications for double precision, and for |R| close to 1. * * BVND calculates the probability that X > DH and Y > DK. * Note: Prob( X < DH, Y < DK ) = BVND( -DH, -DK, R ). * * Parameters * * DH DOUBLE PRECISION, integration limit * DK DOUBLE PRECISION, integration limit * R DOUBLE PRECISION, correlation coefficient * * Translated from Fortran by igv (July 2014). */ static double bvnd(const double DH, const double DK, const double R) { const double TWOPI = 2.0*M_PI; // Gauss-Legendre points static const double X[11][4] = { {0., 0. , 0. , 0. }, {0., -0.9324695142031522, -0.9815606342467191, -0.9931285991850949}, {0., -0.6612093864662647, -0.9041172563704750, -0.9639719272779138}, {0., -0.2386191860831970, -0.7699026741943050, -0.9122344282513259}, {0., 0. , -0.5873179542866171, -0.8391169718222188}, {0., 0. , -0.3678314989981802, -0.7463319064601508}, {0., 0. , -0.1252334085114692, -0.6360536807265150}, {0., 0. , 0. , -0.5108670019508271}, {0., 0. , 0. , -0.3737060887154196}, {0., 0. , 0. , -0.2277858511416451}, {0., 0. , 0. , -0.7652652113349733} }; // Gauss-Legendre weights static const double W[11][4] = { {0., 0. , 0. , 0. }, {0., 0.1713244923791705, 0.04717533638651177, 0.01761400713915212}, {0., 0.3607615730481384, 0.1069393259953183 , 0.04060142980038694}, {0., 0.4679139345726904, 0.1600783285433464 , 0.06267204833410906}, {0., 0. , 0.2031674267230659 , 0.08327674157670475}, {0., 0. , 0.2334925365383547 , 0.1019301198172404 }, {0., 0. , 0.2491470458134029 , 0.1181945319615184 }, {0., 0. , 0. , 0.1316886384491766 }, {0., 0. , 0. , 0.1420961093183821 }, {0., 0. , 0. , 0.1491729864726037 }, {0., 0. , 0. , 0.1527533871307259 } }; int I, IS, LG, NG; double AS, A, B, C, D, RS, XS, BVN; double SN, ASR, H, K, BS, HS, HK; if ( fabs(R) < 0.3 ) { NG = 1; LG = 3; } else if ( fabs(R) < 0.75 ) { NG = 2; LG = 6; } else { NG = 3; LG = 10; } H = DH; K = DK; HK = H*K; BVN = 0.0; if ( fabs(R) < 0.925 ) { if ( fabs(R) > 0.0 ) { HS = ( H*H + K*K )/2.0; ASR = asin(R); for (I = 1; I <= LG; ++I) { for (IS = -1; IS <= 1; IS += 2) { SN = sin( ASR*( IS*X[I][NG] + 1.0 )/2.0 ); BVN += W[I][NG]*exp( ( SN*HK-HS )/( 1.0-SN*SN ) ); } } BVN = BVN*ASR/( 2.0*TWOPI ); } BVN += Phi(-H)*Phi(-K); } else { if ( R < 0.0 ) { K = -K; HK = -HK; } if ( fabs(R) < 1.0 ) { AS = ( 1 - R )*( 1 + R ); A = sqrt(AS); BS = pow(H - K, 2); C = ( 4 - HK )/8.0; D = ( 12 - HK )/16.0; ASR = -( BS/AS + HK )/2.0; if ( ASR > -100.0 ) BVN = A*exp(ASR)*( 1 - C*( BS - AS )*( 1 - D*BS/5 )/3 + C*D*AS*AS/5 ); if ( -HK < 100.0 ) { B = sqrt(BS); BVN -= exp( -HK/2 )*SQTPI*Phi(-B/A)*B *( 1 - C*BS*( 1 - D*BS/5 )/3 ); } A /= 2.0; for (I = 1; I <= LG; ++I) { for (IS = -1; IS <= 1; IS += 2) { XS = pow( A*( IS*X[I][NG] + 1 ), 2); RS = sqrt( 1 - XS ); ASR = -( BS/XS + HK )/2; if ( ASR > -100 ) BVN += A*W[I][NG]*exp( ASR ) *( exp( -HK*XS/( 2*pow(1 + RS, 2) ) )/RS - ( 1 + C*XS*( 1 + D*XS ) ) ); } } BVN = -BVN/TWOPI; } if ( R > 0.0 ) BVN += Phi( -std::max( H, K ) ); else { BVN = -BVN; if ( K > H ) { if ( H < 0.0 ) BVN += Phi(K) - Phi(H); else BVN += Phi(-H) - Phi(-K); } } } return BVN; } namespace npstat { double inverseGaussCdf(const double x) { if (!(x >= 0.0 && x <= 1.0)) throw std::domain_error( "In npstat::inverseGaussCdf: argument outside of [0, 1] interval"); if (x == 1.0) return GAUSS_MAX_SIGMA; else if (x == 0.0) return -GAUSS_MAX_SIGMA; else return invgauss(x); } double incompleteBeta(const double a, const double b, const double x) { if (!(x >= 0.0 && x <= 1.0)) throw std::domain_error( "In npstat::incompleteBeta: argument outside of [0, 1] interval"); if (!(a > 0.0 && b > 0.0)) throw std::invalid_argument( "In npstat::incompleteBeta: parameters must be positive"); return incbet(a, b, x); } double inverseIncompleteBeta(const double a, const double b, const double x) { if (!(x >= 0.0 && x <= 1.0)) throw std::domain_error( "In npstat::inverseIncompleteBeta: " "argument outside of [0, 1] interval"); if (!(a > 0.0 && b > 0.0)) throw std::invalid_argument( "In npstat::inverseIncompleteBeta: parameters must be positive"); return incbi(a, b, x); } double Gamma(double x) { if (x <= 0.0) throw std::domain_error( "In npstat::Gamma: argument must be positive"); // The Stirling formula overflows for x >= MAXGAM if (x >= MAXGAM) throw std::overflow_error( "In npstat::Gamma: argument is too large"); const unsigned ix = x; if (ix && x == static_cast(ix)) return ldfactorial(ix - 1U); static const double P[] = { 1.60119522476751861407E-4, 1.19135147006586384913E-3, 1.04213797561761569935E-2, 4.76367800457137231464E-2, 2.07448227648435975150E-1, 4.94214826801497100753E-1, 9.99999999999999996796E-1 }; static const double Q[] = { -2.31581873324120129819E-5, 5.39605580493303397842E-4, -4.45641913851797240494E-3, 1.18139785222060435552E-2, 3.58236398605498653373E-2, -2.34591795718243348568E-1, 7.14304917030273074085E-2, 1.00000000000000000320E0 }; double p, z = 1.0, q = x; if (q > 33.0) return stirf(q); while( x >= 3.0 ) { x -= 1.0; z *= x; } while( x < 2.0 ) { if( x < 1.e-9 ) return( z/((1.0 + 0.5772156649015329 * x) * x) ); z /= x; x += 1.0; } if( x == 2.0 ) return(z); x -= 2.0; p = polevl( x, P, 6 ); q = polevl( x, Q, 7 ); return( z * p / q ); } double incompleteGamma(const double a, const double x) { return igam(a, x); } double incompleteGammaC(const double a, const double x) { return igamc(a, x); } double inverseIncompleteGammaC(const double a, const double x) { if (!(x >= 0.0 && x <= 1.0)) throw std::domain_error( "In npstat::inverseIncompleteGammaC: " "argument outside of [0, 1] interval"); if (x <= 0.5) return igami(a, x); else return inverseIncompleteGamma(a, 1.0 - x); } double inverseIncompleteGamma(const double a, const double x) { if (!(x >= 0.0 && x <= 1.0)) throw std::domain_error( "In npstat::inverseIncompleteGamma: " "argument outside of [0, 1] interval"); if (x >= 0.5) return igami(a, 1.0 - x); else { const double targetEps = 8.0*MACHEP; double xmin = 0.0; double xmax = igami(a, 0.5); while ((xmax - xmin)/(xmax + xmin + 1.0) > targetEps) { const double xtry = (xmax + xmin)/2.0; if (igam(a, xtry) > x) xmax = xtry; else xmin = xtry; } return (xmax + xmin)/2.0; } } long double dawsonIntegral(const long double x_in) { const unsigned deg = 20U; const long double x = fabsl(x_in); long double result = 0.0L; if (x == 0.0L) ; else if (x <= 1.0L) { const long double epsilon = LDBL_EPSILON/2.0L; const long double twoxsq = -2.0L*x*x; long double term = 1.0L; long double dfact = 1.0L; long double xn = 1.0L; unsigned y = 0; do { result += term; y += 1U; dfact *= (2U*y + 1U); xn *= twoxsq; term = xn/dfact; } while (fabsl(term) > epsilon * fabsl(result)); result *= x; } else if (x <= 1.75L) { static const long double coeffs[deg + 1U] = { +4.563960711239483142081e-1L, -9.268566100670767619861e-2L, -7.334392170021420220239e-3L, +3.379523740404396755124e-3L, -3.085898448678595090813e-4L, -1.519846724619319512311e-5L, +4.903955822454009397182e-6L, -2.106910538629224721838e-7L, -2.930676220603996192089e-8L, +3.326790071774057337673e-9L, +3.335593457695769191326e-11L, -2.279104036721012221982e-11L, +7.877561633156348806091e-13L, +9.173158167107974472228e-14L, -7.341175636102869400671e-15L, -1.763370444125849029511e-16L, +3.792946298506435014290e-17L, -4.251969162435936250171e-19L, -1.358295820818448686821e-19L, +5.268740962820224108235e-21L, +3.414939674304748094484e-22L }; const long double midpoint = (1.75L + 1.0L)/2.0L; const long double halflen = (1.75L - 1.0L)/2.0L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } else if (x <= 2.5L) { static const long double coeffs[deg + 1U] = { +2.843711194548592808550e-1L, -6.791774139166808940530e-2L, +6.955211059379384327814e-3L, -2.726366582146839486784e-4L, -6.516682485087925163874e-5L, +1.404387911504935155228e-5L, -1.103288540946056915318e-6L, -1.422154597293404846081e-8L, +1.102714664312839585330e-8L, -8.659211557383544255053e-10L, -8.048589443963965285748e-12L, +6.092061709996351761426e-12L, -3.580977611213519234324e-13L, -1.085173558590137965737e-14L, +2.411707924175380740802e-15L, -7.760751294610276598631e-17L, -6.701490147030045891595e-18L, +6.350145841254563572100e-19L, -2.034625734538917052251e-21L, -2.260543651146274653910e-21L, +9.782419961387425633151e-23L }; const long double midpoint = (2.5L + 1.75L)/2.0L; const long double halflen = (2.5L - 1.75L)/2.0L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } else if (x <= 3.25L) { static const long double coeffs[deg + 1U] = { +1.901351274204578126827e-1L, -3.000575522193632460118e-2L, +2.672138524890489432579e-3L, -2.498237548675235150519e-4L, +2.013483163459701593271e-5L, -8.454663603108548182962e-7L, -8.036589636334016432368e-8L, +2.055498509671357933537e-8L, -2.052151324060186596995e-9L, +8.584315967075483822464e-11L, +5.062689357469596748991e-12L, -1.038671167196342609090e-12L, +6.367962851860231236238e-14L, +3.084688422647419767229e-16L, -3.417946142546575188490e-16L, +2.311567730100119302160e-17L, -6.170132546983726244716e-20L, -9.133176920944950460847e-20L, +5.712092431423316128728e-21L, +1.269641078369737220790e-23L, -2.072659711527711312699e-23L }; const long double midpoint = (3.25L + 2.5L)/2.0L; const long double halflen = (3.25L - 2.5L)/2.0L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } else if (x <= 4.25L) { static const long double coeffs[deg + 1U] = { +1.402884974484995678749e-1L, -2.053975371995937033959e-2L, +1.595388628922920119352e-3L, -1.336894584910985998203e-4L, +1.224903774178156286300e-5L, -1.206856028658387948773e-6L, +1.187997233269528945503e-7L, -1.012936061496824448259e-8L, +5.244408240062370605664e-10L, +2.901444759022254846562e-11L, -1.168987502493903926906e-11L, +1.640096995420504465839e-12L, -1.339190668554209618318e-13L, +3.643815972666851044790e-15L, +6.922486581126169160232e-16L, -1.158761251467106749752e-16L, +8.164320395639210093180e-18L, -5.397918405779863087588e-20L, -5.052069908100339242896e-20L, +5.322512674746973445361e-21L, -1.869294542789169825747e-22L }; const long double midpoint = (4.25L + 3.25L)/2.0L; const long double halflen = (4.25L - 3.25L)/2.0L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } else if (x <= 5.5L) { static const long double coeffs[deg + 1U] = { +1.058610209741581514157e-1L, -1.429297757627935191694e-2L, +9.911301703835545472874e-4L, -7.079903107876049846509e-5L, +5.229587914675267516134e-6L, -4.016071345964089296212e-7L, +3.231734714422926453741e-8L, -2.752870944370338482109e-9L, +2.503059741885009530630e-10L, -2.418699000594890423278e-11L, +2.410158905786160001792e-12L, -2.327254341132174000949e-13L, +1.958284411563056492727e-14L, -1.099893145048991004460e-15L, -2.959085292526991317697e-17L, +1.966366179276295203082e-17L, -3.314408783993662492621e-18L, +3.635520318133814622089e-19L, -2.550826919215104648800e-20L, +3.830090587178262542288e-22L, +1.836693763159216122739e-22L }; const long double midpoint = (5.5L + 4.25L)/2.0L; const long double halflen = (5.5L - 4.25L)/2.0L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } else if (x <= 7.25L) { static const long double coeffs[deg + 1U] = { +8.024637207807814739314e-2L, -1.136614891549306029413e-2L, +8.164249750628661856014e-4L, -5.951964778701328943018e-5L, +4.407349502747483429390e-6L, -3.317746826184531133862e-7L, +2.541483569880571680365e-8L, -1.983391157250772649001e-9L, +1.579050614491277335581e-10L, -1.284592098551537518322e-11L, +1.070070857004674207604e-12L, -9.151832297362522251950e-14L, +8.065447314948125338081e-15L, -7.360105847607056315915e-16L, +6.995966000187407197283e-17L, -6.964349343411584120055e-18L, +7.268789359189778223225e-19L, -7.885125241947769024019e-20L, +8.689022564130615225208e-21L, -9.353211304381231554634e-22L +9.218280404899298404756e-23L }; const long double midpoint = (7.25L + 5.5L)/2.0L; const long double halflen = (7.25L - 5.5L)/2.0L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } else { const int nterms = 50; long double term[nterms + 1]; const long double epsilon = LDBL_EPSILON/2.0L; const long double xsq = x*x; const long double twox = 2.0L*x; const long double twoxsq = 2.0L*xsq; long double xn = twoxsq; long double factorial = 1.0L; int n = 2; term[0] = 1.0L; term[1] = 1.0L/xn; for (; n <= nterms; ++n) { xn *= twoxsq; factorial *= (2*n - 1); term[n] = factorial/xn; if (term[n] < epsilon) break; } if (n > nterms) n = nterms; for (; n >= 0; --n) result += term[n]; result /= twox; } if (x_in < 0.0L) result *= -1.0L; return result; } long double inverseExpsqIntegral(const long double x_in) { const long double smallx = sqrtl(LDBL_EPSILON); const long double x = fabsl(x_in); long double result = 0.0L; // At the very beginning, the forward function behaves as x + x^3/3 if (x < smallx) result = x; // Are we below the integral upper limit of 1/2? else if (x < 0.5449871041836222236624201L) { static const long double coeffs[] = { 0.2579094020165510878448047L, 0.2509616007789004064162848L, -0.008037441097819308338583395L, -0.0009645198158737501684728918L, 0.0001298963476188684667851777L, 2.83303393956705648285913e-6L, -1.879031091504004970753414e-6L, 8.89033414496187588540693e-8L, 2.190598580039940838615595e-8L, -2.958785324082878640818902e-9L, -1.393055660359914508583854e-10L, 5.934164278634234871603216e-11L, -2.038318546964313720058418e-12L, -8.718325681035840557959458e-13L, 1.013108738230624177815487e-13L, 7.80188419892238096141015e-15L, -2.410959710954351931245264e-15L, 4.28768887009677018630627e-17L, 4.075320326600366577580808e-17L, -3.925657400094226595458262e-18L, -4.52032825617882087470387e-19L, 1.08172986581534678569357e-19L, 7.619258918621063777867205e-23L, -2.043372106982501109979208e-21L, 1.577535083736277024466253e-22L, 2.659861973971217145144534e-23L }; static const unsigned deg = sizeof(coeffs)/sizeof(coeffs[0]) - 1U; const long double midpoint = 0.2724935520918111118312101L; const long double halflen = midpoint; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } // Are we below the integral upper limit of 1? else if (x < 1.462651745907181608804049L) { static const long double coeffs[] = { 0.7736025700148903044157852L, 0.2483094370727645101600738L, -0.0236047661681907541554642L, 0.001718250905865626423590199L, -3.383723505479159579142638e-6L, -0.00002706856684823251983606182L, 5.564585822987557167324481e-6L, -6.297400054716145999859568e-7L, 1.793804726508794465230363e-8L, 9.974525492565834910508807e-9L, -2.630092849889704961303573e-9L, 3.584362156719671871930982e-10L, -1.848479157680947543606944e-11L, -4.505786185834443101923114e-12L, 1.497497786049082262992397e-12L, -2.345976552921671025035597e-13L, 1.675613965321671764579391e-14L, 2.080526674977816958132577e-15L, -9.222206647061674196609463e-16L, 1.635927767802667303380293e-16L, -1.458481834781130100146786e-17L, -8.477734754986502993701378e-19L, 5.883853143077372799032853e-19L, -1.178516243865358739330946e-19L, 1.244597734746846080964198e-20L, 1.834734630829857873717171e-22L, -3.79697806729613468087694e-22L, 8.640525246432963515541941e-23L, -1.072975125060893395507772e-23L }; static const unsigned deg = sizeof(coeffs)/sizeof(coeffs[0]) - 1U; const long double midpoint = 1.003819425045401916233234L; const long double halflen = 0.4588323208617796925708142L; result = chebyshevSeriesSum(coeffs, deg, (x-midpoint)/halflen); } // Are we below the integral upper limit of 3/2? else if (x < 4.063114058624186262070998L) { static const long double coeffs[] = { 1.237474264255216149826922L, 0.2515376732927722040966062L, 0.01255809348578646846677752L, -0.001555984193633101089113648L, -0.00003266175488449684464535135L, 0.00001861132693718087179337154L, 3.038788024730725430182981e-7L, -3.055587383693358717339231e-7L, 2.195404333387651282823557e-10L, 5.223891523975639259623118e-9L, -8.74172529522376286448331e-11L, -9.289565568740836123320397e-11L, 3.037929091913979095346417e-12L, 1.697081673433425594636218e-12L, -8.376226831184981721677363e-14L, -3.144033465184008496066638e-14L, 2.107731553034119316083898e-15L, 5.865893869124205701237769e-16L, -5.061884931584969465772662e-17L, -1.096017981587409945712935e-17L, 1.182058217010476282761648e-18L, 2.040660311333339268713157e-19L, -2.708723954709022960299832e-20L, -3.767175977793052434129249e-21L, 6.118697880934383791686931e-22L, 6.976999259292224769625502e-23L }; static const unsigned deg = sizeof(coeffs)/sizeof(coeffs[0]) - 1U; const long double midpoint = 0.9003422806239746400989134L; const long double halflen = 0.2836972833794905219117775L; result = chebyshevSeriesSum(coeffs, deg, (sqrtl(logl(x))-midpoint)/halflen); } // Are we below the integral upper limit of 2? else if (x < 16.4526277655072302247364L) { static const long double coeffs[] = { 1.750318675095801084701066L, 0.2503081949080952607456335L, -0.0003683988977021577213047602L, -0.0003053217430231520892093692L, 0.00004984978059867808464521542L, -2.910067173452788393366537e-6L, -1.234138038866370318487004e-7L, 3.698515524218906567339791e-8L, -2.598716127513287147666407e-9L, -8.028249755820916232073154e-11L, 3.38719730193993051252875e-11L, -2.805688393204576054583932e-12L, -4.636644559366711367183408e-14L, 3.430651608591744447051587e-14L, -3.234144260654980846019301e-15L, -1.479571272220063709849795e-17L, 3.663824718035758023592578e-17L, -3.866514289799821609911763e-18L, 2.280975526810673800594954e-20L, 4.039168081650367776987e-20L, -4.73516370177886271900388e-21L, 7.25430770939280219423293e-23L, 4.526173229162084391018431e-23L }; static const unsigned deg = sizeof(coeffs)/sizeof(coeffs[0]) - 1U; const long double midpoint = 1.428752297062238501405865L; const long double halflen = 0.2447127330587733393951737L; result = chebyshevSeriesSum(coeffs, deg, (sqrtl(logl(x))-midpoint)/halflen); } // Are we below the integral upper limit of 3? else if (x < 1444.545122892714154713760L) { static const long double coeffs[] = { 2.50294099696623892772079L, 0.4994357077985400736034567L, -0.00291622179390211582491955L, 0.0005756988194441929049790467L, -0.00002810720142885454120222505L, -0.00001091749276776510487675747L, 3.29093808250468757948933e-6L, -4.888188874032568897843031e-7L, 4.172102909667302951211426e-8L, -4.599575541801073701575834e-10L, -6.053857271985015243693571e-10L, 1.509946372242848844836431e-10L, -2.455367299247951048634247e-11L, 2.672752200212154657233337e-12L, -9.103858678297231810743428e-14L, -3.718269817594061105825935e-14L, 1.069438018054794309965648e-14L, -1.744407969349348513403047e-15L, 1.883808777693390076608131e-16L, -7.095219918870126016632479e-18L, -2.561647683564104671111204e-18L, 7.953272251856515822155603e-19L, -1.363875890657281736925731e-19L, 1.519215969460881098898382e-20L, -6.014662801208487581697383e-22L, -2.017882982607174277058989e-22L, 6.264772567011851344575421e-23L }; static const unsigned deg = sizeof(coeffs)/sizeof(coeffs[0]) - 1U; const long double midpoint = 2.185393865913887927395372L; const long double halflen = 0.5119288357928760865943334L; result = chebyshevSeriesSum(coeffs, deg, (sqrtl(logl(x))-midpoint)/halflen); } else { // At this point, sqrtl(logl(x)) looks more or less like // a straight line. We can do Newton-Raphson efficiently // in this transformed variable. const long double target = sqrtl(logl(x)); long double xtry = target, xold = 0.0L; while (fabsl((xtry - xold)/xtry) > 2.0*LDBL_EPSILON) { xold = xtry; const long double daws = dawsonIntegral(xtry); const long double sqrval = sqrtl(xtry*xtry + logl(daws)); const long double deriv = 0.5L/sqrval/daws; xtry += (target - sqrval)/deriv; } result = xtry; } if (x_in < 0.0L) result *= -1.0L; return result; } long double normalDensityDerivative(const unsigned order, const long double x) { static std::vector > coeffs; const unsigned ncoeffs = coeffs.size(); if (order >= ncoeffs) coeffs.resize(order + 1); unsigned nc = coeffs[order].size(); if (nc == 0U) { buildGaussDeriCoeffs(order, &coeffs[order]); nc = coeffs[order].size(); } assert(nc == order + 1U); return polySeriesSum(&coeffs[order][0], order, x)* expl(-x*x/2.0L)/SQRTTWOPIL; } double bivariateNormalIntegral(const double rho, const double x, const double y) { const double absrho = fabs(rho); if (absrho > 1.0) throw std::domain_error( "In npstat::bivariateNormalIntegral: " "correlation coefficient is outside of the [-1, 1] interval"); if (absrho < 0.925) return Phi2(-x, -y, rho); else return bvnd(x, y, rho); } + + double besselK(const double nu, const double x) + { + if (nu < 0.0 || x <= 0.0) throw std::domain_error( + "In npstat::besselK: arguments outside of the allowed range"); + + const unsigned maxIter = 10000U; + const double fpMin = DBL_MIN/DBL_EPSILON; + const double xmin = 2.0; + + const int nl = std::floor(nu + 0.5); + const double xmu = nu - nl; + assert(std::abs(xmu) <= 0.5); + const double xmu2 = xmu*xmu; + const double xi = 1.0/x; + const double xi2 = 2.0*xi; + double h = nu*xi; + if (h < fpMin) + h = fpMin; + double b = xi2*nu; + double d = 0.0; + double c = h; + + bool converged = false; + for (unsigned i=0; i= 0; --l) + { + const double ritemp = fact*ril + ripl; + fact -= xi; + ripl = fact*ritemp + ril; + ril = ritemp; + } + + double rkmu, rk1; + if (x < xmin) + { + const double xover2 = x/2.0; + const long double ln2overx = -logl(xover2); + const long double sigma = xmu*ln2overx; + const double sqrEps = sqrt(DBL_EPSILON); + const double pimu = M_PI*xmu; + const double prefact = std::abs(pimu) < 0.1*sqrEps ? 1.0 : pimu/sin(pimu); + const long double sqrlEps = sqrtl(LDBL_EPSILON); + const long double fact2 = std::abs(sigma) < 0.1*sqrlEps ? 1.0L : sinhl(sigma)/sigma; + const long double xover2gam = powl(xover2, xmu); + const double gamplus = Gamma(1.0 + xmu); + const double gamminus = Gamma(1.0 - xmu); + const double gam2 = 0.5/gamplus + 0.5/gamminus; + double gam1; + if (std::abs(xmu) > 0.01) + gam1 = (0.5/gamminus - 0.5/gamplus)/xmu; + else + { + static const long double coeffs[] = { + -0.577215664901532860606512L, + 0.0420026350340952355290039L, + 0.0421977345555443367482083L, + 0.00721894324666309954239501L, + 0.000215241674114950972815730L, + 0.0000201348547807882386556894L + }; + gam1 = polySeriesSum(coeffs, 5U, xmu2); + } + long double p = 0.5/xover2gam*gamplus; + long double q = 0.5*xover2gam*gamminus; + long double ff = prefact*(gam1*coshl(sigma) + gam2*fact2*ln2overx); + long double sum = ff; + long double sum1 = p; + long double ck = 1.0L; + const double xover2sq = xover2*xover2; + + converged = false; + for (unsigned i=1; i<=maxIter && !converged; ++i) + { + ff = (i*ff + p + q)/(i*i - xmu2); + ck *= (xover2sq/i); + p /= (i - xmu); + q /= (i + xmu); + const long double del = ck*ff; + sum += del; + const long double del1 = ck*(p - i*ff); + sum1 += del1; + if (std::abs(del) < std::abs(sum)*DBL_EPSILON) + converged = true; + } + if (!converged) throw std::runtime_error( + "In npstat::besselK: series failed to converge"); + rkmu = sum; + rk1 = sum1*xi2; + } + else + { + b = 2.0*(1.0 + x); + d = 1.0/b; + h = d; + double delh = d; + double q1 = 0.0; + double q2 = 1.0; + const double a1 = 0.25 - xmu2; + double q = a1; + c = a1; + double a = -a1; + double s = 1.0 + q*delh; + + converged = false; + for (unsigned i=1; i<=maxIter && !converged; ++i) + { + a -= 2*i; + c = -a*c/(i+1.0); + const double qnew = (q1 - b*q2)/a; + q1 = q2; + q2 = qnew; + q += c*qnew; + b += 2.0; + d = 1.0/(b+a*d); + delh = (b*d-1.0)*delh; + h += delh; + const double dels = q*delh; + s += dels; + if (std::abs(dels/s) <= DBL_EPSILON) + converged = true; + } + if (!converged) throw std::runtime_error( + "In npstat::besselK: failed to converge in cf2"); + + h = a1*h; + rkmu = std::sqrt(M_PI/(2.0*x))*std::exp(-x)/s; + rk1 = rkmu*(xmu + x + 0.5 - h)*xi; + } + + for (int i=1; i<=nl; ++i) + { + const double rktemp = (xmu + i)*xi2*rk1 + rkmu; + rkmu = rk1; + rk1 = rktemp; + } + + return rkmu; + } } Index: trunk/npstat/stat/AbsUnbinnedGOFTest1D.cc =================================================================== --- trunk/npstat/stat/AbsUnbinnedGOFTest1D.cc (revision 741) +++ trunk/npstat/stat/AbsUnbinnedGOFTest1D.cc (revision 742) @@ -1,31 +1,57 @@ #include +#include "npstat/nm/findRootUsingBisections.hh" + #include "npstat/stat/AbsUnbinnedGOFTest1D.hh" namespace npstat { void AbsUnbinnedGOFTest1D::simulateStatistic( AbsRandomGenerator& g, const unsigned long sz, const unsigned nPseudo, std::vector* stats) const { if (!sz) throw std::invalid_argument( "In npstat::AbsUnbinnedGOFTest1D::simulateStatistic: " "sample size must be positive"); if (nPseudo < 2U) std::invalid_argument( "In npstat::AbsUnbinnedGOFTest1D::simulateStatistic: " "insufficient number od pseudo experiments"); assert(stats); stats->clear(); stats->reserve(nPseudo); std::vector sampleVec(sz); double* sample = &sampleVec[0]; for (unsigned ips=0; ipsrandom(g, sample+i); std::sort(sampleVec.begin(), sampleVec.end()); stats->push_back(this->testStatistic(sample, sz, true)); } std::sort(stats->begin(), stats->end()); } + + double AbsUnbinnedGOFTest1D::inverseExceedance( + const double pvalue, const unsigned long sz, + const double smin, const double smax) const + { + static const double tol = 1.0e-14; + + if (!hasAnalyticPValue()) throw std::runtime_error( + "In npstat::AbsUnbinnedGOFTest1D::inverseExceedance: " + "analytic p-value calculation is not implemented by " + "the derived class"); + + if (pvalue <= 0.0 || pvalue >= 1.0) throw std::invalid_argument( + "In npstat::AbsUnbinnedGOFTest1D::inverseExceedance: " + "p-value argument must be inside (0, 1) interval"); + + double q = 0.0; + if (!findRootUsingBisections( + GOFTest1DPVFunctor(*this, sz), pvalue, + smin, smax, tol, &q)) throw std::invalid_argument( + "In npstat::AbsUnbinnedGOFTest1D::inverseExceedance: " + "the initial interval does not bracket the root"); + return q; + } } Index: trunk/npstat/stat/OrthoPolyGOFTest1D.hh =================================================================== --- trunk/npstat/stat/OrthoPolyGOFTest1D.hh (revision 741) +++ trunk/npstat/stat/OrthoPolyGOFTest1D.hh (revision 742) @@ -1,78 +1,83 @@ #ifndef NPSTAT_ORTHOPOLYGOFTEST1D_HH_ #define NPSTAT_ORTHOPOLYGOFTEST1D_HH_ /*! // \file OrthoPolyGOFTest1D.hh // // \brief Orthogonal polynomial goodness-of-fit test // // Author: I. Volobouev // // November 2020 */ #include #include "npstat/nm/Matrix.hh" #include "npstat/nm/AbsClassicalOrthoPoly1D.hh" #include "npstat/stat/AbsUnbinnedGOFTest1D.hh" namespace npstat { class OrthoPolyGOFTest1D : public AbsUnbinnedGOFTest1D { public: template OrthoPolyGOFTest1D(const AbsDistribution1D& distro, const Quadrature& q, const AbsClassicalOrthoPoly1D& poly, unsigned maxdeg); inline virtual ~OrthoPolyGOFTest1D() {} + virtual std::string shortName() const; + inline unsigned maxDegree() const {return maxdeg_;} inline unsigned numDeviations() const {return 2U*maxdeg_;} /** // Get the value of the test statistic with all // terms weighted equally */ inline virtual double testStatistic( const double* data, const unsigned long sz, bool) const {return testStat(data, sz);} inline virtual double testStatistic( const float* data, const unsigned long sz, bool) const {return testStat(data, sz);} inline virtual bool hasAnalyticPValue() const {return true;} virtual double analyticPValue(double stat, unsigned long sz) const; + virtual double inverseExceedance(double pvalue, unsigned long sz, + double smin, double smax) const; + /** // Get individual normalized deviations from which // the test statistic is constructed */ template void normalizedDeviations( const Numeric* data, unsigned long lenData, double* deviations, unsigned lenDeviations) const; private: template void calculateNormalizingMatrix(const Quadrature& q, Matrix* m) const; template double testStat(const Numeric* data, unsigned long lenData) const; std::shared_ptr poly_; Matrix invsqr_; std::vector buf_; std::vector statBuf_; unsigned maxdeg_; }; } #include "npstat/stat/OrthoPolyGOFTest1D.icc" #endif // NPSTAT_ORTHOPOLYGOFTEST1D_HH_ Index: trunk/npstat/stat/UnbinnedGOFTests1D.cc =================================================================== --- trunk/npstat/stat/UnbinnedGOFTests1D.cc (revision 741) +++ trunk/npstat/stat/UnbinnedGOFTests1D.cc (revision 742) @@ -1,21 +1,176 @@ +#include +#include #include +#include #include "kstest/kstest.hh" +#include "npstat/nm/SpecialFunctions.hh" + #include "npstat/stat/UnbinnedGOFTests1D.hh" #include "npstat/stat/StatUtils.hh" +namespace { + double cvm_G(const double y) + { + return -std::exp(-y)*(npstat::besselK(0.25, y) + npstat::besselK(0.75, y)); + } + + double cvm_H(const double y) + { + return std::exp(-y)*(npstat::besselK(1.25, y) - 3.0*npstat::besselK(0.75, y) - + 2.0*npstat::besselK(0.25, y)); + } + + double cvm_Ak(const unsigned k, const double x) + { + const double fourkp1 = 4U*k + 1U; + const double fourkp3 = 4U*k + 3U; + const double fourkp5 = 4U*k + 5U; + const double sixteenx = 16.0*x; + + return 7*std::pow(fourkp1, 1.5)*cvm_G(fourkp1*fourkp1/sixteenx) + + 16*std::pow(fourkp3, 1.5)*cvm_G(fourkp3*fourkp3/sixteenx) + + 7*std::pow(fourkp5, 1.5)*cvm_G(fourkp5*fourkp5/sixteenx); + } + + double cvm_Bk(const unsigned k, const double x) + { + const double fourkp1 = 4U*k + 1U; + return std::pow(fourkp1, 2.5)*cvm_H(fourkp1*fourkp1/(16.0*x)); + } + + double cvm_Ck(const unsigned k, const double x) + { + const double fourkp5 = 4U*k + 5U; + return 24*std::pow(fourkp5, 2.5)*cvm_H(fourkp5*fourkp5/(16.0*x)); + } + + double cvm_phi1_partA(const double x) + { + const unsigned maxIter = 1000; + + long double sum = 0.0L; + long double gammaTerm = 0.886226925452758013649084L; + bool converged = false; + for (unsigned k=0; k 10U && std::abs(delta) < std::abs(sum)*DBL_EPSILON) + converged = true; + } + if (!converged) throw std::runtime_error( + "In cvm_phi1_partA: series did not converge"); + sum /= std::pow(M_PI, 1.5); + sum /= (576*std::pow(x, 1.5)); + return sum; + } + + double cvm_phi1_partB(const double x) + { + const unsigned maxIter = 1000; + + long double sum = 0.0L; + long double gammaTerm = 1.77245385090551602729817L; + bool converged = false; + for (unsigned k=0; k 10U && std::abs(delta) < std::abs(sum)*DBL_EPSILON) + converged = true; + } + if (!converged) throw std::runtime_error( + "In cvm_phi1_partB: series did not converge"); + sum /= std::pow(M_PI, 1.5); + sum /= (2304*std::pow(x, 2.5)); + return sum; + } + + double cvm_phi1_partC(const double x) + { + const unsigned maxIter = 1000; + + long double sum = 0.0L; + long double gammaTerm = 1.32934038817913702047363L; + bool converged = false; + for (unsigned k=0; k 10U && std::abs(delta) < std::abs(sum)*DBL_EPSILON) + converged = true; + } + if (!converged) throw std::runtime_error( + "In cvm_phi1_partC: series did not converge"); + sum /= std::pow(M_PI, 1.5); + sum /= (2304*std::pow(x, 2.5)); + return sum; + } + + double cvm_Vinf(const double x) + { + const unsigned maxIter = 1000; + + long double sum = 0.0L; + long double gammaTerm = 1.77245385090551602729817L; + bool converged = false; + for (unsigned k=0; k 10U && std::abs(delta) < std::abs(sum)*DBL_EPSILON) + converged = true; + } + if (!converged) throw std::runtime_error( + "In cvm_Vinf: series did not converge"); + sum /= (std::sqrt(x)*pow(M_PI, 1.5)); + return sum; + } +} + namespace npstat { double KSTest1D::analyticPValue(const double stat, const unsigned long sz) const { + if (!sz) throw std::invalid_argument( + "In npstat::KSTest1D::analyticPValue: sample size must be positive"); assert(sz <= UINT_MAX); return kstest::KSfbar(sz, stat); } double ADTest1D::analyticPValue(const double stat, const unsigned long sz) const { + if (!sz) throw std::invalid_argument( + "In npstat::ADTest1D::analyticPValue: sample size must be positive"); return 1.0 - AD(sz, stat); } + + // The algorithm for the following comes from S. Csorgo and J.J. Faraway, + // "The Exact and Asymptotic Distributions of Cramer-von Mises Statistics", + // Journal of the Royal Statistical Society, Series B (Methodological), + // Vol. 58, No. 1 (1996), pp. 221-234. + double CvMTest1D::analyticPValue(const double x, + const unsigned long sz) const + { + if (!sz) throw std::invalid_argument( + "In npstat::CvMTest1D::analyticPValue: sample size must be positive"); + const double vinf = cvm_Vinf(x); + const double cdf = vinf + (vinf/12.0 + cvm_phi1_partA(x) + + cvm_phi1_partB(x) + cvm_phi1_partC(x))/sz; + return 1.0 - cdf; + } } Index: trunk/npstat/stat/AbsUnbinnedGOFTest1D.hh =================================================================== --- trunk/npstat/stat/AbsUnbinnedGOFTest1D.hh (revision 741) +++ trunk/npstat/stat/AbsUnbinnedGOFTest1D.hh (revision 742) @@ -1,69 +1,94 @@ #ifndef NPSTAT_ABSUNBINNEDGOFTEST1D_HH_ #define NPSTAT_ABSUNBINNEDGOFTEST1D_HH_ /*! // \file AbsUnbinnedGOFTest1D.hh // // \brief Interface definition for goodness-of-fit tests for 1-d distributions // // Author: I. Volobouev // // November 2020 */ #include #include +#include +#include "npstat/nm/SimpleFunctors.hh" #include "npstat/rng/AbsRandomGenerator.hh" #include "npstat/stat/AbsDistribution1D.hh" namespace npstat { class AbsUnbinnedGOFTest1D { public: inline explicit AbsUnbinnedGOFTest1D(const AbsDistribution1D& d) : distro_(d.clone()) {} inline AbsUnbinnedGOFTest1D(const AbsUnbinnedGOFTest1D& r) : distro_(r.distro_->clone()) {} inline AbsUnbinnedGOFTest1D& operator=(const AbsUnbinnedGOFTest1D& r) { if (this != &r) { delete distro_; distro_ = 0; distro_ = r.distro_->clone(); } return *this; } inline virtual ~AbsUnbinnedGOFTest1D() {delete distro_;} + virtual std::string shortName() const=0; + virtual double testStatistic( const double* data, unsigned long sz, bool isDataSorted) const=0; virtual double testStatistic( const float* data, unsigned long sz, bool isDataSorted) const=0; inline virtual double analyticPValue( double /* stat */, unsigned long /* sz */) const { throw std::runtime_error( "In npstat::AbsUnbinnedGOFTest1D::analyticPValue: " "this function is not implemented by the derived class"); } inline virtual bool hasAnalyticPValue() const {return false;} + /** This is the inverse function to "analyticPValue" */ + virtual double inverseExceedance(double pvalue, unsigned long sz, + double smin, double smax) const; + // The following method should sort the statistic values // in the increasing order virtual void simulateStatistic( AbsRandomGenerator& g, unsigned long sz, unsigned nPseudo, std::vector* stats) const; protected: AbsDistribution1D* distro_; }; + + class GOFTest1DPVFunctor : public Functor1 + { + public: + inline GOFTest1DPVFunctor(const AbsUnbinnedGOFTest1D& goftest, + const unsigned long sz) + : goftest_(goftest), sz_(sz) {} + + inline virtual ~GOFTest1DPVFunctor() {} + + inline virtual double operator()(const double& a) const + {return goftest_.analyticPValue(a, sz_);} + + private: + const AbsUnbinnedGOFTest1D& goftest_; + unsigned long sz_; + }; } #endif // NPSTAT_ABSUNBINNEDGOFTEST1D_HH_ Index: trunk/npstat/stat/UnbinnedGOFTests1D.hh =================================================================== --- trunk/npstat/stat/UnbinnedGOFTests1D.hh (revision 741) +++ trunk/npstat/stat/UnbinnedGOFTests1D.hh (revision 742) @@ -1,87 +1,97 @@ #ifndef NPSTAT_UNBINNEDGOFTESTS1D_HH_ #define NPSTAT_UNBINNEDGOFTESTS1D_HH_ #include "npstat/stat/AbsUnbinnedGOFTest1D.hh" namespace npstat { /** Kolmogorov-Smirnov test */ class KSTest1D : public AbsUnbinnedGOFTest1D { public: inline explicit KSTest1D(const AbsDistribution1D& d) : AbsUnbinnedGOFTest1D(d) {} inline virtual ~KSTest1D() {} + inline virtual std::string shortName() const {return "KS";} + inline virtual double testStatistic( const double* data, const unsigned long sz, const bool b) const {return testStat(data, sz, b);} inline virtual double testStatistic( const float* data, const unsigned long sz, const bool b) const {return testStat(data, sz, b);} inline virtual bool hasAnalyticPValue() const {return true;} virtual double analyticPValue(double stat, unsigned long sz) const; private: template double testStat(const Numeric* data, unsigned long lenData, bool dataIsSorted) const; }; /** Anderson-Darling test */ class ADTest1D : public AbsUnbinnedGOFTest1D { public: inline explicit ADTest1D(const AbsDistribution1D& d) : AbsUnbinnedGOFTest1D(d) {} inline virtual ~ADTest1D() {} + inline virtual std::string shortName() const {return "AD";} + inline virtual double testStatistic( const double* data, const unsigned long sz, const bool b) const {return testStat(data, sz, b);} inline virtual double testStatistic( const float* data, const unsigned long sz, const bool b) const {return testStat(data, sz, b);} inline virtual bool hasAnalyticPValue() const {return true;} virtual double analyticPValue(double stat, unsigned long sz) const; private: template double testStat(const Numeric* data, unsigned long lenData, bool dataIsSorted) const; }; /** Cramer-von Mises test */ class CvMTest1D : public AbsUnbinnedGOFTest1D { public: inline explicit CvMTest1D(const AbsDistribution1D& d) : AbsUnbinnedGOFTest1D(d) {} inline virtual ~CvMTest1D() {} + inline virtual std::string shortName() const {return "CvM";} + inline virtual double testStatistic( const double* data, const unsigned long sz, const bool b) const {return testStat(data, sz, b);} inline virtual double testStatistic( const float* data, const unsigned long sz, const bool b) const {return testStat(data, sz, b);} + inline virtual bool hasAnalyticPValue() const {return true;} + + virtual double analyticPValue(double stat, unsigned long sz) const; + private: template double testStat(const Numeric* data, unsigned long lenData, bool dataIsSorted) const; }; } #include "npstat/stat/UnbinnedGOFTests1D.icc" #endif // NPSTAT_UNBINNEDGOFTESTS1D_HH_ Index: trunk/npstat/stat/OrthoPolyGOFTest1D.cc =================================================================== --- trunk/npstat/stat/OrthoPolyGOFTest1D.cc (revision 741) +++ trunk/npstat/stat/OrthoPolyGOFTest1D.cc (revision 742) @@ -1,12 +1,39 @@ +#include + #include "npstat/stat/OrthoPolyGOFTest1D.hh" #include "npstat/stat/Distributions1D.hh" namespace npstat { double OrthoPolyGOFTest1D::analyticPValue( const double stat, unsigned long /* sz */) const { // ndof = 2*maxdeg_ Gamma1D chisq(0.0, 2.0, maxdeg_); return chisq.exceedance(stat); } + + std::string OrthoPolyGOFTest1D::shortName() const + { + std::ostringstream os; + os << "OP" << maxdeg_; + return os.str(); + } + + double OrthoPolyGOFTest1D::inverseExceedance( + const double pvalue, const unsigned long sz, + const double smin, const double smax) const + { + if (pvalue <= 0.0 || pvalue >= 1.0) throw std::invalid_argument( + "In npstat::OrthoPolyGOFTest1D::inverseExceedance: " + "p-value argument must be inside (0, 1) interval"); + + if (pvalue > 1.0e-10) + { + Gamma1D chisq(0.0, 2.0, maxdeg_); + return chisq.quantile(1.0 - pvalue); + } + else + return AbsUnbinnedGOFTest1D::inverseExceedance( + pvalue, sz, smin, smax); + } } Index: trunk/NEWS =================================================================== --- trunk/NEWS (revision 741) +++ trunk/NEWS (revision 742) @@ -1,823 +1,825 @@ Version 5.3.0 - development * Added some facilities for quadruple precision calculations. They mostly rely on the __float128 type (gcc extension, libquadmath) and on compiling the LAPACK/BLAS libraries with the "-freal-8-real-16" switch. * Added the "GaussLegendreQuadratureQ" class capable of performing Gauss-Legendre quadratures in quadruple precision. * Added the "GaussLegendreQuadrature2D" class capable of performing tensor-product Gauss-Legendre quadratures in two dimensions. * Added OrthoPoly1DDeg functor. * Added class ScalableClassicalOrthoPoly1D. * Added class ClassicalOrthoPoly1DFromWeight. This also necessitated refactoring of the class DensityOrthoPoly1D. * Added classes AbsIntervalQuadrature1D and RectangleQuadrature1D. * Added function "kernelSensitivityMatrix". * Added function "sumOfSquares". * Added a number of .i files. * Added classes AbsUnbinnedGOFTest1D and OrthoPolyGOFTest1D. * Added class PolynomialDistro1D. +* Added a number of unbinned goodness-of-fit tests. + Version 5.2.0 - July 20 2020, by I. Volobouev * Added a number of .i files. * Added classes "DiscreteGauss1DBuilder" and "DiscreteGaussCopulaSmoother". * Added class "MatrixFilter1DBuilder". * Added facilities for saddlepoint approximations: classes AbsCGF1D and SeriesCGF1D, function saddlepointDistribution1D, as well as class ExpTiltedDistribution1D which can be useful in saddlepoint approximation precision studies. * Added class "DensityOrthoPoly1D" for building orthonormal polynomial systems using (almost) arbitrary 1-d densities as weights. * Added "empiricalMoment" method to the AbsDistribution1D class. * Added class "HeatEq1DNeumannBoundary". * Added methods "sampleAverages" and "sampleProductAverages" to the ContOrthoPoly1D class. * Added method "expectation" to the AbsDistribution1D class. Version 5.1.0 - Nov 17 2019, by I. Volobouev * Started using python-specific features (pythonprepend, etc) in the SWIG .i files. * Added function "scannedKSDistance". * Added Python API for persistence of several classes. Version 5.0.0 - Oct 17 2019, by I. Volobouev * C++11 support is now mandatory. * Added more classes and functions to Python API. Overhauled some of the existing interfaces. Updated everything to python3. * Increased the maximum order of "classical" Edgeworth expansions to 14. * Increased the maximum order of conversions between 1-d central moments and cumulants to 20. * Added "histoQuantiles" function. * Added "correctDensityEstimateGHU" function. * Added "randomHistoFill1D" utility function. * Added ComparisonDistribution1D class. * Added classes GaussND, LinTransformedDistroND, and DistributionMixND. * Added a lot of SWIG .i files Version 4.11.0 - July 22 2019, by I. Volobouev * Added support for cumulant calculations for various Wald statistics in the Poisson process model. * Added functions convertCumulantsToCentralMoments and convertCentralMomentsToCumulants (header file cumulantConversion.hh). * Added function "cumulantUncertainties". Version 4.10.0 - July 11 2019, by I. Volobouev * Added SemiInfGaussianQuadrature class. * Added functions arrayMoment, arrayMoments, and arrayCentralMoments. * Added enum EdgeworthSeriesMethod and class EdgeworthSeries1D. * Added DeltaMixture1D class. * Added enum LikelihoodStatisticType. * Added functions "mixtureModelCumulants" and "poissonProcessCumulants" in the header likelihoodStatisticCumulants.hh. Version 4.9.0 - Dec 18 2018, by I. Volobouev * Added "integratePoly" and "jointIntegral" methods to the AbsClassicalOrthoPoly1D class. * Added utility functions "truncatedInverseSqrt" and "matrixIndexPairs". * Added a number of functions to "orthoPoly1DVProducts.hh". * Added classes ChebyshevOrthoPoly1st and ChebyshevOrthoPoly2nd inheriting from AbsClassicalOrthoPoly1D. * Added class HermiteProbOrthoPoly1D. * Added FejerQuadrature class. * Added classe IsoscelesTriangle1D and Logistic1D. * Added classes AbsKDE1DKernel and KDE1DHOSymbetaKernel. * Added static function "optimalDegreeHart" to OSDE1D class. Version 4.8.0 - Jul 9 2018, by I. Volobouev * Added ShiftedLegendreOrthoPoly1D class. * Added Lanczos method to generate recurrence coefficients for the ContOrthoPoly1D class. * Added npstat/stat/orthoPoly1DVProducts.hh file with various utilities for statistical analyis of chaos polynomials. Version 4.7.0 - Jan 13 2018, by I. Volobouev * Added "UniPareto1D" distribution (uniform with Pareto tail to the right). * More coordinate/weight cases for the GaussLegendreQuadrature class. * Added ContOrthoPoly1D class -- continuous orthogonal polynomials with discrete measure. * Added functions "linearLeastSquares" and "tdSymEigen" to the Matrix class. * Added OSDE1D class. * Added classes LocationScaleFamily1D and SinhAsinhTransform1D. * Added new functors (CdfFunctor1D, etc) as AbsDistribution1D helpers. * Small fix in StatAccumulatorArr.cc. Version 4.6.0 - Jan 23 2017, by I. Volobouev * Updated 1-d LOrPE cross validation code (classes AbsBandwidthCV, BandwidthCVLeastSquares1D, BandwidthCVPseudoLogli1D) for use with weighted samples in the case the sample itself is available at the point cross validation is run. Version 4.5.0 - Aug 01 2016, by I. Volobouev * A small fix in OrthoPolyND.icc (switched from cycling over unsigned to unsigned long in the scalar product function). * Implemented Gauss-Hermite quadrature with Gaussian density weight. * Changed the meaning of Quadratic1D and LogQuadratic1D parameters to be consistent with Legendre polynomial coefficients on [-1, 1] (new parameters are 1/2 of old). * Added class MinuitUnbinnedFitFcn1D (to interfaces). * Added function "findRootNewtonRaphson". * Added "statUncertainties" header with various functions. Version 4.4.0 - May 9 2016, by I. Volobouev * Added "timestamp" function. * Improved implementation of BinnedDensity1D::unscaledQuantile function. Certain problems caused by round-off errors are now fixed. * Added the capability to use the single closest parameter cells (thus disabling actual interpolation between parameter values, for speed) to "GridInterpolatedDistribution". Version 4.3.0 - March 19 2016, by I. Volobouev * Put additional assert statements inside OrthoPolyND.icc to prevent phony "array subscript is above array bounds" messages in g++ 4.9.2. * Improved CmdLine.hh. * Additional methods in CopulaInterpolationND and GridInterpolatedDistribution. * Added function "volumeDensityFromBinnedRadial". * Added convenience method "globalFilter" to the OrthoPoly1D class. * Initiated the transition of the Python API from Python 2 to Python 3. Version 4.2.0 - October 29 2015, by I. Volobouev * Added interpolation methods for the marginals to classes "CopulaInterpolationND" and "GridInterpolatedDistribution". * Removed assert on underflow in the "igamc" function. Now in such cases this function simply returns 0.0. Version 4.1.0 - July 27 2015, by I. Volobouev * Made a few additional methods virtual in AbsNtuple. * Declared methods "columnIndices" of AbsNtuple const (as they should be). * Added function "weightedCopulaHisto" to build copulas for samples of weighted points. * Added function "weightedVariableBandwidthSmooth1D" to use variable bandwidth smoothing with weighted histograms. * Added "AbsWeightedMarginalSmoother" interface class for smoothing samples of weighted points. Modified classes ConstantBandwidthSmoother1D, JohnsonKDESmoother, and LOrPEMarginalSmoother to support this interface. * Added class "VariableBandwidthSmoother1D" which implements both AbsMarginalSmoother and AbsWeightedMarginalSmoother interfaces. * Implemented cross-validation for weighted samples. * Implemented "buildInterpolatedCompositeDistroND" for generic construction of multivariate transfer functions. * Implemented "buildInterpolatedDistro1DNP" for generic construction of univariate transfer functions. Version 4.0.1 - June 17 2015, by I. Volobouev * Added "dump_qmc" example executable. Version 4.0.0 - June 10 2015, by I. Volobouev * Complete overhaul of 1-d filter-building code. Addition of new boundary methods is a lot easier now. The user API for choosing a LOrPE boundary method is encapsulated in the new "BoundaryHandling" class. * Implemented a number of new filter builders with different boundary treatments. * Updated the "LocalPolyFilter1D" class so that it holds the local bandwidth factors derived by the filter builders. Version 3.8.0 - June 1 2015, by I. Volobouev * Implemented class ConstSqFilter1DBuilder (declared in the header file npstat/stat/Filter1DBuilders.hh). The "BoundaryMethod" enum has been extended accordingly. Other files using this enum have been updated. * Implemented class FoldingSqFilter1DBuilder. Similar to ConstSqFilter1DBuilder but it also folds the kernel in addition to stretching it. * Added virtual destructors to a number of classes. * Added "externalMemArrayND" with various signatures to allow the use of ArrayND with memory not managed by ArrayND. * Added move constructors and move assignment operators to ArrayND and Matrix classes. Version 3.7.0 - May 10 2015, by I. Volobouev * Better numerical derivative calculation in InterpolatedDistribution1D.cc. * Added class "LocalMultiFilter1D" for fast generation of filters which correspond to each orthogonal polynomial separately. * Added a function calculating the area of n-dimensional sphere. * Added I/O capabilities to the RadialProfileND class. * Added class "LogRatioTransform1D". * Added utility function "multiFill1DHistoWithCDFWeights" (header file histoUtils.hh). * Avoiding underflow of the incomplete gamma in "convertToSphericalRandom". Version 3.6.0 - April 6 2015, by I. Volobouev * Fixed Marsaglia's code calculating the Anderson-Darling statistics (it was breaking down for large values of z). * Added high-level driver function "simpleVariableBandwidthSmooth1D" to automatically build the pilot estimate for "variableBandwidthSmooth1D". * Swithched to log of sigma as Minuit parameter in "minuitFitJohnsonCurves" instead of sigma (linear sigma would sometimes break the fit when Minuit would come up with 0 or negative trial value for it). * Extended "MinuitDensityFitFcn1D" class so that it could be used to fit non-uniformly binned histograms. * Adapted "minuitFitJohnsonCurves" function so that it could be used with histograms templated upon DualHistoAxis. * Added functions "fillArrayCentersPreservingAreas" and "canFillArrayCentersPreservingAreas". * Implemented an interface to monotonous coordinate transformations with the intent of using them in constructing densities. Implemented a number of transforms. * Implemented class "TransformedDistribution1D". * Added class "VerticallyInterpolatedDistro1D1P". * Added utility function "fill1DHistoWithCDFWeights". Version 3.5.0 - February 21 2015, by I. Volobouev * Added "symPDEigenInv" method to the Matrix class. * Added "variableCount" method to unfolding bandwidth scanner classes. * Increased the tolerance parameters in JohnsonSu::initialize and in JohnsonFit constructor. * Bug fix in function "fillHistoFromText". Version 3.4.4 - January 13 2015, by I. Volobouev * Corrected handling of some "assert" statements so that the code compiles correctly with the -DNDEBUG option. Version 3.4.3 - January 5 2015, by I. Volobouev * Implemented class MirroredGauss1D. * Added method "getOracleData" to class UnfoldingBandwidthScanner1D. * Bug fix in FoldingFilter1DBuilder::makeOrthoPoly. Version 3.4.2 - December 15 2014, by I. Volobouev * Implemented InterpolatedDistro1D1P class. Version 3.4.1 - November 07 2014, by I. Volobouev * Implemented "divideTransforms" function for deconvolutions. * Implemented the Moyal distribution. * Added "fillHistoFromText" utility function. * Added "apply_lorpe_1d" example. Version 3.4.0 - October 01 2014, by I. Volobouev * Implemented Hadamard product and Hadamard ratio for matrices. * Bug fix in the "solve_linear_system" lapack interface function. Version 3.3.1 - August 08 2014, by I. Volobouev * Terminate iterative refinement of the unfolding error propagation matrix early in case the solution does not seem to improve. Version 3.3.0 - August 05 2014, by I. Volobouev * Added correction factors to the determination of the number of fitted parameters in various unfolding procedures. Version 3.2.0 - July 25 2014, by I. Volobouev * Added "gaussianResponseMatrix" function for non-uniform binning. * Added Pareto distribution. * Implemented EMS unfolding with sparse matrices. * Added methods "getObservedShape" and "getUnfoldedShape" to the AbsUnfoldND class. * Bug fix in the assignment operator of ProductDistributionND class. Made class ProductDistributionND persistent. * Bug fix in the error propagation for unfolding, in the code which takes into account the extra normalization constant. * Added "productResponseMatrix" function to assist in making sparse response matrices. * Bug fix in the factory constructor of the Cauchy1D class. * Completed implementation of the "RatioOfNormals" class. Version 3.1.0 - June 29 2014, by I. Volobouev * Improved (again) random number generator for the 1-d Gaussian distribution. Something about expectation values of normalized Hermite polynomials over random numbers made by this generator is still not fully understood. The standard deviation of these expectations decreases with the polynomial order (while it should stay constant). It is possible that the numbers of points used are simply insufficient to sample the tails correctly. * Implemented smoothed expectation-maximization (a.k.a. D'Agostini) unfolding for 1-d distributions in classes SmoothedEMUnfold1D and MultiscaleEMUnfold1D. In certain usage scenarios, MultiscaleEMUnfold1D can be more efficient than SmoothedEMUnfold1D. * Implemented smoothed expectation-maximization unfolding for multivariate distributions in a class SmoothedEMUnfoldND. * Added class "UnfoldingBandwidthScanner1D" to study 1-d unfolding behavior as a function of filter bandwidth. * Added class "UnfoldingBandwidthScannerND" to study multivariate unfolding behavior as a function of provided bandwidth values. * Added DummyLocalPolyFilter1D class useful when a filter is needed which does not smooth anything. * Added function "poissonLogLikelihood" (header file npstat/stat/arrayStats.hh). * Added function "pooledDiscreteTabulated1D" (header file npstat/stat/DiscreteDistributions1D.hh). * Implemented class UGaussConvolution1D (convolution of uniform distribution with a Gaussian). * Implemented gamma distribution (class Gamma1D). * Defined interface for comparing binned distributions, AbsBinnedComparison1D. * Implemented several comparison classes for comparing binned distributions: PearsonsChiSquared, BinnedKSTest1D, BinnedADTest1D. Class BinnedKSTest1D pulled in dependence on the "kstest" package. * Made classes LocalPolyFilter1D, LocalPolyFilterND, and SequentialPolyFilterND persistent. * Added code generating dense filter matrices to LocalPolyFilterND and SequentialPolyFilterND (as needed for unfolding). * Made class MemoizingSymbetaFilterProvider persistent. * Implemented function goldenSectionSearchOnAGrid (header file npstat/nm/goldenSectionSearch.hh). * Implemented function parabolicExtremum (header npstat/nm/MathUtils.hh). * Added class DistributionMix1D (header npstat/stat/DistributionMix1D.hh). * Added interface to solving A*X = B, with matrices X and B, to the Matrix class (method "solveLinearSystems"). * Added "reshape" methods to the ArrayND class. * Added "gaussianResponseMatrix" function. * Added a section on unfolding to the documentation. * Added "ems_unfold_1d" example program. Version 3.0.0 - March 14 2014, by I. Volobouev * Added interface to the LAPACK SVD routines. * Added function "lorpeMise1D" to calculate MISE for arbitrary distributions. * Added FoldingFilter1DBuilder class. * Changed interfaces for several high-level functions to use FoldingFilter1DBuilder. The major version number got bumped up. * Split DensityScan1D.hh away from AbsDistribution1D.hh. Version 2.7.0 - March 10 2014, by I. Volobouev * Added code to optimize operations with diagonal matrices. * Added discretizedDistance.hh file for simple L1 and L2 distance calculations with numeric arrays. * Added base class for future unfolding code. * The "reset" method of the Matrix class was renamed into "uninitialize" in order to be consistent with ArrayND. * Added function "multinomialCovariance1D". * Added "polyTimesWeight" method to the OrthoPoly1D class. * Added methods "TtimesThis" and "timesT" to the Matrix class. These methods are more efficient than transpose followed by multiplication. Version 2.6.0 - January 30 2014, by I. Volobouev * Added function "lorpeBackgroundCVDensity1D" which linearizes calculation of the cross validation likelihood in semiparametric fits. Argument "linearizeCrossValidation" was added to MinuitSemiparametricFitFcn1D constructor, "lorpeBackground1D" function, etc. * Added the ability to build filters with center point removed to classes WeightTableFilter1DBuilder and StretchingFilter1DBuilder. The function "symbetaLOrPEFilter1D" now has an appropriate switch. * Added "removeRowAndColumn" method to the Matrix class. * Added CircularBuffer class. * Added various plugin bandwidth functions which work with non-integer polynomial degrees. * Switched to the Legendre polynomial basis for calculating all 1-d orthogonal polynomials (instead of monomial basis). * Added MemoizingSymbetaFilterProvider class. * Added "operator+=" method to the MultivariateSumAccumulator class. * Simplified implementation of the PolyFilterCollection1D class. File PolyFilterCollection1D.icc is removed. * Added "RatioOfNormals" 1-d distribution function. Only the density is currently implemented but not the CDF. * Added ExpMapper1d class. Version 2.5.0 - October 15 2013, by I. Volobouev * Added "getFilterMatrix" method to the LocalPolyFilter1D class. * Added "genEigen" method to the Matrix class (for determination of eigenvalues and eigenvectors of general real matrices). * Refactored the LAPACK interface so that interface functions to floats are automatically generated from interface functions to doubles. See the comment at the end of the "lapack_interface.icc" file for the shell commands to do this. Version 2.4.0 - October 6 2013, by I. Volobouev * Added functions "lorpeBackground1D", "lorpeBgCVPseudoLogli1D", and "lorpeBgLogli1D". * Added minuit interface classes "MinuitLOrPEBgCVFcn1D" and "MinuitSemiparametricFitFcn1D". * Added "ScalableDensityConstructor1D" class for use with Minuit interface functions. * Added classes AbsSymbetaFilterProvider and SymbetaPolyCollection1D. Version 2.3.0 - October 1 2013, by I. Volobouev * Allowed point dimensionality to be larger than the histogram dimensionality in the "empiricalCopulaHisto" function. * Added "keepAllFilters" method to AbsFilter1DBuilder and all derived classes. * Implemented exclusion regions for WeightTableFilter1DBuilder and StretchingFilter1DBuilder. * "symbetaLOrPEFilter1D" function (in the header LocalPolyFilter1D.hh) is updated to take the exclusion mask argument. * Added "continuousDegreeTaper" function which can do something meaningful with the continuous LOrPE degree parameter. Version 2.2.0 - June 30 2013, by I. Volobouev * Added classes DiscreteBernsteinPoly1D and BernsteinFilter1DBuilder. * Added classes DiscreteBeta1D and BetaFilter1DBuilder. * Added BifurcatedGauss1D class to model Gaussian-like distributions with different sigmas on the left and right sides. * Added virtual destructors to the classes declared in the Filter1DBuilders.hh header. * Added a method to the Matrix template to calculate Frobenius norm. * Added methods to the Matrix template to calculate row and column sums. * Added "directSum" method to the Matrix template. * Added constructor from a subrange of another matrix to the Matrix template. * Added code to the LocalPolyFilter1D class that generates a doubly stochastic filter out of an arbitrary filter. * Added "npstat/nm/definiteIntegrals.hh" header and corresponding .cc file for various infrequently used integrals. * Added "betaKernelsBandwidth" function. Version 2.1.0 - June 20 2013, by I. Volobouev * Fixed couple problems which showed up in the robust regression code due to compiler update. * Fixed CensoredQuantileRegressionOnKDTree::process method (needed this-> dereference for some member). Version 2.0.0 - June 15 2013, by I. Volobouev * Updated to use "Geners" version 1.3.0. A few interfaces were changed (API for the string archives was removed because Geners own string archive facilities are now adequate) so the major version number was bumped up. Version 1.6.0 - June 12 2013, by I. Volobouev * Updated some documentation. * Updated fitCompositeJohnson.icc to use simplified histogram constructors. * Bug fix in the "minuitLocalQuantileRegression1D" function. * Changed the "quantileBinFromCdf" function to use unsigned long type for array indices. * Added "weightedLocalQuantileRegression1D" function (namespace npsi) for local regression with single predictor on weighted points. Version 1.5.0 - May 23 2013, by I. Volobouev * Added interfaces to LAPACK routines DSYEVD, DSYEVR, and corresponding single precision versions. * Added the "symPSDefEffectiveRank" method to the Matrix class for calculating effective ranks of symmetric positive semidefinite matrices. * Added converting constructor and assignment operator to the Matrix class. * Run the Gram-Schmidt procedure twice when orthogonal polynomials are derived in order to improve orthogonality. Version 1.4.0 - May 20 2013, by I. Volobouev * Added the "append" method to the AbsNtuple class. Version 1.3.0 - May 10 2013, by I. Volobouev * Added the code for Hermite polynomial series. * Improved random number generator for the 1-d Gaussian distribution. * Added a framework for discrete 1-d distributions as well as two concrete distribution classes (Poisson1D, DiscreteTabulated1D). * Added functions "readCompressedStringArchiveExt" and "writeCompressedStringArchiveExt" which can read/write either compressed or uncompressed string archives, distinguished by file extension. Version 1.2.1 - March 22 2013, by I. Volobouev * Improved CmdLine.hh in the "examples/C++" directory. * Added class QuantileTable1D. * Added classes LeftCensoredDistribution and RightCensoredDistribution. Version 1.2.0 - March 13 2013, by I. Volobouev * Added convenience "fill" methods to work with the ntuples which have small number of columns (up to 10). * Fixed a bug in AbsRandomGenerator for univariate generators making multivariate points. * Added LOrPEMarginalSmoother class. Version 1.1.1 - March 11 2013, by I. Volobouev * Added utility function "symbetaLOrPEFilter1D" which creates 1-d LOrPE filters using kernels from the symmetric beta family (and the Gaussian). * Added high level driver function "lorpeSmooth1D". * Allowed variables with zero variances for calculation of correlation coefficients in "MultivariateSumsqAccumulator". Such variables will have zero correlation coefficients with all other variables. * Added rebinning constructor to the HistoND class. Version 1.1.0 - March 8 2013, by I. Volobouev * Changed NUHistoAxis::fltBinNumber method to produce correct results with interpolation degree 0. It is not yet obvious which method would work best for higher interpolation degrees. * Added functions for converting between StringArchive and python bytearray. They have been placed in a new header: wrap/stringArchiveToBinary.hh. * Added methods "exportMemSlice" and "importMemSlice" to ArrayND. These methods allow for filling array slices from unstructured memory buffers and for exporting array slices to such memory buffers. * Added "simpleColumnNames" function (header file AbsNtuple.hh) to generate trivial column names when ntuple column names are not important. * Added functions "neymanPearsonWindow1D" and "signalToBgMaximum1D". They are declared in a new header npstat/neymanPearsonWindow1D.hh. Version 1.0.5 - December 17 2012, by I. Volobouev * Flush string archives before writing them out in stringArchiveIO.cc. * Added class TruncatedDistribution1D. Version 1.0.4 - November 14 2012, by I. Volobouev * Added utilities for reading/writing Geners string archives to files. * Added BinSummary class. * Doxygen documentation improved. Every header file in stat, nm, rng, and interfaces now has a brief description. Version 1.0.3 - September 27 2012, by I. Volobouev * Fixed some bugs related to moving StorableMultivariateFunctor code from "nm" to "stat". Version 1.0.2 - August 6 2012, by I. Volobouev * Added converting copy constructor to the "LinInterpolatedTableND" class. * Added StorableMultivariateFunctor class (together with the corresponding reader class). * Added StorableInterpolationFunctor class which inherits from the above and can be used with interpolation tables. * Added StorableHistoNDFunctor class which inherits from StorableMultivariateFunctor and can be used to interpolate histogram bins. * Added "transpose" method to HistoND class. * Created DualAxis class. * Created DualHistoAxis class. * Added conversion functions between histogram and grid axes. * Added "mergeTwoHistos" function for smooth merging of two histograms. * Added "ProductSymmetricBetaNDCdf" functor to be used as weight in merging histograms. * Added CoordinateSelector class. Version 1.0.1 - June 29 2012, by I. Volobouev * Implemented class LinInterpolatedTableND with related supporting code. Index: trunk/tests/Makefile =================================================================== --- trunk/tests/Makefile (revision 741) +++ trunk/tests/Makefile (revision 742) @@ -1,124 +1,124 @@ # Set the following two variables correctly UTESTPP_DIR = /home/igv/Code/UnitTest++ FFTW_LIBDIR = /usr/local/lib OFILES_COMMON = test_utils.o OTESTS_FAST = test_ArrayND.o test_OrthoPoly1D.o test_OrthoPolyND.o \ test_Functors.o test_interpolation.o test_LocalPolyFilterND.o \ test_Distributions1D.o test_InterpolatedDistribution1D.o \ test_EmpiricalCopula.o test_BoxNDScanner.o \ test_QuadraticOrthoPolyND.o test_LocalLogisticRegression.o \ test_randomPermutation.o test_gegenbauerSeriesSum.o \ test_GaussHermiteQuadrature.o test_SpecialFunctions.o \ test_GaussLegendreQuadrature.o test_WeightedStatAccumulator.o \ test_RandomSequenceRepeater.o test_logLikelihoodPeak.o \ test_legendreSeriesSum.o test_solveCubic.o test_HistoND.o \ test_SampleAccumulator.o test_GridInterpolatedDistribution.o \ test_GridAxis.o test_Matrix.o test_StatUtils.o \ test_amiseOptimalBandwidth.o test_OrderedPointND.o \ test_CompositeDistributionND.o test_CompositeBuilder.o \ test_ArchiveIO.o test_Ntuples.o test_LocalQuantileRegression.o \ test_rectangleQuadrature.o test_RegularSampler1D.o \ test_Instantiations.o test_chebyshevSeriesSum.o \ test_dawsonIntegral.o test_NMCombinationSequencer.o \ test_CrossCovarianceAccumulator.o test_StatAccumulatorPair.o \ test_findPeak2D.o test_StatAccumulatorArr.o test_SeriesCGF1D.o \ test_goldenSectionSearch.o test_empiricalQuantile.o \ test_StatAccumulator.o test_binomialCoefficient.o \ test_CircularMapper1d.o test_LinInterpolatedTableND.o \ test_StorableMultivariateFunctor.o test_NeymanPearsonWindow1D.o \ test_hermiteSeriesSum.o test_hermiteSeriesSumPhys.o \ test_DiscreteDistributions1D.o test_LocalPolyFilter1D.o \ test_definiteIntegrals.o test_ExpMapper1d.o test_MathUtils.o \ test_MemoizingSymbetaFilterProvider.o test_ResponseMatrix.o \ test_gaussianResponseMatrix.o test_RectangleQuadrature1D.o \ WeightedTestAccumulator.o TestAccumulator.o test_histoUtils.o \ test_LocalMultiFilter1D.o test_InterpolatedCompositeBuilder.o \ test_DensityAveScanND.o test_Distro1DBuilder.o \ test_findRootNewtonRaphson.o test_cumulantUncertainties.o \ test_ClassicalOrthoPolys1D.o test_ContOrthoPoly1D.o \ test_findRootUsingBisections.o test_orthoPoly1DVProducts.o \ test_truncatedInverseSqrt.o test_FejerQuadrature.o \ test_matrixIndexPairs.o test_KDE1DKernel.o test_Interval.o \ test_likelihoodStatisticCumulants.o test_Uncertainties.o \ test_DistributionsND.o EdgeworthSeries1DOld.o test_OSDE1D.o \ test_SequentialPolyFilterND.o test_HeatEq1DNeumannBoundary.o \ test_ScalableClassicalOrthoPoly1D.o LogLogQuadratic1D.o OTESTS_COMPREHENSIVE = test_ArrayND_cmh.o test_DiscreteBernstein.o \ test_KDTree.o test_rescanArray_cmh.o test_lorpeMise1D.o # OPTIMIZE = -std=c++11 -g -ggdb -O0 OPTIMIZE = -std=c++11 -O3 INCLUDES = -I. -I.. -I$(UTESTPP_DIR)/src -I$(FFTW_LIBDIR)/../include CPPFLAGS = $(OPTIMIZE) $(INCLUDES) -Wall -W -Werror -Wno-unused-parameter LIBS = -L../npstat/.libs -lnpstat -L$(UTESTPP_DIR) -lUnitTest++ \ -L$(FFTW_LIBDIR) -lfftw3 -llapack -lblas \ -L/usr/local/lib -lrk -lgeners -lbz2 -lz -lkstest -lm %.o : %.cc g++ -c $(CPPFLAGS) -fPIC -MD $< -o $@ @sed -i 's,\($*\.o\)[:]*\(.*\),$@: $$\(wildcard\2\)\n\1:\2,g' $*.d # Change the "all" target below to suit your development needs. # Useful targets which can be included are: # # fast fast_run # comprehensive comprehensive_run # regression regression_run # # The tests inside the "fast" target are a subset of tests from the # "comprehensive" target all: fast fast_run # all: fast regression fast_run regression_run # all: comprehensive comprehensive_run regression regression_run fast: test_main.o $(OTESTS_FAST) $(OFILES_COMMON) rm -f $@ g++ $(OPTIMIZE) -fPIC -o $@ $^ $(LIBS) fast_run: fast ./fast comprehensive: test_main.o $(OTESTS_FAST) \ $(OTESTS_COMPREHENSIVE) $(OFILES_COMMON) rm -f $@ g++ $(OPTIMIZE) -fPIC -o $@ $^ $(LIBS) comprehensive_run: comprehensive ./comprehensive regression: regression_run: regression @ echo Running regression tests. PROGRAMS = kdtree_speed.cc quantileBinFromCdf.cc incompleteGamma.cc \ convertToSpherical.cc printPermutations.cc showIOTraits.cc \ hugeNtuple.cc hugeNtupleRead.cc histoStdev.cc effectiveDim.cc \ buildInterpolatedCheck.cc legendreRoots.cc jacobiPolyStats.cc \ - jacobiVProducts.cc jacobiEpsTest.cc gauss2DRandom.cc + jacobiVProducts.cc jacobiEpsTest.cc gauss2DRandom.cc besselK.cc PROGRAMS += cpp11Random.cc PROGRAMS += failingTest.cc BINARIES = $(PROGRAMS:.cc=) binaries: $(BINARIES) $(BINARIES): % : %.o $(OFILES_COMMON); g++ $(OPTIMIZE) -fPIC -o $@ $^ $(LIBS) clean: rm -f fast comprehensive $(BINARIES) \ *.out core.* *.o *.d *~ *.gsbd *.gsbmf -include test_main.d -include $(OFILES_COMMON:.o=.d) -include $(OTESTS_FAST:.o=.d) -include $(OTESTS_COMPREHENSIVE:.o=.d) -include $(PROGRAMS:.cc=.d) Index: trunk/tests/besselK.cc =================================================================== --- trunk/tests/besselK.cc (revision 0) +++ trunk/tests/besselK.cc (revision 742) @@ -0,0 +1,35 @@ +#include + +#include "examples/C++/CmdLine.hh" +#include "npstat/nm/SpecialFunctions.hh" + +using namespace npstat; +using namespace std; + +int main(int argc, char const* argv[]) +{ + CmdLine cmdline(argc, argv); + if (argc == 1) + { + cout << "Usage: " << cmdline.progname() << " nu x" << endl; + return 0; + } + + double nu, x; + + try { + cmdline.optend(); + if (cmdline.argc() != 2) + throw CmdLineError("wrong number of command line arguments"); + cmdline >> nu >> x; + } + catch (const CmdLineError& e) { + std::cerr << "Error in " << cmdline.progname() << ": " + << e.str() << std::endl; + return 1; + } + + cout.precision(17); + cout << besselK(nu, x) << endl; + return 0; +} Index: trunk/tests/test_ContOrthoPoly1D.cc =================================================================== --- trunk/tests/test_ContOrthoPoly1D.cc (revision 741) +++ trunk/tests/test_ContOrthoPoly1D.cc (revision 742) @@ -1,399 +1,399 @@ #include "UnitTest++.h" #include "test_utils.hh" #include "npstat/nm/ContOrthoPoly1D.hh" #include "npstat/nm/ClassicalOrthoPolys1D.hh" #include "npstat/nm/FejerQuadrature.hh" #include "npstat/nm/StorablePolySeries1D.hh" #include "npstat/rng/MersenneTwister.hh" #include "npstat/stat/Distributions1D.hh" using namespace npstat; using namespace std; inline static int kdelta(const unsigned i, const unsigned j) { return i == j ? 1 : 0; } namespace { TEST(ContOrthoPoly1D_orthonormalization) { const double eps = 1.0e-15; const OrthoPolyMethod method[] = {OPOLY_STIELTJES, OPOLY_LANCZOS}; const unsigned nMethods = sizeof(method)/sizeof(method[0]); const unsigned npoints = 64; const unsigned maxdeg = 10; const unsigned ntries = 10; std::vector points(2U*npoints); std::vector measure; measure.reserve(npoints); std::vector coeffs(maxdeg + 1); MersenneTwister rng; for (unsigned imeth=0; imeth(kdelta(i, j)), d, eps); } measure.clear(); for (unsigned i=0; i(kdelta(i, 0)), coeffs[i], eps); } } #ifdef USE_LAPACK_QUAD TEST(ContOrthoPoly1D_orthonormalization_quad_STIELTJES) { const lapack_double eps = FLT128_EPSILON*100.0; const OrthoPolyMethod method = OPOLY_STIELTJES; const unsigned npoints = 64; const unsigned maxdeg = 10; const unsigned ntries = 5; const double xmin = -1.0; const double xmax = 1.0; const double range = xmax - xmin; std::vector points(2U*npoints); std::vector measure; measure.reserve(npoints); MersenneTwister rng; for (unsigned itry=0; itry(kdelta(i, j)), d, eps); } } } TEST(ContOrthoPoly1D_orthonormalization_quad_LANCZOS) { const lapack_double eps = FLT128_EPSILON*100.0; const OrthoPolyMethod method = OPOLY_LANCZOS; const unsigned npoints = 64; const unsigned maxdeg = 10; const unsigned ntries = 5; const double xmin = -1.0; const double xmax = 1.0; const double range = xmax - xmin; std::vector points(2U*npoints); std::vector measure; measure.reserve(npoints); MersenneTwister rng; for (unsigned itry=0; itry(kdelta(i, j)), d, eps); } } } #endif TEST(ContOrthoPoly1D_weightCoeffs) { const double eps = 1.0e-7; const unsigned maxdeg = 10; const unsigned ntries = 10; const unsigned npoints = maxdeg + 1U; std::vector points(2U*npoints); std::vector measure; measure.reserve(npoints); std::vector coeffs(maxdeg + 1); MersenneTwister rng; for (unsigned itry=0; itry poly2( poly.makeStorablePolySeries(xmin, xmax)); for (unsigned deg=0; deg<=maxdeg; ++deg) for (unsigned ipow=0; ipow<4; ++ipow) { const double i1 = poly.integratePoly(deg, ipow, xmin, xmax); const double i2 = poly2->integratePoly(deg, ipow); if (fabs(i2) > 1.0e-10) CHECK_CLOSE(1.0, i1/i2, 1.0e-10); } for (unsigned i=0; iseries(&coeffs[0], maxdeg, x); CHECK_CLOSE(w, fvalue, 1000*eps); - CHECK_CLOSE(fvalue, f2, eps); + CHECK_CLOSE(fvalue, f2, 10*eps); } } } TEST(ContOrthoPoly1D_cov8) { const double eps = 1.0e-12; MersenneTwister rng; const unsigned npoints = 64; const unsigned maxdeg = 6; const unsigned ntries = 5; const unsigned degtries = 100; std::vector points(npoints); for (unsigned itry=0; itry points(npoints); Beta1D beta(0.0, 1.0, 2.0, 2.0); GaussLegendreQuadrature glq(nInteg); for (unsigned itry=0; itry& mat = poly.extWeightProductAverages( DensityFunctor1D(beta), glq, maxdeg, 0.0, 1.0); for (unsigned i=0; i<=maxdeg; ++i) for (unsigned j=0; j<=maxdeg; ++j) { const double integ = poly.integrateEWPolyProduct( DensityFunctor1D(beta), i, j, 0.0, 1.0, nInteg); CHECK_CLOSE(integ, mat[i][j], eps); } } } TEST(ContOrthoPoly1D_epsCovariance) { const double eps = 1.0e-15; const unsigned npoints = 64; const unsigned maxdeg = 6; const unsigned ntries = 5; MersenneTwister rng; std::vector points(npoints); for (unsigned itry=0; itry r(nrandom); for (unsigned i=0; i result(maxdeg + 1); poly.sampleAverages(&r[0], r.size(), &result[0], maxdeg); Matrix mat(maxdeg + 1, maxdeg + 1, 0); std::vector acc(maxdeg + 1, 0.0L); std::vector tmp(maxdeg + 1); for (unsigned i=0; i(acc[deg]), result[deg], eps); } // Test "pairwiseSampleAverages" method mat /= nrandom; Matrix averages = poly.sampleProductAverages( &r[0], r.size(), maxdeg); Matrix refmat(mat); CHECK((averages - refmat).maxAbsValue() < eps); } TEST(ContOrthoPoly1D_jacobi_2) { typedef ContOrthoPoly1D::PreciseQuadrature Quadrature; const double eps = 1.0e-13; const unsigned alpha = 3; const unsigned beta = 4; const unsigned maxdeg = 5; const unsigned npoints = 64; JacobiOrthoPoly1D jac(alpha, beta); OrthoPoly1DWeight weight(jac); const unsigned npt = Quadrature::minimalExactRule(2*maxdeg+alpha+beta); Quadrature glq(npt); ContOrthoPoly1D poly(maxdeg, glq.weightedIntegrationPoints(weight, -1.0, 1.0), OPOLY_LANCZOS); for (unsigned i=0; i