Disclaimer

このcourselogにあるコードは、主に学部生・博士課程前期向けの講義・演習・実習で例として提示しているもので、原則直感的に分かりやすいように書いているつもりです。例によってはとても非効率なものもあるでしょうし、「やっつけ」で書いているものあります。また、普段はMATLABを使用していますので、変な癖がでているかもしれません。これらの例を使用・参考にする場合はそれを踏まえてたうえで使用・参考にして下さい。
卒業論文に関する資料:[2015PDF] [word template] [latex template] [表紙] [レポートの書き方] [引用文献など]
A Brief Guide to R (Rのコマンドの手引きおよび例)はこちらをどうぞ

認知情報解析演習 W10

n.var = 10
n.child = 10
child = matrix(sample(0:1,n.var*n.child, replace =T),nrow =n.child)

GA_mutate <- function(child, p){
  mut.mat = matrix(sample(0:1,length(child),replace = T, 
                          prob = c(1-p, p)),nrow=nrow(child))
  return(abs(child-mut.mat))
}

GA_crossover<-function(parent){
  n.parent = nrow(parent)
  rand.order = sample(1:n.parent)
  p1.idx= rand.order[1:(n.parent/2)]
  p2.idx= rand.order[(n.parent/2+1):n.parent]
  co.idx = sample(1:(length(parent)/2),length(parent)/4)
  c1 = parent[p1.idx,]
  c1.copy = c1
  c2 = parent[p2.idx,]
  c1[co.idx] = c2[co.idx]
  c2[co.idx] = c1.copy[co.idx]
  return(child = rbind(c1,c2))
}


GA_survive<-function(parent, child, fitP, fitC){
  nPop=nrow(parent);
  fitT=c(fitP,fitC);
  fitMax=sort(fitT,index.return=TRUE,decreasing = TRUE)
  tempX=rbind(parent,child);
  fittest=tempX[fitMax$ix[1:nPop],]
  return(fittest)
}

GA<-function(parent, nGen,p,X,Y){
  optHist=matrix(0,nrow=nGen,ncol=1)
  fittest.hist = matrix(0,nrow=nGen,ncol = ncol(parent))
  nPop = nrow(parent)
  nVar = ncol(parent)
  fitP = rep(0,nPop)
  fitC = fitP
  for (i.pop in 1:nPop){
    dat.temp = X[,which(parent[i.pop,]==1)]
    fitP[i.pop] = summary(lm(Y~dat.temp))$adj.r.squared
  }
  optHist[1]=max(c(fitP))
  fittest.hist[1,] = parent[which.max(fitP),]
  for (i_gen in 2:nGen) {
    child<-GA_crossover(parent)
    child<-GA_mutate(child,p)
    # assuming same numbers of set of genes
    for (i.pop in 1:nPop){
      dat.temp = X[,which(parent[i.pop,]==1)]
      fitP[i.pop] = summary(lm(Y~dat.temp))$adj.r.squared
      if (sum(child[i.pop,])==0){
        child[i.pop,sample(1:ncol(child),sample(1:nVar,1))]=1
      }
      dat.temp = X[,which(child[i.pop,]==1)]
      fitC[i.pop] = summary(lm(Y~dat.temp))$adj.r.squared
    }
    parent<-GA_survive(parent, child, fitP, fitC)
    optHist[i_gen]=max(c(fitP,fitC))
    fittest.hist[i_gen, ] = parent[1,]
  }
  return(list(fittest = parent, optHist = optHist, fittest.hist = fittest.hist))
}

set.seed(19)
nVar = 10
nData = 50
X=matrix(rnorm(nVar*nData,mean=10,sd=5),ncol=nVar);
y=X%*%c(-2,0,2,0,2,0,-3,0,-2,0)+rnorm(nData,mean=0,sd=8);
summary(lm(y~X[,seq(1,9,2)]))
summary(lm(y~X))
parent = matrix(sample(0:1,300,replace = T),nrow=30)
res<-GA(parent,100,0.1,X,y)
head(res$fittest)

より複雑なモデルへ

C = combn(3,2)
paste("X",C,sep="")
temp.lab =toString(paste("X",C[ , 1], sep = ""))
gsub(",", ":", temp.lab) 

