データ解析基礎論 回帰分析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')

2019 データ解析基礎論A 回帰分析2

# two sample t-test
dat <- read.csv("http://www.matsuka.info/data_folder/hwsk8-17-6.csv")
boxplot(dat$ani, dat$otouto, col = c("skyblue", "orange"),ylab = "Score",
        axnt ="n")
axis(1, at = c(1,2),labels = c("Ani", "Otouto"))
t.test(dat$ani, dat$otouto, var.equal=T)
dat2<-data.frame(score=c(dat$ani,dat$otouto),
                 order=c(rep("ani",10),rep("otouto",10)))

plot(dat2$score~as.numeric(dat2$order), pch=20, xlab="order",
     ylab="score", xlim=c(0.5,2.5), cex=2, xaxt="n")
axis(1, at = c(1,2),labels = c("Ani", "Otouto"))
dat2.lm<-lm(score~order,data=dat2)
abline(dat2.lm,col='red',lwd=3)

# one sample t-test
dat.D = dat$ani - dat$otouto
boxplot(dat.D,col="skyblue",ylab="Difference")
t.test(dat.D)

plot(dat.D~rep(1,10), pch=20, xlab="", ylab="Difference", cex=3, xaxt="n")
dat.D.lm<-lm(dat.D~1)
abline(dat.D.lm,col='red',lwd=3)

# LM with 3 or more categories
dat<-read.csv("http://www.matsuka.info/data_folder/dktb312.csv")
dat2<-data.frame(result=dat$result, method=dat$method,
                 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~c1+c2+c3,data=dat2)
dat.lm <- lm(result~method, dat)

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

# trend analysis
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)
trend.lm2 <- lm(result~duration, contrasts=list(duration = "contr.poly"),data=dat.x)
summary(trend.lm2)

# trend analysis - numeric predictor
set.seed(15)
x = runif(50,0,30)
y = x^0.5 + rnorm(50,9,0.2)
dat.xy = data.frame(y=y,x=x, xsq = x^2)
plot(x,y, xlab = "Hours Studies", ylab="score",
     pch=20, cex=2, col="orange")
summary(lm(y~x+xsq, data=dat.xy))

# regression diagnostic
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)

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)

plot(sort(unif.vars),sort(unif.vars))
plot(sort(norm.vars),sort(unif.vars))

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

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) 

院:認識情報解析 T-test & ANOVA

source("http://www.matsuka.info/univ/course_folder/HDI_revised.txt")
dat <- read.csv("http://peach.l.chiba-u.ac.jp/course_folder/TwoGroupIQ.csv")
dat2sd <- dat[dat$Group=="Smart Drug",]

dataList <- list(y = dat2sd$Score, Ntotal = nrow(dat2sd),
                 meanY = mean(dat2sd$Score), sdY = sd(dat2sd$Score))

model.txt = "
model {
  for (i in 1:Ntotal){
    y[i] ~ dnorm(mu, 1/sigma^2)
  }
  mu ~ dnorm(meanY, 1/(100*sdY)^2)
  sigma ~ dunif(sdY/1000, sdY*1000)
}"
writeLines(model.txt, "model.txt")
parameters = c( "mu" ,  "sigma")
library("rjags")
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)
traceplot(codaSamples)
gelman.plot(codaSamples)
mcmcMat<-as.matrix(codaSamples)

HDI.plot(mcmcMat[,1])
HDI.plot(mcmcMat[,2])
hist(dat2sd$Score, probability = T)
x.temp = seq(40,220,0.1)
n.plot = 100
rand.order = sample(1:nrow(mcmcMat),n.plot)
for (i.plot in 1:n.plot){
  y.temp = dnorm(x.temp, mean = mcmcMat[rand.order[i.plot],1],
                 sd = mcmcMat[rand.order[i.plot],2])
  lines(x.temp, y.temp, type='l',col="orange")
}

#y = c(-2:2,15);Ntotal = length(y);meanY = mean(y); sdY = sd(y)
#dataList = list(y= y, Ntotal=Ntotal, meanY = meanY, sdY = sdY)
dataList <- list(y = dat2sd$Score, Ntotal = nrow(dat2sd),
                 meanY = mean(dat2sd$Score), sdY = sd(dat2sd$Score))

