rational arithmetic — RNG

APL-related discussions - a stream of APL consciousness.
Not sure where to start a discussion ? Here's the place to be
Forum rules
This forum is for discussing APL-related issues. If you think that the subject is off-topic, then the Chat forum is probably a better place for your thoughts !

rational arithmetic — RNG

Postby Roger|Dyalog on Sat Jul 18, 2020 12:48 am

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
Roger|Dyalog
 
Posts: 207
Joined: Thu Jul 28, 2011 10:53 am

Return to APL Chat

Who is online

Users browsing this forum: No registered users and 1 guest