n.comb = ncol(C)
var.labels=c()
for (i.comb in 1: n.comb){
  temp.label = toString(paste("X",C[ , i.comb], sep = ""))
  var.labels = c(var.labels, gsub(",", ":", temp.label) )
}


mk.labels <- function(C){
  var.labels = c()
  n.comb = ncol(C)
  for (i.comb in 1: n.comb){
    temp.label = toString(paste("X",C[ , i.comb], sep = ""))
    var.labels = c(var.labels, gsub(",", ":", temp.label) )
  }
  return(var.labels)
}
                                    
n.var = 3
max.interaction = 3
# single variable
var.labels = paste("X", 1:n.var, sep = "")


# interaction terms
for (i.interaction in 2:max.interaction ){
  combination = combn(n.var, i.interaction)
  var.labels = c(var.labels, mk.labels(combination))
}

model.def = paste("Y ~", gsub(",", "+",   toString(var.labels[c(1,1,1,1,1,0,0) == 1]))) 

set.seed(20)
nVar = 3
nData = 50
X=matrix(rnorm(nVar*nData,mean=10,sd=5),ncol=nVar);
y=X%*%c(1,1,-2)+rnorm(nData,mean=0,sd=3);
dat = data.frame(Y=y,X1=X[,1],X2=X[,2],X3=X[,3])

# example
model.def1 = paste("Y ~", gsub(",", "+",  toString(var.labels[c(1,1,1,0,0,0,0) == 1]))) 
model.lm = lm(model.def1, dat)

Evolutionary Strategy

n.parent = 10
n.theta = 10
parent = list(theta = matrix(rnorm(n.parent*n.theta),nrow=n.parent),
              sigma = matrix(runif(n.parent*n.theta)+1,nrow = n.parent))

ES_crossover<-function(parent){
  n.parent = nrow(parent$theta)
  rand.order = sample(1:n.parent)
  p1.idx= rand.order[1:(n.parent/2)]
  p2.idx= rand.order[(n.parent/2+1):n.parent]
  co.idx = sample(1:(length(parent)/2),length(parent)/4)
  c1 = parent$theta[p1.idx,]
  c1.copy = c1
  c2 = parent$theta[p2.idx,]
  c1[co.idx] = c2[co.idx]
  c2[co.idx] = c1.copy[co.idx]
  sigma = (parent$sigma[p1.idx,]+parent$sigma[p2.idx,])/2
  return(child = list(theta = rbind(c1,c2), sigma = rbind(sigma,sigma)))
}

ES_mutate<-function(child, tau){
  child$sigma = child$sigma*
    exp(tau*matrix(rnorm(length(child$sigma)),nrow=nrow(child$sigma)))
  child$theta = child$theta + 
    matrix(rnorm(length(child$sigma),0,child$sigma),nrow=nrow(child$sigma))
  return(child)
}

set.seed(20); nData = 100
X=matrix(rnorm(9*nData,mean=10,sd=2),ncol=9);X=cbind(rep(1,nData),X)
y=X%*%c(10,2,5,-3,-5,0,0,0,0,0)+rnorm(nData,mean=0,sd=2);
#summary(lm(y~X[,2:10]))

# 評価関数
reg.error<-function(b,X,y){
  yhat<-X%*%b
  return(sum((y-yhat)^2))
}
# 評価の使用例
fitC = apply(child$theta,1,reg.error, X, y)

データ解析基礎論 W10

f=c(24.1,23.9,24.4,24.4,23.5)
m=c(25.6,26.1,25.8,25.9,26)

G.mean=mean(c(f,m))
ss.total=sum((c(f,m)-G.mean)^2)
ss.tr=sum((mean(f)-G.mean)^2)*5+sum((mean(m)-G.mean)^2)*5
ss.error=sum((f-mean(f))^2)+sum((m-mean(m))^2)
df.tr = 2 - 1
df.error = (5-1)*2
ms.tr = ss.tr /df.tr
ms.error = ss.error /df.error

dat<-data.frame(ssize=c(f,m),gender=c(rep("f",5),rep("m",5)))
dat.aov<-aov(ssize ~ gender,data=dat)
summary(dat.aov)

