Limit distribution associated with counts (non-trivial combinatoric problem)

Question

Update on 3/1/2020. All the material below and much more has been incorporated into a comprehensive article on this topic. The question below is discussed in that article, entitled "State-of-the-Art Statistical Science to Tackle Famous Number Theory Conjectures", and available here.

Let $W$ be a word (also called block) consisting $k$ binary digits. Let $S$ be a sequence (also called text or book) consisting of $m$ binary digits, with $k\leq m$. Let $N_{W,S}$ be the number of occurrences of $W$ in $S$. For instance, if $S=010001010011$ and $W=00$, then $N_{W,S} = 3$.

Here $k$ is small and $m$ is large: $k=5$ and $m=20,000$ in my experiment.

For a positive integer $x$, a block $W$ of length $k$ and a random sequence $S$ of length $m$, the number of occurrences of the event $\{N_{W,S} = x\}$ is denoted as $P(N_{W,S} = x)$. So $x$ can be viewed as the realization of a discrete random variable $X$. In particular, $$\mbox{E}(X) = \frac{m-k+1}{2^k}.$$

Let $Z = (X-\mbox{E}(X))/\sqrt{\mbox{Var}(X)}.$

Question: what is the limiting distribution of $Z$, as $m\rightarrow\infty$?

Using simulations, I got a pretty decent approximation. Below is its empirical distribution:

It is perfectly smooth and symmetric at the limit, but the excess kurtosis is $0.63$, so it is not a normal distribution. The empirical percentile distribution of $Z$ is in the table below, maybe there is an almost perfect fit with some known distribution with 0 mean and unit variance.

  P(Z < x)    x
    0.01    -2.36
    0.02    -2.03
    0.03    -1.87
    0.04    -1.71
    0.05    -1.62
    0.06    -1.50
    0.07    -1.46
    0.08    -1.38
    0.09    -1.30
    0.10    -1.26
    0.11    -1.18
    0.12    -1.13
    0.13    -1.09
    0.14    -1.05
    0.15    -1.01
    0.16    -0.97
    0.17    -0.93
    0.18    -0.89
    0.19    -0.85
    0.20    -0.81
    0.21    -0.77
    0.22    -0.77
    0.23    -0.73
    0.24    -0.69
    0.25    -0.65
    0.26    -0.60
    0.27    -0.60
    0.28    -0.56
    0.29    -0.52
    0.30    -0.52
    0.31    -0.48
    0.32    -0.44
    0.33    -0.44
    0.34    -0.40
    0.35    -0.36
    0.36    -0.36
    0.37    -0.32
    0.38    -0.28
    0.39    -0.28
    0.40    -0.24
    0.41    -0.24
    0.42    -0.20
    0.43    -0.16
    0.44    -0.16
    0.45    -0.11
    0.46    -0.11
    0.47    -0.07
    0.48    -0.03
    0.49    -0.03
    0.50    0.01
    0.51    0.01
    0.52    0.05
    0.53    0.05
    0.54    0.09
    0.55    0.13
    0.56    0.13
    0.57    0.17
    0.58    0.17
    0.59    0.21
    0.60    0.25
    0.61    0.25
    0.62    0.29
    0.63    0.33
    0.64    0.33
    0.65    0.37
    0.66    0.37
    0.67    0.42
    0.68    0.46
    0.69    0.46
    0.70    0.50
    0.71    0.54
    0.72    0.54
    0.73    0.58
    0.74    0.62
    0.75    0.66
    0.76    0.66
    0.77    0.70
    0.78    0.74
    0.79    0.78
    0.80    0.82
    0.81    0.82
    0.82    0.86
    0.83    0.91
    0.84    0.95
    0.85    0.99
    0.86    1.03
    0.87    1.11
    0.88    1.15
    0.89    1.19
    0.90    1.23
    0.91    1.31
    0.92    1.39
    0.93    1.44
    0.94    1.52
    0.95    1.64
    0.96    1.72
    0.97    1.88
    0.98    2.09
    0.99    2.46

