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

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) 

library(tidyverse)
library(gridExtra)

epGreedy = function(nTrial,nRep,epsilon) {
  rew.hist = opt.hist = matrix(0,nrow=nTrial,ncol=nRep)
  for (i_rep in 1:nRep){
    Q.true=rnorm(10); opt.ID=which.max(Q.true)
    Q.est = Q.cum = act.count=rep(0,10);
    for (i_trial in 1:nTrial) {
      if (runif(1) < epsilon) {
        action=sample(1:10,1)
      } else {
        action=which.max(Q.est)
      }
      rew.hist[i_trial,i_rep]=rnorm(1)+Q.true[action]
      opt.hist[i_trial,i_rep]=action==opt.ID
      act.count[action]=act.count[action]+1
      Q.cum[action]=Q.cum[action]+rew.hist[i_trial,i_rep]
      Q.est[action]=Q.cum[action]/act.count[action]
    }
  }
  return(list(opt = opt.hist, rew = rew.hist))
}
n.Trial = 1000; n.Rep = 2000
type1=epGreedy(n.Trial, n.Rep, 0.0)
type2=epGreedy(n.Trial, n.Rep, 0.01)
type3=epGreedy(n.Trial, n.Rep, 0.1)

res.all = tibble(trial = rep(1:nTrial, n.Rep * 3),
                 rew = c(as.vector(type1$rew),as.vector(type2$rew),as.vector(type3$rew)),
                 opt = c(as.vector(type1$opt),as.vector(type2$opt),as.vector(type3$opt)),
                 condition = c(rep("0.00", nTrial * n.Rep),
                               rep("0.01", nTrial * n.Rep),
                               rep("0.10", nTrial * n.Rep)))

p1 <- ggplot(res.all) +
  geom_smooth(aes(x = trial, y = rew, color = condition)) +
  ylab("Average Reward")
p2 <- ggplot(res.all) +
  geom_smooth(aes(x = trial, y = opt, color = condition)) +
  ylab("Prop. Optimal Action")
grid.arrange(p1, p2, nrow =2)

softmax = function(nTrial, nRep, tau) {
  rew.hist = opt.hist = matrix(0,nrow=nTrial,ncol=nRep)
  for (i_rep in 1:nRep){
    Q.true=rnorm(10); opt.ID=which.max(Q.true)
    Q.est = Q.cum = act.count=rep(0,10);
    for (i_trial in 1:nTrial) {
      p = exp(Q.est/tau)/sum(exp(Q.est)/tau)
      action = sample(1:10, 1, prob = p)
      rew.hist[i_trial,i_rep]=rnorm(1)+Q.true[action]
      opt.hist[i_trial,i_rep]=action==opt.ID
      act.count[action]=act.count[action]+1
      Q.cum[action]=Q.cum[action]+rew.hist[i_trial,i_rep]
      Q.est[action]=Q.cum[action]/act.count[action]
    }
  }
  return(list(opt = opt.hist, rew = rew.hist))
}

n.Trial = 1000; n.Rep = 1000
eG00=epGreedy(n.Trial, n.Rep, 0.0)
eG01=epGreedy(n.Trial, n.Rep, 0.1)
sm=softmax(n.Trial, n.Rep, 0.1)
res.all = tibble(trial = rep(1:n.Trial, n.Rep * 3),
                 rew = c(as.vector(eG00$rew),as.vector(eG01$rew),as.vector(sm$rew)),
                 opt = c(as.vector(eG00$opt),as.vector(eG01$opt),as.vector(sm$opt)),
                 condition = c(rep("epsilon = 0.0", n.Trial * n.Rep),
                               rep("epsilon = 0.1", n.Trial * n.Rep),
                               rep("softmax", n.Trial * n.Rep)))
p1 <- ggplot(res.all) +
  geom_smooth(aes(x = trial, y = rew, color = condition)) +
  ylab("Average Reward")

p2 <- ggplot(res.all) +
  geom_smooth(aes(x = trial, y = opt, color = condition)) +
  ylab("Prop. Optimal Action")

grid.arrange(p1, p2, nrow =2)


# RL example
V=rep(0,25);                

# defining probability matrix
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

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

