| Re: Looking for the inverse of ERRORF ! [message #7273 is a reply to message #7271] |
Fri, 25 October 1996 00:00  |
mahan
Messages: 17
Registered: October 1994
|
Junior Member |
|
|
David,
I also needed the inverse of the error function. Below your original message
is a listing of my implementation based on Matlab's m-file.
-Steve
____________________________________________________________ _________________
Stephen L. Mahan URL: http://aurora.phys.utk.edu/~mahan
S&E 603 UT Molecular Systems and Voice: 423-974-7850 Home: 423-494-7809
Computational Physics Laboratory Fax: 423-974-7843 Pager: 423-515-3149
Knoxville, TN 37996-1200 E-Mail: smahan@utk.edu
____________________________________________________________ _________________
David van Kuijk (kuijk@mpi.nl) wrote:
: Hi
: For the calculation of confidence intervals it is nice to have a function
: which gives the inverse of the error function. Matlab has this inverse
: function (known as "erfinv"). As far as I know IDL does not have it. Did
: anybody implement it????
: Readya,
: David van Kuijk
: Max-Planck-Institute for Psycholinguistics / Department
: of Language and Speech (University of Nijmegen)
: Snail-mail: P.O. Box 310, 6500 AH Nijmegen, The Netherlands
: E-mail: kuijk@mpi.nl
: Tel: +31 (0)24 3521523
: Fax: +31 (0)24 3521213
function erfinv, y
; Stephen Mahan
; Department of Physics and Astronomy
; September 1996
;
;x=erfinv(y)
l=size(y)
pl=2+l(0)
x=fltarr(l(pl))
a=[0.886226899, -1.645349621, 0.914624893, -0.140543331]
b=[-2.118377725, 1.442710462, -0.329097515, 0.012229801]
c=[-1.970840454,-1.624906493, 3.429567803, 1.641345311]
d=[3.543889200, 1.637067800]
y0=.7
k=where(abs(y) le y0)
if (k(0) ne -1) then begin
z=y(k)*y(k)
x(k)=y(k)*(((a(3)*z+a(2))*z+a(1))*z+a(0))/((((b(3)*z+b(2))*z +b(1))*z+b(0))*z+1)
endif
k=where( (y gt y0) and (y lt 1))
if (k(0) ne -1) then begin
z=sqrt(-alog((1-y(k))/2.))
x(k)=(((c(3)*z+c(2))*z+c(1))*z+c(0))/((d(1)*z+d(0))*z+1)
endif
k=where((y gt -y0) and (y gt -1))
if (k(0) ne -1) then begin
z=sqrt(-alog((1+y(k))/2.))
x(k)=-(((c(3)*z+c(2))*z+c(1))*z+c(0))/((d(1)*z+d(0))*z+1)
endif
x=x-(errorf(x)-y)/(2./sqrt(!dpi) * exp(-x^2))
x=x-(errorf(x)-y)/(2./sqrt(!dpi) * exp(-x^2))
k=where(y eq -1)
if (k(0) ne -1) then x(k)=-k*1./0.
k=where(y eq 1)
if (k(0) ne -1) then x(k)=k*1./0.
k=where(y gt 1)
if (k(0) ne -1) then x(k)=k*0./0.
return,x
end
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
