基于蒙特卡洛模拟复现谢林模型（计算机社会学）

谢林模型引言

在谢林模型中，有一片方形区域，区域又被均匀分成若干小方格；有一群被称为代理的个体，每个代理居住在一个小方格内。这些代理可以分成若干类，所有代理都希望周边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、在阙值比较大时，未必就不会发生隔离。在多次模拟后发现，有小概率在迭代次数很高时突然出现一个某类代理的小聚居区，并以滚雪球之势迅速成长为一个大聚居区。一旦某类代理出现了聚集现象，在多次迭代后其它类型的代理也会更容易出现聚集。