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.:
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}))$
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)$
- 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")
head(binProbsList)
