Fifty Shades of J/Chapter 28

From J Wiki
Jump to: navigation, search

Table of Contents ... Glossary ... Previous Chapter ... Next Chapter

Have you a weight on your mind?

Principal Topics

? (deal) {: (take) # (copy) I. (interval index) #. (base) random exponential, random Normal, Box-Muller formula, boundary values, grouping, rounding, frequency distributions.

Uniform Random numbers

A standard procedure for generating n random numbers in the range (0,1) with all numbers in that range being equi-probable is

   rnd=: ?@#&0
   rnd 5
0.5377 0.06353 0.7059 0.5188 0.8832

Converting a list such as the above to random values in the range [a,b] is a matter of scaling, conveniently supplied by the base verb #.. This transforms [a,b] into the list [(b-a),a] by

   relist=: (-/ , {:)@|.                       NB. [a,b] -> [(b-a),a]
   run=: dyad : '> (rnd y) #.each <relist x'   NB. x=range [a,b], y=n
   5 11 run 5
5.544 10.69 5.442 8.004 7.305

Converting a rnd n sequence to a series of uniform random integers from say 1 to 20 is a simple exercise in multiplying and rounding down, for example :

   <. 20 * rnd 5
14 15 14 1 16

Weighted Random numbers

If the uniform, that is equi-probable, requirement is abandoned in favour of making some integers more probable than others then weights must applied. These can be specified by a list of values, for example by specifying m numbers (not necessarily integers) which define the relative frequencies of the integers in i.m . Weights are more tractable when they are expressed cumulatively and normalised to 1, hence

   cumwts=: +/\ % +/                 NB. cumulative weights
   cumwts 1 4 3 2
0.1 0.5 0.8 1

The primitive I. (interval index) identifies the appropriate cell into which each item should be placed :

   rndw=: cumwts@[ I. rnd@]          NB. weighted random integers
   1 4 3 2 rndw 10
1 1 2 2 0 3 1 1 3 2

As an aside, the syntactic structure of rndw is a fork with separate transformations applied to the left and right tines.

Interestingly, reversing the arguments of rndw is equivalent to transposing the above binary table (note the left argument to I. needs to be sorted in ascending order to work correctly), and the result is the cumulative frequency distribution of the set of random drawings :

   cdrndw=: sort@rnd@] I. cumwts@[   NB. cumulative distn of above
   1 4 3 2 cdrndw 10
1 5 8 10

which happens in the above drawing to match exactly the pattern of cumulated weights. A more general drawing might be

   1 4 3 2 cdrndw 100
10 52 84 100

still visibly confirming that the weights have been correctly applied.

Random Exponential values

Drawings from a random exponential distribution with a given mean are a simple extension of rnd :

   rne=: -@* ^.@rnd            NB. rand neg exp. x=mean, y=n
   5 rne 8
0.9258 16.82 14.64 3.177 1.994 24.33 4.793 13.53
   (mean=: +/ % #) 5 rne 10
5.985

Such values are commonly used in simulations where, if the probability of an arrival in a time interval is proportional to the length of the interval and successive arrivals are independent of each other, then statistically the negative exponential distribution is that followed by the gaps between successive arrivals.

Random Normal values

Now that rne and rnd are available, random values from the standard Normal distribution are obtainable using a pseudo-random technique known as the Box-Muller formula. This involves two separate drawings of random numbers using rnd, one to furnish rne and the other used directly in the right hand multiplicand :

   rno=: 3 : '(%:2 rne y)*2 o.(rnd y)*o.2'   NB. random Normal(0,1)
   rno 8
_0.6224 0.0922 0.9571 _0.3281 0.02999 1.28 _0.1247 0.02563

As an aside, 1 o. (sin) replacing 2 o. (cos) in the right hand multiplicand also delivers random Normal values.

A random drawing from a Normal distribution with a given mean and standard deviation is

   rnorm=: 4 : '(rno y) #.every <|.x'        NB. rand Normal, x=mean,sd, y=n
   10 3 rnorm 8
11.18 10.37 10.55 12.97 7.146 11.75 10.61 11.83

Grouping and Rounding

When continuous variables are involved, it is often useful to group values by rounding :

   R=: <.@(0.5&+)             NB. round y to nearest integer
   round=: 4 : 'R&.(%&x) y'   NB. round y to nearest x
   0.2 round 5 rne 10
7.2 0.2 0.8 6 2.2 5.2 1.4 0 1 2.6
   5 round 100 15 rnorm 10
110 70 30 145 115 105 105 75 100 115

Frequency distributions of random samples

Interval Index (I.) provides a means of obtaining indexes of items in intervals defined by the left argument which defines boundaries in ascending order, Items below the lowest boundary have index 0, and #x for those which exceed the highest boundary. Boundary items are assigned to the lower interval.

   fr=: dyad : '#/.~ sort x I. y'   NB. frequency distribution
   y=: 2.7 1.8 0.5 3.3 1.4 3.2 3.1 0
   1 1.5 2 2.5 3 fr y
2 1 1 1 3

The plausibility of the rnorm routine can be checked by, for example

   (67.5 + 5 * i.14) fr 5 round 100 15 rnorm 100
3 1 2 9 7 10 10 14 14 13 5 4 4 3 1

This illustration confirms the general bell shape with values more than two standard deviations away from the mean accounting for about 5% of the total.

Code Summary

rnd=: ?@#&0                                   NB. random uniform(0,1)
run=: dyad : '(rnd y) #.each <relist x'       NB. rand uniform
relist=: (-/ , {:)@|.
cumwts=: +/\ % +/                             NB. cumulative weights
rndw=: cumwts@[ I. rnd@]                      NB. weighted random integers
cdrndw=: sort@rnd@] I. cumwts@[               NB. cumulative distn of rndw

rne=: -@* ^.@rnd                              NB. rand neg exp. x=mean, y=n
rno=: 3 : '(%: 2 rne y) * 2 o.(rnd y) * o.2'  NB. Box-Muller formula
rnorm=: dyad : '(rno y) #.every <|.x'         NB. rand Normal, x=mean,sd, y=n
round=: dyad : 'R&.(%&x)y'                    NB. round y to nearest x
R=: <.@(0.5&+)                                NB. round y to nearest integer
fr=: dyad : '+/(=/ ~.) sort x I. y'           NB. frequency distribution

Script

File:Fsojc28.ijs