# calculating state-values iteratively
# fixed number of iterations 
nRep=1000; gamma=0.9;
for (i_rep in 1:nRep) {
  V.old=V
  for (i_state in 1:25) {
    V[i_state]=sum(P[i_state,]*(R[i_state,]+gamma*V.old[trM[i_state,]]))
  }
}
print(matrix(V,nrow=5)[5:1,])

# jisshu 1
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)

# defining Reward & trans. prob.
R=-1;P=0.25;

# policy improvement
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)
R=-1;P=0.25;V=rep(0,15);
delta=1; gamma=1; tol=1e-10; 
bestP=sample(1:4,14,replace=T)
stable=F;counter=0;
while (stable==F){
  counter=counter+1
  # iterative policy evaluation
  while (delta>tol) {
    delta=0;
    for (i_state in 1:14) {
      v=V[i_state]
      V[i_state]=sum(P*(R+gamma*V[trM[i_state,]]))
      delta=max(delta,abs(v-V[i_state]))
    }
  }
  # policy improvement
  stable=F
  for (i_state in 1:14) {
    b=bestP[i_state]
    bestP[i_state]=which.max(V[trM[i_state,]])
    ifelse((bestP[i_state]==b),stable<-T,stable<-F)
  }
}

install.packages("tidyverse")
library(tidyverse)

random.number <- rnorm(1000)
mean(random.number)
mean(random.number <- rnorm(1000))

rnorm(1000) %>% mean()

# CLT
NperSample = 10
SampleSize = 300000

# "traditional"
random.number <- runif(NperSample * SampleSize)
dat <- matrix(random.number, nrow=NperSample)
means <- colMeans(dat)
dens <- density(means)
hist(means, breaks = 100, probability = T, main = "Distributionf of Means")
lines(dens, lwd = 3, col = "orange")


runif(NperSample * SampleSize) %>%
  matrix(nrow=NperSample) %>%
  colMeans() -> means
hist(means, breaks=100,probability = T, main = "Distributionf of Means")
means %>% density() %>% lines(,lwd=3,col='orange')


histWdensity <- function(dat, nbreaks=30, main.txt){
  dens <- density(dat)
  hist(dat, breaks = nbreaks, probability = T, main = main.txt)
  lines(dens, lwd = 3, col = "orange")
}

runif(NperSample * SampleSize) %>%
  matrix(nrow=NperSample) %>%
  colMeans() %>% histWdensity(nbreaks = 100, main.txt = "Distributionf of Means")

dat <- read.csv("http://www.matsuka.info/data_folder/sampleData2013.txt")
dt <- as_tibble(dat)
dt.la <- filter(dt, affil == "LA")

dt.la2 <- filter(dt, affil == "LA" & grade == "SP")
dt.laNot2 <- filter(dt, affil == "LA" & grade != "SP")

dt.GB <- select(dt, grade, nbooks)
dtWOgender <- select(dt, -gender)

dt.arranged <- arrange(dt, affil, grade)

dt.weekly <- mutate(dt,nbooksWeek = nbooks / 52)

dt.atLeastOneBook <- mutate(dt, atleastOneBook = (nbooks/52) >= 1)
dt.atLeastOneBook <- mutate(dt, atleastOneBook = (nbooks/12) >= 1)

dt.BWindex = mutate(dt, nbooksWeek = nbooks / 52, idx = nbooksWeek / (12*7-Hworked))

dt.byGrade <- group_by(dt, grade)
summarize(dt.byGrade, ave.books = mean(nbooks,na.rm = TRUE), ave.Hworked = mean(Hworked, na.rm = TRUE))

dt.byGrAf <- group_by(dt, grade, affil)
dt.summ <- summarize(dt.byGrAf, ave.books = mean(nbooks,na.rm = TRUE), ave.Hworked = mean(Hworked, na.rm = TRUE), N = n())

dt.summ2 <- dt %>%
  group_by(grade, affil) %>%
  summarize(ave.books = mean(nbooks,na.rm = TRUE), ave.Hworked = mean(Hworked, na.rm = TRUE), N = n()) %>% filter(N > 2) %>% arrange(desc(ave.books))

plot(x = dt.summ2$ave.books, y = dt.summ2$ave.Hworked, pch=20, cex = 3,xlab = "Ave. # books read",ylab = "Ave hours worked")

