## 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

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<q1332913   y← n↑(x<q)⌿x   sort←{(⊂⍋⍵)⌷⍵}   ⊢ c← sort {⍺,≢⍵}⌸y0 1006501  997592 1001873  998174  995745 1003766  995397 1005568  996619  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⍴p0.000437      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸n↑(x<q)⌿x←?m⍴p0.000532      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸n↑(x<q)⌿x←?m⍴p0.00063      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸n↑(x<q)⌿x←?m⍴p0.00042      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸?n⍴q0.000418      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸?n⍴q0.000468      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸?n⍴q0.000863      KS 1⌷⍤1 ⊢sort {⍺,≢⍵}⌸?n⍴q0.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 56 1 2 6      irand 1 2 3 4 51 1 1 0 3      irand 1 2 3 4 57 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 10.000586273      KS 1⌷⍤1 ⊢ c←sort {⍺,≢⍵}⌸ {10⊥irand ⍵}¨ n⍴⊂1 90.000851105      KS 1⌷⍤1 ⊢ c←sort {⍺,≢⍵}⌸ {10⊥irand ⍵}¨ n⍴⊂2 00.001239      1.63 ÷ n*0.5  ⍝ ⍺=0.01 critical value0.00163`
Roger|Dyalog

Posts: 207
Joined: Thu Jul 28, 2011 10:53 am