Re: Inverse error function (erf)?
- From: "Vebjorn Ljosa" <vebjorn@xxxxxxxxx>
- Date: 4 Sep 2005 15:02:06 -0700
Harald Hanche-Olsen wrote:
>
> No, but googling for "inverse error function" yields the page
> http://home.online.no/~pjacklam/notes/invnorm/ which describes a good
> algorithm and lists implementations in many languages (including bc
> and Visual Basic!). You should have no difficulty writing your own in
> Common Lisp from the description or translating one of the other
> implementations. Maybe you could contribute it to that page
> afterwards?
I translated Peter Acklam's ltqnorm function to Common Lisp. The
function finds the lower tail quantile for standard normal distribution
function. (This is what I needed the inverse of erf for in the first
place.)
One question for c.l.l before I ask Peter to add the implementation to
his web page: Would it make sense to signal an error
(e.g.,cl:floating-point-overflow) or other condition when the result is
positive or negative infinity? Or is there a portable way to return
IEEE 754 inf or -inf? Currently, I'm just returning
most-positive-double-float or most-negative-double-float.
Thanks,
Vebjorn
(defun ltqnorm (p)
"Lower tail quantile for standard normal distribution function.
This function returns an approximation of the inverse cumulative
standard normal distribution function. I.e., given P, it returns
an approximation to the X satisfying P = Pr{Z <= X} where Z is a
random variable from the standard normal distribution.
The algorithm uses a minimax approximation by rational functions
and the result has a relative error whose absolute value is less
than 1.15e-9.
Reference: http://www.math.uio.no/~jacklam/";
(unless (<= 0 p 1)
(error "P must be a probability, i.e., between 0 and 1."))
;; Coefficients in rational approximations.
(let ((a0 -3.969683028665376e+01) (a1 2.209460984245205e+02)
(a2 -2.759285104469687e+02) (a3 1.383577518672690e+02)
(a4 -3.066479806614716e+01) (a5 2.506628277459239e+00)
(b0 -5.447609879822406e+01) (b1 1.615858368580409e+02)
(b2 -1.556989798598866e+02) (b3 6.680131188771972e+01)
(b4 -1.328068155288572e+01)
(c0 -7.784894002430293e-03) (c1 -3.223964580411365e-01)
(c2 -2.400758277161838e+00) (c3 -2.549732539343734e+00)
(c4 4.374664141464968e+00) (c5 2.938163982698783e+00)
(d0 7.784695709041462e-03) (d1 3.224671290700398e-01)
(d2 2.445134137142996e+00) (d3 3.754408661907416e+00)
(low 0.02425) (high 0.97575))
(cond ((= p 0) most-negative-double-float)
((= p 1) most-positive-double-float)
((< p low)
;; Rational approximation for lower region.
(let ((q (sqrt (- (* 2 (log p))))))
(/ (+ (* (+ (* (+ (* (+ (* (+ (* c0 q)
c1)
q)
c2)
q)
c3)
q)
c4)
q)
c5)
(1+ (* (+ (* (+ (* (+ (* d0 q)
d1)
q)
d2)
q)
d3)
q)))))
((> p high)
;; Rational approximation for upper region.
(let ((q (sqrt (- (* 2 (log (- 1 p)))))))
(- (/ (+ (* (+ (* (+ (* (+ (* (+ (* c0 q) c1)
q)
c2)
q)
c3)
q)
c4)
q)
c5)
(1+ (* (+ (* (+ (* (+ (* d0 q)
d1)
q)
d2)
q)
d3)
q))))))
(t
;; Rational approximation for central region.
(let* ((q (- p 0.5))
(r (* q q)))
(/ (* (+ (* (+ (* (+ (* (+ (* (+ (* a0 r)
a1)
r)
a2)
r)
a3)
r)
a4)
r)
a5)
q)
(1+ (* (+ (* (+ (* (+ (* (+ (* b0 r)
b1)
r)
b2)
r)
b3)
r)
b4)
r))))))))
.
- Follow-Ups:
- Re: Inverse error function (erf)?
- From: Raymond Toy
- Re: Inverse error function (erf)?
- From: Harald Hanche-Olsen
- Re: Inverse error function (erf)?
- From: josephoswaldgg@xxxxxxxxxxx
- Re: Inverse error function (erf)?
- Prev by Date: async interrupts, without-aborts
- Next by Date: Re: GC in Jon's raytracing benchmark
- Previous by thread: async interrupts, without-aborts
- Next by thread: Re: Inverse error function (erf)?
- Index(es):
Relevant Pages
|
