Postprocessing Techniques for Gradient Percolation Predictions on the Square Lattice

John Tencer and Kelsey Meeks Forsberg
Physical Review E (2021)
View at Publisher

Abstract

In this work, we revisit the classic problem of site percolation on a regular square lattice. In particular, we investigate the effect of quantization bias errors on percolation threshold predictions for large probability gradients and propose a mitigation strategy. We demonstrate through extensive computational experiments that the assumption of a linear relationship between probability gradient and percolation threshold used in previous investigations is invalid. Moreover, we demonstrate that, due to skewness in the distribution of occupation probabilities visited the average does not converge monotonically to the true percolation threshold. We identify several alternative metrics which do exhibit monotonic (albeit not linear) convergence and document their observed convergence rates.

Introduction

For statistical mechanics problems dealing with transport properties and particle connectivity, percolation theory is an important resource in predicting composite behavior and dispersion in random media and provides a tool for linking microstructure and macroscopic material properties [1]. It is often described in terms of the critical parameter at which bulk connectivity is established, the percolation threshold $p_c$. Below the percolation threshold, large connected components do not exist.

Percolation is a well-studied physical phenomena because of its broad applicability, including the physical percolation of fluids through rock [2–4], as well as resistor networks [5], disease spread [6], and many problems in material science [7,8]. Studies of these phenomena often focus on either lattice or continuum systems. Lattice percolation is described by regular or irregular networks, where sites or bonds are occupied with some probability p, and occupied sites form connected pathways. Here we will focus on site percolation on a regular square lattice.

For certain lattice systems, such as bond percolation on a square lattice or site percolation on a triangular lattice, the percolation threshold may be determined analytically. However, for many other lattice systems the percolation threshold must be estimated numerically. Many techniques have been developed over the years for evaluating the percolation threshold in lattice systems including hull gradient [9–16], planar crossing [17–21], histogram Monte Carlo [22], invaded cluster algorithms [23,24], toroidal wrapping [4,25,26], cylindrical correlation [27], dynamic programming [28], and transfer matrices [29].

For simulations on finite-sized lattices the results must be extrapolated to an infinitely large lattice. This is typically done by taking advantage of known critical exponents for the given universality class [20,25,30]. For gradient percolation, an analogous problem presents itself in the need to extrapolate results to a gradient-free lattice. Unfortunately, this convergence rate is not known to be related to the known critical exponents.

It has been common practice when estimating the percolation threshold via gradient percolation to evaluate to high statistical precision \(p_c^*\) for decreasing values of $\nabla p$, where \(p_c^*\) is the average value of $p$ sampled during the hull walk. For large values of $\nabla p$ these have been observed to converge approximately linearly to $p_c$ [11,16]. However, they do not converge precisely linearly, and additional unconsidered sources of error can complicate the estimation of the true convergence rate. Consider, for example, the results in Figure 1. On a linear scale, the percolation threshold \(p_c^*\) appears to converge linearly to the known value $p_c$ as $\nabla p \rightarrow 0$. However, nonlinear behavior is occurring for small probability gradients. In fact, when examined more closely (insert), the percolation threshold is observed to overshoot the known value suggesting that data commonly used for extrapolation is outside the asymptotic regime.

In this work, we consider extremely small occupation probability gradients (as low as $655360^{-1}$). Smaller occupation probability gradients should provide more accurate estimates of the percolation threshold if commonly made assumptions regarding convergence hold. Previous work using the gradient hull method to predict percolation thresholds has used occupation probability gradients as low as $320^{-1}$ [11], $4000^{-1}$ [12], $10000^{-1}$ [14], $12800^{-1}$ [10], $16400^{-1}$ [16], and $100000^{-1}$ [15]. A linear extrapolation of $\nabla{p}\rightarrow0$ has been used in all cases. In this work, we consider significantly smaller occupation probability gradients and demonstrate that the assumed linear convergence rate as $\nabla{p}\rightarrow0$ is incorrect.

Methods

Frontier Walk

We calculate percolation thresholds using the gradient percolation method [12]]. In this method, lattice sites in the $x$-$y$ plane are occupied with a probability $p=p(y)$, where $p$ is chosen to vary linearly with $y$. The occupation probability gradient creates an extended frontier between the percolating and non-percolating regions. The resulting frontier is roughly perpendicular to the gradient of $p$ (parallel to the x-axis). The frontier is generated and traversed using a self-avoiding hull-generating walk, according to the rule that an occupied site will be traversed by the walk and a unoccupied site will reflect the walk as in [9]. Periodic boundary conditions are used on the right and left edges of the domain to reduce the size of the computational domain without limiting the length of the walk.

