Disclaimer

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

データ解析基礎論a 分散分析3

source("http://www.matsuka.info/univ/course_folder/cutil.R")

dat<-read.csv("http://www.matsuka.info/data_folder/dktb316.txt")
dat.aov<-aov(words~method+subj+Error(subj:method),data=dat)
dat.aov2<-aov(words~method+Error(subj+subj:method),data=dat)
dat.aov.BTW <-aov(words~method,data=dat)

summary(dat.aov)
RB.qtukey(dat.aov,dat, 0.05)

dat<-read.csv("http://www.matsuka.info/data_folder/dktb3211.txt")
interaction.plot(dat$duration,  # x軸
                 dat$method,    # まとめる変数    
                 dat$result,    # y軸 
                 pch=c(20,20), 
                 col=c("skyblue","orange"),
                 ylab="score",
                 xlab="Duration",
                 lwd=3, 
                 type='b',
                 cex=2, 
                 trace.label="Method")

dat.aov <- aov(result~method*duration+Error(s+s:duration),dat)
TSME<-SPF.tsme(dat.aov,dat,"result")

dat<-read.csv("http://www.matsuka.info/data_folder/dktb3218.txt")
dat.aov<-aov(result~method+duration+method:duration + Error(s+method:s+duration:s+method:duration:s), data=dat)
TSME<-RBF.tsme(dat.aov, dat, "result")

認知情報解析演習 居住区問題

n.circle = 20; n.sharp = 20; size = 10
loc = sample(1:size^2, n.circle+n.sharp)
type = c(rep(1,n.circle),rep(2,n.sharp))
# circle = 1; sharp = 2
people = cbind(type,loc)
state.mat = matrix(0,size,size)
state.mat[people[,2]]=people[,1]

p.move = #1/(2*n.circle)
p.other = 0.1
c.count = 3

idx = cbind(rep(1:size,size),sort(rep(1:size,size)))

#
term = F
while (term == F){
  term = T
  rand.order = sample(1:nrow(people))
  for (i.p in 1:nrow(people)){
    if (people[rand.order[i.p],1]==1){
      # circle
      if (runif(1) < p.move){
        empty = 1:(size^2)
        empty = empty[-sort(people[,2])]
        people[rand.order[i.p],2] = sample(empty,1)
        state.mat = matrix(0,size,size)
        state.mat[people[,2]]=people[,1]
        term = F
      } 
    } else {
      # sharp
      x.min = max(idx[people[rand.order[i.p],2],1]-1,1)
      x.max = min(idx[people[rand.order[i.p],2],1]+1,size)
      y.min = max(idx[people[rand.order[i.p],2],2]-1,1)
      y.max = min(idx[people[rand.order[i.p],2],2]+1,size)
      circle.in = sum(state.mat[x.min:x.max,y.min:y.max]==1)
      if (circle.in >= c.count){
        empty = 1:(size^2)
        empty = empty[-sort(people[,2])]
        people[rand.order[i.p],2] = sample(empty,1)
        state.mat = matrix(0,size,size)
        state.mat[people[,2]]=people[,1]
        term = F
        #print('moved')
      }
    }
  }
  for (i.p in 1:nrow(people)){
    
    if (people[rand.order[i.p],1]==2){
      x.min = max(idx[people[i.p,2],1]-1,1)
      x.max = min(idx[people[i.p,2],1]+1,size)
      y.min = max(idx[people[i.p,2],2]-1,1)
      y.max = min(idx[people[i.p,2],2]+1,size)
      circle.in = sum(state.mat[x.min:x.max,y.min:y.max]==1)
      print(circle.in)
      if (circle.in >= c.count){
        term = F
        break
      }
    }
  }
}
plot(0,0, type= 'n', xlim = c(0,11),ylim=c(0,11))
lab = c("0","#")
text(idx[people[,2],1],idx[people[,2],2],lab[people[,1]],cex=3)
ab = seq(0.5,10.5,1)
for (i in 1:11){
  abline(h=ab[i],col='red')
  abline(v=ab[i],col='red')
}

認知情報解析演習 石とりゲーム(bugあり)

col1 = matrix(c(rep(0,4),c(1,0,0,0),c(1,1,0,0),c(1,1,1,0),rep(1,4)),nrow=4,byrow=T)
col2 = matrix(c(rep(10,4),c(11,10,10,10),c(11,11,10,10),c(11,11,11,10),rep(11,4)),nrow=4,byrow=T)
col3 = matrix(c(rep(100,4),c(101,100,100,100),c(101,101,100,100),c(101,101,101,100),rep(101,4)),nrow=4,byrow=T)
act.list = list()
state.temp = list()
counter = 0
Q1 = list()
Q2 = list()
for (i.c1 in 1:5){
  if (sum(col1[,i.c1])==0){
    act1 = c()
  } else {
    act1 = seq(1,sum(col1[,i.c1]),1)
  }
  for (i.c2 in 1:5){
    if (sum(col2[,i.c2])==40){
      act2 = c()
    } else {
      act2 = seq(11,sum(col2[,i.c2]==11)*11,11)
    }
    for (i.c3 in 1:5){
      if (sum(col3[,i.c3])==400){
        act3 = c()
      } else {
        act3 = seq(101,sum(col3[,i.c3]==101)*101,101)
      }
      counter = counter + 1
      state.temp[[counter]] = cbind(col1[,i.c1],col2[,i.c2],col3[,i.c3])
      act.list[[counter]] = c(act1,act2,act3)
      Q1[[counter]] = rep(0, length(c(act1,act2,act3)))
      Q2[[counter]] = rep(0, length(c(act1,act2,act3)))
    }
  }
}
rm.stone <- function(act, st.shape){
  if (act == -99){s}
  if (act > 100){
    n.remove = act%%100
    n.stone = length(which(st.shape[,3]==101))
    start = (5-n.stone)
    st.shape[(start:(start+n.remove-1)),3] = 100
  } else {
    if (act > 10){
      n.remove = act%%10
      n.stone = length(which(st.shape[,2]==11))
      start = (5-n.stone)
      st.shape[(start:(start+n.remove-1)),2] = 10
    } else {
      n.remove = act
      n.stone = length(which(st.shape[,1]==1))
      start = (5-n.stone)
      st.shape[(start:(start+n.remove-1)),1] = 0
    }
  }
  return(st.shape)
}

id.state <- function(st.shape, state.temp){
  for (i.st in 1:125){
    if  (all(st.shape == state.temp[[i.st]])){
      state.idx = i.st 
      break
    }
  }
  return(state.idx)
}

ck.act <- function(Q, act.vec, eta){
  if (is.null(act.vec)){
    return(list(act = -99, act.idx = -99))
    break
  }
  if (length(act.vec)==1){
    act = act.vec
  } else {
    p = exp(Q[[state.idx]])/sum(exp(Q[[state.idx]]))
    act = sample(act.vec, 1, prob = p)
  }
  act.idx = which(act.vec==act)
  return(list(act = act, act.idx = act.idx))
}