model.txt ="
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dt( mu , 1/sigma^2 , nu )
  }
  mu ~ dnorm( meanY , 1/(100*sdY)^2 )
  sigma ~ dunif( sdY/1000 , sdY*1000 )
  nu <- nuMinusOne+1
  nuMinusOne ~ dexp(1/29)
}"
writeLines(model.txt, "model.txt")
parameters = c( "mu" ,  "sigma", "nu")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=1000 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
#traceplot(codaSamples)
#gelman.plot(codaSamples)
mcmcMat<-as.matrix(codaSamples)
HDI.plot(mcmcMat[,1])
HDI.plot(mcmcMat[,2])

dataList <- list(y = dat$Score, Ntotal = nrow(dat),
                 meanY = mean(dat$Score), sdY = sd(dat$Score),
                 Ngroup = length(unique(dat$Group)),
                 group = dat$Group)

model.txt ="
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dt( mu[group[i]] , 1/sigma[group[i]]^2 , nu )
  }
  for (j in 1:Ngroup){
    mu[j] ~ dnorm( meanY , 1/(100*sdY)^2 )
    sigma[j] ~ dunif( sdY/1000 , sdY*1000 )
  }
  nu <- nuMinusOne+1
  nuMinusOne ~ dexp(1/29)
}"
writeLines(model.txt, "model.txt")
parameters = c( "mu" ,  "sigma", "nu")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=1000 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
#traceplot(codaSamples)
#gelman.plot(codaSamples)
mcmcMat<-as.matrix(codaSamples)
HDI.plot(mcmcMat[,1]-mcmcMat[,2])
HDI.plot(mcmcMat[,2])


dat<-read.csv("http://www.matsuka.info/univ/course_folder/FruitflyDataReduced.csv")
y=dat$Longevity
yMean = mean(y)
ySD = sd(y)
Ntotal = length(y)

MeanSD2gamma <- function( mean, sd ) {
  shape = mean^2 / sd^2
  rate = mean / sd^2
  return(data.frame(shape,rate))
}

ModeSD2gamma <- function( mode, sd ) {
  rate = ( mode + sqrt( mode^2 + 4 * sd^2 ) )/( 2 * sd^2 )
  shape = 1 + mode * rate
  return(data.frame(shape,rate))
}

temp.gamma = ModeSD2gamma( mode=sd(y)/2 , sd=2*sd(y) )
gShape = temp.gamma$shape
gRate = temp.gamma$rate

x=as.numeric(dat$CompanionNumber)
Ngroup = length(unique(x))
xc=dat$Thorax
xcMean=mean(xc)

dataList = list(
  y = y ,
  x = x ,
  Ntotal = Ntotal ,
  Ngroup = Ngroup ,
  yMean = yMean ,
  ySD = ySD ,
  gShape = gShape,
  gRate  = gRate
)

model.txt="
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dnorm( a0 + a[ x[ i_data ] ], 1/ySigma^2)
  }

  ySigma ~ dunif( ySD/100, ySD*10 )
  a0 ~ dnorm( yMean, 1/(ySD*5)^2)

  for ( i_group in 1:Ngroup ) {
    a[ i_group ] ~ dnorm(0.0, 1/aSigma^2)
  }
  aSigma ~ dgamma( gShape, gRate )

  for ( i_group in 1:Ngroup ) { m[ i_group ] <- a0 + a[ i_group ] }
  b0 <- mean( m[ 1:Ngroup ] )
  for (i_data in 1:Ngroup ) { b[ i_data ] <- m[ i_data ] - b0 }
  }
"
writeLines(model.txt, "model.txt")

parameters = c( "b0" ,  "b" ,  "ySigma",  "aSigma" )

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)

# fig 19.3
plot(x,dat$Longevity,xlim=c(0,5.1),ylim=c(0,120))
for (i_c in 1:100) {
  sd=mcmcMat[i_c,"ySigma"]
  lim=qnorm(c(0.025,0.975))
  means=mcmcMat[i_c,2:6]+mcmcMat[i_c,"b0"]
  for (i_group in 1:5){
    lower=means[i_group]+lim[1]*sd
    upper=means[i_group]+lim[2]*sd
    yseq=seq(lower,upper,length.out = 100)
    dens.y=dnorm((yseq-means[i_group])/sd)
    dens.y=dens.y/max(dens.y)*0.75
    lines(-dens.y+i_group,yseq,type='l',col="orange")
  }
}