Each simulation begins with an empty blank lattice. Sites are neither occupied nor unoccupied until encountered by the walk with the exception of sites in the left- and right-most columns which are initialized to insure that the walk proceeds from left to right. The right-most column is initialized as unoccupied while the left-most column is initialized as half occupied and half unoccupied. The walk begins in the middle of the left-most column at the first occupied site facing away from the last unoccupied site.

The walk proceeds by looking to the left and determining if that site is occupied or unoccupied according to a random draw and the occupation probability of that site. The $i^{th}$ row of sites on the lattice has a single occupation probability given by $p_i=p(y_i)=p(y_0) + i\nabla p$. If the adjacent site is occupied, the walk advances to the new site updating its direction. If the adjacent site is unoccupied, the walk direction is rotated 90\textdegree to the right and the process is repeated. Once the walk reaches the right edge of the computational domain, periodic boundary conditions are used to wrap it back to the left side. Sites one column to the right of the current walk position are progressively reset to blank allowing the walk to backtrack up to the width of the computational domain The width is required to be sufficiently large that the walk never backtracks to a column which has already been reset. The walk is determined to have made one pass through the domain when the left-most column is reset.

A single simulation consists of 501 passes through the domain, with the data from the first pass being discarded to eliminate any effect from the initialization procedure. The width of the walk scales as $\nabla{p}^{4/7}$ and the height of the computational domain is adjusted to always be approximately an order of magnitude larger than the expected width of the walk so that the walk never impacts the top or bottom boundary. The width of the computational domain is set to $2048$ for $\nabla{p}\geq 1/40960$, $4096$ for $1/40960>\nabla{p}\geq1/163840$, and $8192$ for smaller gradients. The width of the domain is increased for especially small gradients to guarantee that the width of the computational domain is always greater expected hull width.

The parameters of each simulation are the occupation probability gradient $\nabla p$ and the minimum occupation probability $p_0=p(y_0)$.
Each simulation produces a count, $N_i$ for the number of times each row was visited during the simulation. The quantity $\sum_i \frac{N_i p_i}{N_{total}}$ with $N_{total}=\sum_i N_i$ has been used as an approximation for the percolation threshold [11,12,16].

Quantization Bias Errors

The occupation probability of the sites visited is approximately normal \begin{equation} p \sim \mathcal{N}(\mu,\sigma^2). \end{equation} The average value of $p$ sampled during the walk is an estimate of the percolation threshold $p_c$ [11]. However, it is generally a biased estimate because the discrete nature of the lattice implies that, in general, $\sum_i \frac{N_i p_i}{N_{total}} \nrightarrow \mu$, especially for large probability gradients.

Within a given simulation, only $p \in p_{att}=\lbrace p_0, p_1, …\rbrace$ values of the occupation probability may be observed, where $p_i=p_0 + i\nabla p$. Because only occupation probabilities corresponding to integer lattice rows may be observed, the distribution of the observed occupation probabilities becomes \begin{multline} pdf(p) = w(\mu,\sigma^2,p_i)\delta(p-p_i),
i\in\lbrace 0,1, \cdots , N_{rows}-1\rbrace \end{multline} with \begin{multline} w(\mu,\sigma^2,p)=\int_{p-\frac{1}{2}\nabla p}^{p+\frac{1}{2}\nabla p} pdf(\mu,\sigma^2,\tilde{p})d\tilde{p} \ =\Phi(\mu,\sigma^2,p+\frac{1}{2}\nabla p)-\Phi(\mu,\sigma^2,p-\frac{1}{2}\nabla p), \label{eq:pmf} \end{multline} where $w$ is the mass probability function, $\delta$ is the delta function, and $\Phi$ is the cumulative probability distribution function. Using the discrete probability distribution function, the expected mean $m$ is determined to be \begin{equation} m(\mu,\sigma^2, p_0)=\sum_i w(\mu,\sigma^2,y_i)p_i \approx \frac{1}{N_{total}\nabla p}\sum_i N_i p_i. \label{eq:m_discrete} \end{equation} The bias introduced by the discrete lattice $\lvert \mu - m\rvert$ is a function of the minimum occupation probability $p(y_0)$, with $\lvert \mu - m\rvert = 0 \Leftrightarrow \lbrace \mu, \mu + \nabla p/2\rbrace \cap p_{att} \neq \emptyset$.

