Re: GSL Gaussian: possible bug

[James Theiler <jt@lanl.gov>]

On Thu, 15 Apr 1999, Sven N. Thommesen wrote:

] 
] James,
] 
] I was looking into what it would take to incorporate the GSL distributions
] into Swarm, using your Gaussian as an example. I have a couple of comments,
] and a possible bug report.
] 
] (I hope you take this in the spirit in which it is offered, which is only
] to be helpful -- I don't want to be 'meddling in the affairs of wizards' ...)

Please meddle!

] 
] a) First, the bug report: the next-to-last line in the Ratio method code says
] 
] 	while (x*x > -4.0*log(u));
] 
] My 2nd edition Knuth has this:
] 
] 	while (x*x > -4.0/log(u));		// i.e. division, not multiplication

I don't have Knuth right here, but that's my code, so I'll bet you're
right. Thanks.

In practice, though (at least on my home computer), it was cheaper to use
the polar method and throw one away, so that chunk of code is currently
commented out.  Still, it's better to have it be correct! 

] b) Several places you have the idiom
] 
] 	do {
] 	   u = gsl_rng_uniform(r);
] 	}
] 	while (u==0);
] 
] Isn't that what 'gsl_rng_uniform_pos' does? 

I suppose so, but this avoids a layer of function-call overhead, and one
sees explicitly what it is doing.  To be honest, I didn't know about
gsl_rng_uniform_pos.  If it's implemented as a macro or inline (which I
have lately come to appreciate as more subtle that I once thought), then
it might be a good idea to use it.

] c) You say that you don't use the Polar method to produce two variates at a
] time, since this would require the object to keep a static variable. Ok, I
] can see that it's cleaner if all distributions are state-less, memory-less
] objects. However: (a) heavy users of Gaussian might appreciate not running
] at half speed; and (b) does this mean that the generators, which all have
] static state, 'screw up re-entrant or threaded code' ?

(a) Sure.  That's (yet another) reason for supplying source code, so
one can adjust to taste.  I've been bit by polar gaussians before (give
it the same seed, but get different answers since an odd number of 
random numbers were generated on the last go), but it might not be a bad
idea to package up a static variable version.

(b) The generators all require void pointers to their states, so those
states are not static.  But in practice, I usually put a wrapper around
them which does have a static variable.  Again, a more formal way of
doing that might not be a bad idea.

] d) On the gsl_ran_ugaussian_tail() : you say that if sigma==0, this method
] will NOT produce an ordinary gaussian. Would it be more intuitive to put an
] IF test at the top, and call the normal gaussian in the special case of
] sigma==0?

Good idea.  For that matter, there is probably some value sigma_o, of 
order 1, such that if sigma<sigma_o, it is cheaper to use ordinary
gaussian with rejection.  That begins to fall into the category of 
packaging, and in general I'd prefer to keep the basic routines basic,
but in this case it would make the routine more robust for the users
(also a worthy ideal...)

Again, thanks for the suggestions.
jt


---------------------------------------------
James Theiler                     jt@lanl.gov
MS-D436, NIS-2, LANL        tel: 505/665-5682
Los Alamos, NM 87545        fax: 505/665-4414
----- Space and Remote Sensing Sciences -----