Solving "all-marginals" problem for independent sets on grid Suppose I have a distribution over independent sets on an $n\times n$ grid where the probability of independent set occupying nodes $(i_1,j_1),\ldots,(i_k,j_k)$ is proportional to $\lambda_{i_1,j_1}\cdot \lambda_{i_2,j_2} \cdots\lambda_{i_k,j_k}$. Let $p_v$ be the probability that vertex $v$ in this grid is occupied, how do I compute $p_v$ for every vertex $v$ in the grid?

For instance, when $n=2$, the four probabilities are the ratios below


Motivation: comparing accuracy of "correlation decay" based approximate counting with exact for some not-too-large graphs

Clarification 11/6: This is the same $p_v$ as in the first formula of page 3 of [Weitz, 2006]
Here's a dynamic-programming algorithm that is fast enough to solve a given 18 x 18 grid in about 20 minutes.  The algorithm is a modification of the algorithm from this paper for counting independent sets in a given $m\times n$ grid.

Fix the input: the grid dimension $n \times n$ and the values $\lambda_{ij}$ for $i,j\in\{1,\ldots,n\}$.  Represent each independent set in the grid by its 0/1 incidence matrix $M$ where $M_{ij}$ is 1 if the grid vertex $(i,j)$ is in the independent set, and 0 otherwise.
Let $S_n$ denote the set of 0/1 vectors of dimension $n$ with no two consecutive 1's.  (I think the size of $S_n$ is the $n+2$nd Fibonacci number, so $|S_{18}|=6765$.)

An $n\times n$ 0/1 matrix $M$ represents an independent set if and only if the following two conditions hold:


Each row $M_i$ of $M$ is in $S_n$.
Every two consecutive rows $M_i$ and $M_{i+1}$ are compatible, meaning that there is no $j$ such that $M_{ij}$ and $M_{i+1,j}$ are both 1.


Construct directed acyclic graph $G=(V,E)$ where $V$ and $E$ are defined as follows.  The vertex set $V$ contains source and sink vertices $s$ and $t$, and, for each row $i\in\{1,\ldots,n\}$ and $x\in S_n$, a vertex $v(i, x)$.  The edge set $E$ contains the following edges: for each $x\in S_n$, an edge from $s$ to $v(1, x)$ and an edge from $v(n, x)$ to $t$; for each $i\in\{1,\ldots,n-1\}$, and each $x\in S_n$, edges from $v(i, x)$ to $v(i+1, y)$ for every $y\in S_n$ that is compatible with $x$.

The $s$-$t$ paths in $G$ correspond bijectively to the independent sets in the grid.  Each $s$-$t$ path is of the form $(s, v(1, x_1), v(2, x_2), \ldots, v(n, x_n), t)$, corresponding to the independent set whose incidence matrix $M$ has $M_i = x_i$ (that is, the $i$th row of $M$ is $x_i$).

Assign weights to the vertices of $G$ as follows.  Give $s$ and $t$ weight $w(s)=w(t)=1$.  Give each vertex $v(i, x)$ weight $w(v(i,x)) = \prod_{j:x_j=1} \lambda_j$.  Then, for any $s$-$t$ path in $G$, the product of the weights of the vertices on that path is the product of the $\lambda_{ij}$'s of the vertices in the corresponding independent set in the grid.

Now, for every vertex $v$ in $G$, define $P(v)$ to be the sum, over all paths from $s$ to $v$, of the product of the weights of the vertices on the path.  Define $P'(v)$ to be the sum, over all paths from $v$ to $t$, of the product of the weights of the vertices on the path except for $v$.  Compute these two quantities for every vertex in $G$.  This can be done in time linear in the size of $G$ (assuming constant-time exact arithmetic) using a standard dynamic-programming approach.  E.g. use $P(s) = 1$ and, for $v\ne s$, 
$$P(v) = w(v)\sum_{u: (u, v)\in E} P(u).$$

For each vertex $v$ in $G$, define $Q(v) = P(v)P'(v)$.  This is the sum, over all $s$-$t$ paths in $G$ that go through $v$, of the product of the weights of the edges on the path.  Then $Q(s)$ is the sum, over all independent sets in the grid, of the product of the $\lambda$'s for the vertices in the independent set. For each vertex $v(i, x)$, then, $Q(v(i,x))$ is the sum, over all independent sets $M$ in the grid whose $i$th row $M_i$ is $x$, of the product of the $\lambda$'s for the vertices in the independent set.

Finally, for each vertex $g=(i,j)$ in the grid graph, compute $p_g$ as the sum, over all $x\in S_n$ such that $x_j = 1$, of $Q(v(i, x))$, divided by $Q(s)-1$.  (Subtract 1 as the problem specification doesn't count the empty set in the denominator.)



For the instance with $n=18$ and each $\lambda_{ij}=1$, this Python script computed the answer in about 17 minutes.  (The graph $G$ had about 122K vertices and 159M edges.)  For non-uniform lambdas, in the probability calculations, you'll need to take some care regarding arithmetic precision.  (If you want exact answers, you'll need to use rational arithmetic.  You could compute close approximations using floating-point arithmetic, but then you'll probably have to be careful about the order in which the sums are performed.)

Here's the output of the script for the instance with $n=7$ and each $\lambda_{ij}=1$. Each entry of the 7 x 7 matrix below is the probability for the corresponding vertex of the grid, times the denominator (which for $\lambda_{ij}=1$ is the number of non-empty independent sets in the grid).

N = 7

denominator = 1280128949

402968942 298303063 333041917 316303258 333041917 298303063 402968942
298303063 296843520 279776968 292351405 279776968 296843520 298303063
333041917 279776968 296033698 285917681 296033698 279776968 333041917
316303258 292351405 285917681 293815906 285917681 292351405 316303258
333041917 279776968 296033698 285917681 296033698 279776968 333041917
298303063 296843520 279776968 292351405 279776968 296843520 298303063
402968942 298303063 333041917 316303258 333041917 298303063 402968942