dat<-read.table("http://www.matsuka.info/data_folder/aov01.txt")
dat.aov<-aov(shoesize~club, dat)
summary(dat.aov)

qT<-qtukey(0.95, 3, 67)
HSD<-qT*sqrt((2.243*(1/23+1/24))/2)
means<-tapply(dat$shoesize,dat$club,mean)
abs(outer(means,means,"-"))
abs(outer(means,means,"-"))>HSD

TukeyHSD(dat.aov)

dat<-read.csv("http://www.matsuka.info/data_folder/dktb312.csv")
dat$method=factor(dat$method, levels(dat$method)[c(1,3,4,2)])

dat.aov<-aov(result~method, data=dat)
summary(dat.aov)
dat.means<-tapply(dat$result,dat$method,mean)
new.alpha = 1-(1-0.05)^(1/5)
cont=c(-3,1,1,1)
bunshi=sum(cont*dat.means)
bunbo=sqrt(5.29*(sum((cont^2)/8)))
t.value=bunshi/bunbo
2*(1-pt(t.value,28))

iterative prisoners’ dilemma

pmat = np.array([[-25,0],[-20,-5]])

def iterPD(action1, action2, pmat):
# A = defect & defect (-10); B = defect & coop (1)
# C = coop & defect   (9)  ;   D = coop & coop (11)

    result = action1 + action2
    if result == -10:
        po = [pmat[0,0],pmat[0,0]]
    elif result == -9:
        po = [pmat[1,0],pmat[0,1]]
    elif result == 10:
        po = [pmat[0,1],pmat[1,0]]
    else:
        po = [pmat[1,1],pmat[1,1]]
    
    return(po)

def strCoop(mode):
    if mode == "p1":
        return(1)
    else:
        return(10)
    
def strDef(mode):
    if mode == "p1":
        return(0)
    else:
        return(-10)
    
def strTfT(lastPlay, mode):
    if lastPlay == "D":
        if mode == "p1":
            return(0)
        else:
            return(-10)
    else:
        if mode == "p1":
            return(1)
        else:
            return(10)
Nrep = 10
p1Hist = np.zeros([Nrep,2])
p2Hist = np.zeros([Nrep,2])
  
for i_rep in range(Nrep):
    p1 = strDef("p1")
    p2 = strCoop("p2")
    po = iterPD(p1,p2,pmat)
    p1Hist[i_rep,...] = [p1,po[0]]
    p2Hist[i_rep,...] = [p2,po[1]]

lastPlay = "C"
for i_rep in range(Nrep):
    p1 = strCoop("p1")
    p2 = strTfT(lastPlay, "p2")
    if p1 == 0:
        lastPlay = "D"
    else:
        lastPlay = "C"
    po = iterPD(p1,p2,pmat)
    p1Hist[i_rep,...] = [p1,po[0]]
    p2Hist[i_rep,...] = [p2,po[1]]

データ解析基礎論A W09

dat <- read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
dat.lm1 <- lm(sales~price,data = dat)
summary(dat.lm1)

dat.lm2 <- lm(sales~price+material,data = dat)
summary(dat.lm2)

dat.lm01<-lm(sales~price, data=as.data.frame(scale(dat)))
plot(dat.lm01,which=1)
plot(dat.lm01,which=2)

norm.vars=rnorm(300)
qqnorm(norm.vars)
qqline(norm.vars,col='red',lwd=2)

unif.vars=runif(300)
qqnorm(unif.vars)
qqline(unif.vars,col='green',lwd=2)

par(mfrow=c(2,2))
plot(dat.lm01)

dat<-read.csv("http://www.matsuka.info/data_folder/dktb312.csv")
dat$method=factor(dat$method, levels(dat$method)[c(1,3,4,2)])
dat.lm<-lm(result~method,data=dat)
plot(dat.lm)

dat <- read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
dat2 <- rbind(dat,c(0,0,100,250))
tail(dat2)
plot(dat2[1:50,],col=c(rep("black",50),"red"),pch=20,cex=3)
plot(dat2[1:50,],pch=20,cex=3)
dat.lm<-lm(sales~design,data=dat)
summary(dat.lm)
dat2.lm<-lm(sales~design,data=dat2)
summary(dat2.lm)

dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg02.csv")
dat.lm<-lm(sales~., data=dat) 
install.packages("DAAG")
library(DAAG)
vif(dat.lm)

dat<-read.csv("http://www.matsuka.info/data_folder/waa07_2.csv")
poly.lm<-lm(grade~study+study.sq, data= dat)

dat<-read.table("http://www.matsuka.info/data_folder/tdkPATH01.txt",header=TRUE)
plot(dat,pch=20,cex=2)
dat.lm1<-lm(absence~interest,dat)
dat.lm2<-lm(study~interest+knowledge,dat)
dat.lm3<-lm(grade~knowledge+study+absence,dat)

dat.STD<-as.data.frame(scale(dat))
summary(lm(absence~interest,dat.STD))
summary(lm(study~interest+knowledge,dat.STD))
summary(lm(grade~study+absence+knowledge,dat.STD))

院・認知情報解析

タカ・ハトゲーム

food = 6.0
cost = 8.0
Ngen = 1000
n_hawk = 10.0
n_dove = 100.0
dt = 0.01
p_hawk = np.zeros(Ngen)
p_dove = np.zeros(Ngen)
p_hawk[0] = n_hawk/(n_hawk + n_dove)
p_dove[0] = n_dove/(n_hawk + n_dove)
p = np.array([p_hawk[0], p_dove[0]])
pay_mat = np.array([[(food-cost)/2.0,food],[0.0, food/2.0]])
for i_gen in range((Ngen-1)):
    mean_W = np.sum(np.multiply(pay_mat, np.outer(p,p)))
    hawk_W = np.sum(np.inner(pay_mat[0,...], p))
    dove_W = np.sum(np.inner(pay_mat[1,...], p))
    p_hawk[i_gen + 1] = p_hawk[i_gen] + p_hawk[i_gen]*(hawk_W - mean_W)/mean_W*dt
    p_dove[i_gen + 1] = p_dove[i_gen] + p_dove[i_gen]*(dove_W - mean_W)/mean_W*dt
    p = [p_hawk[i_gen + 1], p_dove[i_gen + 1]]

ケーキの分配ゲーム

dt = 0.02
n_cake = 10 
pCake = np.random.uniform(0,1, n_cake + 1)
pCake = pCake/sum(pCake)
cakeMat = np.zeros((n_cake + 1, n_cake + 1))
for i_cake in range((n_cake)):
    cakeMat[i_cake + 1, 0:(n_cake - i_cake)] = i_cake + 1

p_hist = np.zeros((Ngen,n_cake + 1))
for i_gen in range(Ngen):
    Ws = np.sum(np.multiply(pCake,cakeMat),axis=1)
    Wmean =  np.sum(np.multiply(cakeMat, np.outer(pCake,pCake)))
    pCake = pCake + pCake*(Ws - Wmean)/Wmean*dt
    p_hist[i_gen,...] = pCake

fig = plt.figure(figsize = (9,6))
ax = fig.add_subplot(111, ylim =(0,1.05),)
for i_plot in range(n_cake):
    ax.plot(p_hist[...,i_plot], label="Demand = " + str(i_plot))

ax.set_ylabel('Proportion')
ax.set_xlabel('Time')
ax.legend()

