タカハトゲーム 

タカ・ハトゲーム
タカ戦略とハト戦略のみが存在する世界において進化的に安定な状態を探るゲーム。
### 設定
# pay-offは以下の通り:
タカ vs. ハト、餌(food)を総取り
タカ vs. タカ、餌からコストを引いたものを2分割
ハト vs. タカ、無し(0)
ハト vs. ハト、餌を2分割
# 引数:
個体数の初期値([タカ、ハト])、世代数、コスト、餌
# 個体数の推移
p.hawk[i_gen]=(p.hawk[i_gen-1]*(p.hawk[i_gen-1]*eHawkHawk+p.dove[i_gen-1]*eHawkDove))
 /mean.W;

HDgame<-function(N,generation,cost,food){
  N.hawk=N[1]; N.dove=N[2]; 
  p.hawk=rep(0,generation);p.hawk[1]=N.hawk/(N.dove+N.hawk)
  p.dove=rep(0,generation);p.dove[1]=1-p.hawk[1]
  eHawkHawk=0.5*(food-cost);eHawkDove=food
  eDoveHawk=0;eDoveDove=0.5*food
  for (i_gen in 2:generation){
    mean.W=p.hawk[i_gen-1]^2*eHawkHawk+p.hawk[i_gen-1]*p.dove[i_gen-1]*
     (eHawkDove+eDoveHawk)+p.dove[i_gen-1]^2*eDoveDove
    p.hawk[i_gen]=(p.hawk[i_gen-1]*(p.hawk[i_gen-1]*eHawkHawk
      +p.dove[i_gen-1]*eHawkDove))/mean.W;
    p.dove[i_gen]=1-p.hawk[i_gen]
  }
  plot(1:generation,p.hawk,type='o',col='red',ylim=c(0,1),pch=19,
    main="Result of Hawk-Dove Game",ylab="Proportion",xlab="Generation")
  lines(1:generation, p.dove,type='o',col='blue',pch=17)
  legend("topright",c("Hawk","Dove"),pch=c(19,17),col=c('red','blue'))
}
# 実行例
HDgame(c(20,80),50,10,6)

# 実装例、その2
HDgameDE<-function(N,generation,cost,food){
  N.hawk=N[1]; N.dove=N[2]; 
  tStep=0.01;ts=seq(0,generation,tStep);Nts=length(ts);
  p.hawk=rep(0,Nts);p.hawk[1]=N.hawk/(N.dove+N.hawk)
  p.dove=rep(0,Nts);p.dove[1]=1-p.hawk[1]
  eHawkHawk=0.5*(food-cost);eHawkDove=food
  eDoveHawk=0;eDoveDove=0.5*food
  for (i_gen in 2:Nts){
    mean.W=p.hawk[i_gen-1]^2*eHawkHawk+p.hawk[i_gen-1]*p.dove[i_gen-1]*
     (eHawkDove+eDoveHawk)+p.dove[i_gen-1]^2*eDoveDove
    p.hawk[i_gen]=p.hawk[i_gen-1]+(p.hawk[i_gen-1]*(p.hawk[i_gen-1]*eHawkHawk
     +p.dove[i_gen-1]*eHawkDove-mean.W)/mean.W)*tStep;
    p.dove[i_gen]=1-p.hawk[i_gen]
  }
  plot(1:Nts,p.hawk,type='l',col='red',ylim=c(0,1),lwd=5,
   main="Result of Hawk-Dove Game",ylab="Proportion",xlab="Generation")
  lines(1:Nts, p.dove,type='l',col='blue',lwd=5)
  legend("topright",c("Hawk","Dove"),lwd=5,col=c('red','blue'))
}
# 実行例
HDgameDE(c(20,80),50,10,6)

D_H01

読みやすい?例

