Pair interactions on finite square lattice

Question

I am looking for an exact or approximate solution to a statistical lattice-particle problem:

Given a lattice of size $L\times L$ where $\rho\cdot L^2$ particles are randomly distributed, calculate the probability distribution function for $M$ defined as number of pair interactions, i.e. number of edges between adjacent occupied cells. Interactions can be horizontal or vertical, and system has open boundaries (i.e. particle on the edge only has 3 or 2 neighbors to interact with).

My final aim is to design an $H_0$ statistics, given a single realization of the $L\times L$ lattice, and estimate the probability this is a non-random distributed set of particles.

Answer 1

Ok, some disperse thoughts, which may get you a little closer to an answer...

I started writing a long answer, on how your problem as stated is closely related to the hypergeometric distribution, but if your lattice is large enough, and the density of points small enough, you can approximate it to a geometric distribution which makes life much easier. It is a similar difference to one of those probability problems of picking balls from a bag, without replacement, or with replacement. In the first case, you will probably be stuck with a combinatorial approach to the problem, where you need to consider every possible arrangement of the particles. In the latter, you can do the approximation that any lattice point has probability $\rho$ of holding a particle, makinf things easier.

I'll be more than happy to elaborate on those, but if all you are after is figuring out whether your distribution has been generated by a non-random process, I would go a different way. If you want to give your statistic some formal support, you will want to read about Poisson Spatial Point Processes, here is a good intro.

But from a practical point of view, I would suggest using the autocorrelation as the base for your analysis. You could also go with the power spectrum, since it is closely related to the autocorrelation through the Wiener-Khinchin theorem, but it is easier to make physical sense of the autocorrelation.

So once you have your autocorrelation $C(i,j)$ computed, the value at $(i,j)$ is the probability that there is a particle at a distance $(i,j)$ from another particle. If your distribution is fully random, this probability should be constant an equal to $\rho$ (actually equal to $\frac{\rho L^2-1}{L^2-1}$, the geometric approximation comes in handy here as well), except at $(0,0)$ where it is $1$.

To come up with a single statistic for the whole distribution, you can consider that the value at $(i,j)$ has been generated by looking at each particle, and seeing if there was another at that distance. This is equivalent to flipping a loaded coin with a heads probability of $\rho$, $\rho L^2$ times, counting the number of heads, then dividing by $\rho L^2$. What you get before the last division follows a binomial distribution with parameters $n=\rho L^2$ and $p=\rho$. You can approximate the binomial distribution with a normal distribution with parameters $\mu=np=L^2 \rho^2$ and $\sigma^2 = np(1-p) = L^2 \rho^2(1-\rho))$. After dividing by $\rho L^2$, the distribution will still be normal, but with parameters $\mu= \rho$ and $\sigma^2 = \frac{1-\rho}{L^2}$. You can normalize this normal distribution by substracting $\mu$ and dividing by $\sigma$.

So when you substract $\rho$ from your autocorrelation everywhere, and then divide by $\sqrt{1-\rho}/L$, you than have $\rho L^2 - 1$ values drawn from $N(0,1)$ distributions. Actually only half of those are really independent values, because $C(i,j) = C(-i,-j)$. But if you square each of the $(\rho L^2-1)/2$ independent values and then add them all together, the result should be distributed as a $\chi^2$ with $(\rho L^2-1)/2$ degrees of freedom, which should suit your hypothesis testing needs.

Of course, you want to make sure that the actual values of $\rho$ an $L$ allow to do the "binomial-as-normal" approximation.

Answer 2

The solution, which is essentially exact for large lattices, is just that the distribution function is a Gaussian. This is a universal result for the sum of any locally fluctuating quantity. It is therefore not very useful as a criterion for determining whether a given instantiation of on/off sites is random.

I take a torus, as Jaime suggested, and further, I replace your constraint of a fixed density with a grand-canonical system, where there is a probability p for having a particle at any given site. This is equivalent to what you are doing for large systems, and this is already essentially exact for a 10 by 10 lattice.

To compute the mean and variance of your quantity, you can use partition function methods. The probability of a site being occupied is p, and unoccupied is 1-p, so define a spin variable S on each site which takes the value 0 or 1, according to whether the site is occupied or unoccupied. If you make the weight of being unoccupied equal to 1, the weight of being occupied is $p\over 1-p$, so you can define an energy for being occupied using:

$$ e^J = {p\over 1-p} $$

Then write the partition function

$$ Z = \sum_{S} e^{J_x S_x} $$

Where $J_x$ is a local parameter at each site which you set to J at the end of the day, to reproduce your system, and the sum is over all configurations of spins chosen to be either 0 or 1 at each site

This partition function is independent at every x (this is the advantage of allowing each site an independent probability of occupation, instead of constraining the total number--- it is the same reason people use grand canonical methods in statistical mechanics), so it factorizes:

$$ Z = \prod_x (1 + e^J) = \prod_x {1\over 1-p} $$

Your quantity is

$$ M = \sum_{\langle x,y\rangle} S_x S_y $$

Where the brackets means sum over each nearest neighbor pair once and only once, and you evaluate this by taking the expectation value, which is very simple because it is an independent system:

$$ \langle M \rangle = \sum \langle S_x S_y \rangle = \sum \langle S_x\rangle \langle S_y\rangle = 2L^2 p^2 $$

That's the mean of M. You can also evaluate it formally using derivatives of the partition function with respect to $J_x$, but the final answer is very simple--- it's just an independent system, so the number of pairs is independent.

To find the variance of M, you subtract the mean of M and square, and the result is

$$ \langle M^2 \rangle - \langle M\rangle^2 = ( \sum (S_x - p)(S_y - p) )^2 = \sum_{\langle x,y\rangle} \sum_{\langle x',y'\rangle} (S_x - p)(S_y -p) (S_{x'} - p ) (S_{y'} - p) $$

And the expected value of the thing on the far right is zero unless the pair x,y is the same as the pair x',y', otherwise there is a zero expected value.

The resulting variance is

$$ 2L^2 (\langle S_x^2\rangle - \langle S_x\rangle^2)^2 = (p(1-p))^2$$

Where I have used the fact that the variance of the spin at one site is $p(1-p)$, which is a straightforward computation.

The mean and variance specify the Gaussian distribution uniquely. This answer is slightly wrong, inasmuch as you have fixed boundaries and a global number constraint, but it is exact in the large system limit.

This doesn't work for what you want, namely to determine if a bunch of points is random. Nothing works for this, really, because the notion of "random" is too vague--- you need to specify more about the alternative probability distribution it could be, some prior knowledge. If you are given a single configuration, for all you know, it could have been selected from the probability distribution which is just a delta-function on this one configuration.

But in practice, you can use the correlation function relations

$$ \langle S_{x_1} S_{x_2} S_{x_3} ... S_{x_n} \rangle = p^n $$

which should hold up to errors which are as the Gaussian width of the object on the left hand side, which gives you a good criterion for independent randomness on each site.