Maxime Devos <[email protected]> writes:
> Op 04-10-2023 om 18:14 schreef Keith Wright:
>> From: Zelphir Kaltstahl <[email protected]>
>>
>>> my goal is normal distributed floats (leaving aside the finite
>>> nature of the computer and floats).
>> The following is either provably correct up to round-off,
>> or totally stupid.
> It's not stupid at all, but neither is it correct (the basic approach is
> correct, and holds more generally for other probability distributions as
> well, but some important details are off ...).
Not surprising, I just made it up on the fly...
>> First define the cumulative distribution for the normal distribution:
>>
>> $$f(x)= \frac{1}{\sqrt{\pi}} \int_{-\infty}^{x} e^{-x^2/2} dx $$
>
> Instead of sqrt(pi) you need sqrt(2pi).
Seems plausible.
> Also, this is the pdf, not the cdf. For the cdf, you need integrate this
> expression from -infinity to x.
I was using TeX notation. That's exactly what $\int_{-\infty}^{x}$
means.
>> Computing the inverse of the cumulative normal distribution is left as
>> an exercise, because I don't know how, but it seems possible.
> While the pdf is easy to invert (multiply by constant factor, take log,
> multiply by constant factor), resulting in an expression only involving
> constants, multiplication, logarithms (and, depending on simplification,
> subtraction), the cdf isn't.
I was thinking of numerical integration, but too busy to write the code.
I was feeling guilty about being so sloppy, and thinking of writing
it more carefully, but...
> (the basic approach is correct, and holds more generally for other
> probability distributions as well)
I just saw some pictures in my head and thought it must be a general
theorem, but it's so simple and general that it must be "well known".
Somebody sometime must have already written it better than I could;
but I don't know who. It is Theorem (X?) on page (Y?) of book (Z?).
-- Keith
PS: While we are cleaning up details, substitute "modulo" for "/" in:
From: Maxime Devos <[email protected]>
> So, to generate an (approximately) uniform random number on [0,1), you
> can simply do
>
> (define (random-real)
> (exact->inexact (/ (random N) N)))
>
> for a suitably large choice of integer N>0.
A choice that makes this nice (on a 32 bit machine) is $N = 2^{32}$.
