diff options
| -rw-r--r-- | arith.c | 24 | ||||
| -rw-r--r-- | arith.h | 1 | ||||
| -rw-r--r-- | eval.c | 1 | ||||
| -rw-r--r-- | txr.1 | 19 |
4 files changed, 45 insertions, 0 deletions
@@ -2725,6 +2725,30 @@ val cum_norm_dist(val arg) } } +/* + * Source: + * Odeh and Evans approximation. + * 1974 + */ +val inv_cum_norm(val arg) +{ + val arg_flo = to_float(lit("inv-cum-norm"), arg); + double p = c_flo(arg_flo); + int is_upper_half = (p >= 0.5); + double r = is_upper_half ? 1 - p : p; + if (r < 1E-20) { + return flo(is_upper_half ? 10 : -10); + } else { + double y = sqrt(-2*log(r)); + double a = ((((4.53642210148e-5*y + 0.0204231210245)*y + + 0.342242088547)*y + 1)*y + 0.322232431088); + double b = ((((0.0038560700634*y + 0.10353775285)*y + + 0.531103462366)*y + 0.588581570495)*y + 0.099348462606); + double z = y - a / b; + return flo(is_upper_half ? z : -z); + } +} + static val rising_product(val m, val n) { val acc; @@ -37,6 +37,7 @@ val in_int_ptr_range(val bignum); ucnum c_unum(val num); val unum(ucnum u); val cum_norm_dist(val x); +val inv_cum_norm(val p); val n_choose_k(val n, val k); val n_perm_k(val n, val k); val divides(val d, val n); @@ -5903,6 +5903,7 @@ void eval_init(void) reg_fun(intern(lit("exp"), user_package), func_n1(expo)); reg_fun(intern(lit("sqrt"), user_package), func_n1(sqroot)); reg_fun(intern(lit("cum-norm-dist"), user_package), func_n1(cum_norm_dist)); + reg_fun(intern(lit("inv-cum-norm"), user_package), func_n1(inv_cum_norm)); reg_fun(intern(lit("n-choose-k"), user_package), func_n2(n_choose_k)); reg_fun(intern(lit("n-perm-k"), user_package), func_n2(n_perm_k)); reg_fun(intern(lit("fixnump"), user_package), func_n1(fixnump)); @@ -34099,6 +34099,25 @@ distribution function: the integral, of the normal distribution function, from negative infinity to the .metn argument . +.coNP Function @ inv-cum-norm +.synb +.mets (inv-cum-norm << argument ) +.syne +.desc +The +.code inv-cum-norm +function calculates an approximate to the inverse of the cumulative +normal distribution function. The argument, a value expected to lie +in the range [0, 1], represents the integral of the normal distribution +function from negative infinity to some domain point +.IR p . +The function calculates the approximate value of +.IR p . +The minimum value returned is -10, and the maximum value returned is 10, +regardless of how closely the argument approaches, respectively, +the 0 or 1 integral endpoints. For values less than zero, or exceeding 1, the +values returned, respectively, are -10 and 10. + .coNP Functions @ n-choose-k and @ n-perm-k .synb .mets (n-choose-k < n << k ) |