Test statistics - Week 1 Review quiz/assignment

This is Dariia’s Week1 Review Weekend quiz (“assignment”) on evaluating randomness statistical tests from https://www.bigmarker.com/wolfram-u/wsg21-week-1-review-hash-functions-security-features  (time about 40:50).

I first tried to implement the tests from some (random) Googled site, but it didn’t make much sense. But when read the recommended NIST document “A Statistical Test Suite for Random and Pseudorandom Number Generators for Cryptographic Applications” (https://nvlpubs.nist.gov/nistpubs/Legacy/SP/nistspecialpublication800-22r1a.pdf)   it was much clearer stated and with good examples.

Here is a write up of the first three tests from the NIST document:
* 2.1 Frequency (Monobit) Test
* 2.2 Frequency Test within a Block
* 2.3 Runs test

All the implemented functions returns the p-value of the test.

I should also note that the implementations below follows the descriptions as carefully as possible.

Below, the ‘testSequence’ is the small test sequence from each section in the NIST documentation, and ‘sequence’ is a larger random sequence.

Define the larger random sequence:

SeedRandom[Method->"Congruential"];
seqSize = 100000;
sequence=Table[RandomInteger[], seqSize];

2.1 Frequency (Monobit) Test (page 2-4f)

* Replace 0 -> -1
* Sum the sequence
* Compute the test statistic sObs = |Sn] /Sqrt[n]
* Compute p-value:  pVal = erfc(s_obs/Sqrt[2])
* If pVal < 0.01 then it’s not random, otherwise it’s random

testFrequency[sequence_] := Module[{sObs,pVal},
sObs = Abs[Total[sequence /. {0 ->-1}]] / Sqrt[Length[sequence]] // N;
pVal = Erfc[sObs/Sqrt[2]]
]

We first test with the small sequence in the description to ensure that it' s correct . 
  p - value should be 0.527089 (>> 0.01 i . e . random).

testSequence = testSequence = {1,0,1,1,0,1,0,1,0,1}
testFrequency[testSequence]
If[pVal < 0.01, Print["Not random based on the Frequency test"],Print["It looks random based on the Frequency test"]]

Test on the larger sequence:

testFrequency[sequence]
If[pVal < 0.01, Print["Not random based on the Frequency test"],Print["It looks random based on the Frequency test"]]

2.2 Frequency Test within a Block (page 2 - 4 f)

Input:
* sequence: the sequence
* m: length of each block
* n: length of bit string

* Partition into M = |m/M| non overlapping blocks. Ignore unused bits
* Determine p(i), the portion of ones in each block
* Compute the ChiSq: 4 * m * sum((p(i) - 1/2)^2,i in 1..N )
* Compute p-value:  incomplete gamma(N/2,ChiSq/2). 
  Note: incomplete gamma is defined in section 5.5.3 (page 99) as 
       Gamma[N/2,ChiSq/2] / Gamma[N/2].
   Expected value for test sequence:  0.801252.

testFrequencyBlock[sequence_,m_] := Module[{len,part,p,nn, testStat,pVal},
len = Length[sequence];
part = Partition[sequence,m];
p = Total[#]/m & /@ part;
nn = Floor[len/m];
testStat = 4*m*Total@N@Table[(p[[i]] - 1/2)^2, {i,1,nn}];
pVal = Gamma[nn/2,testStat/2]/ Gamma[nn/2]//N
]

testSequence = {0,1,1,0,0,1,1,0,1,0}
testFrequencyBlock[testSequence,3]
If[% < 0.01, Print["Not random based on the Frequency Block test"],Print["It looks random based on the Frequency Block test"]]

Let’s test for all block sizes:

Table[testFrequencyBlock[testSequence,n],{n,1,5}]

And for the larger sequence with block sizes of 1 to 15. Also check if some block size give a p-value < 0.01, but I’m not sure how that would be interpreted.

testFrequencyBlock[sequence,3]
tests = Table[testFrequencyBlock[sequence,n],{n,1,15}]
{#,tests[[#]] < 0.01} & /@ Range[Length[tests]]
Select[%,#[[2]] == True &]

2.3 Runs test (page 2 - 5 f)

Parameters:
* s: the sequence
* n: length of the string

Test statistic:
* vn(obs): the total number of runs

* Compute the pre-test proportion p of ones in the input sequence
* Determine if the prerequisite Frequency test is passed. If not then Runs test is not needed
* Compute the test statistics vn(obs): 
      sum(r(k),k=1..n-1) + 1
      where r(k) = 0 if s(k) == s(k+1) and r(k) = 1 otherwise 
* Compute p-value pval = erfc(  |vn(obs) - 2*n*p*(1-p)/ (2*sqrt(2*n)*p*(1-p))

testRuns[sequence_] := Module[{len,p,vnObs,pVal},
len = Length[sequence];
p = Total[sequence]/len//N;
If[testFrequency[sequence] < 0.01,0,
vnObs = Total@Table[If[sequence[[i]] == sequence[[i+1]],0,1],{i,1,len-1}] + 1;
pVal = Erfc[Abs[vnObs - 2*len*p*(1-p)]/(2*Sqrt[2*len]*p*(1-p) )]
]
]

testSequence = {1,0,0,1,1,0,1,0,1,1}
testRuns[testSequence]
If[% < 0.01, Print["Not random based on the Runs test"],Print["It looks random based on the Runs test"]]

Test on the larger sequence. Yes, it still looks random.

testRuns[sequence]

Now use these test on a non-random sequence. Only the runs test decides that it’s non-random.

nonRandomSeq = Table[Mod[i,2],{i,1,10}]
testFrequency[nonRandomSeq]
testFrequencyBlock[nonRandomSeq,3]
testRuns[nonRandomSeq]

And on a larger non random sequence, same thing: Both the Frequency test and the Frequency Block test accepts it as random, but the Runs test is not fooled!

nonRandomSeq = Table[Mod[i,2],{i,1,100000}];
testFrequency[nonRandomSeq]
testFrequencyBlock[nonRandomSeq,3]
testRuns[nonRandomSeq]

The above larger tests was on the ”Congruential” random method.
Let’s see if the result above is different using a more advanced method (Rule30CA).

SeedRandom[Method->"Rule30CA"];
seqSize = 100000;
sequence2=Table[RandomInteger[], seqSize];

testFrequency[sequence2]
testFrequencyBlock[sequence2,3]
testRuns[sequence2]