dataList = list(
  y = y ,
  x = x ,
  xc = xc ,
  Ntotal = Ntotal ,
  Ngroup = Ngroup ,
  yMean = yMean ,
  ySD = ySD ,
  xcMean = xcMean ,
  acSD = sd(xc) ,
  gShape = gShape,
  gRate  = gRate
)
model.txt="
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dnorm( mu[i_data], 1/ySigma^2)
    mu [ i_data] <- a0 + a[ x[ i_data ] ] + ac*( xc[ i_data ] - xcMean )
  }

  ySigma ~ dunif( ySD/100, ySD*10 )
  a0 ~ dnorm( yMean, 1/(ySD*5)^2)

  for ( i_group in 1:Ngroup ) {
    a[ i_group ] ~ dnorm(0.0, 1/aSigma^2)
  }
  aSigma ~ dgamma( gShape, gRate )
  ac ~ dnorm(0, 1/(2*ySD/acSD)^2 )
  b0 <- a0 + mean( a[ 1:Ngroup ] ) - ac*xcMean
  for (i_data in 1:Ngroup ) { b[ i_data ] <- a[ i_data ] - mean( a[1:Ngroup]) }
}"
writeLines(model.txt, "model.txt")
parameters = c( "b0" ,  "b" ,  "ySigma",  "aSigma", "ac" )

jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=5000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)

# fig 19.5 only "none0" will be plotted
plot(Longevity~Thorax,xlim=c(0.6,1),ylim=c(0,120),data=dat[dat$CompanionNumber=="None0",],type='n')
Thorax.value=c(0.65,0.8,0.95)
for (i_c in 1:100) {

  # regression lines
  xs=c(0.5,1.1)
  ys=mcmcMat[i_c,"b0"]+mcmcMat[i_c,3]+mcmcMat[i_c,"ac"]*xs
  lines(xs,ys,col="orange")

  # density
  sd=mcmcMat[i_c,"ySigma"]
  lim=qnorm(c(0.025,0.975))
  means=mcmcMat[i_c,"b0"]+mcmcMat[i_c,3]+mcmcMat[i_c,"ac"]*Thorax.value
  for (i_thorax in 1:3){
    lower=means[i_thorax]+lim[1]*sd
    upper=means[i_thorax]+lim[2]*sd
    yseq=seq(lower,upper,length.out = 100)
    dens.y=dnorm((yseq-means[i_thorax])/sd)
    dens.y=dens.y/max(dens.y)*0.05
    lines(-dens.y+Thorax.value[i_thorax],yseq,type='l',col="orange")
  }
}
abline(v=Thorax.value,col='green')
points(Longevity~Thorax,xlim=c(0.6,1),ylim=c(0,120),data=dat[dat$CompanionNumber=="None0",])



dat<-read.csv("http://www.matsuka.info/univ/course_folder/NonhomogVarData.csv")
y=dat$Y
yMean = mean(y)
ySD = sd(y)
Ntotal = length(y)
x=as.numeric(dat$Group)
Ngroup = length(unique(x))

model.txt = "
model {
  for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dt( a0 + a[ x[ i_data ] ], 1/ySigma[x[i_data]]^2, nu)
}
nu <- nuMinusOne + 1
nuMinusOne ~ dexp(1/29)
for (i_group in 1:Ngroup) {
ySigma[i_group]~dgamma(ySigSh,ySigR)
}
ySigSh <- 1+ySigMode * ySigR
ySigR  <- ((ySigMode + sqrt(ySigMode^2+4*ySigSD^2))/(2*ySigSD^2))

ySigMode ~ dgamma(gShape,gRate)
ySigSD ~ dgamma(gShape,gRate)
a0 ~ dnorm( yMean, 1/(ySD*10)^2)

for ( i_group in 1:Ngroup ) {
a[ i_group ] ~ dnorm(0.0, 1/aSigma^2)
}
aSigma ~ dgamma( gShape, gRate )

