Posted on by

セルオートマトン（２次元）

セルオートマトン
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')
}
Posted on by

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

セルオートマトン

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

最適化法の比較

認知情報科学基礎実習 課題03
提出期限：2013.01.31 23:59

勾配法、単純確率的最適化法、焼き鈍し法を比較してみましょう。
目的関数：x^4-16*x^2-5*x+y^4-16*y^2-5*y
初期値：一様分布（-0.5~0.5）に従う乱数
繰り返し数：1000回
各最適化法の結果の平均を可視化して（下の図のようなもの）傾向を考察してみてください。
各最適化は下のものを参考にしてください。
各最適化法のパラメターを変えてみてましょう。

3optims

 #Objective Function 
funZ<-function(x,y) {x^4-16*x^2-5*x+y^4-16*y^2-5*y}
xmin=-5;xmax=5;n_len=100;
x<-seq(xmin, xmax,length.out=n_len);y<-x;z<-outer(x,y,funZ)

 #gradient descent 
minGD<-function(initXY=c(0,0), maxItr=1000, stepSize=0.01){
  x=initXY[1];y=initXY[2];z=funZ(x,y); 
  opHist=matrix(0,nrow=maxItr,ncol=3)
  opHist[1,]=c(x,y,z)
  for (i_loop in 2:maxItr) {
    x=x-stepSize*(4*x^3-32*x-5);
    y=y-stepSize*(4*y^3-32*x-5);
    z=funZ(x,y);
    opHist[counter,]=c(x,y,z);
  }
  return(opHist[1:counter,])
}

res<-minGD(c(0.5,0.5),1000,0.01)
contour(x,y,z,nlevels=30,drawlabel=F)
lines(res[,1:2],col='cyan',type='o',pch=19)

#Simple Stochastic Optimization
minSTOCH<-function(initXY=c(0,0),maxItr=1000,width=0.1) {
  x=initXY[1];y=initXY[2];z=funZ(x,y); 
  opHist=matrix(0,nrow=maxItr,ncol=3)
  opHist[1,]=c(x,y,z)
  for (i_loop in 2:maxItr) {
    xTemp=x+rnorm(1,mean=0,sd=width);
    yTemp=y+rnorm(1,mean=0,sd=width);
    zTemp=funZ(xTemp,yTemp);
    if (z > zTemp) {
  	x=xTemp;y=yTemp;z=zTemp;
    }
    opHist[i_loop,]=c(x,y,z);
  }
  return(opHist)
}

opHist<-minSTOCH(c(0,0),1000,0.1)
contour(x,y,z,nlevels=30,drawlabel=F)
lines(opHist[,1:2],type='o',lwd=2,col='green',pch=20)

#Simulated Annealing
minSA<-function(initXY=c(0,0),maxItr=1000,C=1,eta=0.99,width=10) {
  x=initXY[1];y=initXY[2];z=funZ(x,y); 
  opHist=matrix(0,nrow=maxItr,ncol=3)
  opHist[1,]=c(x,y,z)
  for (i_loop in 2:maxItr) {
    xTemp=x+rnorm(1,mean=0,sd=width)
    yTemp=y+rnorm(1,mean=0,sd=width)
    zTemp=funZ(xTemp, yTemp);
    delta=zTemp-z;
    prob=1/(1+exp(delta/(C*width)))
    if (runif(1) < prob) {
      x=xTemp;y=yTemp;z=zTemp;
     } 
    opHist[i_loop,]=c(x,y,z);
    width=width*eta
  }
  return(opHist)
}

res<-minSA(c(0,0),1000,1,0.99,10)
contour(x,y,z,nlevels=30,drawlabel=F)
lines(res[,1:2],type='o',lwd=2,col='green',pch=20)
Posted on by

数理社会学: なぜ協力するのか

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

数理社会学・なぜタバコがやめられないか

FIGmathsoc_ch1

社会を＜モデル＞でみる：数理社会学への招待
01章：なぜタバコがやめられないか
モデルの説明：
べッカーの依存症モデル
効用関数
仮定１：今日の消費量x(t)は、昨日の消費量x(t-1)に依存する
仮定２：時間が経てば魅力が減る(0<=C<=1) 仮定３：満足度はu*(x(t),C*x(t-1))で表し、依存状態はもっとも高い効用uを得るとする 仮定４：はじめはあまり依存せず、途中にはまりだし、最後には慢性的にはまる つまり、昨日の消費が増えるにしたがって最大の満足を得るには今日もたくさん消費する必要がある。しかし、消費が増えると魅力が減りそれほど増えなくなる。 ＃実装例（ベッカーのモデルを忠実には再現できていない） 今日の消費量を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)