dt <- read_csv("http://www.matsuka.info/data_folder/sampleData2013.txt",col_names = TRUE)

dt.summ3 <- dt %>%
  group_by(grade, gender) %>%
  summarize(ave.books = mean(nbooks,na.rm = TRUE), ave.Hworked = mean(Hworked, na.rm = TRUE))

dt.summ3G <- dt.summ3 %>% gather('ave.books', 'ave.Hworked', key = 'ave.values', value = "BorW")

dt.summ3SformG <- spread(dt.summ3G, key = ave.values, value =BorW)

dt.sumLA <- dt %>% filter(affil=="LA") %>% group_by(grade) %>% summarize(ave.books = mean(nbooks))

toeic <- tibble(
  grade = factor(c("SP", "JR")),
  score = c(800,830),
)

new.dt1 <- dt.sumLA %>% inner_join(toeic, by = "grade")

dt.sumLA <- add_row(dt.sumLA, grade = c("MS"), ave.books = (13))
toeic2 <- tibble(
  grade = factor(c("SP", "JR","DR")),
  score = c(800,830,900),
)
new.dt3 <- full_join(dt.sumLA, toeic2)

new.dt4 <- left_join(dt.sumLA, toeic2)
new.dt5 <- right_join(dt.sumLA, toeic2)

# CLT
NperSample = 10
SampleSize = 300000

runif(NperSample * SampleSize) %>%
  matrix(nrow=NperSample) %>%
  colMeans() %>% tibble(sample.mean = .) -> means

ggplot(means,aes(x = sample.mean, y = ..density..)) +
  geom_histogram(bins=200) +
  geom_density(colour = "orange",size=2)

ggplot(means,aes(x = sample.mean, y = ..density..)) +
  geom_histogram(bins=200) +
  geom_line(stat = "density", colour = "orange",size=2)

runif(NperSample * SampleSize) %>%
  matrix(nrow=NperSample) %>%
  colMeans() %>% tibble(sample.mean = .) %>%
  ggplot(., aes(x = sample.mean, y = ..density..)) +
  geom_histogram(bins=100,colour = "grey20") +
  geom_line(stat = "density", colour = "skyblue",size=2)


dat <- read.csv("http://www.matsuka.info/data_folder/sampleData2013.txt")
dt <- as_tibble(dat)




ggplot(dt, aes(x = Hworked, y = nbooks)) +
  geom_point(size = 3)

ggplot(dt) +
  geom_point(aes(x = Hworked, y = nbooks, color = grade),size = 3)

ggplot(dt) +
  geom_point(aes(x = Hworked, y = nbooks, shape = grade),size = 5)

ggplot(dt) +
  geom_point(aes(x = Hworked, y = nbooks),size = 5) +
  facet_wrap(~ grade, nrow = 1)

ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks))

ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks, linetype = grade))

ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks)) +
  facet_wrap(~ grade, nrow = 4)

ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks)) +
  geom_point(aes(x = Hworked, y = nbooks), size = 4)

ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks), colour = "black", se = FALSE) +
  geom_point(aes(x = Hworked, y = nbooks, color = grade), size = 4)

ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks, color = grade), se = FALSE) +
  geom_point(aes(x = Hworked, y = nbooks, color = grade), size = 4)


plot1 <- ggplot(dt) +
  geom_smooth(aes(x = Hworked, y = nbooks, color = grade), se = FALSE) +
  geom_point(aes(x = Hworked, y = nbooks, color = grade), size = 4)
plot1 + xlab("Hours worked") + ylab("Number of books read")

plot1 + xlab("Hours worked") +  ylab("Number of books read") +
  theme(axis.title.x = element_text(face = "italic",size = 14, colour = "navy"),
        axis.title.y = element_text(face = "bold",size = 10, colour = "darkgreen"))

ggplot(filter(dt, affil == "LA")) +
  geom_point(aes(x = Hworked, y = nbooks, color = grade), size = 4)


dt$grade <- fct_relevel(dt$grade, "FR","SP","JR","SR")
group_by(dt, grade) %>% summarize(ave.books = mean(nbooks, na.rm = T)) %>%
  ggplot() + geom_bar(aes(x = grade, y = ave.books), stat = "identity")


