Ziggurat algorithm help

[ausairman]

Hi Everyone,

I realise this is a long shot asking here, but if someone knows of a working ziggurat algorithm that would be great. Currently I'm using the one from here:

http://people.sc.fsu.edu/~jburkardt/cpp_src/ziggurat_original/ziggurat_original.html

the function r4_uni_value(), which claims to return a normal distribution with variance of 1 and mean of 0, is actually returning a mean of 0 and a variance of around 13 (I measured this with 100k samples). Given that the variance is way off I basically have to assume the whole thing is broken..

I'd really appreciate it if someone could take a look at that code, or point me to an existing version of the ziggurat algorithm, or any normal random number generator that works. I put that last bit in bold because I've come across so many algorithms that people have touted that either require a degree in statistics to operate or blatantly don't work...

[kev82]

Firstly, if you read the page, r4_uni doesn't give you a normal value, r4_norm does.

However, I think you have a programming error, because the range of r4_uni is 1, and the variance can never be larger than the range squared. I can't see how you can get a value of 13.

Anyway, I tested r4_norm with the following code (adapted from their test program), the histogram looks pretty damn normal to me, and the variance comes out at 1.04.

#include <stdio.h>
#include "zig.h"

int main() {
	unsigned long int seed;

	float fn[128];
	int kn[128];
	float wn[128];

	float value;

	r4_nor_setup ( kn, fn, wn );
	seed = 7658291;

	for ( int i = 0; i <= 1000; i++ ) {
		value = r4_nor ( &seed, kn, fn, wn );
		printf("%f\n", value);
	}

	return 0;
}

[ausairman]

What algorithm are you using? I am using the one from the link above, it doesn't take any arguments. I wrote the wrong equation, but I was testing with the r4_nor_value() function, it just returns rubbish (values are all around +/- 3.5).

[ausairman]

Ok I found the version that you were using and it does seem to work fine. Thanks for the help, I now have a working normal random number generator! I suspect there may simply be a bug in the original version (which I was using), which claims to have better performance because it allows function inlining.

[ausairman]

Speaking of the new version, can you tell me what this line is supposed to do: (line 232 of the source file)

if ( x * x <= y + y );
{
break;
}

This seems like a mistake to me; the code will always break...

[kev82]

Can you link me to the source code file you're using, I can't see that line. I'm using the one on the page you gave above.

http://people.sc.fsu.edu/~jburkardt/cpp_src/ziggurat_original/ziggurat_original.C

[ausairman]

That is the original, that I was trying to use but it didn't seem to work (returning a mean of 3.1 and a variance of 13). So I gave up and started using the new one here:

http://people.sc.fsu.edu/~jburkardt/cpp_src/ziggurat/ziggurat.C

This is the one that has that rather suspect line in it, which was the reason I avoided it initially and went with the original (also because it claims to have better performance due to inlining). I have since just settled for it assuming that the creators had some reason for putting that semicolon there that I just don't understand...

[ausairman]

So I'm now trying it with the original version (the one that you linked).

the following code:

seed = time(0);
zigset(seed);
for (int i = 0; i != 20; i++)
  cout << r4_nor_value() <<endl;

returns this:

3.51426
3.64295
-3.90885
-3.59487
3.67149
3.60543
-3.66418
-3.93646
3.72225
-3.61628
3.47616
3.83132
-3.52794
4.0986
-3.542
3.61096
3.51265
-4.0645
3.62416
3.98983

I don't know what you did, but from where I'm standing the algorithm simply does not work...