If instead of one sequence $S$, you consider $n$ random sequences $S_1,\cdots,S_n$ all of same length $m$, and independent from each other, then the variance for the counts $N_{W,S}$, computed across all sequences bundled together, satisfies $$\mbox{Var}(X)\rightarrow\frac{m-k+1}{2^k}\cdot \Big(1-\frac{1}{2^k}\Big) \mbox{ as } n\rightarrow\infty.$$ This result can be used to test if sequences found in actual data sets are both random and independent from each other.

The challenge

The problem is that the successive $m-k+1$ blocks $W$ of length $k$ do overlap in any sequence $S$ of length $m$, resulting in lack of independence between the various counts $N_{W,S}$. If the blocks (and thus their counts) were independent instead, then the counts would follow a multinomial distribution, with each of the $n\cdot (m-k+1)$ probability parameters being $\frac{1}{2^k}$, and $Z$ would be asymptotically normal. Here this is not the case: the excess kurtosis does not converge to zero. There is convergence to smooth, symmetrical distributions as $n$ and $m$ increase, but that limit is never Gaussian. My big question is: what is it then?

That said, for the first two moments (expectation and variance) attached to $N_{W,S}$, we get the same values (at least asymptotically) as those arising from the multinomial model. But not anymore for higher moments.

The following code performs simulations and computes the variances, expectations, kurtosis and all the counts $N_{W,S}$. Note that the variance and kurtosis depend on $S$, but they stabilize as $n$ is increasing. The expectation depends only on $m$ and $k$.

use strict;

my $k;
my $k1;
my $k2;
my $j;
my $rand;
my $m;
my $even;
my $block;
my @digits;
my @ablock;
my @biglist;
my $bigstring;
my $nbigstrings;
my $binary;
my %hash;
my %hlist;
my @blocksum;
my $tweight;
my $sum;
my $sum2;
my $avg;
my $var;
my $kurtosis;
my $num;

my $count;

$rand=500;
$k1=5; # bits of small word
$k2=2**$k1;
$m=7; # bits in big string  # m > k1 otherwise var = 0
$nbigstrings=5000; # number of sampled big strings 

open(OUT2,">collatzr.txt");      

@biglist=();
%hlist=();

for ($bigstring=0; $bigstring<$nbigstrings; $bigstring++) { 

  @digits=();
  @ablock=();
  $binary="'";

  for ($k=0; $k<$m; $k++) { # compute 200 digits
    $rand=(10232193*$rand + 3701101) % 54198451371;
    $even=int(2*$rand/54198451371);     
    @digits[$k]=$even;  
    $binary=$binary."$even";      
  }
  print OUT2 "\n$binary\n";

  for ($k=0; $k<$m-$k1+1; $k++) { ## kmax - 5
    $block="";
    for ($j=0; $j<$k1; $j++) {
      $block+=($digits[$k+$j]* 2**$j);
    }
    $ablock[$block]++;
  }

  if ($bigstring%1000 == 0) { print "iter... $bigstring\n"; select()->flush(); }

  for ($block=0; $block<$k2; $block++) {
    if ($ablock[$block] eq "") { $ablock[$block]=0; }
    $count=$ablock[$block];
    $hash{$count}++;   #{$ablock[$block]}++;  # number of occurences of $count (used as weight in AVG, VAR)
    $blocksum[$block]+=$count;
    $hlist{$block}=$hlist{$block}."\t$count"; # disuse if it uses too much memory
    print OUT2 "$block\t$count\n";
  }
}

close(OUT2);

#-- summary stats

open(OUT,">coll2.txt");

$tweight=0;
$sum=0;
$sum2=0;
$kurtosis=0;

