谢林隔离模型Matlab仿真实现 The Implementation of Schelling Model of Segregation With Matlab

简介

该模型由马里兰大学的经济学家托马斯·谢林建立。主要感兴趣的隔离有两种：一种是种族隔离，另一种是由收入差异导致的隔离。
谢林构建了一个基于代理人的模型。该类模型要点有三个：代理人，规则，以及群体表现。
谢林模型是个关于人们选择在哪居住的模型。他把每个人放置在棋盘上，整个城市看作是一块巨大的棋盘。棋盘上每一个小格子可以住一个人也可以空着。
假设红格子表示住着富人，灰格子表示住了穷人，白格子表示没有主人。对在中间X格子的富人来说，周围7个邻居有3个是同类，他应该留下还是搬走？

谢林设定了规则称作阈值基本准则：每个人都有一个阈值（threshold），根据阈值做决定是留下还是搬走。比如假设X的阈值是33%，则此时他会选择留下。而当他的富人邻居有一位搬走时，由于只剩2/7是同类所以那时会选择搬走。
（上述内容引自此处）

对于谢林隔离模型，我们重点关心在不同参数下其是否收敛，若收敛则形态如何这样的问题。

仿真程序介绍

该程序由Matlab写成，在图形化界面中含有分布图示、参数调整滑动条、各轮满意度走势的统计图和各轮参数统计表四部分。

参数表示

名称 含义 similarity threshold 相似度阈值，低于则需要搬迁 red/blue ratio 红蓝像素比，此处即两种人群的比重 size 像素点个数（正方形） empty ratio 空置率，无人的像素点数

部分仿真与分析

• 以相似度阈值（容忍度）为自变量
  similarity threshold=30%, red/blue ratio=50%/50%,size=30*30,empty ratio=10%

初始化状态如下，此时满意度为92.72%。

经过9轮搬迁，满意度达到100%，可以看出相同颜色像素点有“连成片”的现象。

similarity threshold=50%, red/blue ratio=50%/50%,size=30*30,empty ratio=10%
在相似度阈值设定为50%，初始满意度为57.41%。

18轮之后，满意度达到100%，可以看出非常明显的聚集现象。刚开始搬动频繁，越往后变化越小。

similarity threshold=80%, red/blue ratio=50%/50%,size=30*30,empty ratio=10%

相似度阈值达到80%时，初始满意度只有12.47%。

246轮后，出现了非常明显的红/蓝分区，事实上对于大部分像素点，其相似度阈值已经达到了100%。这里的结果与程序采用的随机算法有关，在随机选择新一轮位置时，待搬迁像素的可选择范围不仅包括当前空白区域，也包括因搬迁而空出的区域。在极端情况下，谢林模型应该会呈现无法收敛的状态。

• 以红/蓝比为自变量（种族优势）
  similarity threshold=50%, red/blue ratio=80%/20%,size=30*30,empty ratio=10%

150轮后，仍然得不到100%满意，当一种颜色占主导地位时，仍然有聚集的特征，但也有明显的特殊情况（混居）。

• 空置率提高
  similarity threshold=50%, red/blue ratio=50%/50%,size=30*30,empty ratio=50%

尽管中间有很多缓冲带，但是仍然会出现聚集情况。非常容易想到，空置率的提高会使开始的收敛速度加快，本次仅13轮。

部分核心代码