for ( i_group in 1:Ngroup ) {
m[ i_group ] <- a0 + a[ i_group ]
}

b0 <- mean( m[ 1:Ngroup ] )

for (i_data in 1:Ngroup ) {
b[ i_data ] <- m[ i_data ] - b0
}
}"
writeLines(model.txt, "model.txt")

dataList = list(
  y = y ,
  x = x ,
  Ntotal = Ntotal ,
  Ngroup = Ngroup ,
  yMean = yMean ,
  ySD = ySD ,
  gShape = gShape,
  gRate  = gRate
)

parameters = c( "b0" ,  "b" ,  "ySigma",  "aSigma" ,"nu","ySigMode","ySigSD")

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)


# plotting
plot(x,y,xlim=c(0,4.1),ylim=c(70,130))
n2plot=100
idx=sample(1:nrow(mcmcMat),n2plot)
for (i_sample in 1:n2plot) {
  i_c = idx[i_sample]
  lim=qnorm(c(0.025,0.975))
  means=mcmcMat[i_c,2:6]+mcmcMat[i_c,"b0"]
  for (i_group in 1:4){
    sd=mcmcMat[i_c,9+i_group]
    lower=means[i_group]+lim[1]*sd
    upper=means[i_group]+lim[2]*sd
    yseq=seq(lower,upper,length.out = 100)
    dens.y=dnorm((yseq-means[i_group])/sd)
    dens.y=dens.y/max(dens.y)*0.75
    lines(-dens.y+i_group,yseq,type='l',col="orange")
  }
}

認知情報解析演習a Monte Carlo 01

north=c(1:3,15,1:10)
east=2:15;east[ c(3,7,11)]=c(3,7,11)
south=c(5:15,12:14)
west=c(15,1:13);west[ c(4,8,12)]=c(4,8,12)
trM=cbind(north,east,south,west)

r=-1;P=rep(0.25,4);V = rep(0,14)max.iter = 10000;
state.count=rep(0,15)
for (i.iter in 1:max.iter){
  state = sample(1:14,1)
  state.seq = state
  while(state!=15){
    action = sample(1:4,1,prob = P)
    state.seq = c(state.seq,trM[state,action])
    state = trM[state,action]  
  }
  uniq.seq = unique(state.seq)
  for (i.uniq in 1:(length(uniq.seq)-1)){
    first.visit = which(state.seq == uniq.seq[i.uniq])[1]
    V[uniq.seq[i.uniq]] = V[uniq.seq[i.uniq]] + r*(length(state.seq)-first.visit-1)
  }
  state.count[uniq.seq] = state.count[uniq.seq] + 1
}
V = matrix(c(0,V/state.count[1:14],0),nrow=4)

2019データ解析基礎論a 回帰分析

dat <- read.csv("http://www.matsuka.info/data_folder/hwsk8-17-6.csv")
plot(dat$otouto,dat$ani,pch=20, cex =3,
     xlab = "score of younger brother",
     ylab = "score of elder brother")

dat.lm <- lm(ani~otouto, data=dat)
summary(dat.lm)
abline(dat.lm, col = 'red',lwd = 2.5)

dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
dat.reg1<-lm(sales~material,data=dat)
dat.reg2<-lm(sales~price,data=dat)
dat.reg3<-lm(sales~design,data=dat)

dat.regALL<-lm(sales~material+price+design,data=dat)
dat.regALL<-lm(sales~.,data=dat)

dat <- read.csv("http://www.matsuka.info/data_folder/hwsk8-17-6.csv")
t.test(dat$ani, dat$otouto, var.equal=T)
dat2<-data.frame(score=c(dat$ani,dat$otouto),order=c(rep("ani",10),rep("otouto",10)))
plot(dat2$score~as.numeric(dat2$order), pch=20, xlab="order",
     ylab="score", xlim=c(0.5,2.5), cex=2,xaxt="n")
axis(1,c(1,2),c("ani","otouto"))
dat2.lm<-lm(score~order,data=dat2)
abline(dat2.lm,col='red',lwd=3)

