# 認知情報解析　課題２

```# initializing Q matrix
Q = P = matrix(1/4,nrow=25,ncol=4) #

# defining deterministic transition matrix
north=c(2:25,25)
north[ c(5,10,15,20,25)]=c(5,10,15,20,25)
east=c(6:25,21:25)
west=c(1:5,1:20)
south=c(1,1:24)
south[ c(1,6,11,16,21)]=c(1,6,11,16,21)
trM=cbind(north,east,south,west)
trM[10,]=6
trM[20,]=18

R=matrix(0,nrow=25,ncol=4)
R[which(trM==1:25)]=-1
R[10,]=10
R[20,]=5

nRep=1000; gamma=0.9; P = 0.25
for (i_rep in 1:nRep) {
Q.old = Q
for (i_state in 1:25) {
for (i_act in 1:4){
Q[i_state, i_act]=R[i_state, i_act]+gamma * P * sum(Q.old[trM[i_state,i_act]])
}
}
}
```

# 2019年度　データ解析基礎論B W06

```x.temp = 0:9
mass= matrix(dbinom(x.temp,9,0.5),nrow=1)
colnames(mass) <- paste(0:9)
barplot(mass)

sum(dbinom(7:9,9,0.5))
2*(0.5-pnorm(750,800,40))
pnorm(700,800,40)
qnorm(c(0.025, 0.975),800,40)
qnorm(0.99,800,40)

zA=(165-150)/(sqrt(15^2/10))
1-pnorm(zA)
(1-pnorm(zA))*2

mean.M <-mean(dat\$h[dat\$gender=="M"])
sigma = 10
n.M = length(dat\$h[dat\$gender=="M"])
z.value=(mean.M-171)/(sqrt(sigma))

ssize = c(24,25,26,23.5,25,27,24,22,27.5,28)
ssize.mean = mean(ssize)
ssize.var = var(ssize)
N = 10
t.value=(ssize.mean-24)/(sqrt(ssize.var/N))
(1-pt(abs(t.value),df=9))*2
t.test(ssize, mu=24)
```

# 2019 データ解析基礎論A DAA04

```dat<- read.csv("http://www.matsuka.info/data_folder/datWA01.txt")
plot(dat\$shoesize, dat\$h, main="Relationship b/w shoesize and height",
xlab = "shoesize", ylab="height", pch=19, col="red")
txt = paste("r =", round(cor(dat\$shoesize,dat\$h), 4))
text(22, 175, txt, cex = 1.5)
abline(h = mean(dat\$h), col='blue');
abline(v = mean(dat\$shoesize), col='darkgreen');
text(x = 21, y = mean(dat\$h)+3,
paste("mean height =", round(mean(dat\$h),2)),

text(x = mean(dat\$shoesize)+0.2, y = 145,
paste("mean shoesize =", round(mean(dat\$shoesize),2)),

abline(lm(dat\$h~dat\$shoesize), lty=2, lwd=2)

plot(dat[dat\$gender=='F',]\$shoesize, dat[dat\$gender=='F',]\$h,
main="Relationship b/w shoesize and height", xlab='shoesize', ylab='height',
cex.lab=1.5, pch=19, col='blue', xlim=c(20,29), ylim=c(140,190))
plot(dat[dat\$gender=='M',]\$shoesize, dat[dat\$gender=='M',]\$h,
main="Relationship b/w shoesize and height", xlab='shoesize', ylab='height',
cex.lab=1.5, pch=15, col='green', xlim=c(20,29), ylim=c(140,190))

plot(dat[dat\$gender=='F',]\$shoesize, dat[dat\$gender=='F',]\$h,
main="Relationship b/w shoesize and height", xlab='shoesize', ylab='height',
cex.lab=1.5, pch=19, col='blue', xlim=c(20,29), ylim=c(140,190))
points(dat[dat\$gender=='M',]\$shoesize,dat[dat\$gender=='M',]\$h,
pch = 15, col = 'green')
legend("topleft", c('Female','Male'), pch =c(19,15),
col = c('blue','green'), cex = 1.5)

plot(dat, pch=20, col='blue')

plot(dat.pca, pch = rownames(dat.pca), cex = 1.7, col = "blue")

summary(dat)
dat.meter = dat[,1:2]*0.01
dat.h.ext5 = dat\$h+5
hANDshoe = dat\$h+dat\$shoesize
dat.h.meter.ext5 = dat\$h*0.01+0.05

```

# 院：認知情報解析学

```source("http://www.matsuka.info/univ/course_folder/cuUtil02.R")
dat.fa1 <- factanal(dat,1)
dat.fa2 <- factanal(dat,2)
dat.fa3 <- factanal(dat,3)
dat.fa4 <- factanal(dat,4)

install.packages('sem')
library(sem)

model01=cfa(reference.indicator=FALSE)
F1:extrovert,cheerful, leadership, antisocial, talkative, motivated, hesitance, popularity

mod1<-sem(model01, data = dat)

model02=cfa(reference.indicator=FALSE)
F1: extrovert, leadership, motivated, hesitance
F2: cheerful, antisocial, talkative, popularity

mod2<-sem(model02, data = dat)

```

