library(tidyverse) dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/logisticReg01.txt ") ggplot(dat) + geom_point(aes(x = study, y = pass), size =3) + xlab("Hours studied") + ylab("Passing") dat.lm <- lm(pass~study,dat) ggplot(dat,aes(x = study, y = pass)) + geom_point(size =3) + xlab("Hours studied") + ylab("Passing") + geom_abline(slope = dat.lm$coefficients[2], intercept = dat.lm$coefficients[1]+1, color = "red",size = 2) par(mfrow=c(2,2)) plot(dat.lm) real.p = data.frame( real.p = table(dat$pass, dat$study)[2,] / colSums(table(dat$pass, dat$study)), x = 0:30) ggplot(real.p,aes(x = x, y = real.p)) + geom_point(size =3) + xlab("Hours studied") + ylab("Passing (actual probability)") dat.lr <- glm(pass~study,family=binomial, data=dat) summary(dat.lr) coef = coefficients(dat.lr) temp.x = seq(0,30,0.1) pred.pass.p = data.frame(pred.p = 1/(1+exp(-(coef[1]+coef[2]*temp.x))), x = temp.x) ggplot(dat, aes(x = study,y = pass)) + geom_point(size=3) + geom_line(aes(x =x, y= pred.p + 1), data = pred.pass.p, color = 'red',size = 1)+ xlab("Hours studied") + ylab("Passing") ## pred.pass.p = 1/(1+exp(-(coef[1]+coef[2]*c(10:15)))) odds=pred.pass.p/(1-pred.pass.p) exp(coef[2]) odds[2:6]/odds[1:5] dat<-read.csv("http://www.matsuka.info/data_folder/datWA01.txt") dat.lr<-glm(gender~shoesize,family=binomial,data=dat) anova(dat.lr, test ="Chisq") dat.lr0<-glm(gender~1,family="binomial",data=dat) dat.lrS<-glm(gender~shoesize,family=binomial,data=dat) dat.lrh<-glm(gender~h,family="binomial",data=dat) M=matrix(c(52,48,8,42),nrow=2) rownames(M)<-c("present", "absent") colnames(M)<-c("smoker",'non-smoker') dat<-as.data.frame((as.table(M))) colnames(dat)<-c("cancer","smoking","freq") dat=dat[rep(1:nrow(dat),dat$freq),1:2] rownames(dat)<-c() dat.glm<-glm(cancer~smoking,family=binomial,data=dat) dat<-read.csv("http://www.matsuka.info/data_folder/cda7-16.csv") dat.glm<-glm(survival~age, family=binomial,data=dat) dat.glm2<-glm(survival~Ncigarettes,family=binomial,data=dat) dat.glm3<-glm(survival~NdaysGESTATION,family=binomial,data=dat) dat.glmAllAdd=glm(survival~age+Ncigarettes+NdaysGESTATION,family=binomial,data=dat) dat.glmAllMult=glm(survival~age*Ncigarettes*NdaysGESTATION,family=binomial,data=dat) library(MASS) stepAIC(dat.glmAllMult) dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/poissonReg01.txt ") ggplot(dat) + geom_histogram(aes(eye.count), fill='red') + xlab("Number of times looking at eyes") ggplot(dat, aes(x=attr, y = eye.count)) + geom_point(size =2) + ylab("Number of times looking at eyes")+ xlab("Attractiveness") + geom_abline(slope = dat.lm$coefficients[2], intercept = dat.lm$coefficients[1], color = "red",size = 2) dat.lm <- lm(eye.count~attr,data = dat) dat.pr<-glm(eye.count~attr,family = "poisson",data=dat) cf = coefficients(dat.pr) x.temp <- seq(0,10,0.1) pred.c = data.frame(x=x.temp, pred = exp(cf[1]+cf[2]*x.temp)) ggplot(dat, aes(x=attr, y = eye.count)) + geom_point(size =2) + ylab("Number of times looking at eyes")+ xlab("Attractiveness") + geom_abline(slope = dat.lm$coefficients[2], intercept = dat.lm$coefficients[1], color = "red",size = 2)+ geom_line(aes(x = x.temp, y= pred),data = pred.c, color="blue", size=2)
Monthly Archives: November 2019
データ解析基礎論B 多変量解析
install.packages("ggfortify") install.packages("ggdendro") library(ggfortify) library(ggdendro) # pca dat<-read.table("http://www.matsuka.info/data_folder/tdkPCA01.txt") dat.pca<-princomp(dat) autoplot(dat.pca, label = TRUE, label.size = 6, loadings = TRUE, loadings.colour = 'red', loadings.label = TRUE, loadings.label.size = 5) autoplot(dat.pca, shape = FALSE, label.size = 6, loadings = TRUE, loadings.colour = 'red', loadings.label = TRUE, loadings.label.size = 5) cldata<-data.frame(var1=c(4,1,5,1,5), var2=c(1,5,4,3,1)) rownames(cldata) = c("A","B","C","D","E") autoplot(dist(cldata)) cldata.cluster=hclust(dist(cldata),method="average") ggdendrogram(cldata.cluster, rotate = T, theme_dendro = FALSE)+ xlab("Individual") dat<-read.csv("http://matsuka.info/data_folder/tdkClust.csv", header=TRUE, row.names=1) autoplot(dist(dat)) dat.cluster=hclust(dist(dat)) ggdendrogram(dat.cluster, rotate = T, theme_dendro = FALSE)+ xlab("Occupation") dat.pca = princomp(dat) autoplot(dat.pca, label = TRUE, shape = FALSE, label.size = 4, loadings = TRUE, loadings.colour = 'red', loadings.label = TRUE, loadings.label.size = 5) dat.HC.S=hclust(dist(dat), method = "single") dat.HC.C=hclust(dist(dat), method = "complete") dat.HC.A=hclust(dist(dat), method = "average") dat.HC.W=hclust(dist(dat), method = "ward.D") ggdendrogram(dat.HC.S, rotate = T, theme_dendro = FALSE)+ xlab("Occupation")+ggtitle("Method = Single") ggdendrogram(dat.HC.C, rotate = T, theme_dendro = FALSE)+ xlab("Occupation")+ggtitle("Method = Complete") ggdendrogram(dat.HC.A, rotate = T, theme_dendro = FALSE)+ xlab("Occupation")+ggtitle("Method = Average") ggdendrogram(dat.HC.W, rotate = T, theme_dendro = FALSE)+ xlab("Occupation")+ggtitle("Method = Ward's MV") dat.kmeans=kmeans(dat, centers=3, nstart=10) pairs(dat, main = "Clustering Occupations", pch = 21, bg = c("red", "blue", "green") [unclass(dat.kmeans$cluster)]) autoplot(dat.kmeans, dat, size = 3, label = TRUE, label.size = 5) source("http://www.matsuka.info/univ/course_folder/cuUtil02.R") res<-cu.KMC.rep(dat,10,100) autoplot(dat.kmeans, dat, frame = TRUE, frame.type = 'norm') + ylim(-0.7,0.7)+xlim(-1.2,0.7) autoplot(dat.kmeans, dat, frame = TRUE)+ ylim(-0.7,0.7)+xlim(-1.2,0.7) dat<-data.frame(writing=c(68,85,50,54,66,35,56,25,43,70), interview=c(65,80,95,70,75,55,65,75,50,40), cl=c(rep("A",5),rep("N",5))) library(MASS) dat.lda<-lda(cl~.,data=dat) intcpt = (dat.lda$scaling[1]*dat.lda$means[1,1]+dat.lda$scaling[2]*dat.lda$means[1,2]+ dat.lda$scaling[1]*dat.lda$means[2,1]+dat.lda$scaling[2]*dat.lda$means[2,2])/2 new.dim.slope = dat.lda$scaling[1]/dat.lda$scaling[2] disc.intcpt = intcpt / dat.lda$scaling[2] disc.slope = -dat.lda$scaling[1] / dat.lda$scaling[2] ggplot(dat, aes(x = writing, y= interview, color = cl)) + geom_point(size = 4) + geom_abline(aes(intercept = intcpt, slope = new.dim.slope )) + geom_abline(aes(intercept = disc.intcpt, slope = disc.slope ),color = "red") + xlim(30,100)+ylim(30,100) dat<-read.csv("http://matsuka.info/data_folder/tdkDA01.csv", header=T) dat.lda<-lda(class~.,dat) lda.pred<-predict(dat.lda,dat) table(lda.pred$class, dat$class) dat.ldaCV<-lda(class~.,dat, CV=T) dat<-read.csv("http://matsuka.info/data_folder/tdkDA02.csv",header=T) dat.lda=lda(class~.,data=dat) lda.pred <- predict(dat.lda)$x %>% as.data.frame %>% cbind(class = dat$class) ggplot(lda.pred) + geom_point(aes(x=LD1, y=LD2, color = class), size = 2.5) dat<-data.frame(p1=c(4,1,5,1,5),p2=c(1,5,4,3,1)) rownames(dat)<-c('a','b','c','d','e') dat.mds<-cmdscale(dist(dat),2) ggplot(dat.mds, aes(x = dat.mds[,1],y = dat.mds[,2])) + geom_text(aes(label = row.names(dat.mds)), size = 6)
dynamic programming 実装例
V=c(rep(0,100),1);V.hist=c() p=c(0.4,0.6); gamma=1;delta=1; tol=1e-20 max.a=rep(0,101) while (delta>tol) { delta=0; for (i_state in 1:99) { v=V[i_state+1] temp=matrix(0,nrow=1,ncol=i_state) for (i_action in 1:i_state) { temp[i_action]=sum(p*(gamma*c(V[(min(i_state+i_action,100)+1)], V[(max(i_state-i_action,0)+1)]))) } V[i_state+1]=max(temp) max.a[i_state+1]=which.max(round(temp,8)) delta=max(delta,abs(v-V[i_state+1])) } V.hist=rbind(V.hist,V) } # plotting results par(mfrow=c(1,2)) plot(V.hist[1,],type='l',lwd=2,xlab="Capital",ylab="Value Estimates",col='red') lines(V.hist[2,],lwd=2,col='blue') lines(V.hist[3,],lwd=2,col='green') lines(V.hist[32,],lwd=2,col='black') legend("topleft",c("sweep 1","sweep 2","sweep 3", "sweep 32"), col=c("red","blue","green","black"),lwd=2) barplot(max.a,xlab="Capital",ylab="Final Policy",col="white")
広域システム特別講義II 教師あり学習1a
# with THRESHOLD (theta) AND.gate <- function(x1, x2){ w1 = 0.5 w2 = 0.5 theta = 0.7 y.temp = w1*x1 + w2*x2 if (y.temp <= theta){ y = 0 } else { y = 1 } return(y) } AND.gate <- function(x1, x2){ w1 = 0.5; w2 = 0.5; theta = 0.7 return(as.numeric(w1*x1 + w2*x2 > theta)) } NAND.gate <- function(x1, x2){ w1 = -0.5; w2 = -0.5; theta = -0.7 return(as.numeric(w1*x1 + w2*x2 > theta)) } OR.gate <- function(x1, x2){ w1 = 0.5; w2 = 0.5; theta = 0.3 return(as.numeric(w1*x1 + w2*x2 > theta)) } # with BIAS (b) AND.gate <- function(x1, x2){ w1 = 0.5 w2 = 0.5 b = -0.7 y.temp = w1*x1 + w2*x2 + b if (y.temp <= 0){ y = 0 } else { y = 1 } return(y) } AND.gate <- function(x1, x2){ w1 = 0.5; w2 = 0.5; b = -0.7 return(as.numeric(w1*x1 + w2*x2 + b > 0)) } NAND.gate <- function(x1, x2){ w1 = -0.5; w2 = -0.5; b = 0.7 return(as.numeric(w1*x1 + w2*x2 + b > 0)) } OR.gate <- function(x1, x2){ w1 = 0.5; w2 = 0.5; b = -0.3 return(as.numeric(w1*x1 + w2*x2 + b > 0)) } NOR.gate <- function(x1, x2){ w1 = -0.5; w2 = -0.5; b = 0.3 return(as.numeric(w1*x1 + w2*x2 + b > 0)) } plot.logic <- function(logic.oper){ x1 = c(0,0,1,1); x2 = c(0,1,0,1); if (logic.oper == "and") { w1 = 0.5; w2 = 0.5; theta = 0.7; true.point = AND.gate(x1,x2) } else if (logic.oper == "or") { w1 = 0.5; w2 = 0.5; theta = 0.3; true.point = OR.gate(x1,x2) } else if (logic.oper == "nand") { w1 = -0.5; w2 = -0.5; theta = -0.7; true.point = NAND.gate(x1,x2) } else if (logic.oper == "nor"){ w1 = -0.5; w2 = -0.5; theta = -0.3; true.point = NOR.gate(x1,x2) } else {warning("incompatible operator");stop() } plot(c(0,0,1,1),c(0,1,0,1),xlim = c(-0.5, 1.5), ylim = c(-0.5, 1.5), pch = 20, cex= 2, col = true.point+1) abline(a = theta/w1, b = -w1/w2, lwd = 3) } XOR.gate <- function(x1, x2){ gate1 <- NAND.gate(x1,x2) gate2 <- OR.gate(x1,x2) y <- AND.gate(gate1,gate2) return(y) } plot.XOR <- function(){ x1 = c(0,0,1,1); x2 = c(0,1,0,1); w11 = -0.5; w21 = -0.5; theta1 = -0.7 w12 = 0.5; w22 = 0.5; theta2 = 0.3 true.point = XOR.gate(x1, x2) plot(c(0,0,1,1),c(0,1,0,1),xlim = c(-0.5, 1.5), ylim = c(-0.5, 1.5), pch = 20, cex= 2, col = true.point+1) abline(a = theta1/w11, b = -w11/w21, lwd = 3) abline(a = theta2/w12, b = -w12/w22, lwd = 3) } multi.forwd <- function(x,y){ return(x*y) } multi.bckwd <- function(x, y, dout){ dx = dout * y dy = dout * x return(list(dx = dx, dy = dy)) } apple = 100; n.apple = 2; tax = 1.1 apple.pre.tax = multi.forwd(apple, n.apple) apple.post.tax = multi.forwd(apple.pre.tax, tax) dprice = 1 d.apple.post.tax = multi.bckwd(apple.pre.tax, tax, dprice) d.apple = multi.bckwd(apple, n.apple, d.apple.post.tax$dx)$dx d.n.apple = multi.bckwd(apple, n.apple, d.apple.post.tax$dx)$dy add.forwd <- function(x,y){ return(x + y) } add.bckwd <- function(x, y, dout){ dx = dout dy = dout return(list(dx = dx, dy = dy)) } # network step.func <- function(x){ return(as.numeric(x > 0)) } x = seq(-5, 5, 0.1) y = step.func(x) plot(x,y, ylab = 'y', xlab = 'a', type ="l", lwd =2) sigmoid.func <- function(x){ return(1/(1+exp(-x))) } y = sigmoid.func(x) plot(x,y, ylab = 'y', xlab = 'a', type ="l", lwd =2) y.step = step.func(x) y.sigm = sigmoid.func(x) plot(x,y.step, ylab = 'y', xlab = 'a', type ="l", lwd =2) lines(x,y.sigm, lwd =2, lty = 2) relu.func <- function(x){ return(pmax(0,x)) } y.relu = relu.func(x) plot(x,y.relu, ylab = 'y', xlab = 'a', type ="l", lwd =2) A = matrix(1:4, nrow = 2, byrow = T) B = matrix(5:8, nrow = 2, byrow = T) A = matrix(1:6, nrow = 3, byrow = T) B = matrix(7:8, nrow = 2, byrow = T) x = c(1,0.5) W1 = matrix((1:6)*0.1, nrow = 2) B1 = (1:3)*0.1 A1 = x%*%W1 + B1 Z1 = sigmoid.func(A1) W2 = matrix((1:6)*0.1, nrow = 3) B2 = c(0.1, 0.2) A2 = Z1%*%W2 + B2 Z2 = sigmoid.func(A2) W3 = matrix((1:4)*0.1, nrow = 2) B3 = c(0.1, 0.2) A3 = Z2%*%W3+ B3 Z3 = A3 # function to initialize 3L network init.3L.network <- function(){ W1 = matrix((1:6)*0.1, nrow = 2) B1 = (1:3)*0.1 W2 = matrix((1:6)*0.1, nrow = 3) B2 = c(0.1, 0.2) W3 = matrix((1:4)*0.1, nrow = 2) B3 = c(0.1, 0.2) return(list(W1 = W1, B1 = B1, W2 = W2, B2 = B2, W3 = W3, B3 = B3)) } # feedforward process forward.3L <- function(network, x){ A1 = x%*%network$W1 + network$B1 Z1 = sigmoid.func(A1) A2 = Z1%*%network$W2 + network$B2 Z2 = sigmoid.func(A2) A3 = Z2%*%network$W3 + network$B3 Z3 = sigmoid.func(A3) A3 = Z3 return(A3) } network<-init.3L.network() y = forward.3L(network, c(1, 0.5)) a = c(1010,1000,990) exp(a)/sum(exp(a)) softmax.func <- function(x){ max.x = max(x) return(exp(x-max.x)/sum(exp(x-max.x))) } train <- read.csv('http://www.Matsuka.info/univ/course_folder/MNSTtrain.csv', header=TRUE) train <- data.matrix(train) train.x <- train[,-1] train.y <- train[,1] train.x <- t(train.x/255) download.file("http://www.matsuka.info/univ/course_folder/trNetwork.Rdata", "trNetwork.Rdata") load("trNetwork.Rdata") network=trNetwork n.train = ncol(train.x) correct.cl = 0 conf.matrix = matrix(0,10,10) for (i.loop in 1:n.train){ y = forward.3L(network,train.x[,i.loop]) max.y = max.col(y) conf.matrix[max.y, (train.y[i.loop]+1)] = conf.matrix[max.y, (train.y[i.loop]+1)] + 1 } accuracy = sum(diag(conf.matrix))/n.train # learning apple = 100; n.apple = 2; tax = 1.1 orange = 150; n.orange = 3; apple.price = multi.forwd(apple, n.apple) orange.price = multi.forwd(orange, n.orange) all.price = add.forwd(apple.price, orange.price) price = multi.forwd(all.price, tax) dprice = 1 d.all.price = multi.bckwd(all.price, tax, dprice) d.apple.price = add.bckwd(apple.price, orange.price, d.all.price$dx)$dx d.orange.price = add.bckwd(orange, n.orange.price, d.all.price$dx)$dy d.apple = multi.bckwd(apple, n.apple, d.apple.price)$dx d.n.apple = multi.bckwd(apple, n.apple, d.apple.price)$dy d.orange = multi.bckwd(orange, n.orange, d.orange.price)$dx d.n.orange = multi.bckwd(orange, n.orange, d.orange.price)$dy relu.forwd <- function(x){ return(pmax(x,0)) } relu.bckwd <- function(x, dout){ dout[which(x <= 0)] = 0 return(dout) } sigmoid.forwd <- function(x){ return(1/(1+exp(-x))) } sigmoid.bckwd <- function(x, dout){ y = sigmoid.forwd(x) return(dout*(1-y)*y) } affine.forwd <- function(x, W, b){ return(x%*%W + matrix(1, nrow = nrow(x), ncol = 1)%*%b) } affine.bckwd <- function(x, W, b, dout){ dx = dout%*%t(W) dW = t(x)%*%dout db = colSums(dout) return(list(dx = dx, dW = dW, db = db)) } softmax.forwd <- function(x, target){ max.x = apply(x,1,max) C = ncol(x) x = x - max.x%*%matrix(1,nrow=1,ncol=C) y = exp(x)/rowSums(exp(x)) delta = 1e-7; R = nrow(as.matrix(y)) return(-sum(target*log(y + delta))/R) } softmax.bckwd <- function(x, target, dout = 1){ max.x = apply(x, 1, max) R = nrow(x) C = ncol(x) x = x - max.x%*%matrix(1,nrow=1,ncol=C) y = exp(x)/rowSums(exp(x)) return((y-target)/R) } init.network <- function(n.neurons){ n.layer = length(n.neurons) W = list(); b = list() for (i.layer in 1:(n.layer-1)){ W[[i.layer]] = matrix(rnorm(n.neurons[i.layer]*n.neurons[(i.layer+1)],sd = 0.1), nrow=n.neurons[i.layer]) b[[i.layer]] = matrix(rnorm(n.neurons[(i.layer+1)],sd = 0.1), nrow = 1) } return(list(W = W,b = b)) } sigmoid.func <- function(x){ return(1/(1+exp(-x))) } relu.func <- function(x){ y = apply(x,2,function(x) pmax(x,0)) return(y) } activation <- function(A, actFun){ if (actFun == "sigmoid"){ return(sigmoid.func(A)) } if (actFun == "relu"){ return(relu.func(A)) } if (actFun == "softmax"){ return(softmax(A)) } } feedforward <- function(network, x, actFun) { n.layer <- length(network$W) batch.size = nrow(x) for (i.layer in 1:n.layer){ A = x%*%network$W[[i.layer]] + matrix(1,nrow=batch.size,ncol = 1)%*%network$b[[i.layer]] x = activation(A, actFun[i.layer]) } return(x) } cross.entropy = function(y, target){ delta = 1e-7; R = nrow(as.matrix(y)) return(-sum(target*log(y + delta))/R) } loss.network = function(params, x, t, actFun){ y = feedforward(params,x,actFun) return(cross.entropy(y, t)) } softmax <- function(x, target){ max.x = apply(x,1,max) C = ncol(x) x = x - max.x%*%matrix(1,nrow=1,ncol=C) y = exp(x)/rowSums(exp(x)) return(y) } train.x = as.matrix(iris[,1:4]) train.y.temp = as.numeric(iris[,5]) train.y = matrix(0,nrow = nrow(train.x), ncol =3) train.y[which(train.y.temp==1), 1]=1 train.y[which(train.y.temp==2), 2]=1 train.y[which(train.y.temp==3), 3]=1 params = init.network(c(4,15,3)) batch_size = 10; n.iter =5000; lambda =0.05 n.train = nrow(train.x) params = init.network(c(4,30,3)) batch_size = 10; n.iter =5000; lambda =0.01 n.train = nrow(train.x) loss = rep(0,n.iter) for (i.iter in 1:n.iter){ batch_mask = sample(1:n.train, batch_size) x.batch = train.x[batch_mask,] t.batch = train.y[batch_mask,] a1 = affine.forwd(x.batch,params$W[[1]],params$b[[1]]) z1 = sigmoid.forwd(a1) a2 = affine.forwd(z1,params$W[[2]],params$b[[2]]) z2 = softmax.forwd(a2,t.batch) dwSM = softmax.bckwd(a2, t.batch, 1) dwA2 = affine.bckwd(a1,params$W[[2]],params$b[[2]],dwSM) dwSG = sigmoid.bckwd(a1,dwA2$dx) dwA1 = affine.bckwd(x.batch,params$W[[1]],params$b[[1]],dwSG) params$W[[2]] = params$W[[2]] - lambda*dwA2$dW params$b[[2]] = params$b[[2]] - lambda*dwA2$db params$W[[1]] = params$W[[1]] - lambda*dwA1$dW params$b[[1]] = params$b[[1]] - lambda*dwA1$db loss[i.iter] = loss.network(params,x.batch,t.batch,c("sigmoid","softmax")) } plot(loss,type='l', xlab = "trial")
広域システム特別講義II 強化学習1B
mk_MC_seq <- function(){ # defining transitino matrix & probs. 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) P=rep(0.25,4) # creating sequence 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] } return(state.seq = state.seq) } # update Vs cal.cumV.MC <- function(cumV, state.seq,state.count){ r = -1 uniq.seq = unique(state.seq) for (i.uniq in 1:(length(uniq.seq)-1)){ first.visit = which(state.seq == uniq.seq[i.uniq])[1] cumV[uniq.seq[i.uniq]] = cumV[uniq.seq[i.uniq]] + r*(length(state.seq)-first.visit) } state.count[uniq.seq] = state.count[uniq.seq] + 1 return(list(cumV = cumV, state.count = state.count)) } # main script max.iter = 10000; state.count=rep(0,15); cumV = rep(0,14) for (i.iter in 1:max.iter){ state.seq = mk_MC_seq() updates = cal.cumV.MC(cumV, state.seq, state.count) state.count = updates$state.count cumV = updates$cumV } V = matrix(c(0,cumV/state.count[1:14],0),nrow=4) # function to calc. card values card.value<-function(adj.cards) { sum.cards=sum(adj.cards) if (any(adj.cards==1) & sum.cards<=11) { sum.cards=sum.cards+10; usableA=1 #true } else {usableA=2} #false return(c(sum.cards,usableA)) } # function to calc. reward calc.reward<-function(p.val,d.val) { if (p.val>21) { reward=-1 } else { if (d.val>21) { reward=1 } else { if (p.val==d.val) { reward=0 } else { reward=ifelse(p.val>d.val,1,-1) } } } } # playing a single game play.BJ <- function(policy){ cards=sample(rep(1:13,4)) player=cards[1:2]; adj.player=pmin(player,10) dealer=cards[3:4]; adj.dealer=pmin(dealer,10) cards=cards[-(1:4)] d.val=card.value(adj.dealer) p.val=card.value(adj.player) state.hist=c(adj.dealer[1],p.val[1],p.val[2]) while (p.val[1] < policy) { player=c(player,cards[1]); adj.player=pmin(player,10) cards=cards[-1] p.val=card.value(adj.player) state.hist=rbind(state.hist,c(adj.dealer[1],p.val[1],p.val[2])) } while (d.val[1] < 17) { dealer=c(dealer,cards[1]); adj.dealer=pmin(dealer,10) cards=cards[-1] d.val=card.value(adj.dealer) } return(list(p.val = p.val, d.val = d.val, state.hist = state.hist)) } # main function BJ_MC_fixedPolicy<-function(policy=20,maxIter=1e6){ rew.sum=array(0,dim=c(10,10,2)) rew.count=array(0,dim=c(10,10,2)) for (i_play in 1:maxIter) { result <- play.BJ(policy) p.val = result$p.val d.val = result$d.val state.hist = result$state.hist rew=calc.reward(p.val[1],d.val[1]) n.state=nrow(state.hist) if (is.null(n.state)) { n.state=1 state.hist=t(as.matrix(state.hist)) } for (i_state in 1:n.state) { if (state.hist[i_state,2] > 11 & state.hist[i_state,2] < 22) { rew.sum[state.hist[i_state,1],(state.hist[i_state,2]-11),state.hist[i_state,3]]= rew.sum[state.hist[i_state,1],(state.hist[i_state,2]-11),state.hist[i_state,3]]+rew rew.count[state.hist[i_state,1],(state.hist[i_state,2]-11),state.hist[i_state,3]]=rew.count[state.hist[i_state,1],(state.hist[i_state,2]-11),state.hist[i_state,3]]+1 } } } return(rew.sum/rew.count) } play.BJ2 <- function(policy){ cards=sample(c(rep(1:10,4),rep(10,12))) player=cards[1:2]; dealer=cards[3:4]; cards=cards[-(1:4)] d.val=card.value(dealer); p.val=card.value(player) while( p.val[1] < 12 ) { player=c(player,cards[1]) cards=cards[-1] p.val=card.value(player) } # initial random action action=sample(0:1,1) state.hist=c(dealer[1],p.val[1],p.val[2],(action+1)) # player's action while (action==1 & p.val[1]<22) { player=c(player,cards[1]) cards=cards[-1] p.val=card.value(player) state.hist=rbind(state.hist,c(dealer[1],p.val[1],p.val[2],(action+1))) if (p.val[1]<22) { action=policy[dealer[1],(p.val[1]-11),p.val[2]] } } # dealer's action while (d.val[1]<17) { dealer=c(dealer,cards[1]) cards=cards[-1] d.val=card.value(dealer) } return(list(p.val = p.val, d.val = d.val, state.hist = state.hist)) } BJ_MC<-function(maxIter=1e6){ rew.sum=array(0,dim=c(10,10,2,2)) rew.count=array(1,dim=c(10,10,2,2)) Q=array(0,dim=c(10,10,2)) V=array(0,dim=c(10,10,2)) policy=array(sample(0:1,10*10*2,replace=T),dim=c(10,10,2)) # policy: 1 = hit, 0 = stay for (i_play in 1:maxIter) { result = play.BJ2(policy) p.val = result$p.val d.val = result$d.val state.hist = result$state.hist rew=calc.reward(p.val[1],d.val[1]) n.state=nrow(state.hist) if (is.null(n.state)) { n.state=1 state.hist=t(as.matrix(state.hist)) } for (i_state in 1:n.state) { if (state.hist[i_state,2]>11 & state.hist[i_state,2]<22) { ind=state.hist[i_state,]-c(0,11,0,0) rew.sum[ind[1],ind[2],ind[3],ind[4]]= rew.sum[ind[1],ind[2],ind[3],ind[4]]+rew rew.count[ind[1],ind[2],ind[3],ind[4]]=rew.count[ind[1],ind[2],ind[3],ind[4]]+1 Q=rew.sum/rew.count; policy[,,1]=Q[,,1,1] < Q[,,1,2] policy[,,2]=Q[,,2,1] < Q[,,2,2] } } } V[,,1]=(rew.sum[,,1,1]+rew.sum[,,1,2])/(rew.count[,,1,1]+rew.count[,,1,2]) V[,,2]=(rew.sum[,,2,1]+rew.sum[,,2,2])/(rew.count[,,2,1]+rew.count[,,2,2]) return(list(policy,V,Q)) } # TD random walk TD0.ex1<-function(maxItr,alpha,gamma) { V=c(0,rep(0.5,5),0) V.hist=matrix(0,nrow=maxItr+1,ncol=5) V.hist[1,]=V[2:6] P.act=matrix(0.5,ncol=7,nrow=2) for (i_rep in 1:maxItr) { state=5 while (state!=1 & state!=7) { action=sample(c(-1,1),1,prob=P.act[,state]) state.old=state state=state+action r=ifelse(state==7,1,0) V[state.old]=V[state.old]+alpha*(r+gamma*V[state]-V[state.old]) } V.hist[(i_rep+1),]=V[2:6] } return(V.hist) } # constant step-size Monte Carlo constMC.ex1<-function(maxItr,alpha) { V=c(0,rep(0.5,5),0) V.hist=matrix(0,nrow=maxItr+1,5) V.hist[1,]=V[2:6] P.act=matrix(0.5,ncol=7,nrow=2) for (i_rep in 1:maxItr) { state=5; state.hist=state while (state!=1 & state!=7) { action=sample(c(-1,1),1,prob=P.act[,state]) state=state+action state.hist=cbind(state.hist,state) } R=ifelse(state==7,1,0) n.state=length(state.hist) for (i_state in 1:(n.state-1)) { V[state.hist[i_state]]=V[state.hist[i_state]]+ alpha*(R-V[state.hist[i_state]]) } V.hist[(i_rep+1),]=V[2:6] } return(V.hist) } # (re)creating Fig 6.7 alphaTD=c(0.05,0.075,0.1,0.15) alphaMC=c(0.01,0.02,0.03,0.04) n.alphas=length(alphaTD) pchs=0:(0+n.alphas) true.V=1:5*(1/6) n_rep=100 sqs=seq(1,101,2) plot(0,0,type='n',xlim=c(0,100),ylim=c(0,0.25)) for (i_alpha in 1:n.alphas) { rmsTD=matrix(0,101,n_rep) rmsMC=matrix(0,101,n_rep) for (i_rep in 1:n_rep) { resTD=TD0.ex1(100,alphaTD[i_alpha],1) resMC=constMC.ex1(100,alphaTD[i_alpha]) for (i_gen in 1:101) { rmsTD[i_gen,i_rep]=sqrt(mean((resTD[i_gen,]-true.V)^2)) rmsMC[i_gen,i_rep]=sqrt(mean((resMC[i_gen,]-true.V)^2)) } } mTD=rowMeans(rmsTD) mMC=rowMeans(rmsMC) lines(mTD,col='red') lines(mMC,col='blue') lines(sqs,mTD[sqs],col='red',pch=pchs[i_alpha],type='p') lines(sqs,mMC[sqs],col='blue',pch=pchs[i_alpha],type='p') } labs=c("MC, alpha=0.01", "MC, alpha=0.02", "MC, alpha=0.03", "MC, alpha=0.04", "TD, alpha=0.05", "TD, alpha=0.075", "TD, alpha=0.10", "TD, alpha=0.15") legend('topright',labs,col=c(rep('blue',4),rep('red',4)),pch=rep(0:3,2),lwd=1.5) sarsa.ex6.5<-function(maxItr,alpha,gamma,epsilon) { # field size: 7row x 10column # horizontal move -> COLUMN # vertical move -> ROW # effect of wind -> ROW # actions: 1-up, 2-right, 3-down, 4-left act.V=matrix(c(1,0,0,1,-1,0,0,-1),nrow=4,byrow=T) wind=matrix(c(0,0,0,0,0,0,1,0,1,0,1,0,2,0,2,0,1,0,0,0),byrow=T,nrow=10) goal=c(4,8) Qs=array(0,dim=c(7,10,4)) for (i_rep in 1:maxItr) { state=c(4,1) # start if (runif(1) > epsilon) { move=which.max(Qs[state[1],state[2],]) } else { move=sample(1:4,1)} while (!all(state==goal)) { st.old=state mv.old=move state=state+act.V[move,]+wind[state[2],] if (state[1]<1) {state[1]=1} if (state[1]>7) {state[1]=7} if (state[2]<1) {state[2]=1} if (state[2]>10) {state[2]=10} if (runif(1) > epsilon) { move=which.max(Qs[state[1],state[2],]) } else { move=sample(1:4,1)} rew=ifelse(all(state==goal),0,-1) Qs[st.old[1],st.old[2],mv.old]=Qs[st.old[1],st.old[2],mv.old] +alpha*(rew+gamma* Qs[state[1],state[2],move] -Qs[st.old[1],st.old[2],mv.old]) } } return(Qs) } # running example Qs=sarsa.ex6.5(5e6,0.1,1,0.1) # sim optimal actions state=c(4,1);goal=c(4,8); state.hist=state while (!all(state==goal)) { moveID=which.max(Qs[state[1],state[2],]) state=state+act.V[moveID,]+wind[state[2],] if (state[1]<1) {state[1]=1} if (state[1]>7) {state[1]=7} if (state[2]<1) {state[2]=1} if (state[2]>10) {state[2]=10} state.hist=rbind(state.hist,state) } # plotting results plot(0,0,type='n',xlim=c(0,11),ylim=c(0,8),xlab="",ylab="", main="Learned policies -- Sarsa") lines(1,4,type='p',pch=19,col='red',cex=2) lines(8,4,type='p',pch=19,col='red',cex=2) dirs=c("up","right","down","left" ) for (i_row in 1:7) { for (i_col in 1:10) { best.move=dirs[which.max(Qs[i_row,i_col,])] text(i_col,i_row,best.move) } } lines(state.hist[,2],state.hist[,1],col="red",lwd=2) Qlearn.ex6.5<-function(maxItr,alpha,gamma,epsilon) { # field size: 7row x 10column # horizontal move -> COLUMN # vertical move -> ROW # effect of wind -> ROW # actions: 1-up, 2-right, 3-down, 4-left act.V=matrix(c(1,0,0,1,-1,0,0,-1),nrow=4,byrow=T) wind=matrix(c(0,0,0,0,0,0,1,0,1,0,1,0,2,0,2,0,1,0,0,0),byrow=T,nrow=10) goal=c(4,8) Qs=array(0,dim=c(7,10,4)) for (i_rep in 1:maxItr) { state=c(4,1) # start while (!all(state==goal)) { if (runif(1) > epsilon) { move=which.max(Qs[state[1],state[2],]) } else { move=sample(1:4,1)} sIDX=state state=state+act.V[move,]+wind[state[2],] if (state[1]<1) {state[1]=1} if (state[1]>7) {state[1]=7} if (state[2]<1) {state[2]=1} if (state[2]>10) {state[2]=10} max.Q=max(Qs[state[1],state[2],]) rew=ifelse(all(state==goal),0,-1) Qs[sIDX[1],sIDX[2],move]=Qs[sIDX[1],sIDX[2],move] +alpha*(rew+gamma* max.Q-Qs[sIDX[1],sIDX[2],move]) } } return(Qs) } Qs=Qlearn.ex6.5(1e6,0.05,1,0.1) # sim optimal actions state=c(4,1);goal=c(4,8); state.hist=state while (!all(state==goal)) { moveID=which.max(Qs[state[1],state[2],]) state=state+act.V[moveID,]+wind[state[2],] if (state[1]<1) {state[1]=1} if (state[1]>7) {state[1]=7} if (state[2]<1) {state[2]=1} if (state[2]>10) {state[2]=10} state.hist=rbind(state.hist,state) } # plotting results plot(0,0,type='n',xlim=c(0,11),ylim=c(0,8),xlab="",ylab="", main="Learned policies -- Q-learning") lines(1,4,type='p',pch=19,col='red',cex=2) lines(8,4,type='p',pch=19,col='red',cex=2) dirs=c("up","right","down","left" ) for (i_row in 1:7) { for (i_col in 1:10) { best.move=dirs[which.max(Qs[i_row,i_col,])] text(i_col,i_row,best.move) } } lines(state.hist[,2],state.hist[,1],col="red",lwd=2)
データ解析基礎論B W05 Factor Analysis
chisq.test(c(72,23,16,49),p=rep(40,4),rescale.p=F) chisq.test(c(72,23,16,49),p=rep(0.25,4),rescale.p=F) M=matrix(c(52,48,8,42),nrow=2) chisq.test(M,correct=T) #(abs(52-40)-0.5)^2/40+(abs(48-60)-0.5)^2/60 # +(abs(8-20)-0.5)^2/20+(abs(42-30)-0.5)^2/30 dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/FA01.csv") dat.fa<-factanal(dat,1) dat<-read.table("http://www.matsuka.info/data_folder/tdkPCA01.txt") dat.pca<-princomp(dat) dat.fa<-factanal(dat,1) dat.fa<-factanal(dat,1,score="regression") plot(dat.fa$score~dat.pca$score[,1],pch=20,cex=2,xlab="Component Score", ylab="Factor Score") fa_pca.scores = tibble(fa = dat.fa$scores, pca = dat.pca$scores[,1], total.score = rowSums(dat)) ggplot(fa_pca.scores) + geom_point(aes(x = fa, y = pca), size = 3) + xlab("Factor Score") + ylab("Component Score") cor(dat.fa$score,dat.pca$score) ggplot(fa_pca.scores) + geom_point(aes(x = fa, y = total.score), size = 3) + xlab("Factor Score") + ylab("Total Score") dat<-read.csv("http://www.matsuka.info/data_folder/tdkCFA.csv") dat.faWOR<-factanal(dat,2, rotation="none", score="regression") dat.faWR<-factanal(dat,2, rotation="varimax", score="regression") loadingsWOR <- dat.faWOR$loadings[] %>% as.tibble() %>% add_column(variable = row.names(dat.faWOR$loadings)) %>% gather("Factor1","Factor2", key = "factor", value = "loadings") loadingsWR <- dat.faWR$loadings[] %>% as.tibble() %>% add_column(variable = row.names(dat.faWR$loadings)) %>% gather("Factor1","Factor2", key = "factor", value = "loadings") ggplot(loadingsWOR, aes(variable, abs(loadings), fill=loadings)) + facet_wrap(~ factor, nrow=1) + geom_bar(stat="identity") + coord_flip() + scale_fill_gradient2(name = "loadings", high = "blue", mid = "white", low = "red", midpoint=0, guide=F) + ylab("Loading Strength") ggplot(loadingsWR, aes(variable, abs(loadings), fill=loadings)) + facet_wrap(~ factor, nrow=1) + geom_bar(stat="identity") + coord_flip() + scale_fill_gradient2(name = "loadings", high = "blue", mid = "white", low = "red", midpoint=0, guide=F) + ylab("Loading Strength") loadingsWR2 <- as.data.frame(dat.faWR$loadings[]) ggplot(loadingsWR2, aes(x = Factor1, y = Factor2)) + geom_point(size = 3, color = "red") + geom_vline(xintercept=0) + geom_hline(yintercept=0) + geom_text(aes(label = rownames(loadingsWR2))) + ylim(-1.1, 1.1) + xlim(-1.1, 1.1) dat.model1<-factanal(dat,1) dat.model2<-factanal(dat,2) dat.model3<-factanal(dat,3) dat.model4<-factanal(dat,4) source("http://www.matsuka.info/univ/course_folder/cuUtil02.R") cu.lrtest.csq(dat.model3,dat.model4) cu.AIC.csq(dat.model1) library(sem) model01=cfa(reference.indicator=FALSE) F1:extrovert,cheerful, leadership, antisocial, talkative, motivated, hesitance, popularity cv.mat = cov(dat) mod1<-sem(model01,cv.mat,100) model02=cfa(reference.indicator=FALSE) F1: extrovert, leadership, motivated, hesitance F2: cheerful, antisocial, talkative, popularity mod2<-sem(model02, cov(dat), nrow(dat)) opt <- options(fit.indices = c("RMSEA")) summary(mod2)