Geometrically distributed random numbers on Rabbit 2000.

[The Real Andy]

Here we go again.

So if I was to do this, where would be the best place to ask for
advice considering i know stuff all about math, especially when it
comes to probability and stats.

Name a price in AU dollars and the job is yours....

OH, btw, it has to be based on the Knuth algorithm as that is already
approved by a regulatory body. Source available on request.

a _dot_ pearson _at_ tpg _dot_ com _dot_ au <-- still trying to beat
spam.

 

[The Real Andy]

On Wed, 08 Jun 2005 14:52:25 +1000, The Real Andy
<will_get_back_to_you_on_This> wrote:

Here we go again.

So if I was to do this, where would be the best place to ask for
advice considering i know stuff all about math, especially when it
comes to probability and stats.

Name a price in AU dollars and the job is yours....

OH, btw, it has to be based on the Knuth algorithm as that is already
approved by a regulatory body. Source available on request.

a _dot_ pearson _at_ tpg _dot_ com _dot_ au <-- still trying to beat
spam.

So, me thinks.

Range from 1 - 1000

divide the range by 5. Pick a unifomally distributed number in the new
range. Count how many time it takes before you pick a second uniform
random number that matches the first. This works. The stats say that
it will take an average of 200 iterations to get a match, however in
theroy it could take for ever to get a match...Geometric distribution.

Problem, range after being divided is 1,000,000. This takes on average
440 seconds on the rabbit running at 30MHz.

So me thinks more, i wonder if its possible to somehow converge the
numbers to speed things up.

 

[Jonathan Westhues]

"The Real Andy" <will_get_back_to_you_on_This> wrote in message
news:bi9da15bem8koqutud7nu2sptj8d2dj6oc@4ax.com...

So, me thinks.

Range from 1 - 1000

divide the range by 5. Pick a unifomally distributed number in the new
range. Count how many time it takes before you pick a second uniform
random number that matches the first. This works. The stats say that
it will take an average of 200 iterations to get a match, however in
theroy it could take for ever to get a match...Geometric distribution.

Problem, range after being divided is 1,000,000. This takes on average
440 seconds on the rabbit running at 30MHz.

So me thinks more, i wonder if its possible to somehow converge the
numbers to speed things up.

The Art of Computer Programming volume 2 p. 131 says that:

"A convenient way to generate a variable with this [geometric] distribution
is to set

N <- ceil(ln U / ln (1 - p)) [U is uniformly distributed between 0 and
1]

To check this formula, we observe that ceil(ln U / ln(1 - p)) = n if and
only if n - 1 < ln U / ln(1 - p) <= n, that is, (1 - p)^(n-1) > U >= (1 -
p)^n, and that happens with probability p(1 - p)^(n-1) as required. Note
that ln U can be replaced by -Y, where Y has the exponential distribution
with mean 1."

I don't understand why you would use an algorithm that doesn't execute in
deterministic time. Is there any reason that you couldn't just compute the
quantity above (using fixed-point math of course) and be done with it? I
don't know anything about the Rabbit but even if it doesn't have a hardware
multiplier I don't think you could possibly be looking at more than a few
milliseconds.

Jonathan

 

[Jonathan Kirwan]

On Wed, 08 Jun 2005 18:08:22 +1000, The Real Andy
<will_get_back_to_you_on_This> wrote:

So me thinks more, i wonder if its possible to somehow converge the
numbers to speed things up.

Knuth's algorithm is:

INT( LN( uniform ) / LN( 1 - p ) )

where 'uniform' is a uniform deviate and p is the probability. Knuth
mentions this on page 131 in Seminumerical Algorthms. But you already
know that, I assume, since you mentioned it.

Are you only looking for an algorithm that will iterate geometrically?
If so, Knuth discusses, in exercise 6 on page 133, such an algorithm.

Other than that, I don't have a clue about what you want.

Jon

 

[Fred Bloggs]

Jonathan Kirwan wrote:

Knuth's algorithm is:

INT( LN( uniform ) / LN( 1 - p ) )

