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

```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)
```

# 数理社会学：軍拡モデル

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次元のセルオートマトン

#　ルール：
(1) cell[i,j]が０であり、かつ、周りの８つのセルに「1」の状態であるセルが2つある時に、cell[i,j]の次世代での状態は「1」となる
(2) cell[i,j]が1であり、かつ、周りの８つのセルに「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')
}
```

# セルオートマトン（１次元）

セルオートマトン

```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;
}
```

```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章：なぜタバコがやめられないか
モデルの説明：
べッカーの依存症モデル

```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)
```

```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
}
}