Biocaml 0.4-dev : Biocaml_math.ltqnorm

let ltqnorm p =

(*

  Modified from python code by David Koppstein. Original comments follow below.

  First version was written in perl, by Peter J. Acklam, 2000年07月19日.

  Second version was ported to python, by Dan Field, 2004年05月03日.

*)

  if (p <= 0.) || (p >= 1.) then

   raise (ValueError ("Argument to ltqnorm " ^ (Float.to_string p) ^

    " must be in open interval (0,1)"))

  else

   (* Coefficients in rational approximations. *)

   let a =

    [|-3.969683028665376e+01; 2.209460984245205e+02;

    -2.759285104469687e+02; 1.383577518672690e+02;

    -3.066479806614716e+01; 2.506628277459239e+00|] in

   let b =

    [|-5.447609879822406e+01; 1.615858368580409e+02;

    -1.556989798598866e+02; 6.680131188771972e+01;

    -1.328068155288572e+01|] in

   let c =

    [|-7.784894002430293e-03; -3.223964580411365e-01;

    -2.400758277161838e+00; -2.549732539343734e+00;

    4.374664141464968e+00; 2.938163982698783e+00|] in

   let d =

    [|7.784695709041462e-03; 3.224671290700398e-01;

    2.445134137142996e+00; 3.754408661907416e+00|] in



   (* Define break-points. *)

   let plow = 0.02425 in

   let phigh = 1. -. plow in

   let f q =

    (((((c.(0)*.q+.c.(1))*.q+.c.(2))*.q+.c.(3))*.q+.c.(4))*.q+.c.(5)) /.

     ((((d.(0)*.q+.d.(1))*.q+.d.(2))*.q+.d.(3))*.q+.1.)

   in



  (* Rational approximation for lower region: *)

  if p < plow then

   let q = sqrt ((-2.) *. (log p)) in

   f q



  (* Rational approximation for upper region: *)

  else if phigh < p then

   let q = sqrt ((-2.) *. (log (1. -. p))) in

   f q



  (* Rational approximation for central region: *)

  else

   let q = p -. 0.5 in

   let r = q *. q in

  (((((a.(0)*.r+.a.(1))*.r+.a.(2))*.r+.a.(3))*.r+.a.(4))*.r+.a.(5))*.q /.

   (((((b.(0)*.r+.b.(1))*.r+.b.(2))*.r+.b.(3))*.r+.b.(4))*.r+.1.)