gamma=1;alpha = 0.1;n.rep=10000
for (i.rep in 1:n.rep){
  # first action
  state.idx = 125; counter = 1
  st.shape = state.temp[[state.idx]]
  res.act = ck.act(Q1,act.list[[state.idx]],eta)
  act = res.act$act;act.idx = res.act$act.idx
  state.old = state.idx
  act.old = act.idx
  
  # 2nd to last
  while (state.idx != 1) {
    counter = counter + 1
    st.shape <- rm.stone(act, st.shape)
    state.idx <- id.state(st.shape, state.temp)
    if (counter%%2==1) {
      res.act = ck.act(Q1,act.list[[state.idx]],eta)
    } else {
      res.act = ck.act(Q2,act.list[[state.idx]],eta)
    }
    act = res.act$act; act.idx = res.act$act.idx
    if (state.idx == 1){
      if (counter%%2==1){rew1 = -1; rew2 = 1;} else {rew1 = 1; rew2 = -1;}
      Q1[[state.old]][act.old]=Q1[[state.old]][act.old]
         +alpha*(rew1-Q1[[state.old]][act.old])
      Q2[[state.old]][act.old]=Q2[[state.old]][act.old]
         +alpha*rew2-Q2[[state.old]][act.old])
    } else {
      rew1 = 0; 
      rew2 =0;
      if (counter%%2==1){
        Q1[[state.old]][act.old]=Q1[[state.old]][act.old]
          +alpha*(rew1+gamma* Q1[[state.idx]][act.idx]-Q1[[state.old]][act.old])
      } else {
        Q2[[state.old]][act.old]=Q2[[state.old]][act.old]
          +alpha*(rew2+gamma* Q2[[state.idx]][act.idx]-Q2[[state.old]][act.old])
      }
    }
    
    state.old = state.idx
    act.old = act.idx
  }
}  

2019年度 データ解析基礎論a 分散分析2

source("http://www.matsuka.info/univ/course_folder/cutil.R")
adj.alpha(5,0.05)

# anova
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)

# multiple comparison 
# do not use these command
# use simpler one  - "cu.bonF1F"
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))


# cu.bonF1F
new.alpha<-adj.alpha(5,0.05) 
cu.bonF1F(dat.aov,dat,c(-3,1,1,1),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,1,0,0),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,0,1,0),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,0,0,1),new.alpha)
cu.bonF1F(dat.aov,dat,c(0,-2,1,1),new.alpha)

# cu.scheffe1F
new.f<-adj.f(4,3,28,0.05)
cu.scheffe1F(dat.aov,dat,c(-3,1,1,1))


# 2 way ANOVA 
dat <- read.csv("http://peach.l.chiba-u.ac.jp/course_folder/datW03R.txt")
interaction.plot(dat$gender,
                 dat$affil,
                 dat$shoesize, 
                 pch=c(20,20), 
                 col=c("skyblue","orange"), 
                 xlab="gender", ylab="shoesize", 
                 lwd=3,type='b',cex=2,
                 trace.label="Affiliation")
dat.aov=aov(shoesize~gender*affil, data=dat)
dat.aov.sum=summary(dat.aov)

# testing simple main effect
# do not use these command
# use simpler one - CRF.tsme
means<-tapply(dat$shoesize, list(dat$gender,dat$affil), mean)
SS_gen_CS<- 5*(means[2,1]^2 + means[1,1]^2 -0.5*sum(means[,1])^2) # SS_gender CS
SS_gen_PS<- 5*(means[2,2]^2 + means[1,2]^2 -0.5*sum(means[,2])^2) # SS_gender PS
dat.aov.sum=summary(dat.aov)   # ANOVA table
MSe=dat.aov.sum[[1]][4,3]      # MSE from ANOVA table or MSe=0.62
dfE=dat.aov.sum[[1]][4,1]      # DF for error or dfE=16
dfG=1                          # DF for gender
F_gen_CS=(SS_gen_CS/dfG)/MSe   # F-value for gender effect given CS
F_gen_PS=(SS_gen_PS/dfG)/MSe   # F-value for gender effect given PS
P_gen_CS=1-pf(F_gen_CS,1,dfE)  # p-value for gender effect given CS
P_gen_PS=1-pf(F_gen_PS,1,dfE)  
SS_affil_F<- 5*(means[1,1]^2+means[1,2]^2-0.5*sum(means[1,])^2) #SS_affil | F
SS_affil_M<- 5*(means[2,1]^2+means[2,2]^2-0.5*sum(means[2,])^2) #SS_affil | M
dfA=1				          # DF for affil
F_affil_F=SS_affil_F/dfA/MSe         # F-value for affiliation effect | F
F_affil_M=SS_affil_M/dfA/MSe         # F-value for affiliation effect | M
P_affil_F=1-pf(F_affil_F,1,dfE)      # p-value for affiliation effect | F
P_affil_M=1-pf(F_affil_M,1,dfE)      # p-value for affiliation effect | M

# testint simple main effect w/ CRF.tsme
tsme = CRF.tsme(dat.aov, dat)
             
# another 2-way ANOVA
dat <-read.csv("http://www.matsuka.info/data_folder/dktb321.txt")
interaction.plot(dat$duration,dat$method,dat$result, pch=c(20,20), 
                 col=c("skyblue","orange"), ylab="score", xlab="Duration",
                 lwd=3,type='b',cex=2,trace.label="Method")
mod1=aov(result~method+duration,data=dat)
mod1.sum=print(summary(mod1))

mod2=aov(result~method*duration,data=dat)
mod2.sum=print(summary(mod2))
CRF.tsme(mod2, dat)


# plotting 
dat<-read.csv("http://www.matsuka.info/univ/course_folder/datW03R.txt",header=T)
means<-tapply(dat$shoesize,list(dat$gender, dat$affil),mean)
Ns<-tapply(dat$shoesize,list(dat$gender, dat$affil),length)
sds<-tapply(dat$shoesize,list(dat$gender, dat$affil),sd)
sems<-sds/sqrt(Ns)
plot(c(0,1),means[,1],type='o',col='skyblue', xlim=c(-0.1,1.1), lwd=2, cex=2, pch=20,
     ylim=c(min(means)*0.975, max(means)*1.025), xlab="gender", ylab="shoesize", xaxt="n")
axis(1,c(0,1),c("female","male"))
lines(c(0,1),means[,2],type='o',col='orange',lwd=2,cex=2,pch=20)
lines(c(0,0),c(means[1,1]-sems[1,1],means[1,1]+sems[1,1]),col="skyblue",lwd=2)
lines(c(0,0),c(means[1,2]-sems[1,2],means[1,2]+sems[1,2]),col="orange",lwd=2)
lines(c(1,1),c(means[2,1]-sems[2,1],means[2,1]+sems[2,1]),col="skyblue",lwd=2)
lines(c(1,1),c(means[2,2]-sems[2,2],means[2,2]+sems[2,2]),col="orange",lwd=2)
legend("topleft",c("CogSci","PsySci"),col=c("skyblue","orange"),lwd=2)