dat.D = dat$ani - dat$otouto
boxplot(dat.D,col="skyblue",ylab="Difference")
t.test(dat.D)
plot(dat.D~rep(1,10),pch=20,xlab="",ylab="Difference",cex=3)
dat.D.lm<-lm(dat.D~1)
abline(dat.D.lm,col='red',lwd=3)

院:認知情報解析

y = sample(c(rep(1,15), rep(0,35)))
Ntotal=length(y)
datalist = list(y=y,Ntotal=Ntotal)
source("http://www.matsuka.info/univ/course_folder/HDI_revised.txt")

library(rjags)
txt = "
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dbern( theta )
  }
  theta ~ dbeta( 1, 1 )
}"
writeLines(txt, "~/model.txt")
jagsModel = jags.model(file="~/model.txt",
                       data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)

HDI.plot(mcmcMat)
traceplot(codaSamples)
autocorr.plot(codaSamples,type='l')
gelman.plot(codaSamples)

y1 = sample(c(rep(1,6), rep(0,2)))
y2 = sample(c(rep(1,2), rep(0,5))) 
y = c(y1,y2)
s = c(rep(1,8),rep(2,7))
Ntotal=length(y)
datalist = list(y = y, Ntotal = Ntotal, s = s, Nsubj = 2)

txt = "
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern( theta[s[i_data]] )
}
for ( i_s in 1:Nsubj) {
theta[i_s] ~ dbeta( 2, 2 )
}
}"
writeLines(txt, "~/model.txt")
jagsModel = jags.model(file="~/model.txt",
                       data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)
# checking MCMC
HDI.plot(mcmcMat[,2])
traceplot(codaSamples)
autocorr.plot(codaSamples)
gelman.plot(codaSamples)

x.temp = seq(0,150,0.1)
g1 = dgamma(x.temp, shape =0.01, rate =0.01)
plot(x.temp,g1,type='l')

g2 = dgamma(x.temp, shape =1.56, rate = 0.0312)
plot(x.temp,g2,type='l')

g3 = dgamma(x.temp, shape =6.25, rate = 0.125)
plot(x.temp,g2,type='l')

txt = "
model {
for ( i_data in 1:Ntotal ) {
y[ i_data ] ~ dbern( theta[s[i_data]] )
}
for ( i_s in 1:Nsubj) {
theta[i_s] ~ dbeta( omega*(kappa-2)+1,(1-omega)*(kappa-2)+1)
}
omega ~ dbeta(1,1)
kappa <- kappaMinusTwo + 2
kappaMinusTwo ~ dgamma(0.01, 0.01)
}"
writeLines(txt, "~/model.txt")


dat<-read.csv("http://www.matsuka.info/data_folder/TherapeuticTouchData.csv")
y=dat$y
s=as.numeric(dat$s)
Ntotal=length(dat$y)
Nsubj=length(unique(s))
datalist = list(y=y,s=s,Ntotal=Ntotal,Nsubj=Nsubj)

