🚀 go-pugleaf

In article <1992Oct30.182002.11907@netcom.com> tjm@netcom.com (Thomas J. Merritt) writes:
>In article <BwxsF6.3DF@mentor.cc.purdue.edu> hrubin@pop.stat.purdue.edu (Herman Rubin) writes:

			...................

>Herman,
>	Why don't you just post the real code rather than these uncompilable
>snipettes that don't convey what you are trying to do very well?  Are you
>afraid your ego will get deflated if someone shows you a better way to
>solve your problem?  From the code you have posted I am certain that there is
>a more structured way to accomplish your task with out any loss in performance.

I HAVE posted the real code.  Let me state the problem:  It is desired to
generate exponential random variables with mean 2 in an infinite precision
manner, generating as few bits as practical and using as few bits as 
practical in the procedure.  For various reasons concerned with the
problem, this means generating the integer part and a bit mask for the
fraction, so that the final result is 

	int + (uniform_rv &~ mask).

With this formulation, the integer part and the mask are independent, and
the various bits of the mask are also independent.  One can compute the 
amount of information in each.  It is also desired to do this exactly.

Now one way of doing this is to generate Poisson random variables with
mean .5, the integer part being the number of 0's until a non-zero value
k is obtained, and the fractional part is the minimum of k uniform random
variables.  This also should be done exactly, assuming that perfect random
bits are available.  The algorithm does this in the following manner:

Let q_j, j = 0, ..., be 2^-(j+1)/j!.  The sum of the q_j is .5*exp(.5).
With probability q_j, return j as the Poisson random variable, and there
may be some by-products useful in the minimum phase.  These by-products
are a geometric random variable g4 with p=3/4 and/or a geometric random
variable g16 with p=15/16.  If j >= 5, g16 is produced, and if j >=3,
but j != 5, g4 is produced,  This uses certain properties of base 2
expansions, to be illustrated.  With the remaining probability, nothing
is returned, and the process continued.  A very greedy algorithm is used.

To see what the procedure is, let us write out in hex the q_j for j <= 5.

	0	.8
	1	.4
	2	.1
	3	.02AAAAAAA....
	4	.005555555...
	5	.000888888...

in addition, the sum from 6 to infinity is bounded by .0001111 ...
Now the sum of q_3 and q_4 is .003.  Also, note that 

	q_3 = .02 (1 + 4^-1 + 4^-2 + ...)
and
	q_4 = .004 (1 + 4^-1 + 4^-2 + ...)

So we can start by taking a random variable g with P(g=m) = 2^-m.  Then
we make the following return schedule

	 1		0
	 2		1
	 3		reject
	 4		2
	 5		reject
	 6		reject
	 7		3
	 8		3 or 4
	 9		reject
	10		reject
	11		reject
	12		reject
	13 + 4k		5, g16 = (k+1)
	14 + 4k		reject
	15 + 4k		reject
	16 + 4k		g16 = (k+1), process further

In case 8, we get another geometric, which will again be called g.  Now
cases 7 and 8 together, with each line having half the probability of the
previous one, come out as

	7		3, g4=1
					nonexistent
	8 1		3, g4=2
	8 2		4, g4=1
	8 3		3, g4=3
	8 4		4, g4=2

etc.

In the process further, we again get a new g, and form g4=(g+1)>>1.
Now the probability of 6 is 2/3 of the probability of the situation,
so we return j=6 if g is odd.  Then twice the probability of 7 is 2/7
pf that where g is even, so we reject with probability 5/7.  Then 
we loop, starting with j = 7:

	accept j with probability 1/2
	reject with probability j/(j+1)
	increase j by 1 and loop.

This either rejects, in which case we repeat, returns 0, in which case we 
add 1 to the integer and repeat, or returns something larger than 1, in
which case we exit to the appropriate place in the next stage, which is
to get the mask.  The details given accomplish this, and better can be done,
but we are out in a rare situation already.

The mask stage proceeds in the following manner.  We wish to get the
minimum of j random varibles as a mask to be imposed on a uniform.  One
way to do this is to generate the number k of 0's in j bits, and if 
0 < k < j set this bit 0 (otherwise skip and continue with the same j),
and proceed to get the minimum of k uniforms until we get to k=1.  Now
if j > 5, we do not bother to optimize, and do this crudely.  If j=5,
the position of the first bit where j is reduced is distributed as g16,
which we use.  In this case, the probabilities of 1,2,3, and 4 are 1:2:2:1,
and we do this by taking another geometric random variable g.  With probability
1/2, k is 2 if g is odd and 1 if g is even.  Otherwise, g4 = (g+1)>>1, and
k is 3 if g is odd and 4 if g is even.

For j=4, we get another g.  The first bit corresponding to reducing j is
given by the minimum of g and g4.  If g4 < g, this corresponds to setting
k = 2, and the whole process is done by setting the g4-th and the g-th bits
to the right of the current position.  If g <= g4, this corresponds to 
setting the g-th bit, with probability 1/2 setting k=1 (exit), and otherwise
getting a new g4 = g4 - g + 1, and setting k = 3.  The actual procedure 
accomplishes this with a little less work.

For j=3, the first bit to be set is given by g4.  Then with probability
1/2 we should exit with k=1, and otherwise k=2.  A procedure just as bit
efficient and easier to manage is to move the bit-pointer back one place
and continue with the j=2 procedure, which consists of getting another
g and shifting the pointer right g positions before inserting it into the
mask.

Okay, that is the whole procedure.  It should be fairly clear from this
that if the procedure to getting a geometric (1/2) random variable g is
a fast hardware procedure, the whole process will be much faster.
-- 
Herman Rubin, Dept. of Statistics, Purdue Univ., West Lafayette IN47907-1399
Phone: (317)494-6054
hrubin@pop.stat.purdue.edu (Internet, bitnet)  
{purdue,pur-ee}!pop.stat!hrubin(UUCP)