# 数理社会学・第10章：喧嘩しても分かれない二人

10章：なぜ好きなのにケンカするのか
x:xさんの行為、正の値：２人にとって良い行為;負の値:２人もしくはyさんにとって悪い行為
y:yさんの行為、正の値：２人にとって良い行為;負の値:２人もしくはxさんにとって悪い行為
x(t+1)=x(t)+a*x(t)+b*y(t)
y(t+1)=y(t)+c*x(t)+d*y(t)

```# なぜ好きなのにけんかをするのか
fighting_couple<-function(a,b,c,d) {
timeSep=0.05;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]+b*y[i_gen-1])*timeSep
y[i_gen]=y[i_gen-1]+(c*x[i_gen-1]+d*y[i_gen-1])*timeSep
}
if (i_ini==1) {
plot(x,y,xlim=c(-50,50),ylim=c(-50,50),col=4,type='l',lwd=2,
xlab="X's Action",ylab="Y's Action")
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)
}
}
}

par(mfrow=c(2,3))
fighting_couple(-1,0.0,0.5,-1)
fighting_couple(-1,2,-1,-1)
fighting_couple(0,-1,1,0)
fighting_couple(1,-2,2,0)
fighting_couple(1,0,0.5,1)
fighting_couple(1,-4,-4,0)
```
```fun1<-function(a,b,c,d,N,initX,initY) {
dt=0.05;
x=matrix(0,nrow=N,ncol=1);y=matrix(0,nrow=N,ncol=1);
x[1]=initX;y[1]=initY;
for (i_rep in 2:N) {
x[i_rep]=x[i_rep-1]+(a*x[i_rep-1]+b*y[i_rep-1])*dt
y[i_rep]=y[i_rep-1]+(c*x[i_rep-1]+d*y[i_rep-1])*dt
}
return(data.frame(x,y))
}
fun2<-function(a,b,c,d,N,initX, initY) {
res<-fun1(a,b,c,d,N,initX[1],initY[1]);
nXYs=length(initX)
plot(res\$x,res\$y,xlim=c(-50,100),ylim=c(-50,100),type='l')
for (i_loop in 2:nXYs){
res<-fun1(a,b,c,d,N,initX[i_loop],initY[i_loop]);
lines(res\$x,res\$y,type='l')
}
}
# 実行例
initX=rep(seq(10,50,10),5)
initY=sort(rep(seq(10,50,10),5))
fun2(0,-1,1,0,500,initX,initY)
```

# CLT, LLN, GCD

```# CLT
nSample=10;nRep=10^5;
x=matrix(runif(nSample*nRep),nrow=nSample);
x.means<-colMeans(x)
hist(x.means,50,main='Distribution of Means of Uniform Distribution',
xlab='Means', probability=T)
x.means.dens<-density(x.means)
lines(x.means.dens,lwd=3,lty=1,col='red')
x2=seq(0,1,0.001);CLT=dnorm(x2,mean=0.5,sd=(sqrt(1/12))/(sqrt(nSample)))
lines(x2,CLT,lwd=3,lty=3,col='cyan')
legend("topright",c("Density of Actual Means","Normal Distribution"),
lty=c(1,3), col=c('red','cyan'),lwd=3)

# LAW of LARGE NUMBER
dat<-sample(1:6,1000,replace=T)
prob6<-cumsum(dat==6)/(1:1000)
par(mfrow=c(2,1))
plot(dat,type='p',col="red",main="results of 1000 rolls of a die",xlab="rolls",
ylab="outcome", pch=20)
plot(prob6,type='l',main="Cumulative probability that rolls of a die come out SIX",
xlab="rolls", ylab="probability",lwd=3,ylim=c(0.0,0.5))
abline(h=1/6,col='red',lwd=2,lty=2)

# WHILE
r=-99;v1=1633;v2=355
while (r!=0){
r=v1%%v2
print(paste('v1 =',v1,', v2 = ',v2,',remainder = ',r))
v1=v2
v2=r
}

# REPEAT
v1=1633;v2=355;
repeat {
r=v1%%v2
print(paste("v1 =",v1,",v2 = ",v2,",remainder = ",r))
v1=v2;v2=r
if (r==0){ break}
}
```

# 個体群生態学 人口の推移モデル

N_t+1=r*N_t*(1-N_t)

```FUNpop_t<-function(r=2.5,n0=0.5,gen=100) {
nt=n0;x=seq(0, 1, 0.01);y=r*x*(1-x)
maxY=max(y,1);i_gen=1
plot(x,y,type='l',xlim=c(-0.1, 1.1),ylim=c(0,maxY),xlab="Prop. @ time t",
ylab="Prop @ time t+1");
lines(x,x,type='l')
while ((i_gen < gen) & (nt>0)) {
nt1=r*nt*(1-nt);
lines(c(nt,nt),c(nt, nt1),type='l',col='red')
lines(c(nt,nt1),c(nt1,nt1),type='l',col='red')
nt=nt1;
i_gen=i_gen+1;
}
}
# example
par(mfrow=c(2,2))
FUNpop_t(r=2.7,n0=0.15,gen=100)
FUNpop_t(r=3.1,n0=0.15,gen=100)
FUNpop_t(r=3.8,n0=0.15,gen=100)
FUNpop_t(r=4.2,n0=0.15,gen=100)
```

```FUNpop_t2<-function(gen=100,digit=5,n0=0.5){
r=seq(2.5,4,0.001)
res=matrix(0,ncol=2,nrow=1)
r.length=length(r)
for (i_r in 1:r.length) {
nt=matrix(0,nrow=gen,ncol=1)
nt[1]=n0;
for (i_loop in 2:gen)
{nt[i_loop]=r[i_r]*nt[i_loop-1]*(1-nt[i_loop-1])}
bunnki<-unique(round(nt[duplicated(round(nt,digits=digit))],digits=digit))
bunnki.length<-length(bunnki);
res=rbind(res,cbind(rep(r[i_r],bunnki.length),bunnki))
}
return(res[-1, ])
}
# example
res<-FUNpop_t2(gen=10000, digit=5)
plot(res[,1],res[,2],pch=19,cex=0.1,col="red",xlab="r values",ylab="cyclic points")
```