把格子纸里的格子随机染黑白两色，平均每片色斑有多少格子？

01-16

无限大的格子纸(4连通),每一块等概率染黑或者白
平均连通的每一块有多大?

有很多人说看不懂，还有人说题目也看不懂，这让我很受挫啊。。。。

我再来把题目说一下。

在一个无限大的格子纸上对格子染色。每个格子有 $p$ 的概率被染成黑色。那么问当 $p=0.5$ 时，被染色块（cluster）的平均大小是多少。就上图来说，取这个无限大格子纸的一部分，一共有5个cluster，分别是

大小为8的1块
大小为6的1块
大小为3的1块
大小为1的2块

那么求平均的结果是 $(8 imes 1+6 imes 1+3 imes 1+1 imes 2)cdotfrac{1}{5}=3.8$ 。当然，因为只是取一小块，所以会有涨落（偶然成分）。这个问题很简单，在码农遍地的知乎，随便写个程序模拟一下这个过程，得到平均大小在是7.6左右是很容易的事情。

但问题是，为什么可以做模拟，为什么题目中说无穷大的格子纸，我们拿有限的大小的去模拟是可以得到结果的。这点在其他很多问题上是显然，比如用求和去算积分。但这个上，答案是否定的，题目是运气非常好的说黑白是等概率的。如果有2/3的概率染黑，问黑块的大小。再这样做模拟根本是不对的。因为这时候，最大的那个黑块已经可以和整个空间相比拟（如下图中深蓝色的基团）。

这是我们从无限大的格子纸中取的一块，这一块是这个样子了，那个这一块旁边的情况估计和这一块差不多，也会有一个超大的基团。而这两个超大的黑块显然是相连的，实际上是一个更大的黑块。所以这时候最大的那个黑块的大小趋向无穷，而黑块平均的大小也趋向无穷。这时你用一个有限大小的格子纸是模拟不出无穷大情况的。

==========原答案==========

这个问题是一个经典的统计模型，叫做逾渗（http://www.wikiwand.com/en/Percolation）

模型几乎就是题主描述的（题主好聪明啊）。区别是只讨论一种颜色，比如每个格点有 $p$ 的概率放黑色。当然，这是经典模型，所以现在肯定是被讨论的比较充分了，背景空间，也就是格子纸，不一定是格子纸，也可能格子“块”，三角格子纸等等。这个模型关注的是，随着 $p$ 的变化，什么时候出现一块“超大”的黑色块（更准确的说是出现超大黑色块的概率）。这块连通的黑色有多大呢，无穷大。这个 $p$ 值记为 $p_c$ ，以及 $p$ 在 $p_c$ 附近系统的行为。

出现超大黑色块的概率随着 $p$ 的变化而变化这个过程实际上是一种二级相变。逾渗模型是一种典型的几何相变（这里没有什么热力学量）。

那么针对题主的问题，这里的 $p_capprox 0.593$ （http://www.wikiwand.com/en/Percolation_threshold）当过了这个临界点，就出现连通全空间的超大黑块，这时候平均连通块的大小就会变成无穷（因为空间也是无穷大的）。但比较好，题主的例子里， $p=0.5$ 还没有到达这个临界值。这就说明每一个黑块的大小都是有限的，所以你取多大的格子纸来算这个mean cluster size问题都不大。结果大概是7.6左右。

再上一张图，让我们更好的理解这个过程。这里的 $n$ 表示的是一个 $n imes n$ 的格子纸，横轴是格点染黑的概率，纵轴是平均黑块大小（mean cluster size）。怎么看这张图呢，当 $p<P>时，mean cluster size与<img src=$ 是没有关系的。但当过了这个临界值，会出现覆盖全空间的黑块，这时候覆盖全空间黑块的尺寸显然和 $n$ 有关了，当 $n o infty$ 时，mean cluster size也趋于无穷。

这个问题一般人就不要试图去寻找解析解了，因为即使是用重整化群算 $p_c$ ，也是不准确的。要是做出了mean cluster size的解析表达式，基本也是一个大新闻了。

========7.30更新=========

一觉起来这么多赞（请原谅我没见过世面 2333），就再来说两句。

对于一个系统,最显然的理解就是一些对象之间存在某种联系,把这些对象放在一起作为一个更大的对象进行考虑。而逾渗模型就是对这个过程最直接的一种抽象。相互依赖系统的结构相变、鲁棒性、社团结构划分等等实际问题也与之紧密相连。

