C-faq q13.20

[Thomas Lumley]

One of the examples in q13.20 in the FAQ
(http://c-faq.com/lib/sd16.html, linked from
http://c-faq.com/lib/gaussian.html) seems to be wrong.

The full code is at that link. The part of the code that actually does
the work is

y = (double) rand() / (RAND_MAX + 1.0); /* 0.0 <= y < 1.0 */
bin = (y < 0.5) ? 0 : 1;
y = fabs(y - 1.0); /* 0.0 < y <= 1.0 */
y = std_dev * sqrt((-2.0) * log(y));

return bin ? (mean + y) : (mean - y);

The code is supposed to generate Gaussian random variables, but with
mean=0 and std_dev=1 the smallest number it can generate is
sqrt(-2*log(0.5)), about -1.177. The distribution is two-humped, and
so not even vaguely Gaussian.

In the interests of making this message more c.l.c-conforming I should
probably also add a coding quibble:
bin = (y < 0.5) ? 0 : 1;
should just be
bin = y<0.5;
-thomas

[Spiros Bousbouras]

Thomas Lumley wrote:

One of the examples in q13.20 in the FAQ
(http://c-faq.com/lib/sd16.html, linked from
http://c-faq.com/lib/gaussian.html) seems to be wrong.

The full code is at that link. The part of the code that actually does
the work is

y = (double) rand() / (RAND_MAX + 1.0); /* 0.0 <= y < 1.0 */
bin = (y < 0.5) ? 0 : 1;
y = fabs(y - 1.0); /* 0.0 < y <= 1.0 */
y = std_dev * sqrt((-2.0) * log(y));

return bin ? (mean + y) : (mean - y);

Your first link doesn't contain any code. The code
at the second link does not contain anywhere fabs().
So I have no idea where you got the code you posted
above.

[Registered User]

Spiros Bousbouras wrote:

Thomas Lumley wrote:

One of the examples in q13.20 in the FAQ
(http://c-faq.com/lib/sd16.html, linked from
http://c-faq.com/lib/gaussian.html) seems to be wrong.

The full code is at that link. The part of the code that actually does
the work is

y = (double) rand() / (RAND_MAX + 1.0); /* 0.0 <= y < 1.0 */
bin = (y < 0.5) ? 0 : 1;
y = fabs(y - 1.0); /* 0.0 < y <= 1.0 */
y = std_dev * sqrt((-2.0) * log(y));

return bin ? (mean + y) : (mean - y);

Your first link doesn't contain any code.

Look again. There's a link to the code he's talking about.

[Spiros Bousbouras]

Registered User wrote:

Spiros Bousbouras wrote:

Thomas Lumley wrote:

One of the examples in q13.20 in the FAQ
(http://c-faq.com/lib/sd16.html, linked from
http://c-faq.com/lib/gaussian.html) seems to be wrong.
>
The full code is at that link. The part of the code that actually does
the work is
>
y = (double) rand() / (RAND_MAX + 1.0); /* 0.0 <= y < 1.0 */
bin = (y < 0.5) ? 0 : 1;
y = fabs(y - 1.0); /* 0.0 < y <= 1.0 */
y = std_dev * sqrt((-2.0) * log(y));
>
return bin ? (mean + y) : (mean - y);

Your first link doesn't contain any code.

Look again. There's a link to the code he's talking about.

I looked again ; I can't find it.

[Registered User]

Spiros Bousbouras wrote:

Registered User wrote:

Spiros Bousbouras wrote:

Thomas Lumley wrote:
One of the examples in q13.20 in the FAQ
(http://c-faq.com/lib/sd16.html, linked from
http://c-faq.com/lib/gaussian.html) seems to be wrong.

The full code is at that link. The part of the code that actually does
the work is

y = (double) rand() / (RAND_MAX + 1.0); /* 0.0 <= y < 1.0 */
bin = (y < 0.5) ? 0 : 1;
y = fabs(y - 1.0); /* 0.0 < y <= 1.0 */
y = std_dev * sqrt((-2.0) * log(y));

return bin ? (mean + y) : (mean - y);
>
Your first link doesn't contain any code.

Look again. There's a link to the code he's talking about.

I looked again ; I can't find it.

Ok, I'll be more precise. The page at <http://c-faq.com/lib/sd16.html>
(the first link) says:

<quote>another technique, using explictly the inverse of the Normal
density function, contributed by ``Luben'' </quote>

The phrase "another technique" has got a hyperlink to
<http://c-faq.com/lib/gaussrand.luben.htmlwhich contains the code the
OP's talking about.

[Ben Bacarisse]

"Thomas Lumley" <th****@drizzle.netwrites:

One of the examples in q13.20 in the FAQ
(http://c-faq.com/lib/sd16.html, linked from
http://c-faq.com/lib/gaussian.html) seems to be wrong.

The full code is at that link. The part of the code that actually does
the work is

y = (double) rand() / (RAND_MAX + 1.0); /* 0.0 <= y < 1.0 */
bin = (y < 0.5) ? 0 : 1;
y = fabs(y - 1.0); /* 0.0 < y <= 1.0 */
y = std_dev * sqrt((-2.0) * log(y));

return bin ? (mean + y) : (mean - y);

Indeed. The above will not produce a distribution that is symmetric
about "mean" because the decision as to which side of the mean we
return is correlated to the size of the uniform random number y. One
needs either another random sample:

return rand() RAND_MAX/2 ? mean + y : mean - y;

or (if the your source in random bits is rather precious) the decision
must be based on a bit or bits in the original sample that are not
correlated to the size of the original y:

int r = rand();
y = (double)(r + 1) / (RAND_MAX + 1);
y = std_dev * sqrt(-2 * log(y));
return r & 1 ? mean + y : mean - y;

I can't see the advantage of the fabs. I think the above allows y to
be as small as (randomly) possible whist remaining 0. Of course,
this last version raises the question of the suitability of using the
bottom bits returned by rand() which was well thrashed out about a
year ago.

I should add that even then I don't think the result is normally
distributed (though it will at least be symmetrical).

--
Ben.

[Ben Bacarisse]

Ben Bacarisse <be********@bsb.me.ukwrites:

"Thomas Lumley" <th****@drizzle.netwrites:

>One of the examples in q13.20 in the FAQ
(http://c-faq.com/lib/sd16.html, linked from
http://c-faq.com/lib/gaussian.html) seems to be wrong.

The full code is at that link. The part of the code that actually does
the work is

y = (double) rand() / (RAND_MAX + 1.0); /* 0.0 <= y < 1.0 */
bin = (y < 0.5) ? 0 : 1;
y = fabs(y - 1.0); /* 0.0 < y <= 1.0 */
y = std_dev * sqrt((-2.0) * log(y));

return bin ? (mean + y) : (mean - y);

Indeed. The above will not produce a distribution that is symmetric
about "mean" because the decision as to which side of the mean we
return is correlated to the size of the uniform random number y.

<other stuff of mine snipped>

I should add that even then I don't think the result is normally
distributed (though it will at least be symmetrical).

Drat! Posted by mistake. After writing, I decided that the original
linked code was so broken (this last phrase of mine is far too week --
the result is *definitely* not Gaussian) that it was best to say
nothing!

Anyway, having said more the nothing, the best thing to say is to
suggest the link be removed lest anyone try to use it. (And don't use
the above unless you want some strange, unspecified, symmetric
distribution or numbers!)

--
Ben.

[CBFalconer]

Ben Bacarisse wrote:

>

.... snip ...

>
or (if the your source in random bits is rather precious) the
decision must be based on a bit or bits in the original sample
that are not correlated to the size of the original y:

int r = rand();
y = (double)(r + 1) / (RAND_MAX + 1);
y = std_dev * sqrt(-2 * log(y));
return r & 1 ? mean + y : mean - y;

I can't see the advantage of the fabs. I think the above allows
y to be as small as (randomly) possible whist remaining 0. Of
course, this last version raises the question of the suitability
of using the bottom bits returned by rand() which was well
thrashed out about a year ago.

I should add that even then I don't think the result is normally
distributed (though it will at least be symmetrical).

Cfaq 13.20 covers it.

13.20: How can I generate random numbers with a normal or Gaussian
distribution?

A: Here is one method, recommended by Knuth and due originally
to Marsaglia:

#include <stdlib.h>
#include <math.h>

double gaussrand()
{
static double V1, V2, S;
static int phase = 0;
double X;

if(phase == 0) {
do {
double U1 = (double)rand() / RAND_MAX;
double U2 = (double)rand() / RAND_MAX;

V1 = 2 * U1 - 1;
V2 = 2 * U2 - 1;
S = V1 * V1 + V2 * V2;
} while(S >= 1 || S == 0);

X = V1 * sqrt(-2 * log(S) / S);
} else
X = V2 * sqrt(-2 * log(S) / S);

phase = 1 - phase;

return X;
}

See the extended versions of this list (see question 20.40)
for other ideas.

References: Knuth Sec. 3.4.1 p. 117; Marsaglia and Bray,
"A Convenient Method for Generating Normal Variables";
Press et al., _Numerical Recipes in C_ Sec. 7.2 pp.
288-290.

--
Chuck F (cbfalconer at maineline dot net)
Available for consulting/temporary embedded and systems.
<http://cbfalconer.home.att.net>

[Ben Bacarisse]

CBFalconer <cb********@yahoo.comwrites:

Ben Bacarisse wrote:

>>

... snip ...

<discussion of code linked to from C FAQ snipped>

>>
I should add that even then I don't think the result is normally
distributed (though it will at least be symmetrical).

Cfaq 13.20 covers it.

Well, yes. That's how I (and probably the OP) got to the faulty
linked code in the first place!

The problem is it covers it then links to "another technique" that is
wrong (or maybe is another technique for something else).

--
Ben.

[Thomas Lumley]

Ben Bacarisse wrote:

"Thomas Lumley" <th****@drizzle.netwrites:

One of the examples in q13.20 in the FAQ
(http://c-faq.com/lib/sd16.html, linked from
http://c-faq.com/lib/gaussian.html) seems to be wrong.

The full code is at that link. The part of the code that actually does
the work is

y = (double) rand() / (RAND_MAX + 1.0); /* 0.0 <= y < 1.0 */
bin = (y < 0.5) ? 0 : 1;
y = fabs(y - 1.0); /* 0.0 < y <= 1.0 */
y = std_dev * sqrt((-2.0) * log(y));

return bin ? (mean + y) : (mean - y);

Indeed. The above will not produce a distribution that is symmetric
about "mean" because the decision as to which side of the mean we
return is correlated to the size of the uniform random number y. One
needs either another random sample:

return rand() RAND_MAX/2 ? mean + y : mean - y;

or (if the your source in random bits is rather precious) the decision
must be based on a bit or bits in the original sample that are not
correlated to the size of the original y:

int r = rand();
y = (double)(r + 1) / (RAND_MAX + 1);
y = std_dev * sqrt(-2 * log(y));
return r & 1 ? mean + y : mean - y;

I can't see the advantage of the fabs.

My guess is that
y = fabs(2*y -1.0)
was intended. This makes bin and y independent and is
basically the same as your first suggestion.

I think the above allows y to
be as small as (randomly) possible whist remaining 0. Of course,
this last version raises the question of the suitability of using the
bottom bits returned by rand() which was well thrashed out about a
year ago.

I should add that even then I don't think the result is normally
distributed (though it will at least be symmetrical).

No, it's not remotely normal. It is bimodal, for a start.
In fact, each tail is a Weibull distribution. It also does not have the

claimed standard deviation.

The basic problem is that the original author was trying to invert
the normal density but should have been trying to invert the
normal cumulative distribution function. Inverting the CDF does
give slightly better results (less discreteness in the extreme tails)
than the two valid techniques in q13.20, but it is slower and quite a
lot more complicated.

Fortunately anyone who does any testing of the code will quickly
find out that it is wrong. The link is still unfortunate.
-thomas

[Thad Smith]

Thomas Lumley wrote:

One of the examples in q13.20 in the FAQ
(http://c-faq.com/lib/sd16.html, linked from
http://c-faq.com/lib/gaussian.html) seems to be wrong.

The full code is at that link. The part of the code that actually does
the work is

y = (double) rand() / (RAND_MAX + 1.0); /* 0.0 <= y < 1.0 */
bin = (y < 0.5) ? 0 : 1;
y = fabs(y - 1.0); /* 0.0 < y <= 1.0 */
y = std_dev * sqrt((-2.0) * log(y));

return bin ? (mean + y) : (mean - y);

The code is supposed to generate Gaussian random variables, but with
mean=0 and std_dev=1 the smallest number it can generate is
sqrt(-2*log(0.5)), about -1.177. The distribution is two-humped, and
so not even vaguely Gaussian.

I am not familiar with the formula for Gaussian generators, but suspect
that y = fabs(y - 1.0) should be y = fabs(y - 0.5). That would make it
single humped.

The two algorithms on http://c-faq.com/lib/gaussian.html seem to not
generate independent points. The code attributed to Abramowitz and
Stegun apparently returns identical points for two consecutive calls!
The code from Knuth has correlations between the call with phase==0 and
the following call since S is the same for both.

--
Thad

[Eric Sosman]

Thad Smith wrote On 11/28/06 01:21,:

>
[...]
The two algorithms on http://c-faq.com/lib/gaussian.html seem to not
generate independent points. The code attributed to Abramowitz and
Stegun apparently returns identical points for two consecutive calls!

Yes, whenever sin(2 * PI * V) == cos(2 * PI * V) ...

The code from Knuth has correlations between the call with phase==0 and
the following call since S is the same for both.

No, the values are independent (if you're willing to call
deterministically-generated pseudo-random numbers "independent").
The fact that S is identical is of no consequence; it's just a
scaling factor that normalizes the independent V1,V2. As the
FAQ suggests, "See the references for more information."

--
Er*********@sun.com