院:認識情報解析

library(rjags)
source("http://peach.l.chiba-u.ac.jp/course_folder/HDI_revised.txt")

dat<-read.csv("http://www.matsuka.info/data_folder/HtWtData110.csv")
library(plot3D)
w = seq(80,360,length.out=100)
h = seq(50, 75, length.out=100)
M <- mesh(w,h)
P.male = 1/(1+exp(-1*(0.018*M$x+0.7*M$y-50)))

scatter3D(dat$weight, dat$height, dat$mal, pch = 19, cex = 2,
          theta = 30, phi = 45, ticktype = "detailed", zlim=c(-0.1,1),ylim=c(50,78),xlim=c(80,360),
          xlab = "weight", ylab = "height", zlab = "P(male)",
          surf = list(x = M$x, y = M$y, z = P.male,facets = NA))


y = dat$male; x = dat$weight; Ntotal = length(y)
dataList = list(y = y, x = x, Ntotal = Ntotal)

model.txt = "
data {
  xm <- mean(x)
  xsd <- sd(x)
  for (i in 1:Ntotal){
    zx[i] = (x[i] - xm)/xsd
  }
}
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dbern(ilogit( zbeta0 + zbeta * zx[i_data]))
  }
  zbeta0 ~ dnorm(0, 1/2^2)
  zbeta ~ dnorm(0, 1/2^2)

  beta <- zbeta / xsd
  beta0 <- zbeta0 - zbeta * xm/xsd
}"
writeLines(model.txt, "model1.txt")
parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model1.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

plot(dat$weight,dat$male,xlim=c(90,280),yaxt="n",ylab="Male / Female",
     xlab="Weight", cex=2.5)
axis(2,at = 0:1,labels=c("Femal","Male"))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(90,280,length.out = 100)
for (i_sample in 1:n2plot) {
  temp.y = 1/(1+exp(-1*(mcmcMat[idx[i_sample],2] + mcmcMat[idx[i_sample],1]*temp.x)))
  lines(temp.x, temp.y, col='orange', lwd=2)
}


x = cbind(dat$weight,dat$height);Nx = ncol(x)
dataList = list(y = y, x = x, Ntotal = Ntotal, Nx = Nx)

model.txt = "
data {
  for (j in 1:Nx){
    xm[j] <- mean(x[,j])
    xsd[j] <- sd(x[,j])
    for (i in 1:Ntotal){
      zx[i,j] = (x[i,j] - xm[j])/xsd[j]
    }
  }
}
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dbern(ilogit( zbeta0 + sum(zbeta[1:Nx] * zx[i_data, 1:Nx ])))
  }
  zbeta0 ~ dnorm(0, 1/2^2)
  for (j in 1:Nx){
    zbeta[j] ~ dnorm(0, 1/2^2)
  }
  beta[1:Nx] <- zbeta[1:Nx] / xsd[1:Nx]
  beta0 <- zbeta0 -sum(zbeta[1:Nx] * xm[1:Nx]/xsd[1:Nx])
}"
writeLines(model.txt, "model.txt")

parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)
par(mfrow=c(1,3))
HDI.plot(mcmcMat[,3],xlabel='intercept')
HDI.plot(mcmcMat[,1],xlabel='weight')
HDI.plot(mcmcMat[,2],xlabel='height')

par(mfrow=c(1,1))
plot(dat$weight,dat$height,xlab="Weight", ylab="Height", type="n")
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
for (i_sample in 1:n2plot) {
  abline(a=-1*mcmcMat[idx[i_sample],3]/mcmcMat[idx[i_sample],2],
         b=-1*mcmcMat[idx[i_sample],1]/mcmcMat[idx[i_sample],2],col="orange")

}
points(dat$weight,dat$height,pch=paste(dat$male), cex=1.5)

# un-even data
x = rnorm(300)
pr = 1/(1+exp(2*x))
y = pr < runif(300)
plot(x,y)

remove.id = sample(which(y == 0),120)

Ntotal = length(y[-remove.id])
dataList = list(y = y[-remove.id], x = x[-remove.id], Ntotal = Ntotal)

model.txt = "
data {
xm <- mean(x)
xsd <- sd(x)
for (i in 1:Ntotal){
zx[i] = (x[i] - xm)/xsd
}
}
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern(ilogit( zbeta0 + zbeta * zx[i_data]))
}
zbeta0 ~ dnorm(0, 1/2^2)
zbeta ~ dnorm(0, 1/2^2)

beta <- zbeta / xsd
beta0 <- zbeta0 - zbeta * xm/xsd
}"
writeLines(model.txt, "model1.txt")
parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model1.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

plot(x[-remove.id],y[-remove.id],xlim=c(-3,3))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(-3,3,length.out = 100)
for (i_sample in 1:n2plot) {
  temp.y = 1/(1+exp(-1*(mcmcMat[idx[i_sample],2] + mcmcMat[idx[i_sample],1]*temp.x)))
  lines(temp.x, temp.y, col='orange', lwd=2)
}


x1 = rnorm(150)
x2 = x1*0.9+rnorm(150,0,0.5)
pr = 1/(1+exp(x1+x2))
y = pr < runif(150)
Ntotal = length(y)
dataList = list(y = y, x = cbind(x1,x2), Ntotal = Ntotal, Nx = 2)


model.txt = "
data {
for (j in 1:Nx){
xm[j] <- mean(x[,j])
xsd[j] <- sd(x[,j])
for (i in 1:Ntotal){
zx[i,j] = (x[i,j] - xm[j])/xsd[j]
}
}
}
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern(ilogit( zbeta0 + sum(zbeta[1:Nx] * zx[i_data, 1:Nx ])))
}
zbeta0 ~ dnorm(0, 1/2^2)
for (j in 1:Nx){
zbeta[j] ~ dnorm(0, 1/2^2)
}
beta[1:Nx] <- zbeta[1:Nx] / xsd[1:Nx]
beta0 <- zbeta0 -sum(zbeta[1:Nx] * xm[1:Nx]/xsd[1:Nx])
}"
writeLines(model.txt, "model.txt")

parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)
plot(x1,x2,xlab="x1", ylab="x2", type="n")
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
for (i_sample in 1:n2plot) {
  abline(a=-1*mcmcMat[idx[i_sample],3]/mcmcMat[idx[i_sample],2],
         b=-1*mcmcMat[idx[i_sample],1]/mcmcMat[idx[i_sample],2],col="orange")

}
points(x1,x2,pch=paste(y), cex=1.5)

# guessing
y = dat$male; x = dat$weight; Ntotal = length(y)
dataList = list(y = y, x = x, Ntotal = Ntotal)

