/*************************************************************************\ * * File: betasyminv.c * Environment: ANSI C * Author: Richard Simard * Institution: IRO, UniversitÃ© de MontrÃ©al * Date: 1 September 2005 * * The function * * double BetaSymInv (double alpha, double u) * * computes the inverse of the cumulative distribution function (the * percentage point) of the symmetrical beta distribution with shape * parameters (alpha, alpha) at argument u. * * For 0.05 <= alpha <= 90000, the relative error should be smaller than * 0.5*10^(-14) everywhere. * * For 90000 < alpha <= 100000, the relative error should be smaller than * 0.5*10^(-13) everywhere. * * For 100000 < alpha, the relative error should be smaller than * 10^(-6) everywhere, and smaller than 2*10^(-9) for all u > 10^(-15). * * * Restrictions: alpha > 0, 0 <= u <= 1 * \*************************************************************************/ #include #include #include double lgamma (double); double log1p (double); #define HAVE_LOG1P /* Define if log1p is available in math.h */ #define EPSILON 1.0e-15 /* Tolerance double precision */ #define MAXI 11 /* Max number of Newton's iterations */ #define MAXIB 50 /* Max number of bisection iterations */ #define MAXJ 2000 /* Max number of terms in series */ #define LOG2 0.6931471805599453 /* Ln(2) */ #define LOG4 1.3862943611198906 /* Ln(4) */ #define PI_2 1.5707963267948966 /* Pi / 2 */ #define SQPI_2 0.88622692545275801 /* Sqrt(Pi) / 2 */ #define LOG_SQPI_2 -0.1207822376352453 /* Ln(Sqrt(Pi) / 2) */ #define RAC2 1.41421356237309504880 /* Sqrt(2) */ #define EPSSINGLE 1.0e-5 /* Limited tolerance to increase speed */ #define EPSBETA 1.0e-10 /* < 0.75 sqrt(DBL_EPSILON) */ #define ALIM1 100000.0 /* Limiting alpha for normal approximation */ /**************************************************************************/ static const double InvNormal1P1[7] = { 0.160304955844066229311E2, -0.90784959262960326650E2, 0.18644914861620987391E3, -0.16900142734642382420E3, 0.6545466284794487048E2, -0.864213011587247794E1, 0.1760587821390590 }; static const double InvNormal1Q1[7] = { 0.147806470715138316110E2, -0.91374167024260313396E2, 0.21015790486205317714E3, -0.22210254121855132366E3, 0.10760453916055123830E3, -0.206010730328265443E2, 0.1E1 }; static const double InvNormal1P2[8] = { -0.152389263440726128E-1, 0.3444556924136125216, -0.29344398672542478687E1, 0.11763505705217827302E2, -0.22655292823101104193E2, 0.19121334396580330163E2, -0.5478927619598318769E1, 0.237516689024448000 }; static const double InvNormal1Q2[8] = { -0.108465169602059954E-1, 0.2610628885843078511, -0.24068318104393757995E1, 0.10695129973387014469E2, -0.23716715521596581025E2, 0.24640158943917284883E2, -0.10014376349783070835E2, 0.1E1 }; static const double InvNormal1P3[11] = { 0.56451977709864482298E-4, 0.53504147487893013765E-2, 0.12969550099727352403, 0.10426158549298266122E1, 0.28302677901754489974E1, 0.26255672879448072726E1, 0.20789742630174917228E1, 0.72718806231556811306, 0.66816807711804989575E-1, -0.17791004575111759979E-1, 0.22419563223346345828E-2 }; static const double InvNormal1Q3[9] = { 0.56451699862760651514E-4, 0.53505587067930653953E-2, 0.12986615416911646934, 0.10542932232626491195E1, 0.30379331173522206237E1, 0.37631168536405028901E1, 0.38782858277042011263E1, 0.20372431817412177929E1, 0.1E1 }; double inverse_Normal (double x) /* * Returns the inverse of the cdf of the normal distribution. * Rational approximations giving 16 decimals of precision. * J.M. Blair, C.A. Edwards, J.H. Johnson, "Rational Chebyshev * approximations for the Inverse of the Error Function", in * Mathematics of Computation, Vol. 30, 136, pp 827, (1976) */ { int i, negatif; double z, v, numer, denom; const double XBIG = 1000.0; if (x < 0.0) { fprintf (stderr, "***** ERROR inverse_Normal: u < 0\n"); return DBL_MAX; } if (x > 1.0) { fprintf (stderr, "***** ERROR inverse_Normal: u > 1\n"); return DBL_MAX; } if (x <= 0.0) { fprintf (stderr, "WARNING inverse_Normal: u = 0\n"); return -XBIG; } /* Transform x as argument of InvErf */ z = x; x = 2.0 * x - 1.0; if (x >= 1.0) { fprintf (stderr, "WARNING inverse_Normal: u = 1\n"); return XBIG; } if (x < 0.0) { x = -x; negatif = 1; } else { negatif = 0; } if (x <= 0.75) { v = x * x - 0.5625; numer = denom = 0.0; /* Evaluation of the 2 polynomials by Horner */ for (i = 6; i >= 0; i--) { numer = numer * v + InvNormal1P1[i]; denom = denom * v + InvNormal1Q1[i]; } z = x * numer / denom; } else if (x <= 0.9375) { v = x * x - 0.87890625; numer = denom = 0.0; for (i = 7; i >= 0; i--) { numer = numer * v + InvNormal1P2[i]; denom = denom * v + InvNormal1Q2[i]; } z = x * numer / denom; } else { if (z > 1.0E-1) v = 1.0 / sqrt (-log (1.0 - x)); else v = 1.0 / sqrt (-log (2.0*z)); numer = denom = 0.0; for (i = 10; i >= 0; i--) numer = numer * v + InvNormal1P3[i]; for (i = 8; i >= 0; i--) denom = denom * v + InvNormal1Q3[i]; z = (1.0 / v) * numer / denom; } if (negatif) return -(z * RAC2); else return z * RAC2; } /*=========================================================================*/ static double mylog1p (double x) { /* returns a value equivalent to log (1 + x) accurate also for small x. */ if (fabs(x) > 0.1) { return log (1.0 + x); } else { double term = x; double sum = x; int s = 2; while (fabs(term) > EPSILON*fabs(sum) && s < MAXIB) { term *= -x; sum += term/s; s++; } return sum; } } /*=========================================================================*/ static double inverse1 ( double alpha, /* Shape parameter */ double bu /* u * Beta(alpha, alpha) */ ) /* * This method is used for alpha < 1 and x close to 0. */ { int i, j; double x, xnew; double poc, sum, term; x = pow (bu * alpha, 1.0 / alpha); /* First term of series */ /* If T1/T0 is very small, neglect all terms of series except T0 */ term = alpha * (1.0 - alpha) * x / (1.0 + alpha); if (term < EPSILON) return x; x = bu * alpha / (1.0 + term); xnew = pow (x, 1.0 / alpha); /* Starting point of Newton's iterates */ i = 0; do { ++i; x = xnew; /* Compute the series for F(x) */ poc = 1.0; sum = 1.0 / alpha; j = 1; do { poc *= x * (j - alpha) / j; term = poc / (j + alpha); sum += term; ++j; } while ((term > sum * EPSILON) && (j < MAXJ)); sum *= pow (x, alpha); /* Newton's method */ term = (sum - bu) * pow (x * (1.0 - x), 1.0 - alpha); xnew = x - term; } while ((fabs (term) > EPSILON) && (i <= MAXI)); return xnew; } /*------------------------------------------------------------------------*/ static double inverse2 ( double alpha, /* Shape parameter */ double w /* (0.5 - u)B/pow(4, 1 - alpha) */ ) /* * This method is used for alpha < 1 and x close to 1/2. */ { int i, j; double term, sum, y, ynew, z; double poc; term = (1.0 - alpha) * w * w * 4.0 / 3.0; /* If T1/T0 is very small, neglect all terms of series except T0 */ if (term < EPSILON) return 0.5 - w; ynew = w / (1 + term); /* Starting point of Newton's iterates */ i = 0; do { ++i; y = ynew; z = 4.0 * y * y; /* Compute the series for G(y) */ poc = sum = 1.0; j = 1; do { poc *= z * (j - alpha) / j; term = poc / (2 * j + 1); sum += term; ++j; } while ((term > sum * EPSILON) && (j < MAXJ)); sum *= y; /* Newton's method */ term = (sum - w) * pow (1.0 - z, 1.0 - alpha); ynew = y - term; } while ((fabs (term) > EPSILON) && (i <= MAXI)); return 0.5 - ynew; } /*------------------------------------------------------------------------*/ static double bisect ( double alpha, /* Shape parameter */ double logBua, /* Ln(alpha * u * Beta(alpha, alpha)) */ double a, /* x is presumed in [a, b] */ double b) /* * Bisection method to find a solution x. This method is used for alpha > 1 and * u close to 0. It will very rarely be called. */ { int i, j; double z, sum, term; double x, xprev; if (a >= 0.5 || a > b) { a = 0.0; b = 0.5; } x = 0.5 * (a + b); i = 0; do { ++i; z = -x / (1 - x); /* Compute the series for F(x) */ sum = term = 1.0; j = 1; do { term *= z * (j - alpha) / (j + alpha); sum += term; ++j; } while ((fabs (term/sum) > EPSILON) && (j < MAXJ)); sum *= x; /* Bisection method */ term = log (x * (1.0 - x)); z = logBua - (alpha - 1.0) * term; if (z > log(sum)) a = x; else b = x; xprev = x; x = 0.5 * (a + b); } while ((fabs(xprev - x) > EPSILON) && (i < MAXIB)); return x; } /*------------------------------------------------------------------------*/ static double inverse3 ( double alpha, /* Shape parameter */ double logBua /* Ln(alpha * u * Beta(alpha, alpha)) */ ) /* * This method is used for alpha > 1 and x close to 0. */ { int i, j; double z, x, w, xnew, sum, term; double eps = EPSSINGLE; /* For alpha <= 100000 and u < 1.0/(2.5 + 2.25*sqrt(alpha)), X0 is always to the right of the solution, so Newton is certain to converge. */ const double X0 = 0.497; /* Compute starting point of Newton's iterates */ w = logBua / alpha; x = exp (w); #ifdef HAVE_LOG1P term = (log1p(-x) + logBua) / alpha; #else term = (mylog1p(-x) + logBua) / alpha; #endif z = exp (term); if (z >= 0.25) xnew = X0; else if (z > 1.0e-6) xnew = (1.0 - sqrt(1.0 - 4.0*z)) / 2.0; else xnew = z; i = 0; do { ++i; if (xnew >= 0.5) xnew = X0; x = xnew; sum = log (x * (1.0 - x)); w = logBua - (alpha - 1.0) * sum; if (fabs (w) >= DBL_MAX_EXP * LOG2) { xnew = X0; continue; } w = exp (w); z = -x / (1 - x); /* Compute the series for F(x) */ sum = term = 1.0; j = 1; do { term *= z * (j - alpha) / (j + alpha); sum += term; ++j; } while ((fabs (term/sum) > eps) && (j < MAXJ)); sum *= x; /* Newton's method */ term = (sum - w) / alpha; xnew = x - term; if (fabs(term) < 32.0*EPSSINGLE) eps = EPSILON; } while ((fabs (xnew - x) > EPSILON) && (fabs (xnew - x) > sum * EPSILON) && (i <= MAXI)); /* If Newton has not converged with enough precision, call bisection method. It is very slow, but will be called very rarely. */ if (i >= MAXI && fabs (xnew - x) > 10.0 * EPSILON) return bisect (alpha, logBua, 0.0, 0.5); return xnew; } /*------------------------------------------------------------------------*/ static double inverse4 ( double alpha, /* Shape parameter */ double logBva /* Ln(B) + Ln(1/2 - u) + (alpha - 1)*Ln(4) */ ) /* * This method is used for alpha > 1 and x close to 1/2. */ { int i, j; double term, sum, y, ynew, z; double eps = EPSSINGLE; ynew = exp (logBva); /* Starting point of Newton's iterates */ i = 0; do { ++i; y = ynew; /* Compute the series for G(y) */ z = 4.0 * y * y; term = sum = 1.0; j = 1; do { term *= z * (j + alpha - 0.5) / (0.5 + j); sum += term; ++j; } while ((term > sum * eps) && (j < MAXJ)); sum *= y * (1.0 - z); /* Newton's method */ #ifdef HAVE_LOG1P term = log1p (-z); #else term = mylog1p (-z); #endif term = sum - exp (logBva - (alpha - 1.0) * term); ynew = y - term; if (fabs(term) < 32.0*EPSSINGLE) eps = EPSILON; } while ((fabs (ynew - y) > EPSILON) && (fabs (ynew - y) > sum*EPSILON) && (i <= MAXI)); return 0.5 - ynew; } /*------------------------------------------------------------------------*/ static double belog (double x) /* * Used in PeizerInverse. This is the function * * (1 - x^2 + 2x*Ln(x)) / ((1 - x)^2) */ { const double y = 1.0 - x; double ypow, sum, term; int j; if (x < 1.0E-20) return 1.0; if (x < 0.9 || x > 1.1) return (1.0 - x * x + 2.0 * x * log (x)) / (y * y); if (x == 1.0) return 0.0; /* For x near 1, use a series expansion to avoid loss of precision. */ sum = 0.0; ypow = 1.0; j = 2; do { ypow *= y; term = ypow / (j * (j + 1)); sum += term; j++; } while (fabs (term / sum) > EPSBETA); return 2.0 * sum; } /*-------------------------------------------------------------------------*/ static double PeizerInverse (double alpha, double u) { /* Inverse of the normal approximation of Peizer and Pratt */ double t1, t3, xprev; const double C2 = alpha - 1.0 / 3.0 + 0.025 / alpha; const double z = inverse_Normal (u); double x = 0.5; double y = 1.0 - x; int i = 0; do { i++; t1 = (2.0 * alpha - 5.0 / 6.0) * x * y; t3 = 1.0 - y * belog (2.0 * x) - x * belog (2.0 * y); xprev = x; x = 0.5 + 0.5 * z * sqrt(t1 / t3) / C2; y = 1.0 - x; } while (i <= MAXI && fabs (x - xprev) > EPSBETA); return x; } /*------------------------------------------------------------------------*/ static double CalcB4 (double alpha) { /* Convergent series for Gamma(x + 0.5) / Gamma(x) */ double term = 1.0; double sum = 1.0; int i = 1; while (term > EPSILON*sum) { term *= (i - 1.5)*(i - 1.5) /(i*(alpha + i - 1.5)); sum += term; i++; } return SQPI_2 / sqrt ((alpha - 0.5)*sum); } /*------------------------------------------------------------------------*/ double BetaSymInv (double alpha, double u) /* * Compute the inverse of the symmetrical beta distribution. * Returns a negative value on error, otherwise returns x in [0, 1]. */ { double temp, x, y0; int isUpper; /* True if u > 0.5 */ double B; /* Beta(alpha, alpha) */ double C; /* 4^(alpha-1) * Beta(alpha, alpha) */ double logB; /* Ln(Beta(alpha, alpha)) */ double logC; /* Ln(4^(alpha-1)*Beta(alpha, alpha)) */ double logBua; /* Ln(alpha * u * Beta(alpha, alpha)) */ double logBva; /* Ln(4^(alpha-1)*Beta(alpha,alpha)*(0.5 - u) */ if (alpha <= 0.0) { fprintf (stderr, "***** ERROR BetaSymInv: alpha <= 0\n"); return -1.0; } if (u < 0.0 || u > 1.0) { fprintf (stderr, "***** ERROR BetaSymInv: u is not in [0, 1]\n"); return -2.0; } if (u == 0.0) return 0.0; if (u == 1.0) return 1.0; if (u == 0.5) return 0.5; if (alpha == 1.0) return u; /* alpha = 1 is the uniform law */ if (alpha == 0.5) { /* alpha = 1/2 is the arcsin law */ temp = sin (u * PI_2); return temp * temp; } if (alpha > ALIM1) return PeizerInverse (alpha, u); if (u > 0.5) { u = 1.0 - u; isUpper = 1; } else isUpper = 0; /* Compute Beta(alpha, alpha) and Beta(alpha, alpha)*4^(alpha-1) */ if (alpha <= EPSBETA) { /* For a -> 0, B(a,a) = (2/a)*(1 - 1.645*a^2 + O(a^3)) */ B = 2.0 / alpha; C = B / (4.0*(1.0 - alpha*LOG4)); } else if (alpha <= 1.0) { logB = 2.0 * lgamma(alpha) - lgamma(2.0 * alpha); logC = logB + (alpha - 1.0)*LOG4; B = exp(logB); C = exp(logC); } else if (alpha <= 10.0) { logC = LOG_SQPI_2 + lgamma(alpha) - lgamma(0.5 + alpha); logB = logC - (alpha - 1.0)*LOG4; } else if (alpha <= 200.0) { temp = CalcB4 (alpha); logC = log(temp); logB = logC - (alpha - 1.0)*LOG4; } else { /* Asymptotic series for Gamma(x + 0.5) / (Gamma(x) * Sqrt(x)) */ x = 1.0 / (8.0*alpha); temp = 1.0 + x*(-1.0 + x*(0.5 + x*(2.5 - x*(2.625 + 49.875*x)))); /* This is 4^(alpha - 1)*B(alpha, alpha) */ B = SQPI_2 / (sqrt(alpha) * temp); logC = log(B); logB = logC - (alpha - 1.0)*LOG4; } if (alpha <= 1.0) { /* First term of integrated series around 1/2 */ y0 = C * (0.5 - u); if (y0 > 0.25) x = inverse1 (alpha, B * u); else x = inverse2 (alpha, y0); } else { if (u < 1.0 / (2.5 + 2.25*sqrt(alpha))) { logBua = logB + log (u * alpha); x = inverse3 (alpha, logBua); } else { #ifdef HAVE_LOG1P logBva = logC - LOG2 + log1p (-2.0*u); #else logBva = logC - LOG2 + mylog1p (-2.0*u); #endif x = inverse4 (alpha, logBva); } } if (isUpper) return 1.0 - x; else return x; }