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")
Monthly Archives: July 2019
認知情報解析演習 居住区問題
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"