The bias error introduced by the use of a discrete lattice is analgous to bias errors due to finite resolution measurements [31]. The results of the simulations are more rightly viewed as samples from this discrete distribution than the underlying continuous distribution for $p$ (although for sufficiently small values of $\nabla{p}$, this distinction becomes unimportant).

Dithering

Dither is the intentional application of noise (with the intent of randomizing quantization error) commonly used in the processing of digital images [32,33] and audio recordings [34,35]. Here, we apply dither by randomly varying $p_0$ in order to randomize the bias errors introduced by the use of a discrete lattice with nonzero $\nabla p$. Without dither, $p_0$ would be fixed, i.e., $p_0=\bar{p_0}$. With dither $p_0 \sim \mathcal{U}(\bar{p_0},\bar{p_1})$.

Because the expected value of the discrete distribution, $m$ is a function of $p_0$, it may be possible to eliminate the quantization bias error through dither. Again, the $j^{th}$ simulation produces a set of data points $(p_i^j, N_i^j)$. Without dither, $p_i^j=p_i=\bar{p_0}+i\nabla{p}$ for all $j$. With dither, each $p_i^j$ is different. For each simulation, the expected value $m_j$ may be computed through Equation 4. The variation in $m$ is now due to both the dither and the stochastic nature of the underlying hull walk process.

Distribution Fitting

With or without dither, determining the central tendency $\mu$ of the sample comes down to a question of distribution fitting. We will examine 7 options which we will denote $\mu_{MoM}$, $\mu_{MLE}$, $\mu_{MD1}$, $\mu_{MD2}$, $\mu_{MD\infty}$, $\mu_{Med}$, and $\mu_{Mode}$. The first, $\mu_{MoM}$ uses the approach which has been used in all prior studies.

Method of Moments

In all previous investigations of percolation thresholds using gradient percolation, quantization bias errors have been ignored. The samples generated from the discrete distribution have been assumed to have been drawn from the true distribution. The expected value of the discrete probability distribution is used as the estimate of the percolation threshold, which is equivalent to fitting a normal distribution to the data through the method of moments (MoM) [36].

\begin{equation} \mu_{MoM,j}=\frac{1}{N_{total}^j}\sum_i N_i^j p_i^j \approx \frac{N_{occ}^j}{N_{total}^j} \end{equation}

\[\sigma^2_{MoM,j}=\frac{1}{N_{total}^j} \sum_i (N_i^j p_i^j - \mu_{MoM,j})^2.\]

Each simulation yields an expected value $\mu_{MoM,j}$ and many simulations can be combined by simple averaging, since variations in $N_{total}$ are negligible and not correlated with the variations in $\mu_{MoM,j}$, e.g. \begin{equation} \mu_{MoM}=\frac{1}{N_{sims}} \sum_{j=1}^{N_{sims}} \mu_{MoM,j}, \label{eq:averagesims_expval} \end{equation} where the width of the walk taken is proportional to the standard deviation $\sigma$ and may similarly be found by averaging the computed $\sigma_j$s.

Maximum Likelihood

An alternative method for distribution fitting is maximum likelihood estimation (MLE) [37]. In this approach, the values for the parameters $\mu_j$ and $\sigma_j$ are sought that maximize the likelihood (or equivalently the log-likelihood) of the observed data for each simulation $j$. For the continuous normal distribution, the log-likelihood is given by \begin{equation} \begin{aligned} &\log\Big(\mathcal{L}(\mu_j,\sigma_j)\Big) =
&-\frac{N_{total}^j}{2}\log(2\pi\sigma_j^2) - \frac{1}{N_{total}^j 2\sigma_j^2}\sum_{i} (N_i^j p_i^j-\mu_j)^2. \end{aligned} \end{equation} Differentiating and equating to zero yields \begin{equation} \begin{aligned} \mu_{MLE,j}&=\sum_{i} \frac{N_i^j p_i^j}{N_{total}^j}
\sigma_{MLE,j}^2&=\frac{1}{N_{total}^j} \sum_{i} (N_i^j p_i^j - \mu_{MLE,j})^2, \end{aligned} \end{equation} which is equivalent to the estimates obtained by MoM.

