谢林模型引言
在谢林模型中，有一片方形区域，区域又被均匀分成若干小方格；有一群被称为代理的个体，每个代理居住在一个小方格内。这些代理可以分成若干类，所有代理都希望周边8个格子尽可能住有较多的同类代理，若居住区域不满足代理的居住要求，则该代理会搬到另一区域居住。
如图1-1所示，在5×5的区域内有13名代理，不同的颜色表示不同的代理类型。整个区域还存在着一些无人居住的空白方格，代理们可以搬到这些空白方格去居住。
图1-1
假设代理5希望自己周围8格内同类代理的数量至少比异类代理的数量多2，则该代理会搬到如下图所示的位置。
图1-2
假如代理1也有同样的诉求，那么在图1-1中他的居住要求是能得到满足的。而在图1-2中，由于代理5搬走，代理1的居住要求不再被满足，此时他也会想搬走。当一个代理搬走时，一方面会使得该代理周边同类代理的居住满意度下降，另一方面也会留下空白方格，从而允许其它类型的代理搬入。由此不难推出，这是一种级联现象，微观上单个代理的行为有可能引发整个群体的连锁反应，从而导致宏观上各类型代理在居住分布上的同类大片连续、异类相互隔离的现象。
蒙特卡洛模拟法
基本思想
蒙特卡罗方法是由冯诺依曼和乌拉姆等人发明的，“蒙特卡罗”这个名字是出自摩纳哥的蒙特卡罗赌场，这个方法是一类基于概率的方法的统称，不是特指一种方法。
蒙特卡罗方法也成统计模拟方法，是指使用随机数（或者更常见的伪随机数）来解决很多计算问题的方法。他的工作原理就是两件事：不断抽样、逐渐逼近。
本题的蒙特卡洛思想
实际生活中，每一个代理（人家）搬家一定会提前选择好搬家的地址，也就是说搬家一次后基本不会发生变动（别人变动有概率让他变动），但是在程序中，我们的搬家并不能提前知道要搬到什么地方。于是，我们就让其随便搬家，如果不满足条件，则继续搬家，直到符合要求。
通过大量的模拟，代理一定会达到其应该在的位置并且不再搬家。
总结
通过大量模拟，让计算机模拟的模型（具体算法、图片都是模型）以概率（随着模拟量的增加，概率逐渐趋近于1，但始终不会是1）逼近真实的模型
代码验证
# 谢灵模型的属性：规格、阙值、矩阵、种群数量(最多为5种）、人口 import random import matplotlib.pyplot as plt import numpy as np class schellingModel: def __init__(self, size, threshold, kindNum): # 初始化谢林模型所需的参数 self.size = size self.threshold = threshold self.matrix = np.zeros((size, size)) self.kindNum = kindNum self.population = 0 if kindNum == 2: self.red = 1 self.white = 0 self.blue = 2 elif kindNum == 3: self.red = 1 self.white = 0 self.blue = 2 self.yellow = 3 elif kindNum == 4: self.red = 1 self.white = 0 self.blue = 2 self.yellow = 3 self.purple = 4 else: print("输入的种群数量要小于4,大于2") # 初始化谢林模型参数——种群矩阵、人口数 for i in range(size): for j in range(size): self.matrix[i][j] = random.randint(0, kindNum) if self.matrix[i][j] != 0: self.population += 1 # 判断（i，j）点是否满意————判断（i，j）点周围有多少个相似的邻居 def is_satisfied(self, i, j, kind): same_neighbour = -1 # 要把自己除去 neighbour_matrix = self.matrix[i - 1 if i - 1 > 0 else 0:i + 2 if i + 2 <= self.size else self.size, j - 1 if j - 1 > 0 else 0:j + 2 if j + 2 <= self.size else self.size] for x in range(len(neighbour_matrix)): for y in range(len(neighbour_matrix[x])): if neighbour_matrix[x][y] == kind: same_neighbour += 1 if same_neighbour >= self.threshold: return True else: return False # 随机获取一个位置作为搬家地点 def random_find_place(self): i = random.randint(0, self.size - 1) j = random.randint(0, self.size - 1) while self.matrix[i][j] != 0: i = random.randint(0, self.size - 1) j = random.randint(0, self.size - 1) return (i, j) # 返回i、j构成的元组 # 集体完成一次搬家 def move(self): for i in range(self.size): for j in range(self.size): kind = self.matrix[i][j] judge = self.is_satisfied(i, j, kind) if judge == 0: # 不满意 (x, y) = self.random_find_place() self.matrix[x][y] = kind self.matrix[i][j] = self.white # 把原本的家庭住址变为空 # 计算当前的满意家庭数 def satisfy_num(self): satisfy = 0 for i in range(self.size): for j in range(self.size): if self.matrix[i][j] != 0: kind = self.matrix[i][j] judge = self.is_satisfied(i, j, kind) if judge == 1: satisfy += 1 return satisfy # 画图谢林模型图 def draw(self, time, satisfy): if self.kindNum == 2: redx = [] bluex = [] redy = [] bluey = [] for i in range(self.size): for j in range(self.size): if self.matrix[i][j] == self.blue: bluex.append(i) bluey.append(j) elif self.matrix[i][j] == self.red: redx.append(i) redy.append(j) plt.scatter(redx, redy, c='r', marker='.', linewidths=0) plt.scatter(bluex, bluey, c='b', marker='.', linewidths=0) if satisfy == 0: plt.title('Initial') plt.show() return 0 else: title = str(time) + ' times satisfy:' + str(float(satisfy) / float(self.population)*100) + "%" satisfy_rate = float(satisfy) / float(self.population)*100 plt.title(title) plt.show() return satisfy_rate elif self.kindNum == 3: redx = [] bluex = [] redy = [] bluey = [] yellowx = [] yellowy = [] for i in range(self.size): for j in range(self.size): if self.matrix[i][j] == self.blue: bluex.append(i) bluey.append(j) elif self.matrix[i][j] == self.red: redx.append(i) redy.append(j) elif self.matrix[i][j] == self.yellow: yellowx.append(i) yellowy.append(j) plt.scatter(redx, redy, c='r', marker='.', linewidths=0) plt.scatter(bluex, bluey, c='b', marker='.', linewidths=0) plt.scatter(yellowx, yellowy, c='y', marker='.', linewidths=0) if satisfy == 0: plt.title('Initial') plt.show() return 0 else: title = str(time) + ' times satisfy:' + str(float(satisfy) / float(self.population)*100) + "%" satisfy_rate = float(satisfy) / float(self.population)*100 plt.title(title) plt.show() return satisfy_rate elif self.kindNum == 4: redx = [] bluex = [] redy = [] bluey = [] yellowx = [] yellowy = [] purplex = [] purpley = [] for i in range(self.size): for j in range(self.size): if self.matrix[i][j] == self.blue: bluex.append(i) bluey.append(j) elif self.matrix[i][j] == self.red: redx.append(i) redy.append(j) elif self.matrix[i][j] == self.yellow: yellowx.append(i) yellowy.append(j) elif self.matrix[i][j] == self.purple: purplex.append(i) purpley.append(j) plt.scatter(redx, redy, c='r', marker='.', linewidths=0) plt.scatter(bluex, bluey, c='b', marker='.', linewidths=0) plt.scatter(yellowx, yellowy, c='y', marker='.', linewidths=0) plt.scatter(purplex, purpley, c='p', marker='.', linewidths=0) if satisfy == 0: plt.title('Initial') plt.show() return 0 else: title = str(time) + ' times satisfy:' + str(float(satisfy) / float(self.population)*100) + "%" satisfy_rate = float(satisfy) / float(self.population)*100 plt.title(title) plt.show() return satisfy_rate if __name__ == '__main__': s = schellingModel(150, 4, 2) # 150户人家、阙值点为4，两类人群 s.draw(0, 0) # 画出初始图 i = 1 satisfy_num = 0 while(1): # 两次move计算聚合层度变化率 s.move() satisfy_num1 = s.satisfy_num() satisfy_rate1 = s.draw(i, satisfy_num1) i += 1 s.move() satisfy_num2 = s.satisfy_num() satisfy_rate2 = s.draw(i, satisfy_num2) cluster_rate = (satisfy_rate2-satisfy_rate1)/satisfy_rate1 if cluster_rate < 0.000001: # 什么时候停止算法迭代 break i += 1
模型结果分析
设置一片500×500的区域。设区域中居住有两类代理，每类代理各有100000个，初始时随机分布。设两类代理的居住要求都是周边8各中有至少4名同类代理。使用计算机程序进行模拟迭代，每次迭代时，按随机顺序判断每个代理的居住要求是否得到满足，不满足者搬往随机空白方格，程序模拟运行结果如下：
可以看到，在迭代20次后，整个区域就基本稳定下来（如下图）。
二十万个代理，只需20次迭代就可以从随机混沌的状态演变成同类聚居、异类隔离的格局，由此可见同质性的强大威力。观察稳态的图片不难发现，在不同类的同质块间存在着真空区域，将不同类代理聚居区域隔离开来。
改变模型参数，取20×20的区域，设置两类各15000名代理，设每个代理要求周边有至少5名同类代理，则可以得到如下结果：
（2000次后）
（3000次后）
（5000次后）
（15000次后）
在增强了居住要求之后，同类聚集区域变大了许多，隔离现象也变得非常明显。观察迭代过程不难发现，一开始时会存在多个小聚居区，后来这些小聚居区会逐渐合并成一个大的区域，最终只剩下每类代理各一个超大聚居区。由于居住要求过强，小聚居区边缘出的代理的居住要求往往得不到满足，因此这些代理会向外搬离，并级联地带动聚居区内部的代理也向外搬离，最终汇聚到一个大聚居区中，这很好的验证了谢林模型的迁移连锁效应。不过，同样也是因为居住要求过强，大聚居区边缘区域处于动态平衡状态，模拟过程将会无法收敛，聚居区的大小也存在上界。
取300×300的区域，设两类代理各35000名，设两类代理的居住要求皆为周边8格内有x名同类代理。当x取不同的值时，稳态情况如下：
显然在x取0、1时不发生隔离，在x取2时有很轻微的局部隔离，在x取3时有不明显的隔离，在x取4时有明显的隔离，在x取5时有非常明显的聚集和隔离（难以收敛），在x大于或等于6时不发生隔离。显然，当x过小时，几乎所有代理都满足现状，同质性无法积累导致不发生隔离；而当x过大时，几乎所有代理的居住要求都无法得到满足，由于无法定居而不发生隔离。
总结
1、在要求同类代理占周边8格的约50%时，可以比较容易观察到隔离现象且隔离收敛容易。
2、当阙值要求大于4，模型收敛难度大大提高
3、在阙值比较大时，未必就不会发生隔离。在多次模拟后发现，有小概率在迭代次数很高时突然出现一个某类代理的小聚居区，并以滚雪球之势迅速成长为一个大聚居区。一旦某类代理出现了聚集现象，在多次迭代后其它类型的代理也会更容易出现聚集。