model.txt = "
data {
xm <- mean(x)
xsd <- sd(x)
for (i in 1:Ntotal){
zx[i] = (x[i] - xm)/xsd
}
}
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern(mu[i_data])
mu[i_data] <- (guess*0.5 + (1-guess)*ilogit( zbeta0 + zbeta * zx[i_data]))
}
zbeta0 ~ dnorm(0, 1/2^2)
zbeta ~ dnorm(0, 1/2^2)
guess ~ dbeta(1,9)
beta <- zbeta / xsd
beta0 <- zbeta0 - zbeta * xm/xsd
}"
writeLines(model.txt, "model1.txt")
parameters = c( "beta0" ,  "beta", "guess")
jagsModel = jags.model( "model1.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

plot(x[-remove.id],y[-remove.id],xlim=c(-3,3))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(-3,3,length.out = 100)
for (i_sample in 1:n2plot) {
  temp.y = 1/(1+exp(-1*(mcmcMat[idx[i_sample],2] + mcmcMat[idx[i_sample],1]*temp.x)))
  lines(temp.x, temp.y, col='orange', lwd=2)
}
par(mfrow=c(1,3))
HDI.plot(mcmcMat[,2],xlabel='intercept')
HDI.plot(mcmcMat[,1],xlabel='weight')
HDI.plot(mcmcMat[,3],xlabel='guessing')

par(mfrow=c(1,1))
plot(dat$weight,dat$male,xlim=c(90,280),yaxt="n",ylab="Male / Female",
     xlab="Weight", cex=2.5)
axis(2,at = 0:1,labels=c("Femal","Male"))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(90,280,length.out = 100)
for (i_sample in 1:n2plot) {
  temp.y = mcmcMat[idx[i_sample],3]/2+(1-mcmcMat[idx[i_sample],3])*1/(1+exp(-1*(mcmcMat[idx[i_sample],2] + mcmcMat[idx[i_sample],1]*temp.x)))
  lines(temp.x, temp.y, col='orange', lwd=2)
}

# nomial predictors
model.txt = "
model {
  for ( i.data in 1:Ntotal ) {
    y[ i.data ] ~ dbin(mu[i.data],N[i.data])
    mu[i.data] ~ dbeta(omega[x[i.data]]*(kappa-2)+1,(1-omega[x[i.data]])*(kappa-2)+1)
  }
  for (i.pos in 1:Npos){
    omega[i.pos] <- ilogit(a0+a[i.pos])
    a[i.pos] ~ dnorm(0.0, 1/aSigma^2)
  }
  a0 ~  dnorm(0,1/2^2)
  aSigma ~ dgamma(1.64, 0.32)
  kappa <- kappaMinusTwo +2
  kappaMinusTwo ~ dgamma(0.01,0.01)
  for (i.pos in 1:Npos){
    m[i.pos] <- a0+a[i.pos]
  }
  b0 <- mean(m[1:Npos])
  for (i.pos in 1:Npos){
    b[i.pos] <- m[i.pos] - b0
  }
}"
writeLines(model.txt, "model.txt")

dat<-read.csv("http://www.matsuka.info/data_folder/BattingAverage.csv")
y = dat$Hits
N = dat$AtBats
x = dat$PriPosNumber
Ntotal = length(y)
Npos = length(unique(x))
dataList = list(y = y, x = x, N = N, Ntotal = Ntotal, Npos = Npos)
parameters = c( "b0" ,  "b", "omega")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

par(mfrow=c(3,3))
for (i.pos in 1:9){
  HDI.plot(mcmcMat[,i.pos+10])
}

par(mfrow=c(2,2))
HDI.plot(mcmcMat[,1]-mcmcMat[,2])
HDI.plot(mcmcMat[,2]-mcmcMat[,3])
HDI.plot(mcmcMat[,11]-mcmcMat[,12])
HDI.plot(mcmcMat[,12]-mcmcMat[,13])

# softmax regression
x1 = runif(500, min=-2, max = 2)
x2 = runif(500, min=-2, max = 2)
b0 = c(0,-3,-4,-5)
b1 = c(0,-5,-1,10)
b2 = c(0,-5,10,-1)
l1 = b0[1]+b1[1]*x1+b2[1]*x2
l2 = b0[2]+b1[2]*x1+b2[2]*x2
l3 = b0[3]+b1[3]*x1+b2[3]*x2
l4 = b0[4]+b1[4]*x1+b2[4]*x2
p1 = exp(l1)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p2 = exp(l2)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p3 = exp(l3)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p4 = exp(l4)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
ps = cbind(p1,p2,p3,p4)
y = apply(ps,1,which.max)
plot(x1,x2,pch=y,col=y)

b0 = c(0,-4,-1,-1)
b1 = c(0,-5,1,3)
b2 = c(0,0,-5,3)
l1 = b0[1]+b1[1]*x1+b2[1]*x2
l2 = b0[2]+b1[2]*x1+b2[2]*x2
l3 = b0[3]+b1[3]*x1+b2[3]*x2
l4 = b0[4]+b1[4]*x1+b2[4]*x2
p1 = exp(l1)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p2 = exp(l2)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p3 = exp(l3)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p4 = exp(l4)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
ps = cbind(p1,p2,p3,p4)
y = apply(ps,1,which.max)
plot(x1,x2,pch=y,col=y)

p1 = exp(l1)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p2 = exp(l2)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p3 = exp(l3)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
p4 = exp(l4)/sum(exp(l1)+exp(l2)+exp(l3)+exp(l4))
ps = cbind(p1,p2,p3,p4)
p12 = pmax(p1,p2)
p34 = pmax(p3,p4)
y12vs34 = apply(cbind(p1,p2),1,which.max)
plot(x1,x2,pch=y12vs34,col=y12vs34)
y1vs2 = apply(cbind(p1,p3),1,which.max)
points(x1,x2,pch=y1vs2+2,col=y1vs2+2)
y3vs4 = apply(cbind(p1,p4),1,which.max)
points(x1,x2,pch=y3vs4+6,col=y3vs4+6)



model.txt = "
data {
  for ( j in 1:Nx ) {
    xm[j]  <- mean(x[,j])
    xsd[j] <-   sd(x[,j])
    for ( i in 1:Ntotal ) {
      zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
    }
  }
}
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dcat(mu[1:Nout,i])
    mu[1:Nout,i] <- explambda[1:Nout,i]/sum(explambda[1:Nout,i])
    for (k in 1:Nout){
      explambda[k,i]=exp(zbeta0[k] + sum(zbeta[k,1:Nx] * zx[i, 1:Nx ]))
    }
  }
  zbeta0[1] = 0
  for (j in 1:Nx){
    zbeta[1,j] <- 0
  }
  for (k in 2:Nout){
    zbeta0[k] ~ dnorm(0, 1/2^2)
    for (j in 1:Nx){
      zbeta[k,j]~dnorm(0, 1/2^2)
    }
  }
  for ( k in 1:Nout ) {
    beta[k,1:Nx] <- zbeta[k,1:Nx] / xsd[1:Nx]
    beta0[k] <- zbeta0[k] - sum( zbeta[k,1:Nx] * xm[1:Nx] / xsd[1:Nx] )
  }
}"
writeLines(model.txt, "model.txt")