## a simple version ##
HDgame<-function(N,food,cost,n_gen){
  p.hawk=N[1]/sum(N);p.dove=N[2]/sum(N)
  e.HH=1/2*(food-cost);e.HD=food;e.DH=0;e.DD=1/2*(food)
  p.vec=c(p.hawk,p.dove);p.mat<-outer(p.vec,p.vec)
  payMat<-matrix(c(e.HH,e.DH,e.HD,e.DD),ncol=2)
  hist=matrix(0,n_gen,2);hist[1,]=p.vec
  for (i_gen in 2:n_gen){
    w.H=sum(p.vec*payMat[1,])
    w.D=sum(p.vec*payMat[2,])
    w.mean=sum(p.mat*payMat)
    p.vec=c(p.vec[1]*w.H/w.mean,p.vec[2]*w.D/w.mean)
    p.mat<-outer(p.vec,p.vec)
    hist[i_gen,]=p.vec
  }
  plot(1:n_gen,hist[,1],pch=20,type='o',lwd=2,cex=2,col='red',
   ylim=c(0,1),ylim=c(0,1.25),xlab="time",ylab="proportion",cex.lab=2)
  lines(1:n_gen,hist[,2],pch=20,type='o',lwd=2,cex=2,col='blue')
  legend("topright",c("Hawk","Dove"),col=c('red','blue'),lwd=2,pch=20)
}

## DE version ##
HDgameDE<-function(N,food,cost,n_gen){
  p.hawk=N[1]/sum(N);p.dove=N[2]/sum(N)
  e.HH=1/2*(food-cost);e.HD=food;e.DH=0;e.DD=1/2*(food)
  p.vec=c(p.hawk,p.dove);p.mat<-outer(p.vec,p.vec)
  payMat<-matrix(c(e.HH,e.DH,e.HD,e.DD),ncol=2)
  ts=0.01;Nts=length(seq(1,n_gen,ts));
  hist=matrix(0,Nts,2);hist[1,]=p.vec
  for (i_gen in 2:Nts){
    w.H=sum(p.vec*payMat[1,])
    w.D=sum(p.vec*payMat[2,])
    w.mean=sum(p.mat*payMat)
    p.vec[1]=p.vec[1]+ts*p.vec[1]*(w.H-w.mean)/w.mean
    p.vec[2]=p.vec[2]+ts*p.vec[2]*(w.D-w.mean)/w.mean
    p.mat<-outer(p.vec,p.vec)
    hist[i_gen,]=p.vec
  }
  plot(1:Nts,hist[,1],type='l',lwd=4,cex=2,col='red',ylim=c(0,1.25),
    xlab="time",ylab="proportion",cex.lab=2)
  lines(1:Nts,hist[,2],type='l',lwd=4,cex=2,col='blue')
  legend("topright",c("Hawk","Dove"),col=c('red','blue'),lwd=4)
  return(hist)
}

認知情報解析 デマの拡散モデル

Dema_WO_cr<-function(N,ps,ave.followers,n_rep){
  # N - numbers of S, Ig, Is
  # ps - probs for S->I, Ig->I
  dt=0.01;ts=seq(1,n_rep,dt);n_ts=length(ts);
  S=matrix(0,nrow=n_ts,ncol=1);S[1]=N[1];
  Ig=matrix(0,nrow=n_ts,ncol=1);Ig[1]=N[2];
  I=matrix(0,nrow=n_ts,ncol=1);I[1]=N[3];
  N=S[1]+Ig[1]+I[1];
  F=ave.followers;
  p.S2I=ps[1];p.Ig2I=ps[2]

  # main 
  for (i_time in 1:(n_ts-1)) {
    S[i_time+1]=S[i_time]-F/N*I[i_time]*S[i_time]*dt
    Ig[i_time+1]=Ig[i_time]+((1-p.S2I)*F/N*I[i_time]*S[i_time]-p.Ig2I*F/N*Ig[i_time]*I[i_time])*dt
    I[i_time+1]=I[i_time]+(p.S2I*F/N*I[i_time]*S[i_time]+p.Ig2I*F/N*Ig[i_time]*I[i_time])*dt
  }
  # plotting results
  plot(ts,S,type="l",lwd=4,col="red", main="Distribution of rumors in Twitter "
    ,xlab="time", ylab="Proportions of People",cex.lab=1.3,ylim=c(0,S[1]))
  lines(ts,Ig,type="l",lwd=4,col="blue")
  lines(ts,I,type="l",lwd=4,col="green")
  legend("right", c("Not seen rumors nor corrections","Seen rumors", "tweeted rumours"), 
  col=c("red","blue","green"),lty=rep(1,3),lwd=4)
}

