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

## No comments:

Post a Comment