Identically distributed and correlated systems of random variables I have a rectangular grid of $n \times m$ lightbulbs on a torus, each situated at a grid point. That is, the lightbulbs' coordinates range from $(0,0)$ to $(n-1,m-1)$ and the lightbulb $(0,0)$ has neighbors $(0,1)$, $(1,0)$, $(0,m-1)$, and $(n-1,0)$.
For each edge, I flip a coin that gives heads with probability $p$. Any lightbulb that has 1 or 3 of its edges heads is on; otherwise it is off.
Let $X_{i,j}$ be a random variable that is 1 if the lightbulb at $(i,j)$ is on, and 0 otherwise. Clearly the $X_{i,j}$ are identically distributed. However, they are not independent. I would like to describe the distribution of the system of $X_{i,j}$ (whatever that means; the individual $X_{i,j}$ are tame enough). Assume that $n=m$ if that helps. Formally, what I am asking is: if $X$ is the matrix whose $(i,j)$th entry is $X_{i,j}$, how is $X$ distributed? It is not the same as flipping a coin for each entry.
I am rusty enough in my probability theory that I don't really know where to begin. I notice that this has a similarity to the square lattice Ising model, but I believe it is only a resemblance -- there is no value of $p$, for instance, that leads to nearly all of the lightbulbs being on. If this is too difficult to answer, I would appreciate references to toy systems of this nature so that I could analyze mine by analogy.
 A: I'd suggest the following, construct two edge matrices, one for the horizontal edges which we'll call $H_{i,j}$ and one for the vertical edges which we'll call $V_{i,j}$. Then the state of bulb $X_{i,j}$ can be determined as
$$X_{i,j}=H_{i,j}+H_{i-1,j}+V_{i,j}+V_{i,j-1}$$
in which the states of the edges have to be added modulo 2, whereas the indices operate modulo whatever the size of your torus is (modulo n horizontally, modulo m vertically). The edge variables are behaving as independent Bernoulli variables with probability $p$. 
So, with this you can at least generate such random matrices easily with the computer. I see you are reading about random matrix ensembles. I don't know about such ensembles for the particular situation you are describing, but it is an interesting next step in the analysis. However, the point of such ensembles is most of the time to find a measure on the space of matrices. This measure will depend only on the eigenvalues of the matrices because of symmetry properties. 
First thing to note however is that your matrices are not square, so eigenvalues are not necessarily defined. We could however possibly work with singular values and this would reduce to eigenvalues in the square matrix case.
Here's some more thoughts I've put into the question:
First, since we are working on a torus, we can introduce matrices to represent shifts in horizontal or vertical direction. Our measure on the space of matrices should in the end be invariant under such shifts because it is immaterial where we choose our origin. Let's call the horizontal shift operator $S_H$ and the vertical one $S_V$. Both are square matrices of respectively $n \times n$ and $m \times m$ dimensions. It is also clear that horizontal shifts happen through left multiplication whereas vertical ones through right multiplication, so that
$$X=H+S_H H+V+V S_V =  (\mathbb{I}_H+S_H) H+V(\mathbb{I}_V + S_V)$$
where I have introduced the $n \times n$ and $m \times m$ identity matrices.
Then, it is also easy to construct the probability measures for the edge matrices which are
$$\mathbb{P}(H)=\prod_{i,j}p^{H_{ij}}(1-p)^{1-H_{ij}}=p^{\sum_{i,j}H_{ij}}(1-p)^{nm-\sum_{i,j}H_{ij}}$$
and 
$$\mathbb{P}(V)=\prod_{i,j}p^{V_{ij}}(1-p)^{1-V_{ij}}=p^{\sum_{i,j}V_{ij}}(1-p)^{nm-\sum_{i,j}V_{ij}} \; .$$
Note that indeed $\mathbb{P}(S_H H)=\mathbb{P}(H S_V)=\mathbb{P}(H)$ and $\mathbb{P}(S_H V)=\mathbb{P}(V S_V)=\mathbb{P}(V)$.
Formally, we then have that
$$\mathbb{P}(X)=\sum_{X}\mathbb{P}(H,V)=\sum_{X=H+S_H H+V+V S_V}\mathbb{P}(H)\mathbb{P}(V) \; .$$
Now comes the hardest part. We have to count how many configurations of the edge states lead to the same $X$ and for each of these configurations we have to find out the corresponding probability $\mathbb{P}(H,V)$. So, let's play a game I'll call the flipping game.
Say we flip a horizontal edge, this alters the state of the two adjacent bulbs so that we have to flip more edges to compensate. If we choose to flip their other adjacent horizontal edges, we then need to correct for the bulbs further along, etc... In the end, if we flip all the edges along a horizontal line, we have done nothing to the bulbs. Same for flipping edges along a vertical line.
$$\cdots \circ - \circ - \circ - \circ - \circ - \circ \cdots$$
If instead of correcting by flipping more horizontal edges, we had flipped the adjacent upper vertical edges, we could have closed the flipping by doing one more flip of the upper horizontal edge, which gives a little square pattern. So flipping squares doesn't change the overal state of the bulb matrix. 
$$\begin{matrix} & \vdots &  & \vdots &\\
\cdots & \circ & - & \circ & \cdots \\
 & | &  & | &  \\
\cdots & \circ & - & \circ & \cdots \\ 
 & \vdots &  & \vdots & \end{matrix}$$
The general rule is that flipping along closed loops keeps bulb states invariant. This still makes the counting problem quite hard. Maybe this could be best solved numerically.
