Posted on by

認知情報解析　課題２

# 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]])
    }
  }
}
Posted on by

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

dat<-read.csv("http://www.matsuka.info/data_folder/datWA01.txt")
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)
Posted on by

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)),
        col="blue",adj = 0)

text(x = mean(dat$shoesize)+0.2, y = 145,
        paste("mean shoesize =", round(mean(dat$shoesize),2)),
        col="darkgreen",adj = 0)

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)

dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
plot(dat, pch=20, col='blue')

dat.pca<-read.table("http://www.matsuka.info/data_folder/tdkPCA01.txt")
plot(dat.pca, pch = rownames(dat.pca), cex = 1.7, col = "blue")

dat<-read.table("http://www.matsuka.info/data_folder/aov01.txt",header=T)
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
Posted on by

院：認知情報解析学

source("http://www.matsuka.info/univ/course_folder/cuUtil02.R")
dat<-read.csv('http://www.matsuka.info/data_folder/tdkCFA.csv')
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)
Posted on by

強化学習　方策の比較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))
  }
Posted on by

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
dat<- read.csv("http://www.matsuka.info/data_folder/datWA01.txt")
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)

dat<-read.table("http://www.matsuka.info/data_folder/aov01.txt")
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)
Posted on by

2019 認知情報解析学演習　RL01

set.seed(111)
n.trial = 1000; N = 10; sigma  = 1
Q.star = runif(N); Q = rep(0, N)
count = rep(0,N); Q.cum = rep(0, N)
rew.earned = rep(0,n.trial)
### playing slot-machine
for (i.trial in 1:n.trial){
  max.a = max(Q)
  max.idx = which(Q == max.a)
  if (length(max.idx)>1){
    max.idx = sample(max.idx, 1)
  }
  r.t = rnorm(1, Q.star[max.idx], sd = sigma)
  Q.cum[max.idx] = Q.cum[max.idx] + r.t
  count[max.idx] = count[max.idx] + 1
  Q[max.idx] = Q.cum[max.idx] / count[max.idx]
  rew.earned[i.trial] = r.t
}
plot(rew.earned,type='l')
Q
Posted on by

2019 データ解析基礎論a DAA02

x<-matrix(1:8, nrow=2)
x<-matrix(1:8, nrow=2,byrow=T)
data01<-data.frame(score = c(2,4,3,4),     
                   dose = c(rep(10,2),rep(100,2)),  
                   condition = rep(c('exp','control'),2))

dat01<-read.csv("http://www.matsuka.info/data_folder/temp_data01.txt",  
                header=T)
dat02<-read.csv("http://www.matsuka.info/data_folder/temp_data02.txt",    
                header=T, row.name=1)
dat03<-read.table("http://www.matsuka.info/data_folder/temp_data03.txt",
                  header=T, row.name=4)
dat<-read.csv("http://www.matsuka.info/data_folder/datWA01.txt",   
              header=T);
mean(dat$shoesize[dat$gender == "M"])
mean(dat$shoesize[dat$gender == "F"])
mean(dat$shoesize[dat$h > 180])

v1 = seq(-3,3,0.1)
v2 = v1^2
plot(x = v1, y = v2)
plot(v1, v2, col = 'red')

plot(v1, v2, main = "THIS IS THE TITLE", cex.lab = 1.5,
     xlab = "Label for X-axis",ylab = "Label for Y-axis")

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

dat<- read.csv("http://www.matsuka.info/data_folder/datWA01.txt")
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)

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)

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

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) 

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)

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='green')

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)
Posted on by

2019 データ解析基礎論A DAA01

dat<-data.frame(score=c(78,70,66,76,78,76,88, 76, 76,72,60,72,70,72,84,70),
                cond=c(rep('low',8), rep('high',8)))
boxplot(score~cond, col = c("skyblue",'skyblue4'),data=dat)
summary(aov(score ~ cond, data = dat))


dat <- read.csv("http://www.matsuka.info/data_folder/hwsk8-17-6.csv")
plot(ani~otouto, data=dat,pch=20,cex=3,xlab ="score of Otouto", ylab = "score of Ani")
dat.lm <- lm(ani~otouto, data=dat)
abline(dat.lm, col = 'red',lwd = 2.5)

dat<-read.csv("http://www.matsuka.info/data_folder/datWA01.txt")
dat.glm <- glm(gender~shoesize,family="binomial",data=dat)
plot(as.numeric(gender)-1~shoesize,data=dat,pch=20,cex=3,ylab="P(Male)")
cf = coef(dat.glm)
temp.x = seq(20,30,0.1)
y = 1/(1+exp(-1*(cf[1]+temp.x*cf[2])))
lines(temp.x,y,col='cyan',lwd=2)

dat <- read.table("http://www.matsuka.info/data_folder/tdkPCA01.txt")
dat.pca <- princomp(dat)
biplot(dat.pca)

dat<-read.csv("http://matsuka.info/data_folder/tdkClust.csv", header=TRUE, row.names=1)
dat.cluster=hclust(dist(dat),method="average")
plot(dat.cluster,cex=1.5)

data01<-data.frame(score = c(2,4,3,4),
                   dose = c(rep(10,2),rep(100,2)),
                   condition = rep(c('exp','control'),2))

dat01<-read.csv("http://www.matsuka.info/data_folder/temp_data01.txt",
                header=T)
dat02<-read.csv("http://www.matsuka.info/data_folder/temp_data02.txt",
                header=T, row.name=1)
dat03<-read.table("http://www.matsuka.info/data_folder/temp_data03.txt",
                  header=T, row.name=4)