where 'uniform' is a uniform deviate and p is the probability. Knuth
mentions this on page 131 in Seminumerical Algorthms. But you already
know that, I assume, since you mentioned it.

Are you only looking for an algorithm that will iterate geometrically?
If so, Knuth discusses, in exercise 6 on page 133, such an algorithm.

Other than that, I don't have a clue about what you want.

Jon

It's obviously the distribution of number of plays for hitting a jackpot
on some kind of gaming machine where the jackpot sequence is fixed and
randomly re-selected after each payout.

 

[The Real Andy]

On Wed, 08 Jun 2005 10:34:48 GMT, Fred Bloggs <nospam@nospam.com>
wrote:

Jonathan Kirwan wrote:


Knuth's algorithm is:

INT( LN( uniform ) / LN( 1 - p ) )

where 'uniform' is a uniform deviate and p is the probability. Knuth
mentions this on page 131 in Seminumerical Algorthms. But you already
know that, I assume, since you mentioned it.

Are you only looking for an algorithm that will iterate geometrically?
If so, Knuth discusses, in exercise 6 on page 133, such an algorithm.

Other than that, I don't have a clue about what you want.

Jon

It's obviously the distribution of number of plays for hitting a jackpot
on some kind of gaming machine where the jackpot sequence is fixed and
randomly re-selected after each payout.

Do I know you mr bloggs?? your surname is not johnson by any chance?

 

[The Real Andy]

On Wed, 08 Jun 2005 09:30:11 GMT, Jonathan Kirwan
<jkirwan@easystreet.com> wrote:

On Wed, 08 Jun 2005 18:08:22 +1000, The Real Andy
will_get_back_to_you_on_This> wrote:

So me thinks more, i wonder if its possible to somehow converge the
numbers to speed things up.

Knuth's algorithm is:

INT( LN( uniform ) / LN( 1 - p ) )

where 'uniform' is a uniform deviate and p is the probability. Knuth
mentions this on page 131 in Seminumerical Algorthms. But you already
know that, I assume, since you mentioned it.

Are you only looking for an algorithm that will iterate geometrically?
If so, Knuth discusses, in exercise 6 on page 133, such an algorithm.

Other than that, I don't have a clue about what you want.

Jon

I have not got that book. I must go buy it.

The algorithm you mention above I had just located about half an hour
ago on the internet. I am currently running tests of it (not on the
rabbit however) to see if it suits my needs

The knuth algorith I have mentioned is for a uniformly distributed
RNG. It was already approved by the government regulator for use in
gaming applications, hence the reason for using it.

Many thanks.

 

[The Real Andy]

On Wed, 8 Jun 2005 04:55:31 -0400, "Jonathan Westhues"
<google-for-my-name@nospam.com> wrote:

"The Real Andy" <will_get_back_to_you_on_This> wrote in message
news:bi9da15bem8koqutud7nu2sptj8d2dj6oc@4ax.com...
So, me thinks.

Range from 1 - 1000

divide the range by 5. Pick a unifomally distributed number in the new
range. Count how many time it takes before you pick a second uniform
random number that matches the first. This works. The stats say that
it will take an average of 200 iterations to get a match, however in
theroy it could take for ever to get a match...Geometric distribution.

Problem, range after being divided is 1,000,000. This takes on average
440 seconds on the rabbit running at 30MHz.

So me thinks more, i wonder if its possible to somehow converge the
numbers to speed things up.

The Art of Computer Programming volume 2 p. 131 says that:

"A convenient way to generate a variable with this [geometric] distribution
is to set

N <- ceil(ln U / ln (1 - p)) [U is uniformly distributed between 0 and
1]

To check this formula, we observe that ceil(ln U / ln(1 - p)) = n if and
only if n - 1 < ln U / ln(1 - p) <= n, that is, (1 - p)^(n-1) > U >= (1 -
p)^n, and that happens with probability p(1 - p)^(n-1) as required. Note
that ln U can be replaced by -Y, where Y has the exponential distribution
with mean 1."

I don't understand why you would use an algorithm that doesn't execute in
deterministic time. Is there any reason that you couldn't just compute the
quantity above (using fixed-point math of course) and be done with it? I
don't know anything about the Rabbit but even if it doesn't have a hardware
multiplier I don't think you could possibly be looking at more than a few
milliseconds.