jagsModel = jags.model(file="~/model.txt",data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta","omega","kappa"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)

dat<-read.csv("http://www.matsuka.info/data_folder/BattingAverage.csv")
z = dat$Hits
N = dat$AtBats
s = dat$PlayerNumber
c = dat$PriPosNumber
Ncat = 9
Nsubj = 948

datalist = list(z = z, N=N, s=s, c=c, Ncat=Ncat, Nsubj =Nsubj )
txt = "
model {
 for (i_s in 1:Nsubj) {
   z[i_s] ~ dbin( theta[i_s], N[i_s])
   theta[i_s] ~ dbeta(omega[c[i_s]]*(kappa[c[i_s]]-2)+1,
      (1-omega[c[i_s]])*(kappa[c[i_s]]-2)+1)
 }
 for (i_c in 1:Ncat) {
   omega[i_c] ~ dbeta(omegaO*(kappaO-2)+1,(1-omegaO)*(kappaO-2)+1)
   kappa[i_c] <- kappaMinusTwo[i_c]+2
   kappaMinusTwo[i_c] ~ dgamma(0.01,0.01)
 }
 omegaO ~ dbeta(1,1)
 kappaO <- kappaMinusTwoO+2
 kappaMinusTwoO ~ dgamma(0.01,0.01)
}"
writeLines(txt, "~/model.txt")
jagsModel = jags.model(file="~/model.txt",data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta","omega","kappa"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)
HDI.plot(mcmcMat[," omega[1]"]) #pitcher
HDI.plot(mcmcMat[," omega[2]"]) #catcher
HDI.plot(mcmcMat[," omega[1]"]) #1st base

2019 基礎実習MA01

# random number generators
x=rnorm(n=1,mean=100,sd=15)
y=runif(n=3,min=1,max=10)

N = 10000
# N = 1000
random.data = rnorm(N, mean=0, sd=1)
hist(random.data, nclass = 50, col = "navy", xlab = "Data",
  probability = T, main = "Histogram of Random Data")

# density of generated data
dens = density(random.data)
lines(dens, col = "orange", lwd = 4)

# theoretical density
x = seq(-4,4,0.1)  
true.norm = dnorm(x, mean = 0, sd = 1)
lines(x,true.norm, col = "green", lty = 3, lwd = 4)
legend("topleft",c("empirical", "theoretical"), lty = c(1,3),
  col = c('orange','green'),lwd=4)

# random sampling
sample(1:10,3)   
sample(c(“gu”,“choki”,“pa”),1)
sample(1:10)
sample(0:1, 10, replace=T)
sample(c("Head","Tail"), 10, replace=T)
sample(c("Head","Tail"), 10, replace=T, prob=c(0.9,0.1)) 

# flow control
for (i_loop in 1:5){print(i_loop)}
counter <- 1

while(counter<=10){
  print(counter)
  counter<-counter+1
}

counter <- 1
while(counter^2 <= 10){
  print(c(counter, counter^2))
  counter<-counter+1
}

affil<-"cogsci"
if (affil=="cogsci") {
    print("you are wonderful")
}

affil<-"phil"
if (affil=="cogsci") {
   print("you are wonderful")
 } else {
   print("still, you are wonderful")
}

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

counter=6
repeat{ 
  print(counter)
  counter = counter + 1
  if(counter>5){break}
}

counter=6
repeat{
  if(counter>5){break}
  print(counter)
  counter+counter+1
}

認知情報解析06/05の演習問題

演習で書いたプログラムで問題はありませんでした。
以下、結果を可視化するコマンドを書き加えました:

# initalization
p.win =0.4; p.lose = 0.6; P = c(p.lose, p.win);
R = c(rep(0,100), 1); V = rep(0,101);
gamma = 1; tol = 1e-10; counter=0
cols = c("red","skyblue",'orange','black')
par(mfrow=c(1,2))

# value iteration
repeat{
  counter = counter+1
  delta = 0
  for (i.state in 2:100) {
    v <- V[i.state]
    temp.v = rep(0, (i.state - 1))
    for (i.bet in 1:(i.state - 1)) {
      lower.B = i.state - i.bet
      upper.B = min(i.state + i.bet, 101)
      temp.v[i.bet] = sum(P * (R[c(lower.B, upper.B)] + gamma * V[c(lower.B, upper.B)]))
      V[i.state] = max(temp.v)
    }
    delta <- max(abs(v-V[i.state]), delta)
  }
  # plotting results
  if (counter==1){
    plot(V[2:100], type='l', col=cols[1], lwd=2, xlab="capital", ylab="value")
  } else {
    if (counter<4){
      lines(V[2:100], type='l', col=cols[counter], lwd=2)
    }
  }
  if (delta < tol){break}
}
# end of value iteration

lines(V[2:100],type='l', col=cols[4], lwd=2)
legend("topleft", c("1st sweep","2nd sweep", "3rd sweep","32nd sweep") ,col=cols, lwd=2)

# identifying optimal action
policy = rep(0,101)
for (i.state in 2:100) {
  temp.v = rep(0, (i.state - 1))
  for (i.bet in 1:(i.state - 1)) {
    lower.B = i.state - i.bet
    upper.B = min(i.state + i.bet, 101)
    temp.v[i.bet] = sum(P * (R[c(lower.B, upper.B)] + gamma * V[c(lower.B, upper.B)]))
    policy[i.state] = which.max(round(temp.v,4))
  }
}
barplot(policy,xlab="capital",ylab="Optimal action")