# データ解析基礎論a 分散分析３

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

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)

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.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 分散分析２

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

# anova
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
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
cu.scheffe1F(dat.aov,dat,c(-3,1,1,1))

# 2 way ANOVA
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
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
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")

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

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

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

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 分散分析１

```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.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.aov<-aov(result~method, data=dat)
summary(dat.aov)
source("http://www.matsuka.info/univ/course_folder/cutil.R")

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

# 基礎実習　MA0２

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