ggplot(dt) + geom_bar(aes(x = grade))
dt %>% count(grade)

group_by(dt, grade) %>% summarize(ave.books = mean(nbooks, na.rm = T),
                                  se = sd(nbooks, na.rm =T)/sqrt(n())) %>%
  ggplot(aes(x = grade, y = ave.books)) +
  geom_bar(stat = "identity", fill = "grey70") +
  geom_errorbar(aes(ymin = ave.books - se, ymax = ave.books +se), width = 0.2) +
  ylab("Average # books read")

ggplot(dt) +
  geom_boxplot(aes(x = grade, y = nbooks))
ggplot(dt) +
  geom_boxplot(aes(x = grade, y = nbooks)) +
  coord_flip()


ggplot(dt,aes(x = Hworked, y = nbooks)) +
  stat_density2d(aes(colour =..level..)) +
  geom_point()

ggplot(dt,aes(x = Hworked, y = nbooks)) +
  stat_density2d(aes(alpha =..density..), geom="tile",contour=F) +
  geom_point(alpha =0.4)



ggplot(dt) +
  stat_summary(aes(x = grade, y = nbooks),
               fun.y = mean,
               fun.ymin = function(x) mean(x) - sd(x),
               fun.ymax = function(x) mean(x) + sd(x))


dat <- read.csv("http://www.matsuka.info/data_folder/datWA01.txt")
dt <- as_tibble(dat)
dt.lm <- lm(h~shoesize, dt)
cfs <- coef(dt.lm)
ggplot(dt, aes(x = shoesize, y = h)) +
  geom_point() +
  geom_abline(intercept = cfs[1], slope = cfs[2], col = "red") +
  geom_text( x= 22, y =175, aes(label = paste("r^2  =", round(summary(dt.lm)$r.squared,3))))

10.11ABの講義は家族が急病の為、休講とさせてください。直前の連絡で申し訳ありません。
次回の講義・実習は10.25日です。
補講日は次回、皆さんの都合で調整しましょう。
補講日の候補は11.01、12.27、 および、東大が指定している補講日である、12.25（午後）、01.17（午前）、01.20、01.21です。

課題1.1 勾配法を用いた最小化

を最小化するxとyを慣性を含む勾配法を用いて求めてください。

fun01 = function(x,y){
  z = 1/20*x^2 + y^2
  return(z)
}

結果の可視化は以下のようなものでcontourを作図して:

install.packages("plot3D")
library(plot3D)
x = seq(-10,10,0.2)
y = seq(-5,5,0.2)
M = mesh(x,y)
Z = as.vector(1/20*M$x^2)+as.vector(M$y^2) 
Z.mesh = matrix(Z,nrow(M$x))
contour(x, y, Z.mesh, drawlabels = F, nlevels=40)

その後linesで更新履歴を追加してみると良いかもしれません。

# gdは勾配法でのx、yの更新履歴
lines(gd, type='o', col = 'green', pch=20)

# gdmは慣性ありの勾配法でのx、yの更新履歴
lines(gdm, type='o', col = 'blue', pch=20)

実装例

fun01 = function(x,y){
  z = 1/20*x^2 + y^2
  return(z)
}

fun01.grad <- function(x,y){
  grad.x = 1/10*x
  grad.y = 2*y
  return(list(gx=grad.x, gy = grad.y))
}

# script
n.loop = 100; lambda = 0.9
x = y = rep(0,n.loop)
x[1] = -7; y[1] = 2
for (i.loop in 2:n.loop){
  g = fun01.grad(x[i.loop - 1],y[i.loop - 1])
  x[i.loop] = x[i.loop - 1] - lambda*g$gx
  y[i.loop] = y[i.loop - 1] - lambda*g$gy
}

# function
min.fun01 <- function(init.x, init.y, n.loop, labmda){
  x = y = rep(0,n.loop)
  x[1] = init.x; y[1] = init.y
  for (i.loop in 2:n.loop){
    g = fun01.grad(x[i.loop - 1],y[i.loop - 1])
    x[i.loop] = x[i.loop - 1] - lambda*g$gx
    y[i.loop] = y[i.loop - 1] - lambda*g$gy
  }
  return(list(x = x, y = y))
}

