# Time component in the distribution of N elements randomly falling in M bins

I have a system where N elements randomly fall into M bins, so that each bin contains X elements on average. I assume that the distribution of bin occupancies is Poisson(lambda = X).

I want to understand how does the probability of the event [some bin reaches occupancy K] behave, in the course of (or after) L extra actions: a random element (among those currently present in the bins) is revoked, then a new element is added to a random bin (so that the total number of elements in all bins is kept N).

For example, N is 100 and M is 100. Given a static random sample of elements, the probability of any given bin to have occupancy of at least K = 10 could probably be approximated as A = Poisson(lambda = 1, k >= 10) ~= 10-7. Then, the probability that some bin has occupancy of at least 10 is 1 - (1 - A)100 ~= 0.00001, i. e. it is rather low. Then, we perform L = 1 000 000 pairs of actions: 1) revoke one element (selected uniformly at random among the 100 currently present elements) 2) add one new element, and after each such action we check if some bin has occupancy 10. Intuitively, after a million of such actions the probability that at some point we have seen some bin having 10 elements is rather high, probably approaches 1. I want to understand what could be the generic way to express this probability for L = 1, L = 50, L = 1000, etc.

• Your question is unclear, particularly concerning the meaning of "in the course of (or after) ... extra actions." Could you perhaps illustrate it with a small example of what you have in mind? – whuber Dec 29 '18 at 18:37
• @whuber added an example. – leventov Dec 30 '18 at 4:25
• Thank you -- that helps immensely. Could you add just one more detail? Exactly how are elements selected to be "revoked"? In your example there are $N=100$ elements but they are likely to occupy only around $62$ of the $M=100$ bins. So: do you choose occupied bins uniformly at random, or do you choose any bin uniformly at random, or do you choose individual elements uniformly at random? – whuber Dec 30 '18 at 16:12
• @whuber I choose individual elements uniformly at random. – leventov Dec 30 '18 at 20:39
• This, then, is a random walk on a graph. It is solved using Markov chain methods. Analytical solutions are available for relatively small versions of the problem and asymptotic solutions when $L$ is large. It's unclear what you might have in mind with a "generic way to express this probability." – whuber Jan 3 '19 at 20:38

This is neither the most elegant nor the most efficient way, but it should get job done, I think:

Let $$p_0(bin=0),\: p_0(b=1), \ldots,\: p_0(b = N)$$ be the initial probabilities that a bin contains $$1, 2, \ldots,N$$ elements, respectively. By "initial probabilities" I mean the probabilities you get after you distributed your $$N$$ elements for the first time into your $$M$$ bins. Then, the initial probability of having any bin with at least 10 elements is $$1-(1-p_0(b\geq 10))^M$$.

Now, assuming that you didn't find any bin with at least ten elements, you start with your $$L$$ "extra actions". Hence, you want to know the updated probabilities $$p_1(b=0),\: p_1(b=1), \ldots,\: p_1(b = N)$$.