再说一下计算。程序中用到了一步取余。这是为了去掉边界效应对模拟的影响，简单的理解就是没有边界了，或者说有一个周期性边界。如下图所示（这里每次添的是边，而不是点，不过这两个模型基本是一样的。）注意这个方法并不是用有限来模拟无限，因为整个系统还是有限的。举个极端的情况，如果系统本来的大小是 $2 imes 2$ ，也可以用这样的方法卷起来，但显然是模拟不出7.6这个结果的。所以 @杜昊的回答中说可以通过周期性来模拟无限大系统是不太准确的。

另外对于题主的问题，mean cluster size的量级在10，所以一般随便给个尺度，比如 $n=100$ 就足够了。但当接近 $p_c$ ，上面说mean cluster size会趋于无穷，但图中的曲线显然是连续的。这个当然是好理解的，格子纸是有限的，当然不会出现无穷。这被称为有限尺寸效应。为什么说这个呢，因为很多 $p_c$ 附近的相变行为就可以研究不同尺寸下物理的不同表现来研究。很简单的理解就是这里的超大黑块的尺寸能和空间尺寸相比拟，就可以研究超大黑块关于 $n$ 的行为来知道系统的一些行为。这种方法叫做有限尺寸标度。

import itertools import random import numpy as np import matplotlib.pylab as plt


def simulate(l=100, p=0.5):

    nodes={}

    for i,j in itertools.product(range(l), range(l)):

        nodes[(i,j)]=[(i,j),0]
    def find_father(node):

        if not nodes[node][0]==node:

            nodes[node]=find_father(nodes[node][0])

        return nodes[node]
    def add_to(set_a, set_b):

        f_set_a=find_father(set_a)

        f_set_b=find_father(set_b)

        if not f_set_a==f_set_b:

            f_set_a[1]+=f_set_b[1]

            f_set_b[0]=f_set_a[0]
    for i,j in itertools.product(range(l), range(l)):

        if random.random()&<=p:
            nodes[(i,j)][1]=1
            for di, dj in [(0,1), (0,-1), (1,0), (-1,0)]:
                if nodes[((i+di) % l, (j+dj) % l)][1]&>0:

                    add_to(((i+di) % l, (j+dj) % l),(i,j))
    return np.mean([v[1] for n, v in nodes.items() if n==v[0] and v[1]&>0])

#print(simulate(1000,0.5)) for l, times in [(10, 1000), (20,400), (100, 30)]: print(l) p_all=np.arange(0.4,0.8,0.01) mean_size=[] for p in p_all: print(p) mean_size.append(np.mean([simulate(l,p) for _ in range(times)])) plt.plot(p_all, mean_size,label="n=%d" % l) plt.xlabel("p") plt.ylabel("mean cluster size") plt.legend() plt.show()

上面答案太高大上了，没有看懂，自己就随手写个小程序试一试。

上面答案是多少，真没看懂（囧）。

先报答案：约为7.6

假定有限空间，随机一个 N*N的数组，每个数值为0或1，算出有多少个封闭块，就可以得出平均联通区域大小了。

结果是这样的：

001101001110101100111010101100

111111110010010011000000111110

100100111010100110100010111100

001011001101111011010100100000

100000100110110000011110010101

000110000100111010000010011010

101001011000101110000110011101

111010010001110011010010011000

001110001110010011110100100110

000111111011101001001110100111

011001011011101100111001011100

000101001010000010010011011100

110101111101011010111100011000

111011100010001010111110111111

100010111010100101010101111000

000100011101111111000100001111

101011001111101110101001000111

001111100111110001100111100000

000010111001010011001010011010

101011000000001101110010110100

011110110100101001100101110100

011000001001011011101000111110

110110101001001010000111000101

001100100000010000011000111100

010001110011001000011001101101

111100100010000001111011011101

100111111100010111010111110011

011100111101011000011111000001

011101111110010101111101011011

100101111110100001000110100100

30*30区域 ;144个封闭块;平均每块大小6.25

下面就不打印数组了，有兴趣可以自己试试，我的代码贴在下面。

100*100区域 ;1466个封闭块;平均每块大小6.8212824010914055

500*500区域 ;33344个封闭块;平均每块大小7.497600767754319

1000*1000区域 ;132331个封闭块;平均每块大小7.55680830644369

5000*5000区域 ;3291136个封闭块;平均每块大小7.596161325451151

10000*10000区域 ;13161184个封闭块;平均每块大小7.598100596420505

15000*15000区域 ;29618321个封闭块;平均每块大小7.596649384683217

再大就成这样了,果然不能递归

Exception in thread "main" java.lang.OutOfMemoryError: Java heap space

at Test.&(Test.java:7)

at Test.main(Test.java:10)

非常简单的代码相信会点的人都能看懂，写得很水，欢迎指教

import java.util.Random;