For the discrete distribution described in section IIB, the log-likelihood is given by \begin{equation} \log\Big(\mathcal{L}(\mu,\sigma)\Big) = \sum_{i} N_i\log( w(\mu,\sigma^2,p_i)), \end{equation} where $w$ is given by Equation 3. Similar to MoM, variations in $N_{total}$ are negligible and not correlated with the variations in $\mu$, such that \begin{equation} \label{eq:MLE} \mu_{MLE}=\frac{1}{N_{sims}} \sum_{j=1}^{N_{sims}} \mu_{MLE,j}. \end{equation} For the results and discussion to follow, we will refer to the minimum likelihood estimate for the discrete distribution given by Equation 10 as $\mu_{MLE}$.

Minimum Discrepancy

Another alternate estimate for $\mu$ may be recovered by minimizing the discrepancy $\lVert{\frac{N_i}{N_{total}\nabla p}-w(\mu,\sigma^2,p_i)}\rVert$, where $\lVert{\cdot}\rVert$ is an arbitrary vector norm. From Equation 3, the probability mass function for the $j^{th}$ simulation is given by $P(p=p_i) = w(\mu,\sigma^2,p_i^j)$. The empirical probability mass function is given by $\hat{P}(p=p_i) = \frac{N_i^j}{N_{total}^j\nabla{p}}$.

The discrepancy between the probability mass function and empirical probability mass function is given by \begin{equation} \mathcal{D}(\mu,\sigma) = \lVert{P-\hat{P}}\rVert \end{equation} where $\lVert{\cdot}\rVert$ is an arbitrary vector norm. The parameters $\mu$ and $\sigma^2$ may be determined as the values which minimize the discrepancy. For example, using the Euclidean norm, the objective function is defined as \(\mathcal{D}(\mu_j, \sigma^2_j) = \lVert{P-\hat{P}}\rVert_2 = \sqrt{\sum_i \left(\frac{N_i^j}{N_{total}^j\nabla p} - w(\mu,\sigma^2,p_i^j) \right)^2}.\)

Minimizing the objective function gives an estimate for $\mu_j$ and \(\sigma^2_j\), and this process is repeated for each simulation. The final estimate for $\mu$ is found by averaging \begin{equation} \mu_{MDn}=\frac{1}{N_{sims}} \sum_{j=1}^{N_{sims}} \mu_{MDn,j}, \label{eq:averagesims} \end{equation} where $n$ is either $1$, $2$, or $\infty$ depending on the norm used.

Median and Mode

Due to the presence of skewness in the data, the classical statistical methods described in sections IID1 and IID2 consistently underestimate the location of the distribution peak, resulting in some distinctly undesirable consequences, namely competing convergence rates that combine to yield the nonmonotonic convergence behavior seen in Figure 1. Nonmonotonic convergence behavior may be avoided by using the median or mode as the measure of central tendency (rather than the mean) and the root mean square deviation about the median or mode as a corresponding approximation of the hull width. In this way, two alternative defintions of $\mu$ and $\sigma$ are defined, with

\[\mu_{Med,j}=\langle M \rangle\] \[M = \left\lbrace p: \sum_{i\in\lbrace\xi: p_\xi^j<p\rbrace} N_i^j = \sum_{i\in\lbrace\xi: p_\xi^j>p\rbrace} N_i^j \right\rbrace \\ \sigma_{Med,j}^2=\frac{1}{N_{total}^j} \sum_{i=1}^{N_{total}^j} (p_i^j - \mu_{Med,j})^2\]

and

\[\mu_{Mode,j}=p_\xi^j, \quad\xi=argmax_i N_i^j\] \[\sigma_{Mode,j}^2=\frac{1}{N_{total}^j} \sum_{i=1}^{N_{total}^j} (p_i^j - \mu_{Mode,j})^2.\]

In Equation 14, we use $\langle\cdot\rangle$ to denote the average over the set.

Error Estimation

The uncertainty in $\mu$ is determined by examining the variance of $\mu_j$. The samples $\mu_j$ are observed to be normally distributed about $\mu$ regardless of the method used to compute $\mu$. The uncertainty in $\mu$ is given by its standard error \begin{equation} SE_\mu = \sqrt{\frac{1}{N_{sims}(N_{sims}-1)}\sum_{j=1}^{N_{sims}} \lvert \mu_j - \mu \rvert^2}. \end{equation} When using the method described in Section IID3, the uncertainty in $\mu$ is observed to scale approximately with $\nabla p^{0.136(2)}$, as shown in Figure 2. This exponent is not known to be related to any of the known critical exponents for this lattice geometry.