Jonathan

There is so many variants on the above algorithm is unbelievable. As
stated in another post, i must go and by that book, in fact i will
order it tomorrow.

I will take the algorithm you mention and run some tests over it and
check the results.

Many thanks again.

 

[Jonathan Kirwan]

On Wed, 08 Jun 2005 23:18:30 +1000, The Real Andy
<will_get_back_to_you_on_This> wrote:

snip
The knuth algorith I have mentioned is for a uniformly distributed
RNG. It was already approved by the government regulator for use in
gaming applications, hence the reason for using it.

If you test the algorithm that Westhues and I mention, note that I did
NOT copy it correctly from the book and he did. I had glanced
casually at the marks and misread the bars there. So use the CEIL
function he mentions, though it may work similarly either way (I've
not thought it through.)

Otherwise, the algorithm mentioned in the exercise I pointed out is
quite simple. You imagine a circle centered inside a square and sized
so that it just touches the sides of the square. Then select only
quadrant 1 to work with and ignore the rest. Then imagine two uniform
deviates from 0 to 1, with one of them taken as the x-axis and the
other as the y-axis. Then:

(1) get two random uniform deviates U and V,
(2) if U^2+V^2 >= 1 then goto 1, else X is set to U.

The number of executions of step 2 has the geometric distribution.,
with a minimum of 1 time, an average of 4/PI times, a maximum of
infinity, and a deviation of (4/PI)*(SQRT(1-PI/4).

Get the book.

Jon

 

[Rich Grise]

On Wed, 08 Jun 2005 14:52:25 +1000, The Real Andy wrote:

Here we go again.

So if I was to do this, where would be the best place to ask for
advice considering i know stuff all about math, especially when it
comes to probability and stats.

Name a price in AU dollars and the job is yours....

OH, btw, it has to be based on the Knuth algorithm as that is already
approved by a regulatory body. Source available on request.

a _dot_ pearson _at_ tpg _dot_ com _dot_ au <-- still trying to beat
spam.

So, what is it you're trying to accomplish? A Fourier transform of
pink noise?

Thanks,
Rich

 

[Jonathan Westhues]

"Rich Grise" <richgrise@example.net> wrote in message
newsan.2005.06.08.22.06.37.768913@example.net...

So, what is it you're trying to accomplish? A Fourier transform of
pink noise?

Someone said jackpots, and that would be my bet too.

But the spectrum of a random variable is a different thing than the
distribution; if you let X = ln(U) / ln(1 - p) then I think that X will be
white if U is white.

Jonathan

 

[The Real Andy]

On Wed, 8 Jun 2005 19:13:44 -0400, "Jonathan Westhues"
<google-for-my-name@nospam.com> wrote:

"Rich Grise" <richgrise@example.net> wrote in message
newsan.2005.06.08.22.06.37.768913@example.net...
So, what is it you're trying to accomplish? A Fourier transform of
pink noise?

Someone said jackpots, and that would be my bet too.

But the spectrum of a random variable is a different thing than the
distribution; if you let X = ln(U) / ln(1 - p) then I think that X will be
white if U is white.

Jonathan

Spot on Jonathan. The knuth algoritm works a treat. I do enjoy
learning new stuff!! His book has been ordered. I now declare Knuth as
the God of random numbers!!

Anyhow, after many trials, i have decided to go with a knuth variant,
which is pretty much exactly as described in previous posts. I have
had about 5 schooners (425ml ea) of beer so let me try to explain.
Please excuse if i miss a bracket or 2.

U = uniform distributed random number between 0.0 and 1.0
p = probability of loosing (should it be q?)

say, one has a 1 in 200 chance of winning, so
p = (1.0 - (1/200)

the kntuh variant.

floor( (log(1.0 - U)) / (log(p)) ) + 1.0

So this works a treat with a small range from 1 to 1000 with single
precision.

However, take a range from 1 to 100m for example, then you get a
problem, which I am going to call stepping. I would love to know the
mathematical term for this. It is hard to describe, but as the range
goes onto infinity, it seems to miss numbers, when you plot the
frequency distrbution, it looks as if it has steps.

This problem occours when passing a 16bit random number into the
algorithm. Going to a 30bit (the best knuth rng i have at the moment)
helps this problem. Moving to double precision also helps, but only if
you have the larger random number. So i am guessing that there is two
problems here, one being precision, and the other with the the random
number size.

BTW rich, i had thought of your idea too, it is still in the back of
my mind.

 

[Fred Bloggs]

However, take a range from 1 to 100m for example, then you get a
problem, which I am going to call stepping. I would love to know the
mathematical term for this. It is hard to describe, but as the range
goes onto infinity, it seems to miss numbers, when you plot the
frequency distrbution, it looks as if it has steps.

This problem occours when passing a 16bit random number into the
algorithm. Going to a 30bit (the best knuth rng i have at the moment)
helps this problem. Moving to double precision also helps, but only if
you have the larger random number. So i am guessing that there is two
problems here, one being precision, and the other with the the random
number size.

That's easily taken care by going to a multi-dimensional distibution,
but I'll wait to see what they can extract from Knuth on this.

 

[Jonathan Westhues]

"The Real Andy" <will_get_back_to_you_on_This> wrote in message
news:4i0ga114v7ob3qhn8drajkiajlpfl6jecj@4ax.com...

However, take a range from 1 to 100m for example, then you get a
problem, which I am going to call stepping. I would love to know the
mathematical term for this. It is hard to describe, but as the range
goes onto infinity, it seems to miss numbers, when you plot the
frequency distrbution, it looks as if it has steps.

This problem occours when passing a 16bit random number into the
algorithm. Going to a 30bit (the best knuth rng i have at the moment)
helps this problem. Moving to double precision also helps, but only if
you have the larger random number. So i am guessing that there is two
problems here, one being precision, and the other with the the random
number size.

For a given input U the X = ceil(ln(U)/...) algorithm will always produce
the same output X. That means that if you have 2^16 possible inputs--which
you do, if you started with a 16 bit random number--then you have no more
than 2^16 possible outputs. Clearly, then, it is not possible that the
output would cover every integer from 0 to 100m (which I assume is
100*10^3?).

Actually it is worse than that though; that function will map many values of
U to the same value of X. It has to to change the distribution. The variable
U is uniformly distributed. Assume that the most likely `bin' of the output
occurs with probability p1; that is the probability that you win on the
first try, or

p1 = (1/200) = 5e-3

Then the least likely possibility is that you win after exactly 100e3 tries,
or

p2 = (1/200)*(1 - 1/200)^(100e3 - 1) = 1.0215e-220

So that means that if exactly one input code (value of U) corresponds to X =
100e3, then you need to have p1/p2 codes map on to p1... So that breaks down
pretty badly. You can sum up the geometric series to find out how many bits
of uniformly distributed randomness you need to theoretically cover every
code, and that depends on p.

This is before considering rounding error (`precision'). That can only make
things worse.

Do you really need to go up to a hundred thousand? That is not a large
probability, unless I screwed something up (though if I screwed something up
then p2/p1 is smaller and it's not so bad numerically).

Have you ever seen the Slovak beer? 500 mL cans, 12%.

Jonathan

 

[budgie]

On Thu, 09 Jun 2005 19:14:16 +1000, The Real Andy <will_get_back_to_you_on_This>
wrote:

On Wed, 8 Jun 2005 19:13:44 -0400, "Jonathan Westhues"
google-for-my-name@nospam.com> wrote:
(snip)

Someone said jackpots, and that would be my bet too.

(snip)

Spot on Jonathan.

Still up at JC?

 

[The Real Andy]

On Thu, 9 Jun 2005 15:13:20 -0400, "Jonathan Westhues"
<google-for-my-name@nospam.com> wrote:

"The Real Andy" <will_get_back_to_you_on_This> wrote in message
news:4i0ga114v7ob3qhn8drajkiajlpfl6jecj@4ax.com...
However, take a range from 1 to 100m for example, then you get a
problem, which I am going to call stepping. I would love to know the
mathematical term for this. It is hard to describe, but as the range
goes onto infinity, it seems to miss numbers, when you plot the
frequency distrbution, it looks as if it has steps.

This problem occours when passing a 16bit random number into the
algorithm. Going to a 30bit (the best knuth rng i have at the moment)
helps this problem. Moving to double precision also helps, but only if
you have the larger random number. So i am guessing that there is two
problems here, one being precision, and the other with the the random
number size.

For a given input U the X = ceil(ln(U)/...) algorithm will always produce
the same output X. That means that if you have 2^16 possible inputs--which
you do, if you started with a 16 bit random number--then you have no more
than 2^16 possible outputs. Clearly, then, it is not possible that the
output would cover every integer from 0 to 100m (which I assume is
100*10^3?).

Actually it is worse than that though; that function will map many values of
U to the same value of X. It has to to change the distribution. The variable
U is uniformly distributed. Assume that the most likely `bin' of the output
occurs with probability p1; that is the probability that you win on the
first try, or

p1 = (1/200) = 5e-3

Then the least likely possibility is that you win after exactly 100e3 tries,
or

p2 = (1/200)*(1 - 1/200)^(100e3 - 1) = 1.0215e-220

So that means that if exactly one input code (value of U) corresponds to X =
100e3, then you need to have p1/p2 codes map on to p1... So that breaks down
pretty badly. You can sum up the geometric series to find out how many bits
of uniformly distributed randomness you need to theoretically cover every
code, and that depends on p.

This is before considering rounding error (`precision'). That can only make
things worse.

Do you really need to go up to a hundred thousand? That is not a large
probability, unless I screwed something up (though if I screwed something up
then p2/p1 is smaller and it's not so bad numerically).

Actually, the largest range we will be using is about 5,000,000

Now I have tested the theory, I guess I have to start the search for a
double precision library in C that I can easily port to the rabbit.
There is a guy on the rabbit forums advertising a dpfp lib for $75 but
he wont answer any of my emails.

Have you ever seen the Slovak beer? 500 mL cans, 12%.

Sounds scary, a bit like going back to C after using c++ for so long
I love oo now. I must go and find some of those beers, they must
taste horrible! I kinda get a thrill from doing dangerous stuff like
that. I did once drink 3 of those 1 litre german beers once.

Jonathan

 

[The Real Andy]

On Thu, 09 Jun 2005 13:05:34 GMT, Fred Bloggs <nospam@nospam.com>
wrote:

However, take a range from 1 to 100m for example, then you get a
problem, which I am going to call stepping. I would love to know the
mathematical term for this. It is hard to describe, but as the range
goes onto infinity, it seems to miss numbers, when you plot the
frequency distrbution, it looks as if it has steps.

This problem occours when passing a 16bit random number into the
algorithm. Going to a 30bit (the best knuth rng i have at the moment)
helps this problem. Moving to double precision also helps, but only if
you have the larger random number. So i am guessing that there is two
problems here, one being precision, and the other with the the random
number size.


That's easily taken care by going to a multi-dimensional distibution,
but I'll wait to see what they can extract from Knuth on this.

I think if use the MT prng, combine two 32bit random numbers, and use
double precision i should be fine. I am keen to hear about your method
though, please tell!

 

[The Real Andy]

On Fri, 10 Jun 2005 08:42:57 +0800, budgie <me@privacy.net> wrote:

On Thu, 09 Jun 2005 19:14:16 +1000, The Real Andy <will_get_back_to_you_on_This
wrote:

On Wed, 8 Jun 2005 19:13:44 -0400, "Jonathan Westhues"
google-for-my-name@nospam.com> wrote:
(snip)

Someone said jackpots, and that would be my bet too.

(snip)

Spot on Jonathan.

Still up at JC?

If you mean Jupiters, no. I was made redundant with about 10 others
about 3 year ago when TabCorp purchased Jupiters. I saw the redundancy
package a could not stop smiling for a week! Crazy thing is, the are
now rehiring everyone that they sacked, but my theory is that you
always go fowards, and never look back. I am going to try and get out
of the gaming industry soon.