dat<-read.csv( "http://www.matsuka.info/data_folder/CondLogistRegData1.csv" )
y = dat$Y
x = cbind(dat[,1],dat[,2])
Ntotal = length(y)
Nout = length(unique(y))
dataList = list(y = y, x = x, Nx = 2, Ntotal = Ntotal, Nout = Nout)
parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)
par(mfrow=c(1,3))
HDI.plot(mcmcMat[,7+0],xlab='intercept')
HDI.plot(mcmcMat[,1+0],xlab='b1')
HDI.plot(mcmcMat[,4+0],xlab='b2')

model = "
data {
  for ( j in 1:Nx ) {
    xm[j]  <- mean(x[,j])
    xsd[j] <-   sd(x[,j])
    for ( i in 1:Ntotal ) {
      zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
    }
  }
}
# Specify the model for standardized data:
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dcat( mu[1:Nout,i] )
    mu[1,i] <- phi[1,i]
    mu[2,i] <- phi[2,i] * (1-phi[1,i])
    mu[3,i] <- phi[3,i] * (1-phi[2,i]) * (1-phi[1,i])
    mu[4,i] <- (1-phi[3,i]) * (1-phi[2,i]) * (1-phi[1,i])
    for ( r in 1:(Nout-1) ) {
      phi[r,i] <- ilogit( zbeta0[r] + sum( zbeta[r,1:Nx] * zx[i,1:Nx] ) )
    }
  }
  for ( r in 1:(Nout-1) ) {
    zbeta0[r] ~ dnorm( 0 , 1/20^2 )
    for ( j in 1:Nx ) {
      zbeta[r,j] ~ dnorm( 0 , 1/20^2 )
    }
  }
  for ( r in 1:(Nout-1) ) {
    beta[r,1:Nx] <- zbeta[r,1:Nx] / xsd[1:Nx]
    beta0[r] <- zbeta0[r] - sum( zbeta[r,1:Nx] * xm[1:Nx] / xsd[1:Nx] )
  }
}
"
writeLines( model , "model.txt" )

dat<-read.csv( "http://www.matsuka.info/data_folder/SoftmaxRegData2.csv" )
y = dat$Y
x = cbind(dat[,1],dat[,2])
Ntotal = length(y)
Nout = length(unique(y))
dataList = list(y = y, x = x, Nx = 2, Ntotal = Ntotal, Nout = Nout)
parameters = c( "beta0" ,  "beta")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)
par(mfrow=c(1,1))
plot(x[,1],x[,2],col=y)
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
temp.x = seq(-3,3,length.out = 100)
for (i.cat in 0:2){
  for (i_sample in 1:n2plot) {
    abline(a=-1*mcmcMat[idx[i_sample],7+i.cat]/mcmcMat[idx[i_sample],4+i.cat],
           b=-1*mcmcMat[idx[i_sample],1+i.cat]/mcmcMat[idx[i_sample],4+i.cat],col="orange")
  }
}
model2 = "
data {
  for ( j in 1:Nx ) {
    xm[j]  <- mean(x[,j])
    xsd[j] <-   sd(x[,j])
      for ( i in 1:Ntotal ) {
        zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]
      }
    }
  }
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dcat( mu[1:Nout,i] )
    mu[1,i] <- phi[2,i] * phi[1,i]
    mu[2,i] <- (1-phi[2,i]) * phi[1,i]
    mu[3,i] <- phi[3,i] * (1-phi[1,i])
    mu[4,i] <- (1-phi[3,i]) * (1-phi[1,i])
    for ( r in 1:(Nout-1) ) {
      phi[r,i] <- ilogit( zbeta0[r] + sum( zbeta[r,1:Nx] * zx[i,1:Nx] ) )
    }
  }
  for ( r in 1:(Nout-1) ) {
    zbeta0[r] ~ dnorm( 0 , 1/20^2 )
    for ( j in 1:Nx ) {
      zbeta[r,j] ~ dnorm( 0 , 1/20^2 )
    }
  }
  for ( r in 1:(Nout-1) ) {
    beta[r,1:Nx] <- zbeta[r,1:Nx] / xsd[1:Nx]
    beta0[r] <- zbeta0[r] - sum( zbeta[r,1:Nx] * xm[1:Nx] / xsd[1:Nx] )
  }
}"
writeLines( modelString , con="TEMPmodel.txt" )

2019 データ解析基礎論A 分散分析1

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

# ANOVA w/o aov
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)
ss.tr+ss.error
df.tr = 2 - 1
df.error = (4-1)*2
ms.tr = ss.tr /df.tr
ms.error = ss.error /df.error
1-pf(f.value,1,8)

# ANOVA w/ aov
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)

# TukeyHSD w/o TukeyHSD command
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
TukeyHSD(dat.aov)


dat<-read.csv("http://www.matsuka.info/data_folder/dktb312.csv")
dat.aov<-aov(result~method, data=dat)
summary(dat.aov)
source("http://www.matsuka.info/univ/course_folder/cutil.R")
adj.alpha(5,0.05)

# multiple t-test
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)) > new.alpha

# bonferroni w/ cUTIL
cu.bonF1F(dat.aov,dat,c(-3,1,1,1),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,1,0,0),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,0,1,0),new.alpha)
cu.bonF1F(dat.aov,dat,c(-1,0,0,1),new.alpha)
cu.bonF1F(dat.aov,dat,c(0,-2,1,1),new.alpha)

# scheffe w/ cUTIL
new.f<-adj.f(4,3,28,0.05)
cu.scheffe1F(dat.aov,dat,c(-3,1,1,1))
cu.scheffe1F(dat.aov,dat,c(-1,0,1,0))
cu.scheffe1F(dat.aov,dat,c(-1,0,0,1))

基礎実習 MA02

select.act <- function(state){
  poss.act = which(state=="N")
  if (length(poss.act)==1){
    act = poss.act
  } else {
    act = sample(poss.act, 1)
  }
  return(act)
}

ck.term <- function(state) {
  term = F
  result = "undecided"
  term.idx = matrix(c(1,2,3,4,5,6,7,8,9,1,4,7,
                      2,5,8,3,6,9,1,5,9,3,5,7),ncol=3,byrow=T)
  if (sum(state == "N")==0) {
    term=T
    result = "tie"
  }                      
  for (i.ck in 1:8) {
    O.won = all(state[term.idx[i.ck,]]=="O")
    X.won = all(state[term.idx[i.ck,]]=="X")
    if (O.won ==1 ) {
      term= T
      result = "O won"
      break
    }
    if (X.won ==1){
      term=T
      result = "X won"
      break
    }
  }      
  return(list(term = term, result=result))
}