Results

As a case study, consider the canonical example of site percolation on a two-dimensional square lattice with nearest neighbor (von Neumann) connectivity. This system has been studied extensively in the literature [12,25,29,38–40] and the percolation threshold is known with high precision to be $p_c\approx0.5927460507921$ [29]. Each of our simulations consists of a frontier walk traversing the domain 501 times. We perform a large number of simulations with occupation probability gradients ranging from $8^{-1}$ to $655360^{-1}$. For each simulation, the number of sites visited per pass through the computational domain is proportional to the width of the domain $W$ and related to the probability gradient through the fractal dimension of the hull walk [10] \begin{equation} \frac{N_{total}}{W} \propto \nabla{p}^{-\alpha_N} = \nabla{p}^{-\frac{\nu}{1+\nu}D}. \label{eq:hullwidth} \end{equation} In the simulations, we use the pcg64\textunderscore k1024 random number generator [41].

Equation 17 relates the total number of sites generated per lattice width to the probability gradient $\nabla p$, where $\nu=4/3$ describes the divergence of the correlation length and $D=1+1/\nu$ is the fractal dimension of the frontier [42]. In the present numerical investigation, the exponent $\alpha_N$ is experimentally determined to be approximately $0.429$, which is consistent with the theoretical value of $3/7$.

The width of the percolation hull has been shown to scale as $\nabla{p}^{-4/7}$ based on correlation length arguments [10]. Note that this width is in lattice units. In concentration units, the width ($p_{max}-p_{min}$) scales as $\nabla{p}^{3/7}$, as shown in Figure 3. The scaling of the hull width with $\nabla {p}$ remains consistent regardless of which distribution fitting procedure is utilized.

Figure 4 shows the effect that quantization bias errors can have for large probability gradients and the effectiveness of using dither to randomize the effect and produce a reasonably normal output. Since the lattice offset is an arbitary simulation parameter, it is desirable for the resulting distribution to be independent of that choice (no horizontal variation). Unfortunately, for large probability gradients this is not the case. For smaller probability gradients (see Figure 5), the effect is significantly reduced. The dependence of $\mu_{MoM,j}$ on $p_0$ decays very rapidly as $\nabla{p} \to 0$. While dithering is an effective technique for eliminating the bias associated with quantization, a prudent alternative approach is likely to restrict simulations to such sufficiently small probability gradients that the effect is negligible. If large probability gradients are unavoidable, the minimum discrepancy distribution fitting procedure successfully eliminates the quantization bias effect even for large $\nabla p$ as shown in Figure 6.

While not pictured here, the results using the MLE distribution fitting procedure are indistinguishable from those using MoM. Both approaches define the central tendency using the mean occupation probability of the sites visited. This is in contrast to the median and mode approaches as well as the minimum discrepancy approach which appear to be more robust to to skewness in the distribution than the mean-based MoM and MLE approaches.

Despite the fact that the results of each simulation consistently pass the Anderson-Darling normality test [43,44], there is notable negative skewness \begin{equation} \gamma_j = \frac{\sum_i \left( N_i^j p_i^j -\mu_{MoM,j} \right)^3}{\left(N_{total}^j -1\right) \sigma_{MoM,j}^3} \end{equation} in the data. For large probability gradients, this skewness is due to the walk encountering sites with an occupation probability of 1 more often than sites with an occupation probability of 0 (since $p_c > 1/2$) as shown in Figure 7. For smaller probability gradients, the walk does not impact these hard boundaries, yet the skewness remains as seen in Figure 8. Note that the skewness decays proportionally to the hull width and that the shift $\mu_{Med} - \mu_{MoM}$ converges consistently across the entire range of $\nabla p$ (Figure 9), suggesting that both are related to the asymmetric boundary effect.

It is well-known that of the proposed measures of central tendency, the mean is most impacted by skewness. In this case, the negative skewness of the data causes the mean to be smaller than either the median or the mode. Crucially, for this lattice geometry, the mean is sufficiently shifted by the skewness to be \textit{below} $p_c$ for even fairly small values of $\nabla p$ while the other measures are \textit{above} $p_c$.

