/* cdf/inverse_normal.c * * Copyright (C) 2002 Przemyslaw Sliwa and Jason H. Stover. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or (at * your option) any later version. * * This program is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA. */ /* * Computes the inverse normal cumulative distribution function * according to the algorithm shown in * * Wichura, M.J. (1988). * Algorithm AS 241: The Percentage Points of the Normal Distribution. * Applied Statistics, 37, 477-484. */ #include #include #include "error.h" int gsl_cdf_inv_normal_e(double p, gsl_cdf_result * result) { double q, r, ret_val; q = p - .5; if (fabs (q) <= .425) { r = .180625 - q * q; ret_val = q * (((((((r * 2509.0809287301226727 + 33430.575583588128105) * r + 67265.770927008700853) * r + 45921.953931549871457) * r + 13731.693765509461125) * r + 1971.5909503065514427) * r + 133.14166789178437745) * r + 3.387132872796366608) / (((((((r * 5226.495278852854561 + 28729.085735721942674) * r + 39307.89580009271061) * r + 21213.794301586595867) * r + 5394.1960214247511077) * r + 687.1870074920579083) * r + 42.313330701600911252) * r + 1.); result->val = ret_val; return GSL_SUCCESS; } else { if (q < 0.) { r = p; } else { r = 1.0 - p; } if (r <= 0.0) { DOMAIN_ERROR(result); } r = sqrt (-log(r)); if (r <= 5.0) { r += -1.6; ret_val = (((((((r * 7.7454501427834140764e-4 + .0227238449892691845833) * r + .24178072517745061177) * r + 1.27045825245236838258) * r + 3.64784832476320460504) * r + 5.7694972214606914055) * r + 4.6303378461565452959) * r + 1.42343711074968357734) / (((((((r * 1.05075007164441684324e-9 + 5.475938084995344946e-4) * r + .0151986665636164571966) * r + .14810397642748007459) * r + .68976733498510000455) * r + 1.6763848301838038494) * r + 2.05319162663775882187) * r + 1.0); } else { r += -5.0; ret_val = (((((((r * 2.01033439929228813265e-7 + 2.71155556874348757815e-5) * r + .0012426609473880784386) * r + .026532189526576123093) * r + .29656057182850489123) * r + 1.7848265399172913358) * r + 5.4637849111641143699) * r + 6.6579046435011037772) / (((((((r * 2.04426310338993978564e-15 + 1.4215117583164458887e-7) * r + 1.8463183175100546818e-5) * r + 7.868691311456132591e-4) * r + .0148753612908506148525) * r + .13692988092273580531) * r + .59983220655588793769) * r + 1.0); } if (q < 0.0) { ret_val = -ret_val; } result->val = ret_val; return GSL_SUCCESS; } result->val = ret_val; return GSL_SUCCESS; } /* * No error code. */ #include "eval.h" double gsl_cdf_inv_normal ( double p ) { EVAL_RESULT ( gsl_cdf_inv_normal_e ( p, &result)); } /* int main(void) */ /* { */ /* double quantile, p; */ /* int ifail; */ /* for(p=1.0E-15; p<1.0; p+=1.0E-5) */ /* { */ /* quantile = gsl_normal_quantile(p, &ifail); */ /* printf("VALUE: %.19f\tIFAIL %d\n", quantile, ifail); */ /* } */ /* return 0; */ /* } */