# 強化学習　方策の比較1

```# Qが最大のactionを選択
max.Q <- function(Q){
max.a = max(Q)
max.idx = which(Q == max.a)
if (length(max.idx)>1){
max.idx = sample(max.idx, 1)
}
return(max.idx)
}

# greedy方策
greedy <- function(n.trial, Q.star, N){
Q = Q.cum = count = rep(0, N)
rew.earned = rep(0, n.trial)
for (i.trial in 1:n.trial){
act.idx = max.Q(Q)
r.t = rnorm(1, mean = Q.star[act.idx], sd = 1)
Q.cum[act.idx] = Q.cum[act.idx] + r.t
count[act.idx] = count[act.idx] + 1
Q[act.idx] = Q.cum[act.idx] / count[act.idx]
rew.earned[i.trial] = r.t
}
return(list(Q = Q, rew.earned = rew.earned))
}

# epsilon greedy方策
# epsilon = 0の場合はgreedy方策と同等
eps.greedy <- function(n.trial, Q.star, N, epsilon){
Q = Q.cum = count = rep(0, N)
rew.earned = rep(0, n.trial)
for (i.trial in 1:n.trial){
if (runif(1) < epsilon){
act.idx = sample(1:N, 1)
} else {
act.idx = max.Q(Q)
}
r.t = rnorm(1, mean = Q.star[act.idx], sd = 1)
Q.cum[act.idx] = Q.cum[act.idx] + r.t
count[act.idx] = count[act.idx] + 1
Q[act.idx] = Q.cum[act.idx] / count[act.idx]
rew.earned[i.trial] = r.t
}
return(list(Q = Q, rew.earned = rew.earned))
}

# n.rep回繰り返す関数
comp.eps.greedy <- function(n.trial, n.rep, N, epsilon){
rew.history = matrix(0, nrow = n.trial, ncol = n.rep)
for (i.rep in 1:n.rep){
Q.star = rnorm(N, mean = 0, sd = 1);
res = eps.greedy(n.trial, Q.star, N, epsilon)
rew.history[ , i.rep] = res\$rew.earned
}
return(rew.history)
}

# 実行
EG.000 = comp.eps.greedy(1000, 1000, 10, 0.000)
EG.001 = comp.eps.greedy(1000, 1000, 10, 0.001)
EG.010 = comp.eps.greedy(1000, 1000, 10, 0.010)
EG.100 = comp.eps.greedy(1000, 1000, 10, 0.100)
EG.150 = comp.eps.greedy(1000, 1000, 10, 0.150)

# 結果の可視化
plot(rowMeans(EG.000), type="l", ylab="Average Reward", xlab="Trial",
ylim = c(0,2))
lines(rowMeans(EG.001),col=2)
lines(rowMeans(EG.010),col=3)
lines(rowMeans(EG.100),col=4)
lines(rowMeans(EG.150),col=5)
legend("bottomright",
c("epsilon = 0.000",
"epsilon = 0.001",
"epsilon = 0.010",
"epsilon = 0.100",
"epsilon = 0.150"),
col=1:5,lwd=2 )

# softmax
softmax <- function(n.trial, Q.star, N, tau){
Q = Q.cum = count = rep(0, N)
rew.earned = rep(0, n.trial)
for (i.trial in 1:n.trial){
p = exp(Q*tau)/sum(exp(Q*tau))
act.idx = sample(1:N, 1, prob = p)
r.t = rnorm(1, mean = Q.star[act.idx], sd = 1)
Q.cum[act.idx] = Q.cum[act.idx] + r.t
count[act.idx] = count[act.idx] + 1
Q[act.idx] = Q.cum[act.idx] / count[act.idx]
rew.earned[i.trial] = r.t
}
return(list(Q = Q, rew.earned = rew.earned))
}

comp.softmax <- function(n.trial, n.rep, N, tau){
rew.history = matrix(0, nrow = n.trial, ncol = n.rep)
for (i.rep in 1:n.rep){
Q.star = rnorm(N, mean = 0, sd = 1);
res = softmax(n.trial, Q.star, N, tau)
rew.history[ , i.rep] = res\$rew.earned
}
return(rew.history)
}

# 実行
EG.000 = comp.eps.greedy(1000, 1000, 10, 0.000)
EG.010 = comp.eps.greedy(1000, 1000, 10, 0.010)
EG.100 = comp.eps.greedy(1000, 1000, 10, 0.100)
SM.10 = comp.softmax(1000,1000,10,10)
SM.03 = comp.softmax(1000,1000,10,3)
# 結果の可視化
plot(rowMeans(EG.000), type="l", ylab="Average Reward", xlab="Trial",
ylim = c(0,2))
lines(rowMeans(EG.010),col=2)
lines(rowMeans(EG.100),col=3)
lines(rowMeans(SM.10),col=4)
lines(rowMeans(SM.03),col=5)
legend("bottomright",
c("epsilon = 0.000",
"epsilon = 0.010",
"epsilon = 0.100",
"tau = 10",
"tau = 3"),
col=1:5,lwd=2 )

# epsilon greedy (2nd version)
eps.greedy2 <- function(n.trial, Q.star, N, epsilon, lr, init.Q){
Q = rep(init.Q, N)
rew.earned = rep(0, n.trial)
for (i.trial in 1:n.trial){
if (runif(1) < epsilon){
act.idx = sample(1:N, 1)
} else {
act.idx = max.Q(Q)
}
r.t = rnorm(1, mean = Q.star[act.idx], sd = 1)
Q[act.idx] = Q[act.idx] + lr*(r.t - Q[act.idx])
rew.earned[i.trial] = r.t
}
return(list(Q = Q, rew.earned = rew.earned))
}
```

