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)
Category Archives: R companion
数理社会学:軍拡モデル
軍拡モデル
2国間の軍事拡大モデル
### 設定
# parameters
a: X国の自国の戦力・軍事に対する抑制率
k: Y国に戦力に対するX国の戦力の増加率
g: X国の定常的な軍事費(戦力の増加量)
b: Y国の自国の戦力・軍事に対する抑制率
l: X国に戦力に対するY国の戦力の増加率
h: Y国の定常的な軍事費(戦力の増加量)
# 戦力の推移設定
X国:dx/dt=-ax+ky+g
y国:dy/dt=lx-by+h
# 引数:
a,b,k,l,g,h (上記)
dt 時間の推移
mil_exp<-function(a,b,k,l,g,h,dt) { timeSep=dt;ts=seq(1,50,timeSep);n_ts=length(ts) x=matrix(0,nrow=n_ts,ncol=1);y=matrix(0,nrow=n_ts,,ncol=1) initX=c(rep(-40,5),rep(-20,5),rep(0,5),rep(20,5),rep(40,5)) initY=c(rep(c(-40,-20,0,20,40),5)) initX=initX[-13];initY=initY[-13] lengthINI=length(initX) for (i_ini in 1:lengthINI) { x[1]=initX[i_ini];y[1]=initY[i_ini]; for (i_gen in 2:n_ts) { x[i_gen]=x[i_gen-1]+(-a*x[i_gen-1]+k*y[i_gen-1]+g)*timeSep y[i_gen]=y[i_gen-1]+(l*x[i_gen-1]-b*y[i_gen-1]+h)*timeSep } if (i_ini==1) { plot(x,y,xlim=c(-40,40),ylim=c(-40,40),col=4,type='l',lwd=2, xlab="Nation X's Military Force",ylab="Nation X's Military Force") arrows(x[2],y[2],x[3],y[3],col=4,lwd=2,length=0.15) } else { lines(x, y, col=4, lwd=2) arrows(x[2], y[2], x[3], y[3], col=4,lwd=2, length=0.15) } } } mil_exp(1,2,0.5,0.1,10,10,0.05)
セルオートマトン(2次元)
セルオートマトン
2次元のセルオートマトン
# ルール:
(1) cell[i,j]が0であり、かつ、周りの8つのセルに「1」の状態であるセルが2つある時に、cell[i,j]の次世代での状態は「1」となる
(2) cell[i,j]が1であり、かつ、周りの8つのセルに「1」の状態であるセルが1つか2つある時に、cell[i,j]の次世代での状態は「1」のままである
(3) (1)と(2)以外は全て「0」とする
解答例
par(mfrow=c(7,7)); par(mai=c(1/10,1/10,1/10,1/10)); nGen=48;nR=75;nC=75; field=matrix(0,nR,nC);field[37:38,37:38]=1 image(field,xaxt='n',yaxt='n') for (i_gen in 1:nGen){ field2=cbind(field[,nC],field,field[,1]) field3=rbind(field2[nR,],field2,field2[1,]) field.temp=matrix(0,nrow=nR+2,ncol=nC+2); for (i_row in 2:(nR+1)){ for (i_col in 2:(nC+1)){ nOnes=sum(field3[i_row-1,(i_col-1):(i_col+1)])+ (field3[i_row,i_col-1]+field3[i_row,i_col+1])+ sum(field3[i_row+1,(i_col-1):(i_col+1)]) if (nOnes==2) { field.temp[i_row,i_col]=1 } if (nOnes==1 & field3[i_row,i_col]==1) { field.temp[i_row,i_col]=1} } } field=field.temp[2:(nR+1),2:(nC+1)] image(field,xaxt='n',yaxt='n') }
セルオートマトン(1次元)
セルオートマトン
nCell=201;nGen=100; res=matrix(0,nrow=nGen,ncol=nCell) res[1,ceiling(nCell/2)]=1; for (i_gen in 2:nGen) { temp.vec=c(res[i_gen-1,nCell],res[i_gen-1,],res[i_gen-1,1]) res[i_gen,]=(temp.vec[1:nCell]+temp.vec[2:(nCell+1)]+temp.vec[3:(nCell+2)])%%2; } # その1R(直感的に理解しやすい例) for (i_gen in 2:nGen) { for (i_cell in 2:(nCell-1)) { res[i_gen,i_cell]=sum(res[i_gen-1,(i_cell-1):(i_cell+1)])%%2; } res[i_gen,1]=(res[i_gen-1,nCell]+res[i_gen-1,1]+res[i_gen-1,2])%%2; res[i_gen,nCell]=(res[i_gen-1,(nCell-1)]+res[i_gen-1,nCell]+res[i_gen-1,1])%%2; }
解答例その2 256のルールに対応
dec2bin<-function(num, digits=8) { bin=c() if (num==0){ bin=0 } else { while(num!=0){ rem=num%%2 num=num%/%2 bin=c(rem,bin) } } if (length(bin) < digits){ res=matrix(0,nrow=1,ncol=digits) res[(digits-length(bin)+1):digits]=bin } else {res=bin} return(res) } transFUN<-function(st,ruleID){ output=dec2bin(ruleID,8); a=matrix(c(1,1,1,1,1,0,1,0,1,1,0,0,0,1,1,0,1,0,0,0,1,0,0,0),nrow=8,byrow=T) newSt=output[which(apply(a,1,function(x) {all(x==st)}))] return(newSt) } ECA<-function(nCell, nGen,ruleID){ res=matrix(0,nrow=nGen,ncol=nCell) res[1,ceiling(nCell/2)]=1; for (i_gen in 2:nGen) { for (i_cell in 2:(nCell-1)) { res[i_gen,i_cell]=transFUN(res[i_gen-1,(i_cell-1):(i_cell+1)],ruleID) } res[i_gen,1]=transFUN(c(res[i_gen-1,nCell], res[i_gen-1,1], res[i_gen-1,2]),ruleID) res[i_gen,nCell]=transFUN(c(res[i_gen-1,(nCell-1)], res[i_gen-1,nCell], res[i_gen-1,1]),ruleID) } return(res) } # その2R (計算時間の改善はほぼなし) ECA2<-function(nCell, nGen,ruleID){ res=matrix(0,nrow=nGen,ncol=nCell) res[1,ceiling(nCell/2)]=1; for (i_gen in 2:nGen) { temp.mat=cbind(c(res[i_gen-1,nCell],res[i_gen-1,1:(nCell-1)]),res[i_gen-1,], c(res[i_gen-1,2:nCell],res[i_gen-1,1])) res[i_gen,]=apply(temp.mat,1,transFUN,ruleID) } return(res) }
数理社会学: なぜ協力するのか
社会を<モデル>でみる:数理社会学への招待
09章:なぜ裏切ったほうが得なのに協力するのか
繰り返し囚人のジレンマゲーム
payoff matrix =
-20(defect,defect) -25(coop,defect)
0(defect,coop) -5(coop, coop)
# player1 0:defect 1:coop # player2 -10:defect 10:coop # -10:defDef, 10:defCoop # 9:coopDef 11:coopCoop calcPO<-function(points,payMat){ n_play<-length(points) pO<-matrix(0,nrow=2,ncol=n_play) pO[,points==-10]=c(payMat[1,1],payMat[1,1]) pO[,points==10]=c(payMat[2,1],payMat[1,2]) pO[,points==-9]=c(payMat[1,2],payMat[2,1]) pO[,points==11]=c(payMat[2,2],payMat[2,2]) return(t(pO)) } # env set-up n_play=1000000 payMat=matrix(c(-20,0,-25,-5),nrow=2) result=matrix(0,4,4) colnames(result)<-c("random","tit-for-tat","defector","cooperator") rownames(result)<-c("random","tit-for-tat","defector","cooperator") # scenario01 random vs. random p1<-sample(c(0,1),n_play,replace=T) p2<-sample(c(-10,10),n_play,replace=T) pO<-calcPO(rowSums(cbind(p1,p2)),payMat) result[1,1]=mean(colSums(pO)) # scenario02 random vs. tit for tat p1<-sample(c(0,1),n_play,replace=T) p2<-c(10,p1[1:99]*20-10) pO<-calcPO(rowSums(cbind(p1,p2)),payMat) result[2,1]=sum(pO[,1]) result[1,2]=sum(pO[,2]) # scenario03 random vs. defector p1<-sample(c(0,1),n_play,replace=T) p2<-rep(-10,n_play) pO<-calcPO(rowSums(cbind(p1,p2)),payMat) result[3,1]=sum(pO[,1]) result[1,3]=sum(pO[,2]) # scenario04 random vs. cooperator p1<-sample(c(0,1),n_play,replace=T) p2<-rep(10,n_play) pO<-calcPO(rowSums(cbind(p1,p2)),payMat) result[4,1]=sum(pO[,1]) result[1,4]=sum(pO[,2]) # scenario05&07&10 {tit for tat}{cooperator} vs. {tit for tat}{cooperator} pO<-calcPO(rep(11,n_play),payMat) result[2,2]=sum(pO[,1]) result[2,4]=sum(pO[,2]) result[4,2]=sum(pO[,2]) result[4,4]=sum(pO[,2]) # scenario06 tit for tat vs. defector pO<-calcPO(c(-9,rep(-10,n_play-1)),payMat) result[3,2]=sum(pO[,1]) result[2,3]=sum(pO[,2]) # scenario08 defector vs. defector pO<-calcPO(rep(-10,n_play),payMat) result[3,3]=sum(pO[,1]) # scenario09 defector vs. cooperator pO<-calcPO(rep(10,n_play),payMat) result[4,3]=sum(pO[,1]) result[3,4]=sum(pO[,2]) > colMeans(result) random tit-for-tat defector cooperator -12349020 -10660476 -12504475 -12498240
数理社会学・なぜタバコがやめられないか
社会を<モデル>でみる:数理社会学への招待
01章:なぜタバコがやめられないか
モデルの説明:
べッカーの依存症モデル
効用関数
仮定1:今日の消費量x(t)は、昨日の消費量x(t-1)に依存する
仮定2:時間が経てば魅力が減る(0<=C<=1)
仮定3:満足度はu*(x(t),C*x(t-1))で表し、依存状態はもっとも高い効用uを得るとする
仮定4:はじめはあまり依存せず、途中にはまりだし、最後には慢性的にはまる
つまり、昨日の消費が増えるにしたがって最大の満足を得るには今日もたくさん消費する必要がある。しかし、消費が増えると魅力が減りそれほど増えなくなる。
#実装例(ベッカーのモデルを忠実には再現できていない)
今日の消費量をy、昨日の消費量をxとした場合、それらの関係を以下の式で表す。
y=K/(1+exp(r*K*(shift-x)))
rをgainと呼び、今日と昨日の消費量の変化率である(rが大きいと昨日と同じ満足度を得るにはより多くの消費が必要)
Kは最大消費量(あまり意味なない)
shiftは今日の消費量に対する抑制・ペナルティー効果
smoker<-function(K=20,r=0.02,shift=9,nDays=10) { #------ input args ------------------- #K: max # of consumption #r: satisfaction decay #shift: shift in satisfaction function #initial amount of consumption = 10 #----------------------------------------- x.seq=seq(0,K,0.01); y.seq=K/(1+exp(r*K*(shift-x.seq))) plot(x.seq,y.seq,type="l",lwd=2,xlab=c("consumption at T"), ylab=c("consumption at T+1"),xlim=c(-1,21),ylim=c(-1,21)) lines(x.seq,x.seq,type="l",col="blue",lty="dashed",lwd=2) x=10; y=K/(1+exp(r*K*(shift-x))) lines(c(x,x),c(0,y),type="l",col="red",lwd=2) lines(c(x,y),c(y,y),type="l",col="red",lwd=2) for (i_days in 2:nDays) { x.old=x;x=y; y=K/(1+exp(r*K*(shift-x))) lines(c(x,x),c(x,y),type="l",col="red",lwd=2) lines(c(x,y),c(y,y),type="l",col="red",lwd=2) } } # example par(mfrow=c(2,2)) smoker(K=20,r=0.02,shift=9,nDays=10) smoker(K=20,r=0.01,shift=9,nDays=10) smoker(K=20,r=0.02,shift=11,nDays=10) smoker(K=20,r=0.01,shift=10.5,nDays=10)
実装例 その2 (2015年度)
addiction<-function(s,r,x){ # initialization k=100 xs=seq(0,k,0.1); plot(xs,k/(1+exp(k*r*(s-xs))),type='l',lwd=2,xlim=c(-3,k+5),ylim=c(-3,k+5), xlab="consumption @ time T",ylab="consumptino @ time T+1") lines(c(0,k),c(0,k),col='blue',lty=3,lwd=3) # main for (i_day in 1:14) { u=k/(1+exp(k*r*(s-x))) if (i_day ==1){ lines(c(x,x),c(5,u),col='red',lty=2,lwd=2) } else { lines(c(x,x),c(x,u),col='red',lty=2,lwd=2)} lines(c(x,u),c(u,u),col='red',lty=2,lwd=2) lines(x,u,type='p',cex=2.5,pch=20,col="magenta") text(x,0,'|',cex=1.5,pch=18,col="darkgreen",font=2) x=u } } > addiction(20,0.0005,30)