Monthly Archives: June 2015

Posted on by

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

mathsoc_ch11
社会を＜モデル＞でみる：数理社会学への招待
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)
Posted on by

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}
 }
Posted on by

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

人口の推移モデル
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)

population01

周期性の検証

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

population02