認知情報解析演習 08

  # genetic algorithm - identify the best linear model
  GA_recomb<-function(G) {
    nPop=nrow(G);
    nVar=ncol(G);
    child = G;
    G.permuted = G[sample(1:nPop),]
    recomb.idx=which(matrix(sample(0:1,nPop*nVar,replace=T),nrow=nPop)==1)
    child[recomb.idx] = G.permuted[recomb.idx]
    return(child)
  }
  
  GA_mutate = function(child, p){
    n.pop = nrow(child)
    n.var = ncol(child)
    mut.mat= matrix((runif(n.pop*n.var) < p), nrow = n.pop)
    child = abs(child-mut.mat)
    return(child)
  }
  
  GA_survive<-function(G, child, fitG, fitC){
    nPop=nrow(G);
    fitT=c(fitG,fitC);
    fitMax=sort(fitT,index.return=TRUE,decreasing = TRUE)
    tempX=rbind(G,child);
    G=tempX[fitMax$ix[1:nPop],]
    return(G)
  }
  
  GA<-function(G, nGen,p,X,Y){
    optHist=matrix(0,nrow=nGen,ncol=1)
    G.hist = matrix(0,nrow=nGen,ncol = ncol(G))
    nPop = nrow(G)
    nVar = ncol(G)
    fitG = rep(0,nPop)
    fitC = fitG
    for (i.pop in 1:nPop){
      dat.temp = X[,which(G[i.pop,]==1)]
      fitG[i.pop] = summary(lm(Y~dat.temp))$adj.r.squared
    }
    optHist[1]=max(c(fitG))
    G.hist[1,] = G[which.max(fitG),]
    for (i_gen in 2:nGen) {
      child<-GA_recomb(G)
      child<-GA_mutate(child,p)
      # assuming same numbers of set of genes
      for (i.pop in 1:nPop){
        dat.temp = X[,which(G[i.pop,]==1)]
        fitG[i.pop] = summary(lm(Y~dat.temp))$adj.r.squared
        if (sum(child[i.pop,])==0){
          child[i.pop,sample(1:ncol(child.all),sample(1:nVar,1))]=1
        }
        dat.temp = X[,which(child[i.pop,]==1)]
        fitC[i.pop] = summary(lm(Y~dat.temp))$adj.r.squared
      }
      G<-GA_survive(G, child, fitG, fitC)
      optHist[i_gen]=max(c(fitG,fitC))
      G.hist[i_gen, ] = G[1,]
    }
    return(list(G = G, optHist = optHist, G.hist = G.hist))
  }
  
  set.seed(19)
  nVar = 10
  nData = 50
  X=matrix(rnorm(nVar*nData,mean=10,sd=5),ncol=nVar);
  y=X%*%c(-2,0,2,0,2,0,-3,0,-2,0)+rnorm(nData,mean=0,sd=8);
  summary(lm(y~X[,seq(1,9,2)]))
  summary(lm(y~X))
  G = matrix(sample(0:1,300,replace = T),nrow=30)
  res<-GA(G,100,0.1,X,y)
  head(res$G)
  
  # GA - challenging task
  C = combn(3,2)
  temp.lab =toString(paste("X",C[ , i.comb], sep = ""))
  gsub(",", "*", temp.lab) 
  
  n.comb = ncol(C)
  var.labels = c()
  for (i.comb in 1: n.comb){
    temp.label = toString(paste("X",C[ , i.comb], sep = ""))
    var.labels = c(var.labels, gsub(",", "*", temp.label) )
  }
  
  mk.labels <- function(C){
    var.labels = c()
    n.comb = ncol(C)
    for (i.comb in 1: n.comb){
      temp.label = toString(paste("X",C[ , i.comb], sep = ""))
      var.labels = c(var.labels, gsub(",", "*", temp.label) )
    }
    return(var.labels)
  }
  var.labels = paste("X", 1:n.var, sep = "")
  
  # interaction terms
  for (i.interaction  in 2:max.interaction ){
    combination = combn(n.var, i.interaction)
    var.labels = c(var.labels, mk.labels(combination))
  }
  model.def = paste("Y ~", gsub(",", "+", toString(var.labels[c(1,1,1,1,1,0,0) == 1]))) 
  X=matrix(rnorm(nVar*nData,mean=10,sd=5),ncol=nVar);
  y=X%*%c(1,1,-2)+rnorm(nData,mean=0,sd=3);
  dat = data.frame(y=y,X1=X[,1],X2=X[,2],X3=X[,3])
  
  model.lm = lm(model.def, dat)
  

データ解析基礎論 W08

dat <- read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
dat.lm <- lm(sales~.,data = dat)
summary(dat.lm)
anova(dat.lm)
sse = anova(dat.lm)[[2]][4]
ssr = sum(anova(dat.lm)[[2]][1:3])
df.err = anova(dat.lm)[[1]][4]
df.model = sum(anova(dat.lm)[[1]][1:3])
F.val = (ssr/3) / (sse/46)
1 - pf(F.val, df.model, df.err)

dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/dktb312.csv")
tapply(dat$result,dat$method,mean)

dat2<-data.frame(result=dat$result,
                 c1=c(rep(0,8),rep(1,8),rep(0,16)),
                 c2=c(rep(0,16),rep(1,8),rep(0,8)),
                 c3=c(rep(0,24),rep(1,8)))
dat2.lm<-lm(result~.,data=dat2)


dat3<-data.frame(result=dat$result, 
                 c1=c(rep(-3,8), rep(1,24)), 
                 c2=c(rep(0,8),rep(-2,8),rep(1,16)),
                 c3=c(rep(0,16),rep(-1,8),rep(1,8)))
dat3.lm<-lm(result~c1+c2+c3,data=dat3
            
            
dat<-read.csv("http://www.matsuka.info/data_folder/dktb321.txt")
plot(dat$result~dat$duration,data=dat[dat$method=="method.X",])
result<-dat$result[dat$method=="method.X"]
CL=c(rep(-3,5),rep(-1,5),rep(1,5),rep(3,5))
CQ=c(rep(-1,5),rep(1,5),rep(1,5),rep(-1,5))
CC=c(rep(-3,5),rep(1,5),rep(-1,5),rep(3,5))

trend.lm<-lm(result~CL+CQ+CC)
summary(trend.lm)

dat.x = dat[dat$method=="method.X",]
trend.lm2 <- lm(result~duration, contrasts=list(duration = "contr.poly"),data=dat.x)

dat.x$result[16:20]=dat.x$result[16:20]-3
plot(tapply(dat.x$result,dat.x$duration,mean),
     pch=20,col="red",lwd=3, type="o",cex=3,ylab="mean",xlab="time")

trend.lm2 <- lm(result~duration, contrasts=list(duration = "contr.poly"),data=dat.x)

subj.gender = c(rep("male",30),rep("female",30))
pic.gender = rep(c(rep("male",15),rep("female",15)),2)
# 15m-m, 15m-f, 15f-m, 15f-f 
eye.fix = c(round(rnorm(15,50,5)),
            round(rnorm(15,70,5)),
            round(rnorm(15,65,5)),
            round(rnorm(15,50,5)))
datE = data.frame(eye.fix = eye.fix, subj.gender= subj.gender,pic.gender=pic.gender)
interaction.plot(datE$subj.gender, datE$pic.gender, datE$eye.fix, 
                 trace.label = "Gender of Stimuli", xlab = "Gender of Participants",
                 ylab = "Number of fixation",lwd = 3,type = "b",pch = c(17,20),cex =3,
                 col = c('skyblue','coral'),legend =T)
datE.lm <- lm(eye.fix~subj.gender*pic.gender,data=datE)
summary(datE.lm)

dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv") 
dat.lm01<-lm(sales~price, data=as.data.frame(scale(dat)))
plot(dat.lm01,which=1)
plot(dat.lm01,which=2)
par(mfrow=c(2,2))
plot(dat.lm01)

dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/dktb312.csv")
levels(dat$method)
dat$method=factor(dat$method, levels(dat$method)[c(1,3,4,2)])
dat.lm<-lm(result~method,data=dat)
par(mfrow=c(2,2))
plot(dat.lm)

dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg02.csv")
plot(dat)
dat.lm<-lm(sales~., data=dat)

院・認知情報解析

x(t+1) = x(t) + \left(-\alpha x(t) y(t)\right)dt
y(t+1) = y(t) + \left(-\alpha x(t) y(t) + \beta y(t)\right)dt
z(t+1) = z(t) + \beta y(t)dt

import numpy as np
import matplotlib.pyplot as plt
nRep = 1000
x = np.zeros(nRep)
y = np.zeros(nRep)
z = np.zeros(nRep)
x[0]=100000
y[0]=10
z[0]=0
dt = 0.01
alpha = 0.00003
beta = 0.3
for i_loop in range((nRep-1)):
    x[i_loop+1] = x[i_loop]-alpha*x[i_loop]*y[i_loop]*dt
    y[i_loop+1] = y[i_loop]+(alpha*x[i_loop]*y[i_loop]-beta*y[i_loop])*dt
    z[i_loop+1] = z[i_loop]+beta*y[i_loop]*dt

# plotting results
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, label="not affected")
ax.plot(y, label="affected")
ax.plot(z, label="no longer affected")
ax.set_ylabel('Number of People')
ax.set_xlabel('Time')
ax.legend()

