+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: Looking for more algorithms for quasi-random numbers (/thread-14070.html)
Looking for more algorithms for quasi-random numbers - Namir - 11-29-2019 01:06 PM
Hi All Math Lovers,
In another thread of mine, ttw mentions quasi-random numbers. Quasi-random numbers (QRNs) present a better spread over a range of values than pseudo-random numbers (PRNs). On the other hand, QRNs will often fail randomness tests. They true purpose to to cover more uniformly a range of values in one of more dimensions.
This is part of ttw's response in my other thread, where he mentions QRNs:
Quote:The easiest multi-dimensional quasi-random sequence is the Richtmeyer sequence. One uses the fractional part of multiples of the square roots of primes. Sqrt(2), Sqrt(3), etc. It's quick to do these by just setting x(i)=0 updating by x(i)=Frac(x(i)+Sqrt(P(i))). Naturally one just stores the fractional parts of the irrationals and updates. (List mode). The sequence is also called the Kronecker or Weyl sequence at times.
The above text includes the algorithm of setting x(1)=0 updating by x(i)=Frac(x(i)+Sqrt(P(i))). The array of P() represents prime numbers starting with 2. You can change x(1) to had a uniform random number as a seed (to generate different sequences every time you apply the algorithm) or simply set x(1) = sqrt(P(1)) = sqrt(2).
I am curious about other formulas to calculate sequences of quasi-random numbers. You are welcome to use your imagination. My first attempt was something like:
Code:
n = number of x to generate
m = 100*n
Calculate P() for primes in the range of 1 to m
X(1) = rand or Frac(ln(P(3)) * sqrt(P(1))
j = 2
count = 0
for i=2 to n
X(i) = Frac(X(i) + ln(P(j+1)) * sqrt(P(j-1))
j = j + 1
if j > m then
count = count + 1
j = 2 + count
end
endThe above code produces x() with a mean near 0.5 and standard deviation near 0.28. The auto correlations for the first 50 lags are in the orde rof 10^(-2) to 10^(-4).
I am curious about other formulas to calculate sequences of quasi-random numbers. You are welcome to use your imagination. You can even commit math heresy!!! As long as it works, you are fine (and forgiven) :-)
Namir
RE: Looking for more algorithms for quasi-random numbers - SlideRule - 11-29-2019 04:49 PM
Perusal of Quasi-random sequences in art and integration, John D. Cook Consulting, illumes the phenomena with references to more descriptive books; Random Number Generation and Quasi-Monte Carlo Methods & Monte Carlo and Quasi-Monte Carlo Methods, on the same.
BEST!
SlideRule
RE: Looking for more algorithms for quasi-random numbers - mfleming - 11-30-2019 01:36 AM
(11-29-2019 01:06 PM)Namir Wrote: This is part of ttw's response in my other thread, where he mentions QRNs:
Quote:The easiest multi-dimensional quasi-random sequence is the Richtmeyer sequence. One uses the fractional part of multiples of the square roots of primes. Sqrt(2), Sqrt(3), etc. It's quick to do these by just setting x(i)=0 updating by x(i)=Frac(x(i)+Sqrt(P(i))). Naturally one just stores the fractional parts of the irrationals and updates. (List mode). The sequence is also called the Kronecker or Weyl sequence at times.
Using "quote" in place of "code" will autowrap large blocks of text!
RE: Looking for more algorithms for quasi-random numbers - Namir - 11-30-2019 01:29 PM
(11-30-2019 01:36 AM)mfleming Wrote:(11-29-2019 01:06 PM)Namir Wrote: This is part of ttw's response in my other thread, where he mentions QRNs:
Quote:The easiest multi-dimensional quasi-random sequence is the Richtmeyer sequence. One uses the fractional part of multiples of the square roots of primes. Sqrt(2), Sqrt(3), etc. It's quick to do these by just setting x(i)=0 updating by x(i)=Frac(x(i)+Sqrt(P(i))). Naturally one just stores the fractional parts of the irrationals and updates. (List mode). The sequence is also called the Kronecker or Weyl sequence at times.
Using "quote" in place of "code" will autowrap large blocks of text!
I learned that the hard way :-)
RE: Looking for more algorithms for quasi-random numbers - Namir - 11-30-2019 07:52 PM
The few leads I got from the nice folks on this web were able to lead me to methods that generate sequences of quasi-random numbers that are practically perfectly distributed. I got what I was looking for. Thanks!!!
Namir
RE: Looking for more algorithms for quasi-random numbers - ttw - 12-01-2019 05:52 AM
This is one of the sequences from my paper in "Computational investigations of low-discrepancy point sets II" from the 1994 Las Vegas Conference on Monte Carlo and Quasi Monte Carlo Methods. I have made a single ad-hoc change (described below) that improves distribution for small numbers of points.
The Halton Sequence Phi(N,P) (for odd primes, 2 is a special case not considered here) can be described as:
1. Generate the digits of N in base P (for P an odd prime). Call these digits a(1) to a(k) where k is the maximum number of digits needed. (There should be lots of subscripts but I'll treat each prime separately to reduce index management.)
N=Sum from j=1 to k of a(j)*P^(j-1), that is: a(k)a(k-1)...a(2)(a(1).
2. Reverse the digits: a(1),a(2)....a(k-1),a(k) is resulting string.
3. Treat this string as a fraction with a decimal point (p-ary point?) in front.
Example: P=3, N=5: 5(3)=12. Reverse Phi(5,3)=.21(3) or 2/3+1/9 = 7/9
The sequence is very well distributed for example one starts with 0, 1/3, 2/3, 1/9, 4/9, 7/9, 2/9, 5/9, 8/9, 1/27, 10,27...
This sequence is uniformly in the unit d-cube using d different primes. For example in 3 dimensions using primes 3 and 5 gives the points: (skipping 0 which sits on the corner of the cube).
(1/3, 1/5, 1/7)
(2/3, 2/5, 2/7)
(1/9, 3/5, 3/7)
(4/9, 4/5, 4/7)
(7/9,1/25,5/7)
(2/9, 6/25, 6/7)
(5/9, 11/25, 1/49)
(8/9. 16/25, 8/49)
etc.
The process is sometimes termed a Kakatumi-von Neumann odometer.
There is a problem that I noticed about 1967 or so when I started working on quasi-Monte Carlo. For large base, the Halton Sequence produces strongly correlated points until enough points are generated. (This happens with all quasi-random sequences but not as severely.) Take the first few points using bases 101 and 103.
(1/101, 1/103)
(2/101, 2/103)
...
(100/101, 100/103)
(1/10201, 101/103)
(102/10201, 102/103)
(203/10203, 1/10609)
etc.
In 1993-1995 period, I figured out to multiply each numerator by a number (I called a spin) to break this up. Then I used the fact the fractional parts of square roots of primes are independent to do the following seemingly strange rule.
Give a prime P, to find a multiplier S, do the following.
1. Compute the nearest integers to the fractional part of the Sqrt(P), call these H and L (high and low, one is above and one below the number).
2. Compute the continued fraction of H/P and L/P; each generates a string of partial quotients. The multiplier S is the one of these satisfying the following:
3. A. Chose the one for which the sum of the partial quotients is smallest.
B. If tied, chose the one with the smallest maximum partial quotient.
C. If tied, chose whichever H/P or L/P is closest to the fractional part of the square root.
D. Ad Hoc Alert: if P=41, use 16. (To avoid 17/41 being close to 12/29. The only such case in all primes less than 2^32)
4. To generate the modified Phi sequence Phi(N,P,S): generate the digits of N base P as above and reverse. Multiply each digit by S modulo P and sum as above.
Examples: 3 dimensions: P=3, 5, 7 have S= 2, 2,and 5 respectively.
(2/3, 2/5, 5/7)
(1/3, 4/5, 3/7)
(2/9, 1/5, 1/7)
(8/9, 3/4, 6/7)
etc.
For the pathological case 101 and 103, the multipliers are 6 and 16 respectively (not the best but that's another post sometime).
(6/101, 16/103)
(12/101, 32/103)
(18/101, 48/103)
Clearly more spread out than the original Halton Sequence.
I've got some more but Mordechay Levin's paper ArXiv 1806 shows that even the original Halton Sequence hits the theoretical lower bound for dimensions 2 and up so great changes cannot be had by tinkering. I do have a bit better, but it's even longer to compute and I haven't tested the new ideas thoroughly.
HP50g code:(Number, Prime, Spin, Top, Bottom)
<< 0 1 -> N P S T B << WHILE N REPEAT
N S * P MOD T P * + 'T' STO
P 'B' STO* N P IQUOT 'N' STO END
T B />> >>
Not necessarily the fastest but eas to work with. (If I've done this right.)
RE: Looking for more algorithms for quasi-random numbers - Csaba Tizedes - 12-01-2019 11:46 AM
One article about random sampling - maybe helps