Re: 'erf' function in C
From: Deane Yang (yang2w_at_earthlink.net)
Date: 02/25/04
- Next message: Ben Pfaff: "Re: C dangerous?"
- Previous message: Charles Sullivan: "Re: C dangerous?"
- In reply to: George Marsaglia: "Re: 'erf' function in C"
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
Date: Wed, 25 Feb 2004 03:20:45 GMT
Yes, the paper cited below is the one I used to implement
the inverse cumulative normal function. In that paper
you recommended that such a function be used to generate
random normally distributed numbers, instead of the Box-Muller
method.
1) Would you still recommend this?
2) The approach suggested below seems to be aiming more
for accuracy, rather than speed. On the other hand,
the published algorithm is machine-dependent and therefore
a bit of a headache to maintain.
Is there another fast but portable algorithm with sufficient
accuracy available?
George Marsaglia wrote:
> "Deane Yang" <deane_yang@yahoo.com> wrote in message
> news:36ab6f9c.0402221143.13262bfa@posting.google.com...
>
>>And what about the inverse cumulative normal distribution?
>>I used to implement an algorithm you published,
>>but do not like it because it is machine-dependent.
>>What in your opinion is the best algorithm?
>
>
> With the improvement on my method for
> evaluating cPhi(x) as R(x)*phi(x),
> when the Taylor series for R(x) is about zero,
> (suggested by Bill Daly), the simple C function
>
> double Phi(double x)
> { long double t=0,b=1,s=x,pwr=x;
> int i;
> for(i=2;s!=t;i+=2)
> { b/=(i+1); pwr=pwr*x*x; t=s; s+=pwr*b;}
> return .5+s*exp(-.5*x*x-.91893853320467274178L);
> }
>
> should serve very well for solving the equation
> Phi(X)=U
> for positive X, given U>1/2, by Newton's method:
> Start with an initial estimate
>
> x=sqrt(-1.6*ln(1-U^2));
>
> then repeat
>
> x=x-(Phi(x)-U)/phi(x);
>
> until you get no further changes in x.
>
>
> The paper that you refer to must be where I want
> to generate a normal X by means of solving
> Phi(X)=U, given a random U, uniform in [0,1),
> by getting an integer j, formed from certain
> bits of the exponent part of the floating point
> representation of U. Then use j to access tabled
> values A[j] and B[j] to initalize the Taylor series
> for Phi inverse, using the fractional part of U
> as the argument in the series. It is designed for
> speed in generating a normal variate directly
> as Phi^(-1)(U), but only provides accuracy to 6/7 digits,
> (As I recall, the Taylor series was easy and fast only for the first few
> terms.)
>
> It is described in
> ``Normal (Gaussian) random variables in supercomputers'', (1991)
> Journal of Supercomputing}, v 5, 49--55,
> where the method is particularly suited for parallel computation.
>
> George Marsaglia
>
>
>
- Next message: Ben Pfaff: "Re: C dangerous?"
- Previous message: Charles Sullivan: "Re: C dangerous?"
- In reply to: George Marsaglia: "Re: 'erf' function in C"
- Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
Relevant Pages
|