ck.term2 <- function(state) {
  term = F
  result = "undecided"
  term.idx = matrix(c(1,2,3,4,5,6,7,8,9,1,4,7,
                      2,5,8,3,6,9,1,5,9,3,5,7),ncol=3,byrow=T)
  if (sum(state == "N")==0) {
    term=T
    result = "tie"
  }                      
  O.won = apply(term.idx,1,function(x) all(state[x]=="O"))
  X.won = apply(term.idx,1,function(x) all(state[x]=="X"))
  if (any(O.won) == T ) {
    term= T
    result = "O won"
  } else if (any(X.won) ==1){
    term=T
    result = "X won"
  }
  return(list(term = term, result=result))
}

tictactoe <-function(){ 
  coord = matrix(c(1,3,1,2,1,1,2,3,2,2,2,1,3,3,
                   3,2,3,1), nrow = 9, byrow=T)
  plot(0,0, type="n", xlim= c(0.5,3.5), ylim=c(0.5,3.5))
  abline(h = 1.5); abline(h = 2.5)
  abline(v = 1.5); abline(v = 2.5)
  state = rep("N",9); marker = c("O","X")
  repeat { 
    for (i.player in 1:2) {
      act = select.act(state)
      state[act] = marker[i.player]
      text(coord[act,1],coord[act,2], marker[i.player], cex=4)
      res = ck.term2(state)
      if (res$term == T) { break }
      Sys.sleep(0.5)
    }
    if (res$term == T) { 
      print(paste("result:", res$result))
      break 
    }
  }
}

# 実行例
> tictactoe()
[1] "result: X won"

# 対戦版 - 無作為に動くので非常に弱いです。
# 上のselect.actを変更することでマシになると思います。
play.tictactoe.R <-function(){ 
  coord = matrix(c(1,3,1,2,1,1,2,3,2,2,2,1,3,3,
                   3,2,3,1), nrow = 9, byrow=T)
  plot(coord, pch=paste(1:9), col='gray', cex=3, xlim= c(0.5,3.5), ylim=c(0.5,3.5))
  abline(h = 1.5); abline(h = 2.5)
  abline(v = 1.5); abline(v = 2.5)
  state = rep("N",9); 
  repeat { 
    for (i.player in 1:2){
      if (i.player == 1){
        act = as.numeric(readline(prompt="Enter box ID: "))
      } else {act = select.act(state)}
      state[act] = marker[i.player]
      text(coord[act,1],coord[act,2], marker[i.player], cex=4)
      res = ck.term2(state)
      if (res$term == T) { break }
      Sys.sleep(0.5)
    }
    if (res$term == T) { 
      print(paste("result:", res$result))
      break 
    }
  }
}

# 実行例:
> play.tictactoe.R()
Enter box ID: 5
Enter box ID: 1
Enter box ID: 9
[1] "result: O won"

データ解析基礎論 回帰分析3

dat <- read.csv("http://www.matsuka.info/data_folder/hwsk8-17-6.csv")
dat.lm<-lm(ani~otouto, dat)

plot(hatvalues(dat.lm))
dat.lm$residuals
x = dat$otouto
mean.x = mean(x)

h = 1/length(x)+(x-mean.x)^2/sum((x - mean.x)^2)
hatvalues(dat.lm)

e.std = e/(slm$sigma*sqrt(1-h))
rstandard(dat.lm)

(e^2*h)/(2*slm$sigma^2*(1-h)^2)
cooks.distance(dat.lm)


dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg02.csv")
plot(dat)
dat.lm<-lm(sales~., data=dat) 
install.packages("DAAG")
library(DAAG)
vif(dat.lm)
lm.price<-lm(price~material+design+dump, data=dat)
p.rsq = summary(lm.price)$r.squared
vif.p = 1/(1-p.rsq) 
summary(dat.lm)

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)

院 認識情報解析

library(rjags)
source("http://peach.l.chiba-u.ac.jp/course_folder/HDI_revised.txt")
dat = read.csv("http://peach.l.chiba-u.ac.jp/course_folder/HtWtData30.csv")

x = dat$height
y = dat$weight

modelS.txt = "
data{
  Ntotal <- length(y)
  xm <- mean(x)
  ym <- mean(y)
  xsd <- sd(x)
  ysd <- sd(y)
  for (i in 1:length(y)){
    zx[i] <- (x[i] - xm)/xsd
    zy[i] <- (y[i] - ym)/ysd
  }
}
model{
  for (i in 1:Ntotal){
    zy[i] ~ dt(zbeta0 + zbeta1 * zx[i], 1/zsigma^2, nu)
  }
  zbeta0 ~ dnorm(0, 1/10^2)
  zbeta1 ~ dnorm(0, 1/10^2)
  zsigma ~ dunif(1.03E-3, 1.03E+3)
  nu <- nuMinusOne + 1
  nuMinusOne ~ dexp(1/29.0)

  #Transfrom to original scale:
  beta1 <- zbeta1 * ysd/xsd
  beta0 <- zbeta0 * ysd + ym -zbeta1*xm*ysd/xsd
  sigma <- zsigma * ysd
}
"
writeLines(modelS.txt, "modelS.txt")
dataList = list(x = x, y = y)
jagsModel =jags.model("modelS.txt", data = dataList, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("beta1", "beta0", "sigma","zbeta0","zbeta1"), n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)

plot(dat$height,dat$weight,xlab="height",ylab ='weight',pch=20,col="orange",cex=5)
par(mfrow=c(2,2))
HDI.plot(mcmcMat[,1])
plot(mcmcMat[,2],mcmcMat[,1], xlab='B1',ylab="B0",pch=20, col='orange')
plot(mcmcMat[,1],mcmcMat[,2], xlab='B0',ylab="B1",pch=20, col='orange')
HDI.plot(mcmcMat[,2])

par(mfrow=c(2,2))
HDI.plot(mcmcMat[,4])
plot(mcmcMat[,5],mcmcMat[,4], xlab='B1',ylab="B0",pch=20, col='orange')
plot(mcmcMat[,4],mcmcMat[,5], xlab='zB0',ylab="zB1",pch=20, col='orange')
HDI.plot(mcmcMat[,5])

n.mcmc = nrow(mcmcMat)
par(mfrow=c(1,1))
temp.x = c(0,80)
temp.y = mcmcMat[1,1]+mcmcMat[1,2]*temp.x
plot(temp.x,temp.y,type='l',lwd=2,col="orange",xlab="height",ylab='weight')
idx = sample(1:n.mcmc,100)
for (i.plot in 1:length(idx)){
  temp.y = mcmcMat[idx[i.plot],1]+mcmcMat[idx[i.plot],2]*temp.x
  lines(temp.x,temp.y,lwd=2,col="orange")
}
points(dat$height,dat$weight,pch=20,cex=4)