foreach $count (keys(%hash)) {
  $tweight+=$hash{$count};
  $sum+=$count*$hash{$count};
  $sum2+=$count*$count*$hash{$count};
  print "count weight: $count\t$hash{$count}\n";
  print OUT "count\tweight\t$count\t$hash{$count}\n";
}

$avg=$sum/$tweight; 
$var=($sum2/$tweight)- $avg*$avg;

foreach $count (keys(%hash)) {
  $kurtosis+=$hash{$count}*(($count - $avg)/sqrt($var))**4;
}
$kurtosis = -3+$kurtosis/$tweight;
$num = $avg*$k2;

print "($k1 | $m | $nbigstrings) avg ~ sum2| var | excess_kurt | tweight | missing : $avg ~ $sum2 | $var | $kurtosis | $tweight | $hash{0}\n";


for ($block=0; $block<$k2; $block++) {
  # print "block: $block\t$blocksum[$block]\n";
  print OUT "block\tblocklist\t$block\t$hlist{$block}\n";
}

close(OUT);

Context

I am checking if all blocks of $k=5$ binary digits are distributed as expected (that is, randomly) in the first $m$ binary digits of a bunch of quadratic irrational numbers. To test my hypothesis that this is the case, I need to know the exact distribution of the test statistic for the null hypothesis. The exact distribution is the distribution attached to $Z$. More about this project can be found on Math.StackExchange, here.

Answer 1

Correspondence to a waiting time problem

There is an alternative way to look at this. We can switch between the 'number of words' as being the variable and the 'text size' as being the variable.

Imagine assembling the text or book untill you reach some fixed number, $x$, of words. Then consider the length of the text, $m$, as the variable.

We can relate the two situations by:

$$P(X < x|m) = P(M > m| x)$$

The 'probability that a text of fixed size $m$, has less than $x$ words' equals the 'probability that you need a text of more than size $m$ in order to reach the fixed number of $x$ words'. See for a similar use of this relation here: Does exponential waiting time for an event imply that the event is Poisson-process?

In this alternative way we can consider the problem as a waiting time problem. We describe the distribution of text size $M$ for a given number of words $m$ as a sum of waiting times for the individual words. This is analogous to many other distributions:

Distribution of                 Distribution of 
Waiting time between events     number of events

Exponential                     Poisson
Erlang/Gamma                    over/under-dispersed Poisson  
Geometric                       Binomial
Negative Binomial               over/under-dispersed Binomial

• For instance: the simplest case is the number of occurrences of the words 10 or 01. The waiting time until the word appears (or the distance between the words) relates to the event that we observe a the first character and then the second character, that is we are waiting untill we get twice a 'good' flip. In the case of 10 we are first waiting for a 1 and then for a 0. The waiting time is distributed as $P(K=k) = (k-1)0.5^k$ (this is a Negative Binomial distribution, that you may also see as a sum of two geometric distributed variables).

  The corresponding distribution of the number of words for a given text size relates to the binomial distribution. It relates to the number of times the binary string is changing character. E.g. the string 0100101 changed sign 5 times and half of those changes relate to the word 10. The distribution of these changes are binomial. And the number of words 10 will be half of the number of changes (roughly, you need to round of and take into account whether the text starts with a 0 or 1).

Asymptotic normal distribution for waiting time

So this waiting time for $x$ words can be seen as the sum of $x$ independent identical variables of waiting time for $1$ word. (The waiting time for the first word may in some cases be different distributed. But otherwise, when 'waiting' for all the other words, the waiting time is identical distributed. Because every-time you start waiting for the next occurrence of a word has the identical end of the text, namely the word itself).

If the waiting time for a single event has finite variance then the conditions of the CLT are fulfilled and the waiting time for more events will be asymptotically normal distributed. Then we can model the distribution of $m$ as an approximately normal distribution :

$$P(M \leq m \vert x) \approx \Phi \left( \frac{m-x\theta_1}{\sqrt{x}\theta_2} \right)$$