Dema_W_cr<-function(N,ps,ave.followers,n_rep){
  # N - numbers of S, Ig, Is
  # ps - probs for S->I, Ig->I
  dt=0.01;ts=seq(1,n_rep,dt);n_ts=length(ts);
  S=matrix(0,nrow=n_ts,ncol=1);S[1]=N[1];
  Ig=matrix(0,nrow=n_ts,ncol=1);Ig[1]=N[2];
  I=matrix(0,nrow=n_ts,ncol=1);I[1]=N[3];
  Rg=matrix(0,nrow=n_ts,ncol=1);Rg[1]=N[4];
  R=matrix(0,nrow=n_ts,ncol=1);R[1]=N[5];
  N=S[1]+Ig[1]+I[1]+Rg[1]+R[1];
  F=ave.followers;
  p.S2I=ps[1];p.Ig2I=ps[2];p.S2R=ps[3];p.Ig2R=ps[4];p.I2R=ps[5];p.Rg2R=ps[6];

  # main 
  for (i_time in 1:(n_ts-1)) {
    S[i_time+1]=S[i_time]+(-F/N*I[i_time]*S[i_time]-F/N*R[i_time]*S[i_time])*dt
    Ig[i_time+1]=Ig[i_time]+((1-p.S2I)*F/N*I[i_time]*S[i_time]-p.Ig2I*F/N*Ig[i_time]*I[i_time]
      -F/N*Ig[i_time]*R[i_time])*dt
    I[i_time+1]=I[i_time]+(p.S2I*F/N*I[i_time]*S[i_time]+p.Ig2I*F/N*Ig[i_time]*I[i_time]
      -F/N*I[i_time]*R[i_time])*dt
    Rg[i_time+1]=Rg[i_time]+((1-p.S2R)*F/N*R[i_time]*S[i_time]+(1-p.Ig2R)*F/N*Ig[i_time]*R[i_time]
      +(1-p.I2R)*F/N*I[i_time]*R[i_time]-p.Rg2R*F/N*Rg[i_time]*R[i_time])*dt
    R[i_time+1]=R[i_time]+(p.S2R*F/N*R[i_time]*S[i_time]+p.Ig2R*F/N*Ig[i_time]*R[i_time]
      +p.I2R*F/N*I[i_time]*R[i_time]+p.Rg2R*F/N*Rg[i_time]*R[i_time])*dt
  }
  # plotting results
  plot(ts,S,type="l",lwd=4,col="red", main="Distribution of rumors in Twitter "
    ,xlab="time", ylab="Proportions of People",cex.lab=1.3,ylim=c(0,S[1]*1.25))
  lines(ts,Ig,type="l",lwd=4,col="blue")
  lines(ts,I,type="l",lwd=4,col="green")
  lines(ts,Rg,type="l",lwd=4,col="magenta")
  lines(ts,R,type="l",lwd=4,col="cyan")
  legend("topright", c("Not seen rumors nor corrections","Seen rumors", "tweeted rumours",
    "Seen corrections","tweeted corrections"), 
  col=c("red","blue","green","magenta","cyan"),lty=rep(1,5),lwd=4)
  return(data.frame(S,Ig,I,Rg,R))
}

w=Dema_W_cr(c(10000,0,100,0,1),c(0.001,0.01,0.01,0.01,0.01,0.01),20,50)

twitterDema