From: Kent M Pitman
Subject: Re: Inverse error function (erf)?
Date:
Message-ID: <u64to9wmc.fsf@nhplace.com>
"Vebjorn Ljosa" <·······@ljosa.com> writes:
> Does anybody know of a Common Lisp implementation of the inverse of
> erf, the error function?
Oh, here we go again--another person for 'fre' software.
;)
Kent M Pitman wrote:
>>Does anybody know of a Common Lisp implementation of the inverse of
>>erf, the error function?
>
>
> Oh, here we go again--another person for 'fre' software.
>
> ;)
AAAUUUUGGGHHH! :)
Paul
"Vebjorn Ljosa" <·······@ljosa.com> writes:
> Does anybody know of a Common Lisp implementation of the inverse of
> erf, the error function?
You'll have to find an expression of this function. MathWorld
indicates a MacLaurin serie for erf^-1:
http://mathworld.wolfram.com/InverseErf.html
http://mathworld.wolfram.com/MaclaurinSeries.html
So you only need to expand the serie up to the needed precision,
(or express the general term). A good mathematical 1/2 hour exercise.
--
__Pascal Bourguignon__ http://www.informatimago.com/
In a World without Walls and Fences,
who needs Windows and Gates?
+ "Vebjorn Ljosa" <·······@ljosa.com>:
| Does anybody know of a Common Lisp implementation of the inverse of
| erf, the error function?
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?
--
* Harald Hanche-Olsen <URL:http://www.math.ntnu.no/~hanche/>
- Debating gives most of us much more psychological satisfaction
than thinking does: but it deprives us of whatever chance there is
of getting closer to the truth. -- C.P. Snow
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))))))))
Vebjorn Ljosa wrote:
> 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)
Have you compared the output of this function to one of his posted
implementations?
I suspect you will want to replace the E in the constant exponents with
D (or L, in some implementations) to get IEEE double precision values
for these constants, rather than the single precision you are likely to
get with E. You might also declare them as constant, rather than
binding them like variables.
Also, there are Fortran implementations he links to which have higher
quoted precision, such as the Applied Statistics Algorithm 241.
+ ···············@hotmail.com" <············@gmail.com>:
| I suspect you will want to replace the E in the constant exponents with
| D (or L, in some implementations) to get IEEE double precision values
| for these constants, rather than the single precision you are likely to
| get with E.
If the choice is implementation dependent, as this seems to indicate,
maybe you should rather set *read-default-float-format* to
double-float, or optionally long-float depending on the
implementation.
--
* Harald Hanche-Olsen <URL:http://www.math.ntnu.no/~hanche/>
- Debating gives most of us much more psychological satisfaction
than thinking does: but it deprives us of whatever chance there is
of getting closer to the truth. -- C.P. Snow
+ "Vebjorn Ljosa" <·······@ljosa.com>:
| (/ (+ (* (+ (* (+ (* (+ (* (+ (* c0 q)
| c1)
| q)
| c2)
| q)
| c3)
| q)
| c4)
| q)
| c5)
| (1+ (* (+ (* (+ (* (+ (* d0 q)
| d1)
| q)
| d2)
| q)
| d3)
| q)))
Much as I enjoy the visual effect of this, it looks a lot to me like
evaluation of polynomials using Horner's rule. Maybe a macro like the
following would be helpful:
(defmacro horner (var &rest coeffs)
(labels
((hornerize (coeffs)
(cond ((null coeffs) 0) ; not needed, in your case
((null (rest coeffs)) (first coeffs))
(t `(+ ,(first coeffs) (* ,var ,(hornerize (rest coeffs))))))))
(hornerize coeffs)))
Then the above could be written
(/ (horner q c0 c1 c2 c3 c4 c5) (horner q d0 d1 d2 d4 1))
or something like that. (I didn't stop to figure that out carefully.
Caveat coder.)
--
* Harald Hanche-Olsen <URL:http://www.math.ntnu.no/~hanche/>
- Debating gives most of us much more psychological satisfaction
than thinking does: but it deprives us of whatever chance there is
of getting closer to the truth. -- C.P. Snow
>>>>> "Vebjorn" == Vebjorn Ljosa <·······@ljosa.com> writes:
Vebjorn> One question for c.l.l before I ask Peter to add the implementation to
Vebjorn> his web page: Would it make sense to signal an error
Vebjorn> (e.g.,cl:floating-point-overflow) or other condition when the result is
Vebjorn> positive or negative infinity? Or is there a portable way to return
Vebjorn> IEEE 754 inf or -inf? Currently, I'm just returning
To get infinity (or a floating-point overflow), you could (attempt) to
return (* 1d300 1d300) or (* 1d300 -1d300). Or something along those
lines. This depends on your Lisp supporting infinity or at least
letting you control the IEEE 754 flags. And you also need a compiler
that won't fold those constants away.
But you get the idea.
Ray