As $\nabla p \to 0$ the skewness decreases proportional to $\nabla p^{3/7}$ and the resulting leftward shift of the mean is reduced as shown in Figure 9. The mean is observed to approach the median approximately linearly as $\nabla p \to 0$.

Because $p_c$ is known so precisely for this system, the convergence of $\lvert p_c - \mu(\nabla p) \rvert$ may be easily studied without the added complication of determining $p_c$. Figure 10 shows the convergence of $\mu$ to the true percolation threshold $p_c$ as $\nabla p \to 0$ for the various methods for computing the central tendency. Note the steep dip in the error in the neighborhood of $\nabla{p}=1280^{-1}$ for the mean-based MoM and MLE estimates, indicating where the error transitions from negative to positive. To the right of this point, error introduced in the distribution fitting procedure cancels the error from the finite gradient producing unreliably accurate solutions. To the left of this point, a slight divergence and then continued convergence is observed, as the MoM and MLE values become sandwiched between the minimum discrepancy values and $p_c$. Although the mean-based measures produce the most accurate estimates of the percolation threshold for a given $\nabla p$, the non-monotonic convergence behavior makes extrapolating to $\nabla p=0$ using these values, as has been done previously, invalid. Additionally, the mean-based MoM and MLE values agree very well with each other indicating little benefit of treating the distribution as discrete.

The median, mode, and minimum discrepancy values converge monotonically to the true percolation threshold $p_c$ from above according to a power law \begin{equation} \label{eq:power_law} \mu(\nabla p) = p_c + c \nabla p^r, \end{equation} as shown in Figure 10. The mode and minimum discrepancy values all display similar convergence rates $r$ (see Table 1), although the differences between them are statistically significant. The convergence rates are known within the 95\% confidence intervals given in Table 1. The median has a noticeably faster convergence rate, but a larger error for a given probability gradient over the wide range of probability gradients simulated.

While the minimum discrepancy metric using the 1-norm generally provides the most accurate estimate of the percolation threshold, it results in a significantly less well-behaved optimization problem relative to the 2-norm or $\infty$-norm, which manifests itself in the jittery behavior observed at large $\nabla{p}$ in Figure 10.

The mean-based MoM and MLE metrics are not power law convergent in the range of $\nabla p$ frequently encountered in gradient percolation studies. It is possible that they might be power law convergent for significantly smaller $\nabla p$, but it is not possible to determine this from our data. Taken together, the results of Figures 8, 9, and 10 suggest that a fit of the form \begin{equation} \label{eq:mod_power_law} \mu_{MoM} (\nabla p) \approx p_c + c_1 \nabla p^{r_1} + c_2 \nabla p^{r_2} \end{equation} might be appropriate. Unfortunately, the two power law terms are of opposite sign and similar magnitude resulting in loss of significance. As a result, our data is unable to conclusively determine if Equation 20 is accurate.

The 3 constants in Equation 19 may be simultaneously fit to the experimental data. Doing so produces the convergence rates in Table 1 as well as the estimates of $p_c$ in Table 2.

Conclusions

We show that the average value of the occupation probability visited during a gradient percolation simulation is an unreliable surrogate for the percolation threshold on the regular square lattice. While initially suspected to be related to quantization effects, we demonstrate these effects to be small for reasonably small occupation probability gradients. For large probability gradients, dither may be an effective tool for randomizing the quantization effects; however, when metrics beyond the average occupation probability are used, the dependence on lattice offset is essentially eliminated and dither is likely unnecessary.

Unlike the average, the median, mode, and minimum discrepancy metrics are all observed to converge monotonically to the known value of the percolation threshold according to a power law. The minimum discrepancy metric behaves very similarly to the mode, both in terms of convergence rate and absolute error as a function of probability gradient, with the minimum discrepancy metric producing slightly more accurate results but with a slightly slower convergence rate. Due to the higher computational cost of computing the minimum discrepancy metric, the mode is likely preferable for most applications. The median has the fastest observed convergence rate (approximately linear) and the largest observed absolute error for a given probability gradient within the range considered. The net result of this is that all five approaches produce similarly accurate and similarly uncertain predictions for the extrapolated percolation threshold $p_c$ at $\nabla{p}=0$ when fit with a curve of the form $\mu=p_c + c\nabla{p}^r$.

Future work will seek to analytically derive relationships for the observed convergence rates. It is suspected that these rates should be related to the properties of the underlying probability distribution and the known critical exponents. A clue to these relationships may be found be examining how these rates change for different lattice geometries.