Note that the CLT is more precisely $\sqrt{x}(M/x-\theta_1) \sim N(0,\theta_2)$ so this is where the distribution of $M$ and not $\sqrt{x}(M/x-\theta_1)$ may still differ from a normal distribution. But the scaled variable $M/x$ is approaching a normal dsitribution.

and we can convert this into

$$P(X < x \vert m) = P(M > m \vert x) \approx \Phi \left( \frac{x-m/\theta_1}{\sqrt{x}\theta_2/\theta_1} \right) \underbrace{\approx \Phi \left( \frac{x-m/\theta_1}{\theta_2\sqrt{m/\theta_1^3}} \right)}_{\text{Taylor series approximation}} $$

where $\theta_1$ and $\theta_2$ are respectively the mean and standard deviation of the waiting time for a single event.

Computational model and relation with diffusion/chromatography

Below I demonstrate that the waiting time for this diffusion may be modeled (approximately) with a Gamma distribution. This can be converted to a distribution in space (the cumulative distribution of waiting time relates to the cumulative distribution in space) which resembles an over-dispersed Poisson distribution.

(Possibly there might be some non continuous function to model this more exactly, e.g. the negative binomial for the words 10 and 01. However the Gamma distribution does it pretty well)

This diffusion process resembles somewhat a chromatography process (this is what I first thought about when trying to tackle this problem), possibly there might be more detailed descriptions in the literature relating to chromatography.

Example:

For the case where we look for the number of times the word 11 appears in a text of size $k$ we describe two variables:

• $X_k(x)$: the probability that the number of times that 11 appears in a text of size $k$ is $x$ and the last character is 1
• $Y_k(x)$: the probability that the number of times that 11 appears in a text of size $k$ is $x$ and the last character is 0

Then we can describe the evolution iteratively

$$\begin{array}{} X_k(x) &=& 0.5*X_{k-1}(x-1) &+& 0.5*Y_{k-1}(x) \\ Y_k(x) &=& 0.5*X_{k-1}(x) &+& 0.5*Y_{k-1}(x) \end{array}$$

We can model this computationally. In addition, the arrival time seems to be modeled as an Erlang/Gamma distribution as well.

n <- 202-2   #number of itterations (text size-2)
Xk <- matrix(rep(0,n*(n+1)),n+1) 
Yk <- matrix(rep(0,n*(n+1)),n+1)

# start
Xk[1,1] <- 1   #01
Xk[2,1] <- 1   #11
Yk[1,1] <- 2   #00 or 10

# evolution
for (run in 2:n) {
  Xk[1,run] <- Yk[1,run-1] 
  Yk[1,run] <- Yk[1,run-1] + Xk[1,run-1]
  for(place in c(2:(n+1))) {
    Xk[place,run] <- Xk[place-1,run-1]+Yk[place,run-1]
    Yk[place,run] <- Xk[place,run-1]+Yk[place,run-1]
  }
}


# arrival time
arr <- rep(0,n)
xp <- 20  #position
for (i in 1:n) {
  # how many are already in position >=xp 
  arr[i] <- sum(Xk[c((xp+1):(n+1)),i]+Yk[c((xp+1):(n+1)),i])/2^(i+1)
}


# plot arrival
x <- 3:(n+1)
plot(x,diff(arr),log="y")

# erlang/gamma distribution with same mean and variance
mu <- sum(x*diff(arr))
vr <- sum((x-mu)^2*diff(arr))
scale <- vr/mu
shape <- mu/scale
lines(x,dgamma(x,shape=shape, scale=scale),col=3)
shape*scale
shape*scale^2

We can relate this with the waiting time to make an additional step (given that the last character is 1, what is the distribution of the number of characters that we need to add before a new word arises)

For the word 11 we can compute this with reasonable accuracy:

# computing waiting time for 1 step
# mean and variance
x <- rep(0,30)
y <- rep(0,30)
f <- rep(0,30)
x[1] <- 1
for (i in 2:30) {
  y[i] <- x[i-1]*0.5+y[i-1]*0.5
  x[i] <- y[i-1]*0.5
  f[i] <- x[i-1]*0.5
}
plot(f, log="y")
mwt <- sum(t*f)
vwt <- sum((t-mwt)^2*f)

# compare means and variance
mwt*xp
mu
vwt*xp
vr

So the waiting time seems to be Gamma distributed with mean and variance a multiple times the mean and variance of a single step.

Sidenotes:

1. I have not yet a proof of this approximation.
2. This multiple factor seems to be not very exact. I am not sure yet why this is the case.

3. For more complex words it becomes more difficult to model/compute the waiting time, but potentially you could determine it experimentally by fitting with a gamma distribution.

   

   ### converting from waiting time to distribution in space
   
   
   time <- 50
   k <- 1:time
   mean=k*mwt
   var=k*vwt
   Gk <- 1-pnorm(time,mean=mean,sd=sqrt(var))
   Fk <- 1-pgamma(time,shape=mean^2/var, scale = var/mean)
   plot(Xk[,time]/2^time, log="y", xlim=c(0,time),
        ylab = "density, mass", xlab = "occurences")
   lines(k[-1],diff(Fk),col=2)
   lines(k[-1],diff(Gk),col=3)
   
   sc <- sqrt(vwt/mwt^2)
   overk <- sc*(k-time/4)+time/4
   lines(overk,dpois(k-1,time/mwt)/sc,col=4)
   
   
   legend(5,10^-10,c("distribution of occurences word '11' in text of size 52",
                      "Erlang distributed waiting time",
                      "Gaussian distributed waiting time",
                      "Overdispersed Poisson"),cex=0.7,
          pch=c(21,-1,-1,-1),lty=c(0,1,1,1),col=c(1,2,3,4))
   

Regarding your experimental distribution

When I compare your experimental distribution with a normal distribution then I get something that resembles a normal distribution very much:

Z <- c(-2.36,-2.03,-1.87,-1.71,-1.62,-1.5,-1.46,-1.38,-1.3,-1.26,-1.18,-1.13,-1.09,-1.05,-1.01,-0.97,-0.93,-0.89,-0.85,-0.81,-0.77,-0.77,-0.73,-0.69,-0.65,-0.6,-0.6,-0.56,-0.52,-0.52,-0.48,-0.44,-0.44,-0.4,-0.36,-0.36,-0.32,-0.28,-0.28,-0.24,-0.24,-0.2,-0.16,-0.16,-0.11,-0.11,-0.07,-0.03,-0.03,0.01,0.01,0.05,0.05,0.09,0.13,0.13,0.17,0.17,0.21,0.25,0.25,0.29,0.33,0.33,0.37,0.37,0.42,0.46,0.46,0.5,0.54,0.54,0.58,0.62,0.66,0.66,0.7,0.74,0.78,0.82,0.82,0.86,0.91,0.95,0.99,1.03,1.11,1.15,1.19,1.23,1.31,1.39,1.44,1.52,1.64,1.72,1.88,2.09,2.46)
p <- seq(0.01,0.99,0.01)
plot(p,Z, cex=0.5, pch=21,col=1,bg=1)
lines(p,qnorm(p),col=2)

I wonder whether the observed excess Kurtosis is something peculiar (From the image I actually notice a negative excess Kurtosis, but the same is true for the variance. The points are a bit closer to 0 than the red curve which are values for a normal distribution. So maybe the points in the tails, which are not in these 99 points, are important here.).

There is convergence to smooth, symmetrical distributions as $n$ and $m$ increase, but that limit is never Gaussian.

You mention that you do not observe a tendency to approach a normal distribution when you increase $m$. But you will get a sort of (over/under) dispersed Poisson distribution with a mean $m/2^k$. Only when this mean is sufficiently increasing will the distribution approach a normal distribution.