public class Test {

public static final int NUM = 10000;

int[][] num = new int[NUM][NUM];

public static void main(String[] args) {

Test test = new Test();

test.initNum();

test.scanAll();

}

void initNum() {

Random random = new Random();

for (int i = 0; i &< NUM; i++) {

for (int j = 0; j &< NUM; j++) {

num[i][j] = random.nextInt(2);

System.out.print(num[i][j]);

}

System.out.println();

}

void scanAll() {

int areaCount = 0;

for (int i = 0; i &< NUM; i++) {

for (int j = 0; j &< NUM; j++) {

if (num[i][j] &< 2) {

int areaNum = scanOne(i, j, num[i][j]);

// System.out.println(i + ":" + j + ":" + areaNum);

areaCount++;

}

System.out.println(NUM+"*"+NUM+"区域 ;"+areaCount + "个封闭块;平均每块大小" + (double) NUM * NUM / areaCount);

}

private int scanOne(int i, int j, int value) {

if (i &< NUM j &< NUM i &>= 0 j &>= 0 num[i][j] == value) {

num[i][j] = 2;

return 1 + scanOne(i + 1, j, value) + scanOne(i, j + 1, value)

+ scanOne(i - 1, j, value) + scanOne(i, j - 1, value);

}

return 0;

}

这个问题实际上是Percolation network问题的一个具体情况。我这里有些资料，题主可以一看。

MCpercolation.pdf_免费高速下载

里面介绍了一种普适性的算法（Monte Carlo）来计算平均值。

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~分割一下~~~~~~~~~~~

没想到这个问题有人关注，为了我的知乎首赞我要补充一下。

高票答案中说到了用有限的格子无法模拟无限情况，然而实际上是可以的，这个模拟依托于这个模型中的格子具有周期性。这个周期性和热力学中的易辛模型有相同的性质。可能是我理解有误，不过高票答主在模拟中使用了周期性，却在更新中似乎没有意识到周期性的存在？

举例来说，格子alpha如果是一个有限网络A顶部的小格子，alpha的往上两格（实际上不存在于我们已知的有限网格，因为A已经是顶部了），等价于网络A的底部朝上两格（并且在同一列）。

楼主描述的这种情况叫做square site percolation（中文不明）。这里贴上论文中算法的PYTHON版本，以及MATHEMATICA得出的图。

import math

import random

L = 600 # linear dimension

N = L * L

EMPTY = -N-1

f=open("D:phscis hw\size600.txt","w")

ptr = [] # array of N pointers

nn = [] # 4 nearest neighbors of each N

order = [] # occupation order

def initialize():

# create arrays

global ptr, nn, order

ptr = [ 0 for i in range(N) ]

for i in range(N):

nn.append( [ i for i in range(4) ] )

order = [ 0 for i in range(N) ]

def boundaries():

for i in range(N):

nn[i][0] = (i+1)%N

nn[i][1] = (i+N-1)%N

nn[i][2] = (i+L)%N

nn[i][3] = (i+N-L)%N

if i%L == 0:

nn[i][1] = i+L-1

if (i+1)%L == 0:

nn[i][0] = i-L+1

def permutation():

for i in range(N):

order[i] = i

for i in range(N):

j = int(i + (N - i) * random.random())

temp = order[i]

order[i] = order[j]

order[j] = temp

def findroot(i):

if ptr[i] &< 0:

return i

ptr[i] = findroot(ptr[i])

return ptr[i]

def percolate():

big = 0

for i in range(N):

ptr[i] = EMPTY

for i in range(N):

r1 = order[i]

s1 = r1

ptr[s1] = -1

for j in range(4):

s2 = nn[s1][j]

if ptr[s2] != EMPTY:

r2 = findroot(s2)

if r2 != r1:

if ptr[r1] &> ptr[r2]:

ptr[r2] += ptr[r1]

ptr[r1] = r2

r1 = r2

else:

ptr[r1] += ptr[r2]

ptr[r2] = r1

if -ptr[r1] &> big:

big = -ptr[r1]

f.write(str((i+1.0)/(float(N))))

f.write(" ")

f.write(str(big/float(N)))

f.write("
")

f.close()

initialize()

boundaries()

permutation()

percolate()

print("finished")

这个程序用了一个600乘600的网格进行计算，并缓慢增加概率p的数值，来计算最大cluster(连续色块）/网格size的数值，输出如下。

2.77777777778e-06 0.0

5.55555555556e-06 0.0

8.33333333333e-06 0.0

1.11111111111e-05 0.0

1.38888888889e-05 0.0

1.66666666667e-05 0.0

1.94444444444e-05 0.0

2.22222222222e-05 0.0

…………

第一列代表概率，第二列代表比值。由于模拟本身存在随机性，我们可以多使用几个不同size的网格并分别运行这个模拟，最终可以得到一个比较精确的平均值。

这其中我们会发现一个特殊的概率值，p约等于0.59，在这个点，所有size的网络都出现了cluster_size的大幅度跳升。

这里不同颜色的曲线是不同尺寸网格的数值模拟结果，画到一起是为了让这个结果更加显眼(竖线）。

这个概率叫做 critical probability for site percolation, 也是图论领域的一个重要结果，定义就是，

，site percolation存在无穷大的cluster时的概率的最小值。不过这个数值一直都是一个近似值，没有人给出过解析解。

晚上贴一个我自己做的视频（虽然质量堪忧），是学期末的final project，有兴趣可以看看。。欢迎讨论。

FinalProject.mp4_免费高速下载

占坑，我也来试试！

只针对方格！

假设染黑一个方格的概率为p.

染黑的过程等于寻找一个黑色的联通域，从数据结构的角度讲，这个过程可以用一个队列来实现。

过程如下，

1. 找到最左边最上面的第一个黑格子；

2. 按照先左后右，先上后下的顺序，把周围8个格子中的黑格子加入一个队列；

3. 把队列第一个格子的周围未加入队列的格子加到队尾，队列头这个格子为联通联通域的一部分；

4. 重复以上过程，直到队列为空。

直观的说的这就是一个树，有下面的形状，比如

1. 每一个节点表示一个黑格，所有节点的数目为黑色联通域的面积；

2. 每一个节点的下层节点数确定了黑色块的形状；例如面积为5就有如下的形状,0表示没染色

0 1 0

2 3 4

0 5 0

队列顺序为 1 2 3 4 5

或者

1 0 2

3 4 0

0 0 5

队列的顺序为 1 3 4 2 5

假设黑色块面积为 n, 面积为n 的黑色块概率为f(n), 队列的最后几个元素一定是属于某个黑格的联通域，把这个黑格叫做父节点，这些最后元素的个数可能有1, 2,3,4,5,6,7,8个，分别的概率为p1, p2, p3,p4,p5,p6,p7,p8

例如：

有8个元素的情况最简单， p8 = 0, 按照从左到右，从上到下的顺序不可能出现8的情况。

按照全概率公式，

$f(n) = p1*f(n-1)+p2*f(n-2) + .... + p8*f(n-8)$

设 $G(z)$ 为 $f(n)$ 的z变换，有

$G(z) = sum_{1}^{8}{p_{i} *z^{i} } , f(1) = p*(1-p)^{8}$

面积为1 ，则只有一个格染色，周围全没染。

解此方程就得到面积n 的概率分布，有了概率分布就能求期望了。

关键的问题是求概率 p1 - p8,

已知 p8 = 0;

以最后元素为3 为例，按照我们的染色算法，最后元素的左边和上边不能出现除父节点以外的黑格，可能的情形有，1表示黑， 0 表示白

1 0 1

1 1 1

0 0 1 这种的概率为 $p^{6} *(1-p)^{3}$

或者

1 0 0

1 1 1

0 1 1 这个的概率同上，则 $p^{6} *(1-p)^{3}$

还有

1 0 0

0 1 1

0 1 1 $p^{5} *(1-p)^{4}$

所以3个元素的概率为 $p^{6} *(1-p)^{3} + p^{6} *(1-p)^{3} + p^{5} *(1-p)^{4}$

其他的还需要一个个算，等有时间再来搞，先搬砖去吧，O(∩_∩)O！

大家看看有没有啥问题，欢迎指导，总感觉有可能哪里有问题？

当然可以做计算机模拟，但是要做收敛测试。你先取1000x1000比如结果是7.500，再取2000x2000比如结果时7.501，再取4000x4000比如结果是7.501，你基本刻意认定最终结构是7.501。如果结果随着你的格子数变化而变化，那说明时无穷。

看了前面几个大神的答案，不是很懂……于是自己算了一下

设置的30 * 30无边界方格，模拟100000次求平均下来结果是7.8左右，跟大神们的答案差不多，但是也不排除我哪里算错了

生物本科，平时没什么机会写代码，但是觉得挺有意思，就写了一下，代码很烂，求大神拍

基本算法递归，已知一块黑色区域，每次向外扩展一圈，检查是否有黑色块，直到黑色区域外全是白色块为止，求得这一部分连通黑色区域的面积，代码如下（再次求大神拍）：

至于无限大的p值临界点的问题，还在揣摩大神的代码中

取消关注程序员了，还是可以看见程序……