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

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; jblock]++;
}

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\thash{count}\n"; print OUT "count\tweight\tcount\thash{$$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.

• I do not totally understand what you are doing, but have you taken into account that some patterns have a higher variance than others. To take a simple example, if $k=2$ and $m=3$, then the variance of occurrences of 10 is $\frac14$ which is smaller than the variance of occurrences of 00 of $\frac12$ Feb 10, 2020 at 8:53
• I will try to make my question more clear. And yes, you are right about the variance. Among all the 5-bit strings, 00000 and 11111 have the highest variances, followed by 11010 and 01101. Initially, I thought that I discovered a pattern in the binary digits of numbers such as $\sqrt{3}$ but when I run my experiment with random words [i.i.d. Bernouilli$(\frac{1}{2})$], I run into the exact same pattern, confirming (in some ironic way) that there is indeed no pattern in the digits of $\sqrt{3},\sqrt{5}$ and so on, as this is expected in the case of pure randomness. Feb 10, 2020 at 15:06
• The ratio between lowest and highest count variances for 5-bit words seem to be around $\frac{1}{2}$ as in your 2-bit example. Feb 10, 2020 at 15:06
• @Alex: The negative values (2nd column) in my table are the values on the X-axis, while the first column represent the associated percentile. For instance, if you look at the first entry, it reads as $P(Z<-2.36) = 0.01$. Feb 16, 2020 at 7:45
• If you have a moment, could you plot the empirical PDF? Feb 16, 2020 at 23:47

### 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))


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.

• I will do more experiments, but mostly the tails are not the same as for Gaussian. Kurtosis is also highly volatile but always either well above 0, or well below 0, depending on $m$ and $k$. It you use a different seed in the random generator, the mean and variance barely change [the mean is actually constant equal to $(m-k+1)\cdot 2^{-k}$], variance $\approx$ mean $\times (1-2^{-k})$ (before standardization) but the kurtosis is off, volatile, and of course unchanged after standardization. Feb 23, 2020 at 7:27
• There are some real reasons to believe it's not normal: if the counts are strongly correlated. If they were just weakly correlated (or not correlated) you'd end up with Gaussian. Feb 23, 2020 at 7:27
• @VincentGranville it is true that it is not Gaussian, it is more like some over/under dispersed Poisson distribution. But the asymptotic distribution is Gaussian. --- For any counts of events: If the 'distance between events' is independent identical distributed with mean $\mu_d$ and finite variance $\sigma_d^2$, then the distribution of number events $X$ on size of interval $m$ will be asymptotic Gaussian (for increasing $m$), $$X/m \xrightarrow{p} \mathcal{N} \left( 1/\mu_d, \frac{\sigma_d^2}{m\mu_d^3} \right)$$ Feb 23, 2020 at 12:27
• I am not sure whether there is a reputable source for this theorem. My derivation is a bit intuitive and handwaving with not so much rigour, but you can see this step with the Taylor series expansion as sort of analogous to the delta method. Feb 23, 2020 at 12:27