SAS and R

In example 8.1, we considered some simple tests for the randomness of the digits of Pi. Here we develop a different test and implement it.

If each digit appears in each place with equal and independent probability, then the places between recurrences of a digit should be Pr(gap = x) = .9^x * .1-- the probability the digit recurs immediately is 0.1, the probability it recurs after one intervening digit is 0.09, and so forth. Our test will compare the observed distribution of gaps to the null hypothesis described by the equation above. We'll just look at the gaps between the occurrence of the numeral 1.

R

The piby1 object, created in example 8.1, is a vector of length 10,000,000 containing the digits of Pi. We use the grep() function (section 1.4.6) to identify the places in the vector where the string "1" appears. Then we subtract the adjacent locations to find the distance between them, subtracting 1 to count the spaces.

grep("1", piby1)
gap = pos1[2:length(pos1)] - pos1[1:length(pos1) -1] -1

> head(gap)
[1] 1 33 2 8 18 25

Recollecting that Pi begins 3.141, the first value, 1, in the gap vector reflects accurately that there is a single non-1 digit between the first and second occurrences of the numeral 1. There are then 33 non-1 digits before another 1 turns up.

We next need to calculate the probabilities of each gap length under the null. We'll do that with the : shorthand for the seq() function (section 1.11.3). We'll lump all gaps larger than 88 digits together into one bin.

probgap = 0.9^(0:88) * .1
probgap[89] = 1 - sum(probgap[1:88])

How unusual was that 33-digit gap between the second and third 1?

> probgap[34]
[1] 0.003090315

Under the null, that gap should happen only about once in every 300 appearances.

Next, we lump the observed gaps similarly, then run the one-way chi-squared test.

obsgap = c(gaptable[1:88], sum(gaptable[89:length(gaptable)]))
chisq.test(obsgap, p=probgap)

 Chi-squared test for given probabilities
data: obsgap 
X-squared = 82.3844, df = 88, p-value = 0.6488

This test, more stringent than the simple test for equal frequency of appearance, also reveals no violation of the null.


SAS

In SAS, this operation is much more complex, at least when limited to data steps. Possibly it would be easier with proc sql.

We begin by finding the spaces between appearances of the numeral 1, using a retain statement to "remember" the number of observations since the last time a 1 was observed.

data gapout;
set test;
retain gap 0;
 target=1;
 if digit eq target then do;
 output;
 gap=0;
 end;
 else gap = gap + 1;
run;

Next, we tabulate the number of times each gap appeared, using proc freq to save the result in a data set (section 2.3.1).

proc freq data=gapout (firstobs=2) noprint;
 tables gap / out=c1gaps;
run;

Next, we calculate the expected probabilities, simultaneously lumping the larger gaps together. Here we manually summed the cases-- this would be a somewhat more involved procedure to automate. The final if statement (section 1.5.1) prevents outputting lines with gaps larger than 88.

set c1gaps;
retain sumprob 0;
if _n_ lt 89 then do;
 prob = .9**(_n_ - 1) * .1;
 sumprob = sumprob + prob;
 end;
else if _n_ eq 89 then do;
 prob = 1 - sumprob;
 count=93;
 end;
if _n_ le 89;
run;

As involved as the above may seem, the real problem in SAS is to get the null probabilities into proc freq for the test. SAS isn't set up to accept a data set or variable in that space. One option would be to insert it manually, but it may be worthwhile showing how to do it via a global macro variable.

First, we use proc transpose (section 1.5.4) to make a single observation with all of the probabilities in it.

proc transpose data=c1gaps2 out=c1gaps3;
 var prob;
run;

Next we make a global macro variable to store the printed values of the probabilities, and use the call symput function (section A.8.2) to put the values into the variable; the final trick is to use the of syntax (section 1.11.4) to print all of the probabilities at once.

%global explist;

data probs;
set c1gaps3;
call symput ("explist", catx(' ', of col1-col89));
run;

We're finally ready to run the test. The global macro is called in the tables statement to define the expected probabilities. The weight statement tells SAS that there are count observations in each row.

proc freq data=count1gaps2;
 tables gap / chisq testp= (&explist);
 weight count;
run;
 The FREQ Procedure

 Chi-Square 82.3844
 DF 88
 Pr> ChiSq 0.6488

 Sample Size = 999332

Do the digits of Pi appear in a random order? If so, the trillions of digits of Pi calculated can serve as a useful random number generator. This post was inspired by this entry on Matt Asher's blog.

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