dat<-read.csv("http://www.matsuka.info/data_folder/datWA01.txt",
              header=T);
Posted on by

認知情報解析学演習b

init.RNN <- function(n.uniq,size.hidden){
  W.h = matrix(rnorm(size.hidden*size.hidden), nrow = size.hidden)*0.01
  W.x = matrix(rnorm(n.uniq*size.hidden), nrow = n.uniq)*0.01
  b.h = matrix(rnorm(size.hidden), nrow = 1)*0.01
  W.o = matrix(rnorm(n.uniq*size.hidden),nrow = size.hidden)*0.01
  b.o = matrix(rnorm(n.uniq), nrow = 1)*0.01
  return(list(W.h = W.h, W.x= W.x, b.h = b.h, W.o = W.o, b.o = b.o))
}

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.pred <- 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)
}

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

RNN.forward <- function(h.prev, x, network){
  b.size = nrow(x)
  h = h.prev%*%network$W.h + x%*%network$W.x
  hb = h+matrix(1,nrow=b.size,ncol=1)%*%network$b.h
  h.next = tanh(hb)
  return(h.next = h.next)
}

RNN.backward <- function(dh.next, network, x, h.next, h.prev){
  dt = dh.next * (1- h.next^2)
  db = colSums(dt)
  dW.h = t(h.prev)%*%dt
  dh.prev = dt%*%t(network$W.h)
  dW.x = t(x)%*%dt
  dx = dt%*%t(network$W.x)
  return(list(db = db, dW.h = dW.h, dh.prev = dh.prev, dW.x = dW.x, dx=dx))
}

txt = "You say goodbye and I say hello.  you say goodbay and I say hello"

txt2corpus <- function(txt){
  txt = tolower(txt)
  txt = gsub('[.]', ' . sos',txt)
  words = unlist(strsplit(c('sos',txt), " "))
  uniq.words = unique(words)
  n.uniq = length(uniq.words)
  n.words = length(words)
  corpus = rep(0,n.words)
  corpus = match(words,uniq.words)
  return(corpus)
}

corp2contxt1SRNN = function(corpus){
  len.corp = length(corpus)
  # creating target matrix
  idxT = cbind(1:(len.corp-1), corpus[2:len.corp])
  targ1S = matrix(0,nrow=len.corp-1,ncol=length(unique(corpus)))
  targ1S[idxT]=1
  # creating context matrices
  idxC = cbind(1:(len.corp-1),corpus[1:(len.corp-1)])
  contxt = matrix(0,nrow=len.corp-1,ncol=length(unique(corpus)))
  contxt[idxC]=1
  return(list(y=targ1S,x=contxt))
}

corp = txt2corpus(txt)
dat = corp2contxt1SRNN(corp)

size.hidden = 7
network <- init.RNN(8,size.hidden)

n.rep = 100000;lambda = 0.01;batch.size = 3; time = 3;
h.prev =array(0, c(batch.size, size.hidden, time))
h.next = array(0, c(batch.size, size.hidden, (time+1)))
loss = rep(0,n.rep)
for (i.rep in 1:n.rep){
  for (i.t in 1:time){
    idx = seq(i.t, corp.len, time)
    h.next[, , i.t] = RNN.forward(h.prev[, , i.t], dat$x[idx,], network)
    if (i.t < time){
      h.prev[, , (i.t+1)] = h.next[, , i.t]
    }
  }
  dWHs = matrix(0,nrow=nrow(network$W.h),ncol=ncol(network$W.h))
  dWXs = matrix(0,nrow=nrow(network$W.x),ncol=ncol(network$W.x))
  dBHs = matrix(0,nrow=nrow(network$b.h),ncol=ncol(network$b.h))
  dWOs = matrix(0,nrow=nrow(network$W.o),ncol=ncol(network$W.o))
  dBOs = matrix(0,nrow=nrow(network$b.o),ncol=ncol(network$b.o))
  d.prev = matrix(0,nrow=batch.size,ncol=size.hidden)
  L = 0
  for (i.t in time:1){
    idx = idx = seq(i.t, corp.len, time)
    O = affine.forwd(h.next[,,i.t], network$W.o, network$b.o)
    L = L + softmax.forwd(O, dat$y[idx,])
    ds = softmax.bckwd(O, dat$y[idx,], 1)
    dW.o = affine.bckwd(h.next[,,i.t], network$W.o, network$b.o, ds)
    dWOs = dWOs + dW.o$dW
    dBOs = dBOs + dW.o$db
    RNN.d = RNN.backward(dW.o$dx+d.prev, network, dat$x[idx,],h.next[,,(i.t+1)],h.prev[,,i.t])
    dWHs = dWHs + RNN.d$dW.h
    dWXs = dWXs + RNN.d$dW.x
    dBHs = dBHs + RNN.d$db
    d.prev = RNN.d$dh.prev
  }
  loss[i.rep] = L
  network$W.o = network$W.o - lambda*dWOs
  network$b.o = network$b.o - lambda*dBOs
  network$W.h = network$W.h - lambda*dWHs
  network$W.x = network$W.x - lambda*dWXs
  network$b.h = network$b.h - lambda*dBHs
}
plot(loss,type='l')
par(mfrow=c(9,1))
for (i.t in 1:time){
  idx = idx = seq(i.t, corp.len, time)
  O = affine.forwd(h.next[,,i.t], network$W.o, network$b.o)
  print(softmax.pred(O, dat$y[idx,]))
  for (i in 1:3){
    barplot(softmax.pred(O, dat$y[idx,])[i,])
  }
}