谢林模型——R语言实现

谢林模型

参考文章谢林模型

数据准备

• N×N的网格:::一共有pop个格子
  1.红格子，蓝格子 --代表家庭
  2.白格子 --代表搬家的目的地

用N×N的矩阵M表示
红格子：2
蓝格子：1
白格子：0

• 给每个家庭编号：给矩阵M中的每个元素一个编号，方便后面使用

pop.coord， 记录矩阵中每个元素的编号

• 计算每个家庭的周围状态，即左右上下的家庭都是红家庭还是蓝家庭

neighbor.all 记录矩阵中每个元素的周围，即左右上下分别是红格子、蓝格子、白格子

模型演化

• 只针对红格子或者蓝格子进行处理

popLoc 代表了红蓝格子的编号

• 每一轮演化，只针对一部分家庭进行处理，即针对一部分家庭来判断是否搬家

mp = 0.9 ，就是每轮要判断90%的家庭的周围状态，然后判断他们是否搬家
Chosen ，就是那些选出来要做判断的家庭编号，即矩阵元素的编号

• 针对这个家庭的判断，计算他周围跟自己不同类的家庭个数，和同类的家庭个数，前者大，就搬家

M[i]==1 自己是蓝格子
& 并且
length(which(neighbor.all[[i]]==2)) 周围是红格子的数量

length(which(neighbor.all[[i]]==1)) 周围是蓝格子的数量

• 针对这个家庭的判断，计算他周围跟自己同类的家庭个数，然后除以4， 即可得到同类比例，如果小于阈值75%，就搬家

rate <- length(which(neighbor.all[[i]]==M[i])) / 4

• 搬家：自己原来的位置变为空，搬去的目的地变为自己的家庭状态

M[sample(which(M==0),1)]<- 1 随机选一个空格子变为 蓝格子
M[i]<-0 自己原有的位置变为白格子

完整代码

##下载并加载pacman包
if(!require(pacman))install.packages("pacman")
library(pacman)
p_load(plot.matrix,purrr)

#系统参数设置
#网格数
pop<-500
if(sqrt(pop) %% 1 != 0){
  repeat{      
    pop<-pop+1
    if(sqrt(pop) %% 1==0){
      break
    }
  }
}
print(pop)
N<-sqrt(pop)
empty<-0.1   #空格比例为10%
vacant<-as.integer(N*N*empty) ##空格的数量（整数）
vacant

##两类人群比例1：1（蓝色和红色）
B2R <- 1

#设定函数 找到蓝色红色的个数并赋值成1、2，确定空格数并随机排列组合
rand.init<-function(N,B2R,empty){
  blues<- as.integer(pop*(1-empty)/(1+1/B2R))  #强制转换成整数
  reds<- pop-blues
  M<-rep(1,N*N) #蓝色：1
  M[1:reds]<-2 #红色：2
  M[(N*N-vacant+1):(N*N)]<-0 #空格：0
  # set.seed(0)
  M<-sample(M)      #随机排列组合
  return(matrix(M,nrow=N,ncol=N))
}

#随机生成初始网格状态
M<-rand.init(N,B2R,empty)  #N、B2R、empty

init.M<-M

#画图
par(mfrow=c(1,1))                     #一行一列
plot(init.M,breaks=c(0,1,2,3),        #breaks是0、1、2、3四个数字
     col=c("white","blue","red"))     #颜色是蓝色、红色、白色

#确定网格中某个家庭的周围状态， mx是当前的网格
neighbor<-function(row,col,mx){
  if(row==1){U<-mx[sqrt(pop),col]} else{U<-mx[row-1,col]}  #sqrt(pop)=N,如果在第一行的话，上面就是第23行
  if(row==sqrt(pop)){D<-mx[1,col]} else{D<-mx[row+1,col]}
  if(col==1){L<-mx[row,sqrt(pop)]} else{L<-mx[row,col-1]}
  if(col==sqrt(pop)){R<-mx[row,1]} else{R<-mx[row,col+1]}
  c(L,R,U,D) # 左右上下
}

#网格家庭坐标编号
pop.coord<-cbind(rep(c(1:sqrt(pop)),sqrt(pop)),
                 as.vector(matrix(rep(1:sqrt(pop),sqrt(pop)),
                                  nrow=sqrt(pop),byrow=T)))
pop.coord

#确定网格上每个家庭的周围状态  是一个list,  其index：家庭编号， 元素是list， 代表了这个家庭的周围状态。 
neighbor.all<-map2(pop.coord[,1],pop.coord[,2],neighbor,mx=M)


mp<-0.9 #每期移动比例
for(r in 1:100){ #运行100次
  popLoc<-which(M==1 | M==2) #筛除空格（只留下代表家庭的格子）
  # set.seed(r)
  Chosen<-sample(popLoc,mp*length(popLoc))  #打乱家庭的位置，即随机选一定比例的家庭进行移动
  for(i in Chosen){        #每个家庭循环
    rate <- length(which(neighbor.all[[i]]==M[i])) / 4
    if(rate < 0.75){
        M[sample(which(M==0),1)]<- M[i]
        M[i]<-0
    }
    # if(M[i]==1 & length(which(neighbor.all[[i]]==2))>length(which(neighbor.all[[i]]==1))){ # 自己家是蓝色，且周围是红色的家庭大于蓝色的家庭，搬家
    #   M[sample(which(M==0),1)]<-1
    #   M[i]<-0
    # }
    # if(M[i]==2 & length(which(neighbor.all[[i]]==1))>length(which(neighbor.all[[i]]==2))){
    #   M[sample(which(M==0),1)]<-2
    #   M[i]<-0
    # }
   neighbor.all<-map2(pop.coord[,1],pop.coord[,2],neighbor,mx=M)
  }
}

#可视化
par(mfrow=c(1,2))    #一行2列
plot(init.M,breaks=c(0,1,2,3),col=c("white","blue","red"))
plot(M,breaks=c(0,1,2,3),col=c("white","blue","red"))