res <- min.fun01(-7, 2, 100, 0.9)
temp.x = seq(-10,10,0.2);temp.y = seq(-5,5,0.2)
M = mesh(temp.x,temp.y)
Z = as.vector(1/20*M$x^2)+as.vector(M$y^2)
Z.mesh = matrix(Z,nrow(M$x))
contour(temp.x, temp.y, Z.mesh, drawlabels = F, nlevels=40)
lines(res$x, res$y, type='o', col="green")

min.fun01M <- function(init.x, init.y, n.loop, labmda, gamma){
  x = y = rep(0,n.loop)
  x[1] = init.x; y[1] = init.y; v = rep(0,2)
  for (i.loop in 2:n.loop){
    g = fun01.grad(x[i.loop - 1],y[i.loop - 1])
    v = gamma*v - c(lambda*g$gx, lambda*g$gy)
    x[i.loop] = x[i.loop - 1] + v[1]
    y[i.loop] = y[i.loop - 1] + v[2]
  }
  return(list(x = x, y = y))
}

resM <- min.fun01M(-7, 2, 100, 0.9, 0.4)
lines(resM$x, resM$y, type='o', col="blue")

### with ggplot ###
ES <- tibble(x = as.vector(M$x),
             y = as.vector(M$y),
             z = Z)

res <- tibble(x = c(res$x, resM$x),
              y = c(res$y, resM$y),
              type = c(rep("Std",100), rep("moment.",100)))

ggplot(ES,aes(x = x, y = y, z =z )) + stat_contour(bins = 10) +
  geom_line(aes(x = x, y = y, color = type), data = res)

課題1.2 multi-armed bandit problem
wikipediaの記事
何種類かバージョンはありますが、この課題では各スロットマシーンの真のリワードは、mean = 0、sd = 1の正規分布に従うこととしてください。ただし、ある特定の回で実際に得られるリワードは真のリワード、プラス正規分布にしたがうノイズがあるとします（ここでも、mean = 0、sd = 1としましょう）。
真のリワードは不明で、1000回スロットマシーンを引いた時の総リワードを最大化することを目的とします。
1つ目の方策はgreedy法で、リワードの推定値が最も高いスロットマシーンを選択するといったものです。
2つ目の方策は基本的にはgreedy法で、ある確率epsilonでランダムにスロットマシーンを選択するものです。
各方策につき複数回検証しましょう（M回）。

# スロットマシーンの数を５としました（実際には幾つでも結構です）
N = 5

# スロットマシーンを引く回数を1000としました（実際には幾つでも結構です）
nTrial = 1000

# 真のリワードを生成します。
reward.true = rnorm(N, mean=0, sd=1)
> reward.true
[1] -0.2822860  0.5645874 -0.1968128  0.5430834 -0.3696859

# リワードの推定値を初期化します。
reward.est = rep(0, N)

# リワードの累積和を初期化します。
reward.cum = rep(0, N)

# 各スロットマシーンを引いた回数を初期化します。
sm.count = rep(0, N)

# リワードの履歴を初期化します。
reward.hist = rep(0, nTrial)

# greedy法で、どのスロットマシーンを引くか選択します。
# reward.estの最大値が複数個ある場合は、それらから1つをランダムで選択します。
max.est = which(max(reward.est) == reward.est)
if (length(max.est) > 1){
  selected = sample(max.est, 1)
} else {selected = max.est}

# 今回は5が選択されました。
> selected
[1] 5

# スロットマシーン5を引いてみます。
# 真のリワードは
# > reward.true[selected]
# [1] -0.3696859
# ですが、今回実際に得られるのはこれにノイズが乗ります。
reward = reward.true[selected] + rnorm(1, mean = 0, sd =1)
> reward
[1] -1.61256

reward.hist[1] = reward
# 繰り返す場合は、reward.hist[i.trial] = reward など 


# リワードの推定値をアップデートします
reward.cum[selected] = reward.cum[selected] + reward
sm.count[selected] = sm.count[selected] + 1
reward.est[selected] = reward.cum[selected] / sm.count[selected]
> reward.est
[1]  0.00000  0.00000  0.00000  0.00000 -1.61256

# 2回目
max.est = which(max(reward.est) == reward.est)
if (length(max.est) > 1){
   selected = sample(max.est, 1)
 } else {selected = max.est}