The general formula to calculate these probabilities is: \begin{align} p_{l+1}(b=i) &= p_l(b=i) \\ &-\frac{p_l(b=i)*i}{N}((1-p_l(b=i-1)-p_l(b=i)+2*(p_l(b=i)-\frac{1}{M})) &\quad (1) \\ &+\frac{p_l(b=i+1)*(i+1)}{N}*((1-p_l(b=i-1)-p_l(b=i)-\frac{1}{M})+2*p_l(b=i-1)) &\quad (2) \\ &+\frac{N-(p_l(b=i)*i+p_l(b=i+1)*(i+1))*M}{N^2}*(p_l(b=i-1)-p_l(b=i)) &\quad (3) \\ \end{align}

i.e. we get the probability of a bin having $$i$$ elements after the first $$l$$ "extra action" by looking at the potential changes to the the probability of a bin having $$i$$ elements in the last round $$p_l(b=i)$$. The following three cases can change this probability.:

1. We pick an element from a bin with $$i$$ elements ($$\frac{1}{N}*p_l(b=i)*i$$) and distribute that element into any bin other than the one it came from or any bin other than the ones with $$i-1$$ elements or $$i$$ elements, i.e. $$(1-p_l(b=i-1)-p_l(b=i)$$. In these cases, we lose one bin with $$i$$ elements. In addition, if distribute this element into a bin with $$i$$ elements (except for the one were we took the element from), we lose two bins with $$i$$ elements, i.e. $$2*(p_l(b=i)-\frac{1}{M}))$$

2. We select an element from a bin with $$i+1$$ elements $$\frac{1}{N}*p_l(b=i+1)*(i+1)$$. If we distribute this element into any bin other than bins with $$i-1$$ or $$i$$ elements or the bin we took it from, we get an additional bin with $$i$$ elements, i.e. $$(1-p_l(b=i-1)-p_l(b=i)-\frac{1}{M}$$. If we distribute the element into a bin wih $$i-1$$ elements, we get two additional bins with $$i$$ elements: $$2*p_l(b=i-1)$$

3. We select an element from any bin other than bins with $$i$$ or $$i+1$$ elements, i.e. $$\frac{1}{N^2}*(N-(p_l(b=i)*i+p_l(b=i+1)*(i+1))*M)$$. If we then distribute this element into a bin with $$i-1$$ elements, we gain a bin with $$i$$ elements, i.e. $$p_l(b=i-1)$$. If we distribute it to a bin with $$i$$ elements, we lose a bin with $$i$$ elements, i.e. $$-p_l(b=i)$$

Note that after the first "extra action", the probabilities $$p_l(b = 10)$$ and beyond can be disregarded for the calculation, because you stated in the question that you always check for any bin with at least 10 elements after each "extra action" $$l$$. If you haven't found any such bin and initialised the next "extra action" $$l+1$$, we know that $$p_{l}(b >=10)=0$$. Hence, we also need to adjust the other probabilities simply by dividing each probability by the sum of the remaining probabilities, i.e. $$p_l(b=i) = p_l(b=i)/\sum_{i=0}^9 p_l(b=i)$$.

Now, to calculate the probability $$p_{l+1}(b\geq 10)$$, we need to take the inverse probability of transitioning from a bin with 9 elements to a bin with 10 elements in the $$(l+1)^{\text{th}}$$ "extra action" round, i.e. $$1-p_{l+1}(b\geq 10) = 1-p_{l+1}(b = 10) = 1-(\frac{1}{N^2}*(N-p_l(b=9)*9*M)*p_l(b=9) + \frac{1}{N}*p_l(b=9)*9*(p_l(b=9)-\frac{1}{M})$$

We then need to multiply this probability with the respective inverse probabilities of the previous rounds and take the inverse of that:

$$p_L(b\geq 10) = 1-\prod_{l=0}^{L-1} (1-p_{l+1}(b\geq10))$$

Because of this, the probability of $$p_L(bin\geq=10)$$ indeed approaches one, as the following quick-and-dirty R-simulation shows.

# Edit:

I found two error in my calculations. 1. I forgot to adjust the probabilities for the number of elements in each bin. I corrected the calculations above and the code below. 2. I forgot to adjust for the fact that after each extra action, $$p_l(b\geq 10)$$ is set to zero, so we need to adjust the other probabilities for this as well.

M <- 100
N <- 100
lambda <- 1
binContents <- seq(0,10)
L <- 100

calcBinProb <- function(binContent){
prob <- dpois(binContent,lambda)
}
binProbs0 <- sapply(binContents,calcBinProb) #calculate probabilities of bins having a specific value

updateBinProbs <- function(binProbs,M,N){
from_i <- function(i){
j = i-1
prob <- binProbs[i]*j/N*((1-binProbs[i-1]-binProbs[i])+2*(binProbs[i]-1/M))
}
from_i_plus_1 <- function(i){
j = i-1
prob <- binProbs[i+1]*(j+1)/N*((1-binProbs[i-1]-binProbs[i]-binProbs[i+1])+(binProbs[i+1]-1/M)+2*binProbs[i-1])
}
from_rest <- function(i){
j = i-1
prob <- ((1-binProbs[i]*j*M/N-binProbs[i+1]*(j+1)*M/N)/N)*(binProbs[i-1]-binProbs[i])
}

updatedBinProbs <- rep(0,length(binProbs))
binProbs[length(binProbs)] <- 0
for (i in 1:length(binProbs)){
if (i == 1){
updatedBinProbs[i] <- binProbs[i]+binProbs[i+1]*1/N*((1-binProbs[i]-binProbs[i+1])+(binProbs[i+1]-1/M))-((1-binProbs[i+1]*1*M/N)/N)*(binProbs[i])
}
else if (i<length(binProbs)){
updatedBinProbs[i] <- binProbs[i]-from_i(i)+from_i_plus_1(i)+from_rest(i)
}
else if (i==length(binProbs)){
j <- i-1
updatedBinProbs[i] <-((1-binProbs[i-1]*(j-1)*M/N)/N)*(binProbs[i-1]) + binProbs[i-1]*(j-1)/N*(binProbs[i-1]-1/M)
}
}
return(updatedBinProbs/sum(updatedBinProbs[1:length(binProbs)-1]))
}

len <- length(binProbs0)
p_geq_10 <- rep(0,L)
binProbsList <- vector("list",N)
binProbsList[[1]] <- binProbs0
p_geq_10[1] <- 1-sum(binProbs0[1:len-1])
for (i in 2:L){
binProbs <- updateBinProbs(binProbs,M,N)
binProbsList[[i]] <- binProbs
p_geq_10[i] <-1-prod(1-p_geq_10[1:i-1])*(1-binProbs[len])
}
p_any_geq_10 <- 1-(1-p_geq_10)^M

plot(1:L,p_geq_10,xlab="L",ylab="1-(1-p(bin>=10))^M",type="l",ylim=c(0,1))
lines(1:L,p_any_geq_10,col="red")