# データ解析基礎論B GLM

```library(tidyverse)
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.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.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.glmAllMult=glm(survival~age*Ncigarettes*NdaysGESTATION,family=binomial,data=dat)
library(MASS)
stepAIC(dat.glmAllMult)

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

# データ解析基礎論B 多変量解析

```install.packages("ggfortify")
install.packages("ggdendro")
library(ggfortify)
library(ggdendro)

# pca
dat.pca<-princomp(dat)
autoplot(dat.pca, label = TRUE, label.size = 6,
autoplot(dat.pca, shape = FALSE, label.size = 6,

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

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,

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.lda<-lda(class~.,dat)
lda.pred<-predict(dat.lda,dat)
table(lda.pred\$class, dat\$class)
dat.ldaCV<-lda(class~.,dat, CV=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

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 <- data.matrix(train)
train.x <- train[,-1]
train.y <- train[,1]
train.x <- t(train.x/255)
"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)
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")

```
Posted in UT

# 広域システム特別講義II 強化学習１B

```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)
}
V = matrix(c(0,cumV/state.count[1:14],0),nrow=4)

# function to calc. card values
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))
cards=cards[-(1:4)]
while (p.val[1] < policy) {
cards=cards[-1]
}
while (d.val[1] < 17) {
cards=cards[-1]
}
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)
```
Posted in UT

# データ解析基礎論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.fa<-factanal(dat,1)

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.faWOR<-factanal(dat,2, rotation="none", score="regression")
dat.faWR<-factanal(dat,2, rotation="varimax", score="regression")

gather("Factor1","Factor2", key = "factor", value = "loadings")

gather("Factor1","Factor2", key = "factor", value = "loadings")

facet_wrap(~ factor, nrow=1) +
geom_bar(stat="identity") +
coord_flip() +
high = "blue", mid = "white", low = "red",
midpoint=0, guide=F) +

facet_wrap(~ factor, nrow=1) +
geom_bar(stat="identity") +
coord_flip() +
high = "blue", mid = "white", low = "red",
midpoint=0, guide=F) +

ggplot(loadingsWR2, aes(x = Factor1, y = Factor2)) +
geom_point(size = 3, color = "red") +
geom_vline(xintercept=0) +
geom_hline(yintercept=0) +
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)
```