# Algorithm to estimate the probability that a "0 and 1 matrix" fills up following the bootstrap percolation rules

Presentation of the model: we consider the regular lattice created from $\mathbb{Z}^2$ (It's no more, no less a square lattice).

• At $t=0$, each site is said "active" independently with a probability $p$, "inactive" otherwise.
• At $t$, if a site is inactive, it becomes active if at least $2$ of its neighbours are active.

We call neighbour of $x$ each element of $V(x)=\{y\in\mathbb{Z}^2, ||x-y||=1\}$ (or, if I have a matrix M in python, a neighbour of $M[i][j]$ is $M[i\pm1][j]$ or $M[i][j\pm 1]$).

My problem: I consider a $n\times n$ matrix filled with 0s and 1s, (0 = inactive ; 1 = active), and I want to estimate the probability $p_c$ that a matrix, following the rules above, fills up (at $t\rightarrow \infty$, but this happen really fast in reality). The probability $p_c$ depend on $n$. I want to test it for $n=2$ to $n=100$ at least.

What I already have: I have a program, matrix_creation(n,p) that creates a $n\times n$ matrix with the probability $p$ for each coefficient to be 1, and another, test(M), to test if "M" fills up.

• In order to estimate the probability of the event empirically, simply generate $N$ matrices for large $N$, calculate how many of them "fill up", say $M$, and estimate the probability by $M/N$. You can even estimate the error in this estimate (i.e. calculate confidence intervals). Commented Sep 17, 2017 at 12:05
• The problem is that the generation of the matrix depends on $p$ I choose before. But maybe I can generate a random $p$ and test as you said? I'll, but it seems not to be really accurate, is it? Commented Sep 17, 2017 at 12:11
• Perhaps you can calculate the probability, but it seems you want to estimate it empirically. Using the method I described, you can get a $\pm\epsilon$ approximation using roughly $1/\epsilon^2$ queries. Commented Sep 17, 2017 at 12:14
• I did what you said, it works but the complexity is so high, and it's not really accurate (I mean, you generate random things, and test a number of them, so if I want it to be accurate, I have to generate a great number of them, at least 1000, and it's really long... Maybe my first algorithms are bad too...). Maybe using dichotomy? But idk how to do it with random things, what criteria to stop etc... Commented Sep 17, 2017 at 12:32
• Random simulations are very useful in practice. Sometimes there are tricks that speed up the estimation. I suggest optimizing your program, and using a fast language like C rather than a slow language like python. Commented Sep 17, 2017 at 12:37

Here is a reasonably fast approach to estimate $p_c$. We will perform $N$ different trials, count the number of successes $M$, and estimate the result as $M/N \pm 1/\sqrt{N}$ (better estimates are possible if $p_c$ is close to 0 or to 1).
For each trial, you initialize an $N \times N$ matrix with the initial filling. You then go over the array, calculating the number of active neighbors of each non-filled point in the matrix, and whenever that number reaches 2, you put it in a stack. You then repeatedly pop elements from the stack, updating the matrix of active neighbors, and whenever that number reaches 2 for some non-filled point, you add it to the stack. When the stack is empty, the process is complete, and you can check whether all points have been activated.
There are many optimizations possible: for example, you can use an $(N+2)\times(N+2)$ matrix to help deal with border cases, you can keep track of the number of active points, and so on.
Using this naive approach, I managed to perform $10^6$ trials for $n = 100$, $p = 0.0535$ in about two minutes on my PC. The answer was $M/N \approx 0.5008$, so the critical probability is pretty close to $0.0535$.