par(mfrow=c(1,1))
n2plot = 100
temp.x = c(0,80)
mean.set = seq(50,80,5)
idx=sample(1:nrow(mcmcMat),n2plot)
plot(x,y,xlim=c(50,80),ylim=c(0,400))
for (i_sample in 1:n2plot) {
  temp.y = mcmcMat[idx[i_sample],1]+mcmcMat[idx[i_sample],2]*temp.x
  lines(temp.x,temp.y,lwd=2,col="orange")
  for (i.means in 1:length(mean.set)){
    means = mcmcMat[idx[i_sample],1]+mcmcMat[idx[i_sample],2]*mean.set[i.means]
    y.seq = seq(0,400,length.out=1000)
    dens.y=dt((y.seq-means)/mcmcMat[idx[i_sample],3],29)
    dens.y=dens.y/max(dens.y)
    lines(dens.y,y.seq,type='l',col="orange")
    lines(-dens.y+mean.set[i.means],y.seq,type='l',col="orange")
  }
}
points(x,y,pch=20,cex=3)


model.txt = "
data{
  Ntotal <- length(y)
  xm <- mean(x)
  ym <- mean(y)
  xsd <- sd(x)
  ysd <- sd(y)
  for (i in 1:length(y)){
    zx[i] <- (x[i] - xm) / xsd
    zy[i] <- (y[i] - ym) / ysd
  }
}
model{
  for (i in 1:Ntotal){
    zy[i] ~ dt( zbeta0[s[i]] + zbeta1[s[i]] * zx[i], 1 / zsigma^2, nu)
  }
  for (j in 1:Nsubj){
    zbeta0[j] ~ dnorm( zbeta0mu, 1/(zbeta0sigma)^2)
    zbeta1[j] ~ dnorm( zbeta1mu, 1/(zbeta1sigma)^2)
  }

  zbeta0mu ~ dnorm(0, 1/(10)^2)
  zbeta1mu ~ dnorm(0, 1/(10)^2)
  zsigma ~ dunif(1.0E-3, 1.0E+3)
  zbeta0sigma ~ dunif(1.0E-3, 1.0E+3)
  zbeta1sigma ~ dunif(1.0E-3, 1.0E+3)
  nu <- nuMinusOne + 1
  nuMinusOne ~ dexp(1/29.0)

  for (j in 1:Nsubj){
    beta1[j] <- zbeta1[j] * ysd / xsd
    beta0[j] <- zbeta0[j] * ysd + ym - zbeta1[j] * xm * ysd / xsd
  }
  beta1mu <- zbeta1mu * ysd / xsd
  beta0mu <- zbeta0mu * ysd + ym -zbeta1mu * xm * ysd / xsd
  sigma <- zsigma * ysd
}
"
writeLines(model.txt, "model.txt")

dat = read.csv( file="http://peach.l.chiba-u.ac.jp/course_folder/HierLinRegressData.csv" )
x = dat$X
y = dat$Y
s = dat$Subj

dataList = list(x = x ,  y = y ,s = s , Nsubj = max(s))

jagsModel =jags.model("model.txt", data = dataList, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("beta1mu", "beta0mu", "sigma",  "beta0","beta1"),
                           n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)

par(mfrow=c(1,1))

temp.x = c(40,90)
temp.y = mcmcMat[1,"beta0mu"]+mcmcMat[1,"beta1mu"]*temp.x
plot(temp.x,temp.y,type='l',lwd=2,col="orange",xlab="height",ylab='weight',xlim=c(40,90),ylim=c(40,260))
idx = sample(1:nrow(mcmcMat),100)
for (i.plot in 1:length(idx)){
  temp.y = mcmcMat[idx[i.plot],"beta0mu"]+mcmcMat[idx[i.plot],"beta1mu"]*temp.x
  lines(temp.x,temp.y,lwd=2,col="orange")
}
points(dat$X,dat$Y,pch=20,cex=4)


par(mfrow=c(5,5))
temp.x = c(40,90)
for (i.s in 1:25){
  temp.y = mcmcMat[1,i.s]+mcmcMat[1,i.s+26]*temp.x
  plot(temp.x,temp.y,type='l',lwd=2,col="orange",xlab="height",ylab='weight',xlim=c(40,90),ylim=c(40,260))
  idx = sample(1:nrow(mcmcMat),100)
  for (i.plot in 1:length(idx)){
    temp.y = mcmcMat[idx[i.plot],i.s]+mcmcMat[idx[i.plot],i.s+26]*temp.x
    lines(temp.x,temp.y,lwd=2,col="orange")
  }
  points(dat$X[which(dat$Subj==i.s)],dat$Y[which(dat$Subj==i.s)],pch=20,cex=2)
}

# ch 18
dat = read.csv( file="httP://www.matsuka.info/univ/course_folder/Guber1999data.csv" )
par(mfrow=c(1,1))
plot(dat[,c(2,5,8)],pch=20,cex=3)
y = dat$SATT
x = as.matrix(cbind(dat$Spend,dat$PrcntTake))


dataList = list(x = x ,y = y, Nx = dim(x)[2], Ntotal = dim(x)[1])

model.txt = "
data{
  ym <- mean(y)
  ysd <- sd(y)
  for (i in 1:Ntotal){
    zy[i] <- (y[i] - ym) / ysd
  }
  for (j in 1:Nx){ 
    xm[j] <- mean(x[,j]) 
    xsd[j] <- sd(x[,j])
    for (i in 1:Ntotal){
      zx[i,j] <- (x[i,j] - xm[j])/xsd[j]
    }
  }
}
model{
  for( i in 1:Ntotal){
    zy[i] ~ dt( zbeta0 + sum( zbeta[1:Nx] * zx[i, 1:Nx]), 1/zsigma^2, nu)
  }

  zbeta0 ~ dnorm(0, 1/2^2)
  for ( j in 1:Nx){
    zbeta[j] ~ dnorm(0, 1/2^2)
  }
  zsigma ~ dunif(1.0E-5, 1.0E+1)
  nu <- nuMinusOne + 1
  nuMinusOne ~ dexp(1/29.0)

  beta[1:Nx] <- ( zbeta[1:Nx] / xsd[1:Nx] ) * ysd
  beta0 <- zbeta0 * ysd + ym - sum( zbeta[1:Nx] * xm[1:Nx] / xsd[1:Nx]) * ysd
  sigma <- zsigma * ysd
}
"
writeLines(model.txt, "model.txt")

jagsModel =jags.model("model.txt", data = dataList, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("beta[1]", "beta[2]", "beta0", "sigma"), n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)

plot(as.data.frame(mcmcMat))
par(mfrow=c(1,2))
HDI.plot(mcmcMat[,2])
HDI.plot(mcmcMat[,3])

x = as.matrix(cbind(dat$Spend,dat$PrcntTake,dat$Spend*dat$PrcntTake))

dataList = list(x = x ,y = y, Nx = dim(x)[2], Ntotal = dim(x)[1])
jagsModel =jags.model("model.txt", data = dataList, n.chains = 3, n.adapt = 500)
update(jagsModel, 500)
codaSamples = coda.samples(jagsModel, variable.names = c("beta[1]", "beta[2]", "beta[3]", "beta0", "sigma"), n.iter = ((10000*1)/1), n.adapt = 500)
mcmcMat<-as.matrix(codaSamples)
par(mfrow=c(1,3))
HDI.plot(mcmcMat[,2])
HDI.plot(mcmcMat[,3])
HDI.plot(mcmcMat[,4])

