help with using Wichmann-Hill Pseudo Random Number Generator with box Muller box

[stefman]

this is my attempt of the Wichmann-Hill Pseudo Random Number Generator, if you can find any improvement please help:
Public Static Function RndX() As Double
' Wichmann, B.A. and I.D. Hill, Algorithm AS 183: An Efficient and Portable
' Pseudo-Random Number Generator, Applied Statistics, 31, 188-190, 1982.
' C IX, IY, IZ SHOULD BE SET TO INTEGER VALUES BETWEEN 1 AND 30000 BEFORE FIRST ENTRY
' IX = MOD(171 * IX, 30269)
' IY = MOD(172 * IY, 30307)
' IZ = MOD(170 * IZ, 30323)
' RANDOM = AMOD(FLOAT(IX) / 30269# + FLOAT(IY) / 30307# + FLOAT(IZ) / 30323#, 1#)

' Note: the FORTRAN function AMOD(A,P) = A - (INT(A / P) * P)
' VBA Mod function rounds non-integer values of A,P FIRST, so it won't work in the last line!

Dim ix As Long, iy As Long, iz As Long
Dim F_ix As Single, F_iy As Single, F_iz As Single
Dim sum As Single
Dim seeded As Boolean ' initialized to 'False' during first call to function

If seeded = False Then
ix = 123 ' use a fixed initial seed
iy = 234
iz = 345
seeded = True ' static variable set true so we don't 're-seed' every time
End If

' the main part of the linear congruential algorithm is only three lines long
ix = (171 * ix) Mod 30269
iy = (172 * iy) Mod 30307
iz = (170 * iz) Mod 30323

'CSng convert the integer values to single precision reals for arithmetic
F_ix = CSng(ix) / 30269! ' use ! to ensure denominator is also single precision real
F_iy = CSng(iy) / 30307!
F_iz = CSng(iz) / 30323!
sum = F_ix + F_iy + F_iz

' take the fractional part of the sum as the random number
RndX = sum - Int(sum)
End Function

Then i use the random number to use in the box muller fuction i attempted to make:

Public Function RndBM()

' User function to transform uniform random numbers from randBM (x,y)
' into standard normal variates (Z_1, Z_2) using the Box-Muller transformation
' Box-Muller Ann. Math. Stat. 29, 610-611 (1958)

Dim x As Double, y As Double, s As Double, Z_1 As Double
Const TWO_PI As Double = 2# * 3.14159265358979

x = RndX() ' generate two uniform random numbers, x and y
y = RndX()