> selected 
[1] 2

reward = reward.true[selected] + rnorm(1, mean = 0, sd =1)
> reward
[1] 1.497099

reward.cum[selected] = reward.cum[selected] + reward
sm.count[selected] = sm.count[selected] + 1
reward.est[selected] = reward.cum[selected] / sm.count[selected]

> reward.est
[1]  0.000000  1.497099  0.000000  0.000000 -1.612560

# これをnTrial分繰り返します。

# 2つの方策の良さの検証は、特定のreward.trueの値に依存するべきではないので、
# reward.trueの値を変えてみます。これをM回繰り返してみましょう。

実装例：

## n-armed bandit
nTrial = 1000; epsilon = 0.0;
Q.true = rnorm(10); opt.ID = which.max(Q.true)
Q.est = Q.cum = act.count = rep(0, 10);
rew.hist = opt.hist = rep(0, nTrial);
for (i_trial in 1:nTrial) {
  if (runif(1) < epsilon) {
    action = sample(1:10, 1)
  } else {
    action = which.max(Q.est)
  }
  rew.hist[i_trial] = rnorm(1) + Q.true[action]
  opt.hist[i_trial] = action == opt.ID
  act.count[action] = act.count[action] + 1
  Q.cum[action] = Q.cum[action] + rew.hist[i_trial]
  Q.est[action] = Q.cum[action] / act.count[action]
}
plot(rew.hist, type='l')

epGreedy = function(nTrial,nRep,epsilon) {
  rew.hist = opt.hist = matrix(0,nrow=nTrial,ncol=nRep)
  for (i_rep in 1:nRep){
    Q.true=rnorm(10); opt.ID=which.max(Q.true)
    Q.est = Q.cum = act.count=rep(0,10);
    for (i_trial in 1:nTrial) {
      if (runif(1) < epsilon) {
        action=sample(1:10,1)
      } else {
        action=which.max(Q.est)
      }
      rew.hist[i_trial,i_rep]=rnorm(1)+Q.true[action]
      opt.hist[i_trial,i_rep]=action==opt.ID
      act.count[action]=act.count[action]+1
      Q.cum[action]=Q.cum[action]+rew.hist[i_trial,i_rep]
      Q.est[action]=Q.cum[action]/act.count[action]
    }
  }
  return(list(opt = opt.hist, rew = rew.hist))
}
n.Trial = 1000; n.Rep = 2000
type1=epGreedy(n.Trial, n.Rep, 0.0)
type2=epGreedy(n.Trial, n.Rep, 0.01)
type3=epGreedy(n.Trial, n.Rep, 0.1)

par(mfrow=c(2,1))
plot(colMeans(type3$rew),type='l',xlab="Play",ylab="average reward")
lines(colMeans(type2$rew),type='l',col='red')
lines(colMeans(type1$rew),type='l',col='green')
legend("bottomright",c("epsilon=0.00","epsilon=0.01","epsilon=0.10"),
       col=c("black","red","green"),lty=c(1,1,1))
plot(colMeans(type3$opt),type='l',xlab="Play",ylab="% optimal action")
lines(colMeans(type2$opt),type='l',col='red')
lines(colMeans(type1$opt),type='l',col='green')
legend("bottomright",c("epsilon=0.00","epsilon=0.01","epsilon=0.10"),
       col=c("black","red","green"),lty=c(1,1,1))

### with ggplot ###
res.all = tibble(trial = rep(1:nTrial, n.Rep * 3),
                 rew = c(as.vector(type1$rew), 
                         as.vector(type2$rew),
                         as.vector(type3$rew)),
                 opt = c(as.vector(type1$opt), 
                         as.vector(type2$opt),
                         as.vector(type3$opt)),
                 condition = c(rep("0.00", nTrial * n.Rep),
                               rep("0.01", nTrial * n.Rep),
                               rep("0.10", nTrial * n.Rep)))

p1 <- ggplot(res.all) +
  geom_smooth(aes(x = trial, y = rew, color = condition)) +
  ylab("Average Reward")

p2 <- ggplot(res.all) +
  geom_smooth(aes(x = trial, y = opt, color = condition)) +
  ylab("Prop. Optimal Action")

library(gridExtra)
grid.arrange(p1, p2, nrow =2)
###