Generating pseudo-random numbers is a key piece of much of modern statistical practice, whether for Markov chain Monte Carlo applications or simpler simulations of what raw data would look like under given circumstances.

Generating sufficiently random pseudo-random numbers is not trivial and many methods exist for doing so. Unfortunately, there's no way to prove a series of pseudo-random numbers is indistinguishable from a truly random series of numbers. Instead, there are simply what amount to ad hoc tests, one of which might be failed by insufficiently random series of numbers.

The first 10 million digits of Pi can be downloaded from here. Here, we explore couple of simple ad hoc checks of randomness. We'll try something a little trickier in the next entry.

**SAS**

We start by reading in the data. The file contains one logical line with a length of 10,000,000. We tell SAS to read a digit-long variable and hold the place in the line. For later use, we'll also create a variable with the order of each digit, using the

`_n_`implied variable (section 1.4.15).

data test;

infile "c:\ken\pi-10million.txt" lrecl=10000000;

input digit 1. @@;

order = _n_;

run;

We can do a simple one-way chi-square test for equal probability of each digit in

`proc freq`.

proc freq data = test;

tables digit/ chisq;

run;

Chi-Square 2.7838

DF 9

Pr > ChiSq 0.9723

Sample Size = 10000000

We didn't display the counts for each digit, but none was more than 0.11% away from the expected 1,000,000 occurrences.

Another simple check would be to assess autoregression. We can do this in

`proc autoreg`. The

`dw=2`option calculates the Durbin-Watson statistic for adjacent and alternating residuals. We limit the observations to 1,000,000 digits for compatibility with R.

proc autoreg data=test (obs = 1000000);

model digit = / dw=2 dwprob;

run;

Durbin-Watson Statistics

Order DW Pr < DW Pr > DW

1 2.0028 0.9175 0.0825

2 1.9996 0.4130 0.5870

We might want to replicate this set of tests for series of 4 digits instead. To do this, we just tell the

`data step`to read the line line in 4-digit chunks.

data test4;

infile "c:\ken\pi-10million.txt" lrecl=10000000;

input digit4 4. @@;

order = _n_;

run;

proc freq data = test4;

tables digit4/ chisq;

run;

Chi-Square 9882.9520

DF 9999

Pr > ChiSq 0.7936

Sample Size = 2500000

proc autoreg data=test4 (obs = 1000000);

model digit = / dw=3 dwprob;

run;

Durbin-Watson Statistics

Order DW Pr < DW Pr > DW

1 2.0014 0.7527 0.2473

2 1.9976 0.1181 0.8819

3 2.0007 0.6397 0.3603

So far, we see no evidence of a lack of randomness.

**R**

In R, we use the

`readLines()`function to create a 10,000,000-digit scalar object. In the following line we split the digits using the

`strsplit()`function (as in section 6.4.1). This results in a list object, to which the

`as.numeric()`function (which forces the digit characters to be read as numeric, section 1.4.2) cannot be applied. The

`unlist()`function converts the list into a vector first, so that

`strsplit()`will work. Then the

`chisq.test()`function performs the one-way chi-squared test.

mypi = readLines("c:/ken/pi-10million.txt", warn=FALSE)

piby1 = as.numeric(unlist(strsplit(mypi,"")))

chisq.test(table(piby1), p=rep(0.1, 10))

This generates the following output:

Chi-squared test for given probabilities

data: table(piby1)

X-squared = 2.7838, df = 9, p-value = 0.9723

Alternatively, it's trivial to write a function to automatically test for equal probabilities of all categories.

onewaychi = function(datavector){

datatable = table(datavector)

expect = rep(length(datavector)/length(datatable),length(datatable))

chi2 = sum(((datatable - expect)^2)/expect)

p = 1- pchisq(chi2,length(datatable)-1)

return(p)

}

> onewaychi(piby1)

[1] 0.972252

The Durbin-Watson test can be generated by the

`dwtest`function, from the

`lmtest`package. Using all 10,000,000 digits causes an error, so we use only the first 1,000,000.

> library(lmtest)

> dwtest(lm(piby1[1:1000000] ~ 1))

Durbin-Watson test

data: lm(piby1[1:1e+06] ~ 1)

DW = 2.0028, p-value = 0.9176

alternative hypothesis: true autocorrelation is greater than 0

To examine the digits in groups of 4, we read the digit vector as a matrix with 4 columns, then multiply each digit and add the columns together. Alternatively, we could use the

`paste()`function (section 1.4.5) to glue the digits together as a character string, the use the

`as.numeric()`to convert back to numbers.

> pimat = matrix(piby1, ncol = 4,byrow=TRUE)

> head(pimat)

[,1] [,2] [,3] [,4]

[1,] 1 4 1 5

[2,] 9 2 6 5

[3,] 3 5 8 9

[4,] 7 9 3 2

[5,] 3 8 4 6

[6,] 2 6 4 3

> piby4 = pimat[,1] * 1000 + pimat[,2] * 100 +

+ pimat[,3] * 10 + pimat[,4]

> head(piby4)

[1] 1415 9265 3589 7932 3846 2643

# alternate approach

# piby4_v2 = as.numeric(paste(pimat[,1], pimat[,2],

# pimat[,3], pimat[,4], sep=""))

> onewaychi(piby4)

[1] 0.7936358

> dwtest(lm(piby4[1:1000000] ~ 1))

Durbin-Watson test

data: lm(piby4[1:1e+06] ~ 1)

DW = 2.0014, p-value = 0.753

alternative hypothesis: true autocorrelation is greater than 0

## 2 comments:

Just realized I failed to link to the data! Here it is: http://www.cs.princeton.edu/introcs/data/pi-10million.txt

neat! thanks, just the problem I was looking at recently.

Post a Comment