恋愛モデル
x(t+1) = x(t) + \left( ax(t) + by(t)\right)dt
y(t+1) = y(t) + \left( cx(t) + dy(t)\right)dt

nRep = 15
x0 = np.arange(-5,6,1)
y0 = np.arange(-5,6,1)
dt = 0.02
a = -2
b = 1
c = -1
d = 1
fig = plt.figure(figsize = (10,10))
ax = fig.add_subplot(111, xlim =(-6,6), ylim =(-6,6))
for ix in range(len(x0)):
    for iy in range(len(y0)):
        x = np.zeros(nRep)
        y = np.zeros(nRep)
        x[0] = x0[ix]
        y[0] = y0[iy]
        for i_loop in range((nRep-1)):
            x[i_loop+1] = x[i_loop]+(a*x[i_loop]+b*y[i_loop])*dt
            y[i_loop+1] = y[i_loop]+(c*x[i_loop]+d*y[i_loop])*dt
        ax.plot(x,y,color='b')
        ax.arrow(x[i_loop-1],y[i_loop-1],x[i_loop+1]-x[i_loop],y[i_loop+1]-y[i_loop],width=0.05,color='r')


タバコの依存モデル
U = \frac{K}{1 + \textup{exp}\left(rK(s-x)\right)}

x = np.arange(0,100,1)
r = 0.0007
s = 40
K = 100
u = K/(1+np.exp(r*K*(s-x)))
nRep = 10
x1 = np.zeros(nRep)
x1[0] = 45
fig = plt.figure(figsize = (6,6))
ax = fig.add_subplot(111)
ax.plot(x,u,color='orange')
ax.plot([0,100],[0,100],color='blue')
u_old = 0
for i_loop in range((nRep-1)):
    x1[i_loop+1] = K/(1+np.exp(r*K*(s-x1[i_loop])))
    ax.plot([x1[i_loop],x1[i_loop]],[u_old,x1[i_loop+1]],color='red')
    u_old = x1[i_loop+1]
    ax.plot([x1[i_loop],x1[i_loop+1]],[x1[i_loop+1],x1[i_loop+1]],color='red')   

2018 基礎実習02

append.fun01 <- function(n.rep){
  result = c()
  for (i.loop in 1:n.rep){
    result = c(result,i.loop)
  }
  return(result)
}

append.fun02 <- function(n.rep){
  result = rep(0,n.rep)
  for (i.loop in 1:n.rep){
    result[i.loop] =i.loop
  }
  return(result)
}

system.time(append.fun01(1e5))
system.time(append.fun02(1e5))

流行のモデル スクリプト版

x0 = 100000
y0 = 10
z0 = 0
alpha = 0.00003
beta = 0.3
dt = 0.01
n.Rep = 1000

x = rep(0, n.Rep)
y = rep(0, n.Rep)
z = rep(0, n.Rep)

x[1] = x0
y[1] = y0
z[1] = z0

for (t in 1:(n.Rep-1)){
  x[t+1] = x[t] - alpha*x[t]*y[t]*dt
  y[t+1] = y[t] + (alpha*x[t]*y[t] - beta*y[t])*dt 
  z[t+1] = z[t] + beta*y[t]*dt
}

plot(x,type='l',col='blue',lwd=2)
lines(y,col='red',lwd=2)
lines(z,col='green',lwd=2)

認知情報解析演習 06

Traveling Salesman Problem

TSP.inv <- function(x, n){
  x.length = length(x)
  start.idx = sample(1:x.length,1)
  end.idx = min(x.length, (start.idx + n -1))
  x[start.idx:end.idx]=x[end.idx:start.idx]
  return(x)
}

TSP.switch <- function(x){
  x.length = length(x)
  switch.idx = sample(1:x.length,2)
  x[switch.idx] = x[rev(switch.idx)]
  return(x)
}