% 以下代码省略了控制各个控件的响应函数及起始函数，库来自MATLAB R2016afunction step_btn_Callback(hObject, eventdata, handles)% hObject    handle to step_btn (see GCBO)% eventdata  reserved - to be defined in a future version of MATLAB% handles    structure with handles and user data (see GUIDATA)global size; %正方形区域边长（像素点个数）global empty; %空置率global ratio; %红色（即富人）占比global sim; %相似度容忍阈值global R; %分布矩阵global count; %轮数计算用变量（ROUND count）global map; %涂色用颜色表m = size * size; %计算像素点（居民）数量%根据百分比，计算像素点数量empty_tmp = round(m * empty / 100); ratio_tmp = round((m - empty_tmp) * ratio / 100);sim_tmp = sim / 100;str_tmp = get(handles.round_text,'String'); %获取ROUND NUMBER的字符串not_satis_red=[]; %初始化不满意红色点序列not_satis_blue=[];%初始化不满意蓝色点序列empty_now=[];%初始化空白像素点序列%R(i,j)=1代表空置，2代表红点，3代表蓝点for i=1:size    for j=1:size        switch R(i,j)            case 1                empty_now = [empty_now (j-1)*size+i];            case 2                satis_tmp = 0;                amount_tmp = 0;                %检测不满意情况                    if i ~= 1 && R(i-1,j) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= 1 && j ~= 1 && R(i-1,j-1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= 1 && j ~= size && R(i-1,j+1) == R(i,j) satis_tmp = satis_tmp + 1; end                if j ~= 1 && R(i,j-1) == R(i,j) satis_tmp = satis_tmp + 1; end                if j ~= size && R(i,j+1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= size && j ~= 1 && R(i+1,j-1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= size && R(i+1,j) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= size && j ~= size && R(i+1,j+1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= 1 && j ~= 1 && R(i-1,j-1) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= 1 && R(i-1,j) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= 1 && j ~= size && R(i-1,j+1) ~= 1 amount_tmp = amount_tmp + 1; end                if j ~= 1 && R(i,j-1) ~= 1 amount_tmp = amount_tmp + 1; end                if j ~= size && R(i,j+1) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= size && j ~= 1 && R(i+1,j-1) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= size && R(i+1,j) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= size && j ~= size && R(i+1,j+1) ~= 1 amount_tmp = amount_tmp + 1; end                if satis_tmp<round(sim_tmp*amount_tmp)                    not_satis_red =[not_satis_red (j-1)*size+i];                end            case 3                satis_tmp = 0;                amount_tmp = 0;                if i ~= 1 && R(i-1,j) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= 1 && j ~= 1 && R(i-1,j-1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= 1 && j ~= size && R(i-1,j+1) == R(i,j) satis_tmp = satis_tmp + 1; end                if j ~= 1 && R(i,j-1) == R(i,j) satis_tmp = satis_tmp + 1; end                if j ~= size && R(i,j+1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= size && j ~= 1 && R(i+1,j-1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= size && R(i+1,j) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= size && j ~= size && R(i+1,j+1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= 1 && j ~= 1 && R(i-1,j-1) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= 1 && R(i-1,j) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= 1 && j ~= size && R(i-1,j+1) ~= 1 amount_tmp = amount_tmp + 1; end                if j ~= 1 && R(i,j-1) ~= 1 amount_tmp = amount_tmp + 1; end                if j ~= size && R(i,j+1) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= size && j ~= 1 && R(i+1,j-1) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= size && R(i+1,j) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= size && j ~= size && R(i+1,j+1) ~= 1 amount_tmp = amount_tmp + 1; end                if satis_tmp<round(sim_tmp*amount_tmp)                    not_satis_blue =[not_satis_blue (j-1)*size+i];                end        end    endend%pool是不满意像素点和空白点放在一起的像素点池pool = [not_satis_red not_satis_blue empty_now];index = randperm(length(pool)); %index是与pool数量相同的由1~pool的随机序列pool = pool(index); %pool中元素按照index索引打乱排列，实现随机分配R=reshape(R,1,m);  %将R矩阵转换为1*m矩阵，方便将随机分配后的序列放回for i=1:length(not_satis_red)    R(1,pool(1,i)) = 2;endfor i=1:length(not_satis_blue)    R(1,pool(1,i+length(not_satis_red))) = 3;endfor i=1:length(empty_now)    R(1,pool(1,i+length(not_satis_red)+length(not_satis_blue))) = 1;endR = reshape(R,size,size); %将R矩阵恢复size*size矩阵，为绘图做准备axes(findobj('tag','axes1'));%将R在axes1中绘制出来，颜色表为mapimage(R);colormap(map);%显示轮数ROUNDcount = count + 1;set(handles.round_text,'String', ['ROUND ' int2str(count)]);%检测本轮结束后的满意度satisfied = m - empty_tmp;for i=1:size    for j=1:size        switch R(i,j)            case 2                satis_tmp = 0;                amount_tmp = 0;                if i ~= 1 && R(i-1,j) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= 1 && j ~= 1 && R(i-1,j-1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= 1 && j ~= size && R(i-1,j+1) == R(i,j) satis_tmp = satis_tmp + 1; end                if j ~= 1 && R(i,j-1) == R(i,j) satis_tmp = satis_tmp + 1; end                if j ~= size && R(i,j+1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= size && j ~= 1 && R(i+1,j-1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= size && R(i+1,j) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= size && j ~= size && R(i+1,j+1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= 1 && j ~= 1 && R(i-1,j-1) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= 1 && R(i-1,j) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= 1 && j ~= size && R(i-1,j+1) ~= 1 amount_tmp = amount_tmp + 1; end                if j ~= 1 && R(i,j-1) ~= 1 amount_tmp = amount_tmp + 1; end                if j ~= size && R(i,j+1) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= size && j ~= 1 && R(i+1,j-1) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= size && R(i+1,j) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= size && j ~= size && R(i+1,j+1) ~= 1 amount_tmp = amount_tmp + 1; end                if satis_tmp<round(sim_tmp*amount_tmp)                    satisfied = satisfied - 1;                end            case 3                satis_tmp = 0;                amount_tmp = 0;                if i ~= 1 && R(i-1,j) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= 1 && j ~= 1 && R(i-1,j-1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= 1 && j ~= size && R(i-1,j+1) == R(i,j) satis_tmp = satis_tmp + 1; end                if j ~= 1 && R(i,j-1) == R(i,j) satis_tmp = satis_tmp + 1; end                if j ~= size && R(i,j+1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= size && j ~= 1 && R(i+1,j-1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= size && R(i+1,j) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= size && j ~= size && R(i+1,j+1) == R(i,j) satis_tmp = satis_tmp + 1; end                if i ~= 1 && j ~= 1 && R(i-1,j-1) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= 1 && R(i-1,j) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= 1 && j ~= size && R(i-1,j+1) ~= 1 amount_tmp = amount_tmp + 1; end                if j ~= 1 && R(i,j-1) ~= 1 amount_tmp = amount_tmp + 1; end                if j ~= size && R(i,j+1) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= size && j ~= 1 && R(i+1,j-1) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= size && R(i+1,j) ~= 1 amount_tmp = amount_tmp + 1; end                if i ~= size && j ~= size && R(i+1,j+1) ~= 1 amount_tmp = amount_tmp + 1; end                if satis_tmp<round(sim_tmp*amount_tmp)                   satisfied = satisfied - 1;                end        end    endendset(handles.sat_text,'String',['SATISFIED ' num2str(satisfied/(m-empty_tmp)*100,'%.2f') '%']); %用统计图选择数据global datas;datas = [datas satisfied/(m-empty_tmp)*100];axes(findobj('Tag','axes2'));plot(datas); %存储本轮记录global record;record(end,5) = satisfied/(m-empty_tmp)*100;if satisfied/(m-empty_tmp) == 1    record(end,6) = count;    set(handles.step_btn,'Enable','off');endset(findobj('Tag','record1'),'Data',record);end

GitHub

https://github.com/sqyx008/Schelling-Model-of-Segregation