addiction

Posted on by

認知情報解析　Jo Project beta 2

// openCVを使った顔検出
#include "opencv2/objdetect/objdetect.hpp"
#include "opencv/highgui.h"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv/cv.h"
#include "iostream"
#include "stdlib.h"
#include "fstream"

using namespace cv;
int main(int argc, char **argv) {
  IplImage *src_img=0, *src_gray=0;
  if (argc < 2 || (src_img = cvLoadImage (argv[1], CV_LOAD_IMAGE_COLOR)) == 0)
    return -1;
  src_gray = cvCreateImage (cvGetSize (src_img), IPL_DEPTH_8U, 1);
    
  CvMemStorage *storage = 0;
  storage = cvCreateMemStorage (0);
  cvClearMemStorage (storage);
  cvCvtColor (src_img, src_gray, CV_BGR2GRAY);
  cvEqualizeHist (src_gray, src_gray);
    
  CvMemStorage* FaceStorage = cvCreateMemStorage(0);
  CvHaarClassifierCascade* FaceCascade = (CvHaarClassifierCascade*)
    cvLoad("/opt/local/share/OpenCV/haarcascades/haarcascade_frontalface_alt_tree.xml");
    
  std::ofstream out_file("DetectedFace.out");
  out_file<<"xMin,"<<"yMin,"<<"xMax,"<<"yMax"<<'\n';

  cvClearMemStorage( FaceStorage );
  CvSeq* FaceObjects = cvHaarDetectObjects(src_img,FaceCascade,FaceStorage,1.1,4,0,
    cvSize( 30, 30 ));
  CvRect* rF;
  // detect at most 10 faces
  int n=std::min(10,(FaceObjects ? FaceObjects->total : 0 ));
  for (int i=0;i<n;++i){
    rF = ( CvRect* )cvGetSeqElem( FaceObjects, i );
    out_file<<rF->x<<","<< rF->y<<","<<rF->x+rF->width<<","<<rF->y+rF->height<<'\n';
  }
  cvReleaseImage(&src_img);
  cvReleaseImage(&src_gray); 
  return 0;
}

# 上のopenCVのプログラムをRで使用する関数
JoFaceDetection<-function(file_name) {
  res<-system(paste("./JoFaceDetect ",file_name))
  if (res!=0) {
    warning("incorrect file name");
    return(-99)   	
  } else {
    xys=read.csv("DetectedFace.out");
    if (nrow(xys)!=0) {
      print("successfully detected face(s)");
      return(xys);
    } else {
      print("no face was detected");
      return(-99);
    }
  }
}

実行例

loc<-JoFaceDetection("~/Downloads/IMG_2095.jpg")
[1] "successfully detected face(s)"
> loc
  xMin yMin xMax yMax
1  307  114  521  328
Posted on by

認知情報解析B face detection

徐くんプロジェクト
OpenCVを使った顔と目の検出。

#include "opencv2/objdetect/objdetect.hpp"
#include "opencv/highgui.h"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>