# 2019データ解析基礎論A DAA03可視化２

```v1 = seq(-3,3,0.1)
v2 = v1^2
plot(v1, v2, col = "blue", type = "o", lty = 2, pch = 19,
cex.lab = 1.5, lwd = 3, main = "Y=X*X", xlab = "X",
ylab="X*X", xlim=c(-3.5,3.5), ylim=c(-0.5, 10))

# histogram
hist(dat\$h)
hist(dat\$h, breaks = 20, main = "Histogram of Height",
xlab = "Height", col = 'blue', xlim = c(140, 190))
dens<-density(dat\$h);
hist(dat\$h, main = "Histogram of Height", xlab = "Height",
xlim = c(140,190), probability = T)
lines(dens, lwd = 2, col = "red", lty=2)
plot(v1, v2, col = "blue", type = "l", pch = 19, cex.lab = 1.5,
lwd = 3, xlab = "X", ylab="f(X)", xlim=c(-3.5,3.5), ylim=c(-0.5, 10))
lines(v1, v1^3, col='red',lwd = 3)
legend("bottomright", c("x^2","x^3"), col=c('blue','red'), lwd=2)
legend(-3,0.5, c("x^2","x^3"), col=c('blue','red'), lwd=2)

boxplot(dat\$h,main="Boxplot of Height", ylab="Height", col='cyan', ylim=c(140,190))
boxplot(dat\$h,main="Boxplot of Height", xlab="Height", col="orange", horizontal=T)

boxplot(dat\$h ~ dat\$gender,
main="Distribution of Height by Gender",
ylab="Gender", xlab="Height", col=c('blue','cyan'),
ylim=c(140,190), horizontal=T)

boxplot(dat\$h ~ dat\$gender + dat\$affil,
main="Distribution of  Height by Gender and Affiliation",
ylab="Gender x Affiliation", xlab="Height",
col=c('blue','cyan','red','magenta'), ylim=c(140,190),horizontal=T)

interaction.plot(dat\$gender,
dat\$affil,
dat\$h,
pch=c(20,20),
col=c("skyblue","orange"),
xlab="gender", ylab="height",
lwd=3,type='b',cex=2,
trace.label="Affiliation")

par(mfrow=c(1,2))
hist(dat[dat\$gender=="F",]\$h, main="Dist. of Height for Female Participants",
xlab="Height", xlim=c(140,190), probability=T)

dens.F = density(dat[dat\$gender=='F',]\$h)

lines(dens.F, col='blue',lwd=2)

#２つ目のhistogram
hist(dat[dat\$gender=="M",]\$h, main="Dist. of Height for Male
Participants", xlab="Height", xlim=c(140,190), probability=T,ylim=c(0,0.08))

dens.M = density(dat[dat\$gender=='M',]\$h)

lines(dens.M, col='green', lwd=2)

par(mfrow=c(1,1))
plot(dens.F,col='blue',lwd=2, ylab='density', xlim=c(140,190),
main="Dist. of Height by gender",xlab='Height')
lines(dens.M,col='green',lwd=2)
legend("topleft", c('Female','Male'), col=c('blue','green'), cex=1.5,lwd=2)
text(157.5, 0.04, 'Female', col='blue', cex=2)
text(170, 0.04,'Male', col='green', cex=2)
plot(dat\$shoesize, dat\$h, main="Relationship b/w shoesize and height",
xlab = "shoesize", ylab="height", pch=19, col="red")

txt = paste("r =", round(cor(dat\$shoesize,dat\$h), 3))
text(22, 175, txt, cex = 1.5)

plot(dat[dat\$gender=='F',]\$shoesize, dat[dat\$gender=='F',]\$h,
main="Relationship b/w shoesize and height", xlab='shoesize', ylab='height',
cex.lab=1.5, pch=19, col='blue', xlim=c(20,29), ylim=c(140,190))
lines(dat[dat\$gender=='M',]\$shoesize,dat[dat\$gender=='M',]\$h,
type = 'p', pch = 15, col = 'green')
legend("topleft", c('Female','Male'), pch =c(19,15),
col = c('blue','green'), cex = 1.5)
```