从经济学到代码：用ABM仿真理解谢林隔离现象（Python实战版）

当托马斯·谢林在1971年用硬币和方格纸演示种族隔离现象时，这位经济学家可能没想到，半个世纪后人们会用Python代码重现他的思想实验。这个看似简单的模型揭示了一个深刻悖论：即使每个个体只表现出轻微的偏好，整个群体也会自发形成高度隔离的模式。今天，我们将用不到100行Python代码，让这个经典的社会学模型在计算机中"活"过来。

1. 谢林模型的核心洞察

1.1 微观行为与宏观现象的鸿沟

谢林模型最颠覆性的发现是：群体层面的隔离程度与个体偏好强度并非线性相关。通过调整几个关键参数，我们会看到三种典型状态：

阈值区间	宏观表现	稳定状态	相似度特征
20%-40%	轻度隔离	静态稳定	略高于随机分布
40%-65%	显著隔离	静态稳定	形成明显聚居区
65%以上	动态混沌	持续流动	实际相似度反而下降

注意：当阈值超过65%时，系统会进入"过度挑剔"状态，居民不断搬家导致无法形成稳定社区。

1.2 模型的社会学启示

• 非意图后果：隔离可能是无意识个体行为的自然结果
• 临界阈值：存在使系统性质突变的敏感参数区间
• 政策启示：仅观察宏观隔离无法推断微观歧视程度

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# 初始化参数
GRID_SIZE = 50
THRESHOLD = 0.4  # 尝试调整这个关键参数观察不同现象

2. Python实现完整流程

2.1 智能体世界的构建

我们首先创建一个二维网格世界，包含三类单元格：

• 红色智能体（占比40%）
• 蓝色智能体（占比40%）
• 空白格（占比20%）

def initialize_grid():
    """创建初始随机分布"""
    status = ["red", "blue", "empty"]
    proportions = [0.4, 0.4, 0.2]
    grid = np.random.choice(status, size=(GRID_SIZE, GRID_SIZE), p=proportions)
    return grid

2.2 满意度计算引擎

每个智能体评估周围8个邻居的构成情况：

def calculate_happiness(grid, i, j):
    """计算指定位置的满意度"""
    if grid[i,j] == "empty":
        return None
    
    same = 0
    different = 0
    for x in range(max(0,i-1), min(GRID_SIZE,i+2)):
        for y in range(max(0,j-1), min(GRID_SIZE,j+2)):
            if (x,y) == (i,j) or grid[x,y] == "empty":
                continue
            if grid[x,y] == grid[i,j]:
                same += 1
            else:
                different += 1
                
    total = same + different
    return same / total if total > 0 else 0

2.3 迁移动态系统

不满意智能体的迁移遵循以下规则：

1. 随机选择一个空白格
2. 交换当前位置与空白格的状态
3. 更新满意度矩阵

def update_grid(grid):
    """执行一轮迁移过程"""
    unhappy_positions = []
    empty_positions = list(zip(*np.where(grid == "empty")))
    
    # 找出所有不满意智能体
    for i in range(GRID_SIZE):
        for j in range(GRID_SIZE):
            happiness = calculate_happiness(grid, i, j)
            if happiness and happiness < THRESHOLD:
                unhappy_positions.append((i,j))
    
    # 随机迁移不满意智能体
    np.random.shuffle(unhappy_positions)
    for i,j in unhappy_positions:
        if empty_positions:
            target = empty_positions.pop()
            grid[i,j], grid[target] = "empty", grid[i,j]
    
    return len(unhappy_positions)

3. 可视化与动态分析

3.1 实时动画呈现

使用Matplotlib的动画功能展示演化过程：

def animate_schelling():
    fig, ax = plt.subplots(figsize=(10,8))
    grid = initialize_grid()
    img = ax.imshow(np.where(grid=="red",1,np.where(grid=="blue",2,0)), 
                   cmap="cool", vmin=0, vmax=2)
    
    def update(frame):
        nonlocal grid
        update_grid(grid)
        img.set_array(np.where(grid=="red",1,np.where(grid=="blue",2,0)))
        return img,
    
    ani = FuncAnimation(fig, update, frames=100, interval=200, blit=True)
    plt.title(f"谢林模型仿真 (阈值={THRESHOLD})")
    plt.show()

3.2 关键指标追踪

建立监控系统记录：

• 平均相似度
• 不满意比例
• 迁移次数

def track_metrics(grid):
    metrics = {
        'similarity': [],
        'unhappy_ratio': [],
        'moves': 0
    }
    
    for _ in range(100):  # 模拟100轮
        moves = update_grid(grid)
        metrics['moves'] += moves
        
        # 计算全局指标
        similarities = []
        unhappy = 0
        total = 0
        
        for i in range(GRID_SIZE):
            for j in range(GRID_SIZE):
                if grid[i,j] != "empty":
                    happiness = calculate_happiness(grid, i, j)
                    similarities.append(happiness)
                    if happiness < THRESHOLD:
                        unhappy += 1
                    total += 1
        
        metrics['similarity'].append(np.mean(similarities))
        metrics['unhappy_ratio'].append(unhappy/total if total >0 else 0)
    
    return metrics

4. 进阶实验设计

4.1 多参数对比实验

通过系统性地调整参数，我们可以绘制出系统的相变图：

def parameter_experiment():
    thresholds = np.linspace(0.2, 0.8, 7)
    results = []
    
    for th in thresholds:
        global THRESHOLD
        THRESHOLD = th
        grid = initialize_grid()
        metrics = track_metrics(grid)
        results.append({
            'threshold': th,
            'final_similarity': metrics['similarity'][-1],
            'total_moves': metrics['moves']
        })
    
    # 绘制结果曲线
    plt.figure(figsize=(12,5))
    plt.subplot(121)
    plt.plot([r['threshold'] for r in results], 
             [r['final_similarity'] for r in results], 'o-')
    plt.xlabel('满意度阈值')
    plt.ylabel('最终平均相似度')
    
    plt.subplot(122)
    plt.plot([r['threshold'] for r in results], 
             [r['total_moves'] for r in results], 's-')
    plt.xlabel('满意度阈值')
    plt.ylabel('总迁移次数')
    plt.tight_layout()

4.2 模型变体探索

• 异质智能体：不同群体设置不同阈值
• 移动成本：引入迁移阻力系数
• 空间约束：限制迁移最大距离

class AdvancedAgent:
    def __init__(self, group, mobility=1.0):
        self.group = group
        self.threshold = np.random.normal(0.5, 0.1)  # 个体差异
        self.mobility = mobility  # 迁移倾向
        
    def decide_move(self, neighborhood):
        similar = sum(1 for n in neighborhood if n == self.group)
        ratio = similar / len(neighborhood)
        return ratio < self.threshold and np.random.rand() < self.mobility

5. 从仿真到现实应用

当我们将这个简单模型与现实数据对比时，会发现惊人的一致性。以芝加哥市的人口分布为例，尽管没有明确的隔离政策，城市依然形成了明显的种族聚居区。模型的预测能力使其在以下领域展现出价值：

• 城市规划：预测公共服务需求分布
• 市场分析：模拟消费者聚集效应
• 政策评估：测试不同干预措施效果

在完成基础实现后，建议尝试将这些代码封装成类结构，添加更多现实因素。例如，引入经济梯度、交通网络等要素，会使模型更加贴近真实城市系统。