s = Sqr(-2# * Log(x)) ' transform to normal random numbers, r1 and r2
Z_1 = s * Cos(TWO_PI * y)

RndBM = Z_1 ' return the cosine term
End Function

[stefman]

Wichmann-Hill Pseudo Random Number Generator:

http://links.jstor.org/sici?sici=003...OR-enlargePage

box muller:

http://mathworld.wolfram.com/Box-Mul...formation.html

[jemidiah]

There was a relatively long thread on Wichmann-Hill on these forums a few years ago. Merri gave a few optimization pointers for Ahn's code that may apply to yours. Without being an expert on VB6 optimization myself, yours seems just fine speed-wise. I'm unsure if making the function Static is helpful or not to its speed.

Correctness-wise, your use of ix, iy, and iz to reinitialize themselves on successive calls worries me slightly. Wichmann-Hill suggests ix, iy, and iz be in the range 1 to 30,000, while your usage allows them to be in the range 0 to 30269-1, 0 to 30307-1, or 0 to 30323-1. Conveniently, 30269, 30307, and 30323 are all prime, and the integers modulo primes form an integral domain (actually, a field), so (171 * ix) Mod 30269 etc. will only return 0 when ix (etc.) is itself 0 (mod 30269, etc.). So, reinitializing at least won't make them hit 0. However, they could still become greater than 30,000. This in fact happens: the map ix |-> (171 * ix) Mod 30269 is injective (again, since the integers modulo primes form an integral domain), so is surjective since it's domain and codomain are finite and of the same size. Another potential issue is cycling. The above map has the (thankfully unique) fixed point of ix=0. Setting ix=iy=iz=0 initially, you'll generate the same number again and again forever. Since you never hit 0, this isn't a problem currently.

Again, though, the problem remains of allowing ix, iy, or iz to be > 30,000. If you just try to "mod them down", be careful of the cycling issue, since you'll break the invertibility of the resulting map and may accidentally hit 0 or add small-length cycles.


As for Box-Muller, that seems just fine. My only comment is that I'd prefer your TWO_PI const to be named TAU, but that's purely personal preference. I'd like it if tau were added to the standard mathematical notation for cases like the Fourier transform, or the standard normal PDF, or a more intuitive radians to degrees formula, or a more elegant description of the trig function extrema, or ....


Further edit: I should mention I don't have access to JSTOR right now. This version may be the same as the one you linked.

[jemidiah]

As an explicit example of initial inputs that produce numbers > 30,000, I found (using two online integer-mod-p calculators to double-check the results)

4602 * 171 = 30217 mod 30269
26430 * 172 = 30217 mod 30307
25328 * 170 = 30217 mod 30323

Using ix=4602, iy=26430, iz=25328, after the first iteration each of ix, iy, and iz becomes 30217 which is larger than 30,000. It's unclear to me what the purpose of that requirement is, but the article does say it's necessary. The iterative reinitialization approach produces a period length of 30269*30307*30323 = 27,817,185,604,309 (since the multiplicative order of 171, 172, and 170 modulo 30269, 30307, and 30323 is 30269, 30307, and 30323, from Lagrange's theorem and these modulus's primality). This is apparently the source of the article's sentence starting, "It has a cycle length exceeding 2.78*10^13", so the article's algorithm also uses iterative initialization. (I don't know Fortran enough to deduce this from the code provided.)

However, if initial inputs <= 30,000 are required, using the above initial inputs only delays the problem of ix, iy, and iz being > 30,000 by a single call. It would appear there's a mistake in the writeup or reasoning. I could understand if the authors wanted to keep the code as simple as possible and wrote "should be set to integer values between 1 and 30000" instead of "should be set to integer values between 1 and 30268, 1 and 30306, and 1 and 30322, respectively". That is, it would be reasonable to think there's no statistical reason for ix, iy, and iz being > 30,000, and that the authors wanted to avoid magic constants in their use instructions, expecting people not to look into the wording this carefully. In this case, my objection above is withdrawn (and replaced by a criticism that the paper isn't clear enough).

[stefman]

from looking at the other thread you link to by merri, do you have any insight in the reason he set Const MaxLong As Double = 2 ^ 31 - 1, in the sub RandomizeX, how does the use of it prevent overflow. any help would be usefull and its true that the paper could be more clear, thats why i asked for help.

[jemidiah]

To be blunt, I'll assume most of my last posts went over your head since you didn't mention their real content. To summarize them, I'm now fine with the correctness of your implementations.

Anhn's use of MaxLong, specifically

      d = Abs(Int(Val(Number)))
      If d > MaxLong Then '-- prevent Long overflow
         d = d - Int(d / MaxLong) * MaxLong
      End If
      n = d

is just used so that you can pass a double-precision floating point named Number to RandomizeX. Number is then "modded down" to be in the valid range for Longs. You'll notice Merri's first post doesn't include this logic, since "Number" becomes a Long. I'd say you don't need to worry about it whatsoever.

[stefman]

the help you have been giving has been very useful ,i understood most of what you said after referring back to the paper, it clear the point that the ix, iy, and iz being > 30,000 for me after the first call. the method i used produces numbers less than 0 and greater than 1, while merri version run perfectly have no such problem. this is the main problem i have!

[jemidiah]

the method i used produces numbers less than 0 and greater than 1

I used your code to generate 1,000,000 numbers and none of them were less than 0 or greater than 1. It also doesn't make sense conceptually that you should get that error. Could you give the input ix, iy, and iz which result in this issue?

[stefman]

yea your right it was my fault, i was looking at the wrong output on the excel our sheet, thank you for your help! i actually understand the paper more fully etc...

[jemidiah]

Glad it worked out. I've never looked into PRNGs and this was interesting.