| Subject: | Unneeded voodoo-math in Math::OO::Random::Normal |
| Date: | Wed, 15 May 2013 07:39:45 -0500 |
| To: | bug-Math-Random-OO [...] rt.cpan.org |
| From: | Joe Zbiciak <intvnut [...] gmail.com> |
I tried sending this directly to dagolden, but it bounced. So, here it is
as a bug report.
Show quoted text
---------- Forwarded message ----------
From: Joe Zbiciak <intvnut@gmail.com>
Date: Wed, May 15, 2013 at 7:31 AM
Subject: Unneeded voodoo-math in Math::OO::Random::Normal
David,
I've been reading through the implementation of Math::OO::Random::Normal,
and I think I see some mathematical voodoo in there. I want to see if you
agree it's unneeded voodoo that can be removed.
# _ltqnorm on (0,1) ranges from about -38 to 8, so we'll use
# the bottom half and invert it for the top half
my $norm = ( $rnd <= 0.5 ) ? _ltqnorm($rnd) : -( _ltqnorm(1-$rnd) );
Here's what I think is voodoo here: _ltqnorm can return the range -38 to
-8 only for very small random numbers. By very small, I mean smaller than
about 2^-53. You can see that in Peter's analysis of ltqnorm here:
http://home.online.no/~pjacklam/notes/invnorm/
The comment in your code suggests that 1 - $rnd will somehow generate
outputs in the range 8 .. 38, but I don't see how that's possible. Because
$rnd is a double precision float (at least, I presume it is), for numbers
in the range [0.5,1), the least significant bit of the mantissa (also
called the "units of the last place" or "ulp") is 2^-53. The largest value
of $rnd in [0.5,1) is 1.0 - 2^-53, so the smallest value 1-$rnd could have
is 2^-53.
That is, we never get into that long tail of 8 .. 38 this way.
Depending on how rand() is implemented, we also probably don't get into the
long negative take -38 to -8 *either*. At least, common implementations of
rand() I've seen (I can't comment on Perl's—it may be better than this)
take a uniformly distributed N-bit random integer, where N is on the order
of 53 to 64 bits, typecast this double, and then adjust the mantissa. The
smallest random number greater than zero this could get you is 2^-63, if
I'm not mistaken. That'll get you a bit beyond -8. To get to -38, Peter's
analysis indicates you need 2^-1073... which is rather smaller. :-)
So, what that's saying is for practical random number generators, norm()
won't return a value more than about 8 stdev from the mean, which as
Peter's analysis indicates, is less than about 2^-53 in each tail (2^-52)
total.
All this is a long-winded way of saying that next() could have been written
more simply as:
sub next {
my ($self) = @_;
my $rnd = rand() || 1e-254; # can't have zero for normals
return _ltqnorm( $rnd ) * $self->stdev + $self->mean;
}
And even if you don't buy all my math argument, you might look at the $plow
/ $phigh stuff in _ltqnorm itself, which basically does the "1 - $rnd" for
you anyway. If $rnd < 0.02425, _ltqnorm uses it directly. If $rnd >
1-0.02425, _ltqnorm uses 1-$rnd. So on that basis alone, you should either
remove the $phigh branch of _ltqnorm, or your own 1-$rnd. :-)
(In case you're wondering why I even bothered: I'm writing an optimized
_ltqnorm for a DSP, and I was trying to *really* understand what I was
looking at.)
—Joe
--
We sell spatulas, and that's all!
http://spatula-city.org/~im14u2c/
http://spacepatrol.info/
--
We sell spatulas, and that's all!
http://spatula-city.org/~im14u2c/
http://spacepatrol.info/