TSP.trans <-function(x){
  x.length = length(x)
  trans.idx = sample(1:x.length,2)
  new.x = x[-trans.idx[1]]
  new.x = append(new.x, x[trans.idx[1]],after = (trans.idx[2]-1))
  return(new.x)
}

TSP.init<-function(n_city=10) {
  return(matrix(runif(n_city*2,1,100),nrow=n_city,ncol=2))
}

cities = TSP.init(10)

TSP.totDist<-function(route, cities) {
  route=c(route,route[1]);
  n.cities=nrow(cities);
  totDist=0;
  for (i.cities in 1:n.cities) {
    totDist = totDist + dist(cities[route[i.cities:(i.cities + 1)],])
  }
  return(totDist)
}

TSP.totDist(route,cities)

TSP.demo<-function(n.city=20, maxItr=1000, parm.set = c(1,0.99, 50)) {
  loc=TSP.init(n.city);
  route=1:n.city
  ## param. for simulated annealing - change values if necessary 
  C=parm.set[1];
  eta=parm.set[2];
  TEMP=parm.set[3]; 
  ##
  optDist=TSP.totDist(route,loc)
  optHist=matrix(0,nrow=maxItr,ncol=(length(route)+1))
  optHist[1,]=c(route,optDist)
  for (i_loop in 2:maxItr) {
    rand.op=sample(c('inv','sw','trans'),1,prob=c(0.75,0.125,0.125))
    if (rand.op=='inv') {
      new.route=TSP.inv(route,sample(2:n.city,1))
    } else if (rand.op=='sw') {
      new.route=TSP.switch(route)
      } else { 
        new.route=TSP.trans(route)
    }
    new.dist=TSP.totDist(new.route,loc)
    delta=new.dist-optDist
    prob=1/(1+exp(delta/(C*TEMP)))
    if (runif(1) < prob) {
      route=new.route;optDist=new.dist;
    } 
    optHist[i_loop,]=c(route,optDist);
    TEMP=min(TEMP*eta,0.1)
  }
  par(mfrow=c(1,2)) 
  plot(rbind(loc,loc[1,]),type='o',pch=20,cex=2.5, col='red',
       xlab='location @X',ylab='location @Y',main='Initial route')
  plot(loc[optHist[1000,c(1:n.city,1)],],type='o',pch=20, col='magenta',
       cex=2.5,xlab='location @X',ylab='location @Y',main='Optimized route')
  return(optHist)
}
# regression with SimAnneal.
reg.ERR<-function(b,X,y){
  yhat<-X%*%b
  return(sum((y-yhat)^2))
}
##
SA.reg <- function(maxItr,parm.set=c(C=1,eta=0.995,TEMP=100,sigma=100)){
  set.seed(20)
  nData = 100
  x=matrix(rnorm(9*nData,mean=10,sd=2),ncol=9);X=cbind(rep(1,nData),x)
  y=X%*%c(10,2,5,-3,-5,0,0,0,0,0)+rnorm(nData,mean=0,sd=2);
  lm.sum = summary(lm(y~X[,2:10]))
  C = parm.set[1]
  eta = parm.set[2]
  TEMP= parm.set[3]
  sigma = parm.set[4]
  b = rnorm(10,0,100)
  optERR=reg.ERR(b,X,y)
  optHist=matrix(0,nrow=maxItr,ncol=(10+1))
  optHist[1,]=c(b,optERR)
  for (i_loop in 2:maxItr) {
    new.b = b + rnorm(10,0,sigma*TEMP)
    new.ERR=Fun.reg(new.b, X,y)
    delta=new.ERR-optERR
    prob=1/(1+exp(delta/(C*TEMP)))
    if (runif(1) < prob) {
      b=new.b;
      optERR=new.ERR;
    } 
    optHist[i_loop,]=c(b,optERR);
    TEMP=max(TEMP*eta,0.01)
  }
  print(optHist[maxItr,])
  print(coef(lm.sum)[,1])
  return(list(optHist=optHist, gt = coef(lm.sum)))
}
res=SA.reg(1e5,c(80,0.999,1,1))