empHDI(mcmcMat[,2])

par(mfrow=c(1,1))
perT = seq(0,80,5)
plot(x=c(perT[1],perT[1]),empHDI(mcmcMat[,2]+mcmcMat[,4]*perT[1]),
     xlim=c(-5,85),ylim=c(-15,35),type='l',lwd=3,col="orange",
     xlab = c("value on Percent Take"),ylab="Slope on Spend")
for (i.plot in 2:length(perT)){
  lines(c(perT[i.plot],perT[i.plot]),empHDI(mcmcMat[,2]+mcmcMat[,4]*perT[i.plot]),
        lwd=3,col="orange")
}
abline(h=0,lwd=2,lty=2)

par(mfrow=c(1,1))
perT = seq(0,10,0.5)
plot(x=c(perT[1],perT[1]),empHDI(mcmcMat[,3]+mcmcMat[,4]*perT[1]),
     xlim=c(-0.5,11),ylim=c(-5,1),type='l',lwd=3,col="orange",
     xlab = c("value on Spent"),ylab="Slope on %Take")
for (i.plot in 2:length(perT)){
  lines(c(perT[i.plot],perT[i.plot]),empHDI(mcmcMat[,3]+mcmcMat[,4]*perT[i.plot]),
        lwd=3,col="orange")
}
abline(h=0,lwd=2,lty=2)

認知情報解析 実習問題

# 修正済み
# 時間がかかる場合は、繰り返し回数を減らしてください。
# 効率が悪いと思うので、適宜変更してください。

temp.state = expand.grid(loc1 = 0:2,loc2=0:2,loc3=0:2,
                         loc4 = 0:2,loc5=0:2,loc6=0:2,
                         loc7 = 0:2,loc8=0:2,loc9=0:2)
temp.state = = expand.grid(rep(list(0:2),9))

n.ones = rowSums(temp.state == 1 )
n.twos = rowSums(temp.state == 2 )
omitTwo = which(n.ones < n.twos)
omitOne = which((n.ones-1 ) > n.twos)
omitUniq = unique(c(omitOne, omitTwo))
state = temp.state[-omitUniq,]
poss.act = apply(state, 1, function(x) which(x==0))

temp.win = matrix(1:9,3)
win.idx = matrix(c(temp.win[1,],temp.win[2,],temp.win[3,],
                   temp.win[,1],temp.win[,2],temp.win[,3],
                   diag(temp.win),
                   diag(temp.win[3:1,])),ncol=3,byrow=T)

idx1 = c()
idx2 = c()
idxCM = c()
for (i.win in 1:nrow(win.idx)){
  idx.temp = apply(state, 1, function(x) sum(x[win.idx[i.win,]]==1)==3)
  idx1 = c(idx1, which(idx.temp))
  idxCM.temp = apply(state, 1, function(x) sum(x[win.idx[i.win,]]==1)==2)
  idxCM = c(idxCM, which(idxCM.temp))
  idx.temp = apply(state, 1, function(x) sum(x[win.idx[i.win,]]==2)==3)
  idx2 = c(idx2, which(idx.temp))
}
n0=apply(state,1,function(x) length((which(x==0))))
tie = which(n0==0)

Q = list()
V = list()
rew.sum = list()
rew.count = list()
policy = list()
for (i.state in 1:nrow(state)){
  Q[[i.state]] =  rep(0,length(poss.act[[i.state]]))
  V[[i.state]] = rep(0,length(poss.act[[i.state]]))
  rew.sum[[i.state]] = rep(0,length(poss.act[[i.state]]))
  rew.count[[i.state]] = rep(0,length(poss.act[[i.state]]))
  policy[[i.state]] = rep(1/length(poss.act[[i.state]]),length(poss.act[[i.state]]))
}

R.W  = 10
R.T  = 5
R.L = -10
gamma = 1
epsilon = 0.05
eta = 1

ck.result <- function(st.idx, idx1, idx2, tie){
  term = F
  rew = 0
  result = "not terminal"
  if (match(st.idx ,idx1, nomatch = 0)!=0){
    rew = R.W
    term = T
    result = "win"
  } else if (match(st.idx ,idx2, nomatch = 0)!=0){
    rew = R.L
    term = T
    result = "lose"
  } else if (match(st.idx ,tie, nomatch = 0)!=0){
    rew = R.T
    term = T
    result = "tie"
  }
  return(list(rew = rew, term = term, result = result))
}

n.rep = 10000
game.hist = rep(0,n.rep)

# main loop
for (i.rep in 1:n.rep){
  st.idx = 1
  term = F
  state.temp = rep(0,9)
  state.hist1 = c()
  state.hist2 = c()
  repeat {
    # playing game
    if (length(poss.act[[st.idx]])==1){
      act1 = poss.act[[st.idx]]
    } else{
      p.act = exp(eta*Q[[st.idx]])/sum(exp(eta*Q[[st.idx]]))
      act1 = sample(poss.act[[st.idx]],1, prob = p.act)
    }
    state.hist1 = rbind(state.hist1,c(st.idx, act1))
    state.temp[act1] = 1
    st.idx = which(apply(state, 1, function(x) sum(x==state.temp) )==9)
    res = ck.result(st.idx, idx1, idx2, tie)
    if (res$term == T){
      rew = res$rew
      break
    }
    p.act = exp(eta*Q[[st.idx]])/sum(exp(eta*Q[[st.idx]]))
    act2 = sample(poss.act[[st.idx]],1, prob = policy[[st.idx]])
    state.hist2 = rbind(state.hist2,c(st.idx, act2))
    state.temp[act2] = 2
    st.idx = which(apply(state, 1, function(x) sum(x==state.temp) )==9)
    res = ck.result(st.idx, idx1, idx2, tie)
    if (res$term == T){
      rew = res$rew
      break
    }
  }
  # update Q & policy
  game.hist[i.rep] = res$result!="lose"
  n.st = nrow(state.hist1)
  #print(res$result)
  if (i.rep%%100==0){print(i.rep)}
  for (i.st in 1:n.st){
    act.idx = which(poss.act[[state.hist1[i.st,1]]]==state.hist1[i.st,2])
    rew.sum[[state.hist1[i.st,1]]][act.idx] = rew.sum[[state.hist1[i.st,1]]][act.idx] + rew
    rew.count[[state.hist1[i.st,1]]][act.idx] = rew.count[[state.hist1[i.st,1]]][act.idx] + 1
    Q[[state.hist1[i.st,1]]][act.idx] = rew.sum[[state.hist1[i.st,1]]][act.idx] / rew.count[[state.hist1[i.st,1]]][act.idx]
  }
}

# plotting results
library(pracma)
game.hist.smooth = movavg(game.hist, 400, type="s")
plot(game.hist.smooth,type='l')