int main(int argc, char** argv){
  IplImage* img=cvLoadImage(argv[1]);
  CvMemStorage* FaceStorage = cvCreateMemStorage(0);
  CvMemStorage* EyeStorage = cvCreateMemStorage(0);
  CvHaarClassifierCascade* FaceCascade = (CvHaarClassifierCascade*)  
   cvLoad("/opt/local/share/OpenCV/haarcascades/haarcascade_frontalface_alt_tree.xml");
  CvHaarClassifierCascade* EyeCascade =(CvHaarClassifierCascade*) 
   cvLoad("/opt/local/share/OpenCV/haarcascades/haarcascade_eye_tree_eyeglasses.xml");
  double scale = 1.2;
  
  // Detect Faces
  cvClearMemStorage(FaceStorage);
  CvSeq* FaceObjects = cvHaarDetectObjects(img,FaceCascade,FaceStorage,scale,1.1,0,cvSize(1,1));
  CvRect* rF;
  // Loop through objects and draw boxes
  for (int i=0;i<(FaceObjects ? FaceObjects->total:0);i++ ){
    rF = ( CvRect* )cvGetSeqElem(FaceObjects,i);
    cvRectangle(img,cvPoint(rF->x,rF->y),cvPoint(rF->x+rF->width,rF->y+rF->height ),
      CV_RGB(255,0,0),3,4,0);
  }
  // Detect Eyes
  cvClearMemStorage(EyeStorage);
  CvSeq* EyeObjects = cvHaarDetectObjects(img,EyeCascade,EyeStorage,scale,1.1,0,cvSize(1,1));
  CvRect* rE;
  // Loop through objects and draw boxes
  for (int i=0;i<(EyeObjects ? EyeObjects->total:0);i++){
    rE = ( CvRect* )cvGetSeqElem(EyeObjects,i);
    cvRectangle(img,cvPoint(rE->x,rE->y),cvPoint(rE->x+rE->width,rE->y+rE->height ),
      CV_RGB(0,0,255),3,4,0);
  }
  // Displaying image
  cvNamedWindow("Face&Eye Detection");
  cvShowImage("Face&Eye Detection",img);
  cvWaitKey();
  cvReleaseImage(&img);    
  return 0;
}

opencvFD01input.jpg opencvFD01output.jpg

パラメターを変えないと、うまくは検出されないこともある。口と鼻の検出の精度は低いみたい。

Posted on by

single layer perceptron

勾配法と確率的最適化による学習例：

# ２変数からなる16種類のロジック
x=matrix(c(1,1,1,1,1,0,1,0,1,1,0,0),ncol=4);
y=matrix(c(0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,0,1,1,0,1,0,1,1,0,0,1,0,
  1,1,0,1,0,1,0,1,1,0,0,1,1,1,0,1,1,0,1,1,0,1,1,0,1,1,1,1,1,1,1),byrow=T,ncol=4)

> x　# input matrix
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1   # X0
[2,]    1    1    0    0   # X1
[3,]    1    0    1    0   # X2
> y
      [,1] [,2] [,3] [,4]
 [1,]    0    0    0    0  # none
 [2,]    0    0    0    1  # notX1 & notX2
 [3,]    0    0    1    0  # notX1 & X2
 [4,]    0    1    0    0  # X1 & notX2
 [5,]    1    0    0    0  # X1 & X2
 [6,]    0    0    1    1  # notX1
 [7,]    0    1    0    1  # notX2
 [8,]    1    0    0    1  # {X1 & X2} or {notX1 & notX2}
 [9,]    0    1    1    0  # {X1 & notX2} or {notX1 & X2}
[10,]    1    0    1    0  # X2
[11,]    1    1    0    0  # X1
[12,]    1    1    1    0  # X1 or X2
[13,]    1    1    0    1  # X1 or notX2
[14,]    1    0    1    1  # notX1 or X2
[15,]    0    1    1    1  # notX1 or notX2
[16,]    1    1    1    1  # all

# 勾配法
stepSize=0.1;gamma=5;logicID=2;
w=matrix(rnorm(3),nrow=1);
for (i_tr in 1:2000) {
  for (i_inst in 1:4){
    a=w%*%x[,i_inst,drop=F]
    s=1/(1+exp(-gamma*a))
    w=w+stepSize*(y[logicID,i_inst]-s)*s*(1-s)*gamma*x[,i_inst]
  }
}
for (i_inst in 1:4){
  a=w%*%x[,i_inst,drop=F]
  print(1/(1+exp(-gamma*a)))
}

# 確率的最適化法（焼きなまし法）
gamma=2;logicID=2;
w=matrix(rnorm(3),nrow=1);
err=4;width=1;sTemp=matrix(0,nrow=1,ncol=3);C=0.25
for (i_tr in 1:5000) {
  wTemp=w+rnorm(3,mean=0,sd=width)
  for (i_inst in 1:4) {
    a=wTemp%*%x[,i_inst,drop=F]
    sTemp[i_inst]=1/(1+exp(-gamma*a))
  }
  errTemp=sum((sTemp-y[logicID,])^2)
  delta=errTemp-err;
  prob=1/(1+exp(delta/(C*width)))
  if (runif(1) <  prob) {
    w=wTemp
    err=errTemp
  }
  width=width*0.99
}
for (i_inst in 1:4) {
  a=w%*%x[,i_inst,drop=F]
  print(1/(1+exp(-gamma*a)))
}