rational arithmetic — RNG

[Roger|Dyalog]

The initial implementation of the rational arithmetic facility does not have ?Q ⍵ where ⍵ is a rational number. ⍵ is required to be greater than 0 with a denominator of 1. Therefore, the problem is to generate a random number in the uniform distribution ⍳⍵ where ⍵ is an integer represented as a vector of decimal digits.

The technique is as follows, explained first with ordinary integers. To generate uniform random numbers ?n⍴q, we generate ?m⍴p where p≥q and keep only the numbers which are less than q. m can be as large as we like, as large as necessary to generate n random numbers less than q.

An example:

      p←30
      q←10
      n←1e6

That is, we want to generate 1e6 uniform random numbers less than 10.

      m← n×⌈1.1×p÷q   ⍝ 1.1× to ensure n suitable numbers

   x← ?m⍴p
   +/ x<q
1332913

   y← n↑(x<q)⌿x

   sort←{(⊂⍋⍵)⌷⍵}
   ⊢ c← sort {⍺,≢⍵}⌸y
0 100650
1  99759
2 100187
3  99817
4  99574
5 100376
6  99539
7 100556
8  99661
9  99881

Here, c are the count of the numbers 0,1,2,...,9 in y.

The maximum individual absolute difference between the sample cumulative distribution and the theoretical uniform cumulative distribution is:

      KS← {⌈/|((⊃⌽s)÷⍨s←+\⍵)-(1+⍳m)÷m←≢⍵}
      KS c[;1]
0.00065

The Kolmogorov-Smirnov test specifies that the ⍺=0.01 critical value for n>35 is 1.63÷n*0.5, which is 0.00163, larger than the KS stat of 0.00065 above.

The following experiments contrast the technique versus doing ?n⍴q directly.

      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸n↑(x<q)⌿x←?m⍴p
0.000437
      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸n↑(x<q)⌿x←?m⍴p
0.000532
      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸n↑(x<q)⌿x←?m⍴p
0.00063
      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸n↑(x<q)⌿x←?m⍴p
0.00042

      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸?n⍴q
0.000418
      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸?n⍴q
0.000468
      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸?n⍴q
0.000863
      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸?n⍴q
0.000608

I consulted the J Forum on the above technique, and elicited two pertinent replies. Xiao-Yong Jin, a computational scientist at the U.S. Argonne National Laboratory, points out that the technique is a special case of rejection sampling. User “ethiejiesa” offers a proof of the technique using Bayes’ Theorem (translated from J and using a mix of APL and conventional mathematical notation):

Let k←?p. Then for every integer i in ⍳p

      P(k=i) = 1÷p.

We can then model the selective keeping process with a posterior distribution on this pdf. For every i in ⍳q,

      P(k=i|k<q) = P(k<q|k=i) × P(k=i) ÷ P(k<q)
               === 1 × (÷p) ÷ (q÷p)
               === 1÷q.

Happily, this corresponds to the discrete uniform pdf you want. Of course, the standard provisions apply—the random variables are independent, q≤p, etc.

Now the actual code that we want:

      ⍝
irand←{
 p←((⊃⍵)+0∨.≠1↓⍵),(¯1+≢⍵)⍴10
 {0≥⊃0~⍨⍵-z←?p:∇ ⍵ ⋄ (⍳0∧.=z),(∨\0≠z)⌿z}⍵
}

      irand 1 2 3 4 5
6 1 2 6
      irand 1 2 3 4 5
1 1 1 0 3
      irand 1 2 3 4 5
7 1 6 7

irand executes z←?p until z is less than the argument ⍵. If d←⊃⍵ is the leading digit of ⍵, then p is constructed as (1+d),10,...,10 if the rest of ⍵ is not all zeros, or as d,10,...,10 if it is all zeros.

A few experiments, computing ?n⍴11, ?n⍴19, and ?n⍴20 on arguments represented as vectors of decimal digits:

      KS 1⌷⍤1 ⊢ c←sort {⍺,≢⍵}⌸ {10⊥irand ⍵}¨ n⍴⊂1 1
0.000586273
      KS 1⌷⍤1 ⊢ c←sort {⍺,≢⍵}⌸ {10⊥irand ⍵}¨ n⍴⊂1 9
0.000851105
      KS 1⌷⍤1 ⊢ c←sort {⍺,≢⍵}⌸ {10⊥irand ⍵}¨ n⍴⊂2 0
0.001239
      1.63 ÷ n*0.5  ⍝ ⍺=0.01 critical value
0.00163