## rational arithmetic — RNG

**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 !

1 post
• Page

**1**of**1**### 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:

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

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:

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.

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):

Now the actual code that we want:

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:

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 ⍳pP(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

1 post
• Page

**1**of**1**### Who is online

Users browsing this forum: No registered users and 1 guest

Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group