オンライン実験の案内

オンライン実験の案内です。
様々な都市または国に関する諸課題を実施します。詳細は下のリンクをクリックしてください。
参加者数の人数制限もあり、認知情報科学発展実習・特別実習(および11月25日までは認知情報科学基礎)の受講者のみを対象としています。これらを受講していな方には謝金を支払うことは出来ませんので、ご了承ください。

実験に参加する・内容を確認するにはここをクリックしてください。

質問などは松香まで連絡ください。

Disclaimer

このcourselogにあるコードは、主に学部生・博士課程前期向けの講義・演習・実習で例として提示しているもので、原則直感的に分かりやすいように書いているつもりです。例によってはとても非効率なものもあるでしょうし、「やっつけ」で書いているものあります。また、普段はMATLABを使用していますので、変な癖がでているかもしれません。これらの例を使用・参考にする場合はそれを踏まえてたうえで使用・参考にして下さい。
卒業論文に関する資料:[2015PDF] [word template] [latex template] [表紙] [レポートの書き方] [引用文献など]
A Brief Guide to R (Rのコマンドの手引きおよび例)はこちらをどうぞ

認知情報解析学演習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,])
  }
}

認知情報解析学演習 RNN

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

corp = txt2corpus(txt)

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

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 = matrix(rnorm(size.hidden),nrow=1)*0.01
  return(list(W.h = W.h, W.x= W.x, b = b))
}
MatMult.forwd <- function(x, W){
  return(x%*%W)
}

MatMult.bckwd <- function(x, W, dout){
  dx = dout%*%t(W)
  dW = t(x)%*%dout
  return(list(dx = dx, dW = dW))
}

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

batch.size = 3
time = 3

size.hidden = 7
corp = txt2corpus(txt)
dat<-corp2contxt1SRNN(corp)
network <- init.RNN(ncol(dat$y), size.hidden)
corp.len = nrow(dat$y)
h.prev =array(0, c(batch.size, size.hidden, time))
h.next = array(0, c(batch.size, size.hidden, (time+1)))
W.out = matrix(rnorm(size.hidden * ncol(dat$y)), nrow = size.hidden)
b.out = matrix(rnorm(ncol(dat$y)), nrow = 1)
n.rep = 100000;lambda = 0.25;

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))
  dBs = matrix(0,nrow=nrow(network$b),ncol=ncol(network$b))
  dWOs = matrix(0,nrow=nrow(W.out),ncol=ncol(W.out))
  dBOs = matrix(0,nrow=nrow(b.out),ncol=ncol(b.out))
  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], W.out, b.out)
    L = L + softmax.forwd(O, dat$y[idx,])
    ds = softmax.bckwd(O, dat$y[idx,], 1)
    dW.o = affine.bckwd(h.next[,,i.t], W.out, b.out, 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])
    #RNN.d = RNN.backward(dW.o$dx, network, dat$x[idx,],h.next[,,i.t],h.prev[,,i.t])
    dWHs = dWHs + RNN.d$dW.h
    dWXs = dWXs + RNN.d$dW.x
    dBs = dBs + RNN.d$db
    d.prev = RNN.d$dh.prev
  }
  loss[i.rep] = L
  W.out = W.out - lambda*dWOs
  b.out = b.out - lambda*dBOs
  network$W.h - lambda*dWHs
  network$W.x - lambda*dWXs
  network$b - lambda*dBs
}
plot(loss, type='l')
for (i.t in 1:time){
  idx = idx = seq(i.t, corp.len, time)
  O = affine.forwd(h.next[,,i.t], W.out, b.out)
  print(softmax.pred(O, dat$y[idx,]))
}

院:認知情報解析演習

txt2 = "You say Yes, I say No. You say Stop and I say Go, go, go. Oh no. You say Goodbye and I say Hello, hello, hello. I don't know why you say Goodbye, I say Hello, hello, hello. I don't know why you say goodbye, I say hello. I say High, you say Low. You say Why? And I say I don't know. Oh no. You say Goodbye and I say Hello, hello, hello. I don't know why you say Goodbye, I say Hello, hello, hello. Hello, goodbye, hello, goodbye. Hello, goodbye. I don't know why you say Goodbye, I say Hello. Why, why, why, why, why, why, do you Say Goodbye, goodbye, bye, bye. Oh no. You say Goodbye and I say Hello, hello, hello. I don't know why you say Goodbye, I say Hello, hello, hello. I don't know why you say Goodbye, I say Hello. You say Yes, I say No. I say Yes, but I may mean No. You say Stop, I say Go, go, go. I can stay still it's time to go. Oh, oh no. You say Goodbye and I say Hello, hello, hello. I don't know why you say Goodbye, I say Hello, hello, hello. I don't know why you say Goodbye, I say Hello, hello, hello."
txt2corpus <- function(txt){
  txt = tolower(txt)
  txt = gsub('[.]', ' .',txt)
  txt = gsub('[,]', ' ,',txt)
  words = unlist(strsplit(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(list(corpus = corpus, uniq.words = uniq.words))
}
words.sim <- function(word1, word2, eps=1e-8){
  nw1 = word1/(sqrt(sum(word1^2)) + eps)
  nw2 = word2/(sqrt(sum(word2^2)) + eps)
  return(nw1%*%nw2)
}


most.sim <- function(word, uniq.words, cc, N=5){
  n.uniq = length(uniq.words)
  word2 = cc[which(uniq.words==word),]
  res = data.frame(words = uniq.words, similarity = rep(0,n.uniq))
  top = data.frame(words = rep("",N+1), similarity=rep(0,N+1))
  res$similarity = apply(cc,1, function(x) words.sim(x,word2))
  sort.res = sort(res$similarity, decreasing = T,  index.return=T)
  top$words = res$words[sort.res$ix[1:(N+1)]]
  top$similarity = res$similarity[sort.res$ix[1:(N+1)]]
  self.idx = which(top$words == word)
  return(top[-self.idx,])
}

corp2 = txt2corpus(txt2)
dat2 <- corp2contxt1S(corp2$corpus)
n.word = max(corp2$corpus)
network<-init.W2V(30,10)
res <- w2v.ws1(network,dat2,5,10000)
plot(res$loss,type='l')
most.sim("hello", corp2$uniq.words, res$network$W.in, 5)
embed.forwd <- function(W, idx){
  return(W[idx,])
}

embed.dot.forwd <- function(W, h, idx){
  return(O=rowSums(W[idx,]*h))
}

sigmoidWL.forwd <- function(O,target){
  delta = 1e-7
  y = 1/(1+exp(-O)) - delta
  loss=-sum(target*log(y)+(1-target)*log(1 - y))
  return(loss)
}

sigmoidWL.backwd <- function(O,target,dout=1){
  delta = 1e-5
  y = 1/(1+exp(-O)) - delta
  dW = (y - target)*dout
  return(dW)
}

embed.dot.backwd <- function(W,h,idx,dout){
  dh = dout * W[idx,]
  dW = dout * h
  return(list(dh = dh, dW = dW))
}

embed.backwd <- function(W, idx, dout){
  dW = matrix(0,nrow(W),ncol(W))
  for (i.idx in 1:length(idx)){
    dW[idx[i.idx], ] = dW[idx[i.idx], ] + dout[i.idx, ]
  }
  return(dW)
}

txt2corpus <- function(txt){
  txt = tolower(txt)
  txt = gsub('[.]', ' .',txt)
  words = unlist(strsplit(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)
}
 
corp2contxt <- function(corpus){
  len.corp = length(corpus)
  target = corpus[2:(len.corp-1)]
  col1 = corpus[1:(len.corp-2)]
  col2 = corpus[3:len.corp]
  context = cbind(col1,col2)
  return(list(context=context,target = target))
}

init.W2V <- function(n.uniq,size.hidden){
  W.in = matrix(rnorm(n.uniq*size.hidden),nrow=n.uniq)*0.01
  W.out = matrix(rnorm(n.uniq*size.hidden),nrow=n.uniq)*0.01
  return(list(W.in = W.in, W.out= W.out))
}

network<-init.W2V(8,5)
n.batch = 3;
txt = "You say goodbye and I say hello.."
corp = txt2corpus(txt)
corp = c(8,8,corp)
dat<-corp2contxt(corp)
n.data =  length(dat$target)

n.iter = 3000;
lambda=0.1;
loss = rep(0, n.iter)
for (i.iter in 1:n.iter){
  samp.idx = sample(1:n.data,n.batch)
  batch.c1 = dat$context[samp.idx,1]
  batch.c2 = dat$context[samp.idx,2]
  batch.y = dat$target[samp.idx]
  h1 <- embed.forwd(network$W.in,batch.c1)
  h2 <- embed.forwd(network$W.in,batch.c2)
  h = 0.5 * (h1 + h2)
  # positive only
  s = embed.dot.forwd(network$W.out,h,batch.y)
  z = sigmoidWL.forwd(s,1)
  loss[i.iter] = z
  ds = sigmoidWL.backwd(s, 1, 1)
  dE = embed.dot.backwd(network$W.out,h, batch.y, ds)
  dh = dE$dh*0.5
  dIn1 = embed.backwd(network$W.in,dat$context[samp.idx,1], dh)
  dIn2 = embed.backwd(network$W.in,dat$context[samp.idx,2], dh)
  network$W.out[batch.y,] = network$W.out[batch.y,] - lambda*dE$dW
  network$W.in = network$W.in - lambda*dIn1
  network$W.in = network$W.in - lambda*dIn2
}
plot(loss, type='l')

sigmoid.forwd <- function(O){
  return(1/(1+exp(-O)))
}

samp.idx = c(2:6,8,9,1)
batch.c1 = dat$context[samp.idx,1]
batch.c2 = dat$context[samp.idx,2]
batch.y = dat$target[samp.idx]
h1 <- embed.forwd(network$W.in,batch.c1)
h2 <- embed.forwd(network$W.in,batch.c2)
h = 0.5 * (h1 + h2)
s = MatMult.forwd(h,t(network$W.out))
z = sigmoid.forwd(s)
res=cbind(z,batch.y)
#par(mfrow=c(8,1))
#for (i in 1:8){
#  col.spec = rep("black",8)
#  col.spec[i]="orange"
# barplot(res[i, 1:8],col=col.spec)
#}

par(mfrow=c(4,1))
for (i in 1:4){
  col.spec = rep("black",8)
  col.spec[i]="orange"
  barplot(res[i, 1:8],col=col.spec)
}
par(mfrow=c(4,1))
for (i in 5:8){
  col.spec = rep("black",8)
  col.spec[i]="orange"
  barplot(res[i, 1:8],col=col.spec)
}
unigram.sampler <- function(corpus, target, power, sample.size){
  neg.corpus <- corpus[-which(corpus == target)]
  uniq.word = unique(neg.corpus)
  n.word = length(neg.corpus)
  prob = (as.vector(table(neg.corpus))/n.word)^power
  prob = prob/sum(prob)
  sample.idx = sample(uniq.word, sample.size, prob = prob)
  return(sort(sample.idx))
}
get.neg.sample <- function(corpus, target, power, sample.size){
  sample.id = matrix(0, nrow = length(target), ncol = sample.size)
  for (i.targ in 1:length(target)){
    sample.id[i.targ,] = unigram.sampler(corpus, target[i.targ], power, sample.size)
  }
  return(sample.id)
}
init.W2V <- function(n.uniq,size.hidden){
  W.in = matrix(rnorm(n.uniq*size.hidden),nrow=n.uniq)*0.01
  W.out = matrix(rnorm(n.uniq*size.hidden),nrow=n.uniq)*0.01
  return(list(W.in = W.in, W.out= W.out))
}



network<-init.W2V(8,5)
n.batch = 3;
txt = "You say goodbye and I say hello.."
corp = txt2corpus(txt)
corp = c(8,8,corp)
dat<-corp2contxt(corp)
n.data =  length(dat$target)

n.iter = 10000;
lambda=0.1;
loss = rep(0, n.iter)
sample.size = 2
for (i.iter in 1:n.iter){
  samp.idx = sample(1:n.data,n.batch)
  batch.c1 = dat$context[samp.idx,1]
  batch.c2 = dat$context[samp.idx,2]
  batch.y = dat$target[samp.idx]
  h1 <- embed.forwd(network$W.in,batch.c1)
  h2 <- embed.forwd(network$W.in,batch.c2)
  h = 0.5 * (h1 + h2)
  # positive
  s = embed.dot.forwd(network$W.out,h,batch.y)
  z = sigmoidWL.forwd(s,1)
  loss[i.iter] = z
  ds = sigmoidWL.backwd(s, 1, 1)
  dE = embed.dot.backwd(network$W.out,h, batch.y, ds)
  dh = dE$dh*0.5
  dIn1 = embed.backwd(network$W.in,dat$context[samp.idx,1], dh)
  dIn2 = embed.backwd(network$W.in,dat$context[samp.idx,2], dh)
  network$W.out[batch.y,] = network$W.out[batch.y,] - lambda*dE$dW
  network$W.in = network$W.in - lambda*dIn1
  network$W.in = network$W.in - lambda*dIn2
  # negative
  neg.set <- get.neg.sample(corp,batch.y, 0.75, sample.size)
  for (i.neg in 1:sample.size){
      s = embed.dot.forwd(network$W.out,h,neg.set[,i.neg])
      z = sigmoidWL.forwd(s,0)
      loss[i.iter] = loss[i.iter] + z
      ds = sigmoidWL.backwd(s, 0, dout=1)
      dE = embed.dot.backwd(network$W.out,h,neg.set[,i.neg], ds)
      dh = dE$dh*0.5
      dIn1 = embed.backwd(network$W.in, dat$context[samp.idx,1],dh)
      dIn2 = embed.backwd(network$W.in, dat$context[samp.idx,2],dh)
      network$W.out[neg.set[,i.neg],] = network$W.out[neg.set[,i.neg],] - lambda*dE$dW
      network$W.in = network$W.in - lambda*dIn1
      network$W.in = network$W.in - lambda*dIn2
    }
}
par(mfrow=c(1,1))
plot(loss, type='l')

samp.idx = c(2:6,8,9,1)
batch.c1 = dat$context[samp.idx,1]
batch.c2 = dat$context[samp.idx,2]
batch.y = dat$target[samp.idx]
h1 <- embed.forwd(network$W.in,batch.c1)
h2 <- embed.forwd(network$W.in,batch.c2)
h = 0.5 * (h1 + h2)
s = MatMult.forwd(h,t(network$W.out))
z = sigmoid.forwd(s)
res=cbind(z,batch.y)
#par(mfrow=c(8,1))
#for (i in 1:8){
# col.spec = rep("black",8)
#  col.spec[i]="orange"
#  barplot(res[i, 1:8],col=col.spec)
#}

par(mfrow=c(4,1))
for (i in 1:4){
  col.spec = rep("black",8)
  col.spec[i]="orange"
  barplot(res[i, 1:8],col=col.spec)
}

par(mfrow=c(4,1))
for (i in 5:8){
  col.spec = rep("black",8)
  col.spec[i]="orange"
  barplot(res[i, 1:8],col=col.spec)
}

認知情報解析学演習 DL-2 ch.5

txt = "You say goodbye 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)
}

corp = txt2corpus(txt)

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))
}
dat<-corp2contxt1SRNN(corp)

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 = matrix(rnorm(size.hidden),nrow=1)*0.01
  return(list(W.h = W.h, W.x= W.x, b = b))
}

RNN.forward <- function(h.prev, x, network){
  b.size = nrow(x)
  h = h.prev%*%network$W.h + x.temp%*%network$W.x
  hb = h+matrix(1,nrow=b.size,ncol=1)%*%network$b
  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))
}

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

# checking modules
batch.size = 3; time = 3

corp.len = nrow(dat$y) 
size.hidden = 5
h.prev =array(0, c(batch.size, size.hidden, time))
h.next = array(0, c(batch.size, size.hidden, time))

for (i.t in 1:time){
  idx = seq(i.t, corp.len, batch.size)
  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]
  }
}

W.out = matrix(rnorm(size.hidden * 8), nrow = size.hidden)
b.out = matrix(rnorm(8), nrow = 1)
O = affine.forwd(h.next[,,3], W.out, b.out)
L = softmax.forwd(O, dat$y[c(3,6,9),])
ds = softmax.bckwd(O, dat$y[c(3,6,9),], 1)
dW.o = affine.bckwd(O, W.out, ds)

院:認識情報解析学演習

txt = "You say goodbye and I say hello.";txt = tolower(txt)
txt = gsub('[.]', ' .',txt)
words = unlist(strsplit(txt, " "))
uniq.words = unique(words)

n.uniq = length(uniq.words)
n.words = length(words)

# creating corpus
corpus = rep(0,n.words)
corpus = match(words,uniq.words)

## co-occurance matrix
cc = matrix(0,nrow=n.uniq, ncol=n.uniq)
for (i.c in 1:n.uniq){
  loc = which(corpus==i.c)
  L = which(match(uniq.words,words[pmax(loc-1,1)])!='NA')
  R = which(match(uniq.words,words[pmin(loc+1,n.words)])!='NA')
  cc[i.c, c(L,R)]=cc[i.c, c(L,R)]+1
}
diag(cc)=0

CosSim = function(vec1,vec2){
  eps = 1e-8
  similarity = vec1%*%vec2/((sqrt(sum(vec1))+eps)*(sqrt(sum(vec2))+eps))
  return(similarity)
}

ppmi <- function(cc, eps = 1e-8){
  n.uniq = ncol(cc)
  PPMI = matrix(0, n.uniq, n.uniq)
  N = sum(cc)
  r.sum = rowSums(cc)
  pmi = log2(cc*N/outer(r.sum,r.sum))
  PPMI = pmi
  PPMI[pmi<0]=0
  return(PPMI)
}

txt = "You say goodbye and I say hello.";
txt2corpus <- function(txt){
  txt = tolower(txt)
  txt = gsub('[.]', ' .',txt)
  words = unlist(strsplit(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)
}

corp2contxt1S = function(corpus){
  len.corp = length(corpus)
  # creating target matrix
  idxT = cbind(1:(len.corp-2), corpus[2:(len.corp-1)])
  targ1S = matrix(0,nrow=len.corp-2,ncol=length(unique(corpus)))
  targ1S[idxT]=1
  
  # creating context matrices
  idxC1 = cbind(1:(len.corp-2),corpus[1:(len.corp-2)])
  idxC2 = cbind(1:(len.corp-2),corpus[3:len.corp])
  contxt1 = matrix(0,nrow=len.corp-2,ncol=length(unique(corpus)))
  contxt2 = matrix(0,nrow=len.corp-2,ncol=length(unique(corpus)))
  contxt1[idxC1]=1
  contxt2[idxC2]=1
  return(list(target=targ1S,contxt1=contxt1,contxt2=contxt2))
}

corp = txt2corpus(txt)
dat<-corp2contxt1S(corp)

MatMult.forwd <- function(x, W){
  return(x%*%W)
}

MatMult.bckwd <- function(x, W, dout){
  dx = dout%*%t(W)
  dW = t(x)%*%dout
  return(list(dx = dx, dW = dW))
}

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

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

### 実装例
w2v.ws1 <- function(network, dat, n.batch, n.iter){
  n.data = nrow(dat$target)
  lambda=0.1;
  loss = rep(0, n.iter)
  for (i.iter in 1:n.iter){
    samp.idx = sample(1:n.data,n.batch)
    batch.c1 = dat$contxt1[samp.idx,]
    batch.c2 = dat$contxt2[samp.idx,]
    batch.y = dat$target[samp.idx,]
    h1 <- MatMult.forwd(batch.c1, network$W.in)
    h2 <- MatMult.forwd(batch.c2, network$W.in)
    h = 0.5 * (h1 + h2)
    s = MatMult.forwd(h,network$W.out)
    z = softmax.forwd(s,batch.y)
    loss[i.iter] = z
    ds = softmax.bckwd(s, batch.y, 1)
    da = MatMult.bckwd(h,network$W.out,ds)
    da$dW = da$dW*0.5
    dIn1 = MatMult.bckwd(batch.c1,network$W.in,da$dx)
    dIn2 = MatMult.bckwd(batch.c2,network$W.in,da$dx)
    network$W.out = network$W.out - lambda*da$dW
    network$W.in = network$W.in - lambda*dIn1$dW
    network$W.in = network$W.in - lambda*dIn2$dW
  }
  return(list(network = network, loss = loss))
}

init.W2V <- function(n.uniq,size.hidden){
  W.in = matrix(rnorm(n.uniq*size.hidden),nrow=n.uniq)*0.01
  W.out = matrix(rnorm(n.uniq*size.hidden),nrow=size.hidden)*0.01
  return(list(W.in = W.in, W.out= W.out))
}

txt = "You say goodbye and I say hello.";
corp = txt2corpus(txt)
dat <- corp2contxt1S(corp)
network<-init.W2V(7,5)
res <- w2v.ws1(network,dat,3,1000)
plot(res$loss,type='l')

広域システム特別講義II ベイズ統計

library(rjags)
source("http://peach.l.chiba-u.ac.jp/course_folder/HDI_revised.txt")

island.hopping <- function(n.rep=1e5, init.st=4) {
# example from DBDA 2nd Ed. (Kruschke) ch. 07
# intro to MCMC, island hopping
  
  # initialization
  state = 1:7
  curr.st = init.st
  state.counter = rep(0,7)
  state.counter[curr.st] = 1
  state.history=rep(0, n.rep)
  state.history[1]=curr.st
  # main loop
  for (i.rep in 2:n.rep) {
    move = sample(c(-1, 1), 1)
    prob.move = min((curr.st + move)/curr.st, 1)
    if (runif(1) < prob.move) {
      curr.st = curr.st + move
      if (curr.st > 7){ curr.st = 7 }
    } 
    state.history[i.rep] = curr.st
    state.counter[curr.st] = state.counter[curr.st]+1
  }
  par(mfrow=c(2, 1))
  par(mai=c(0, 1, 1, 1))
  par(mar=c(4, 3, 1, 3))
  barplot(state.counter, xlab="theta", ylab="Frequency")
  plot(state.history, 1:n.rep, type='o', log='y', xlab="theta", ylab="Time Step (log)")
}

island.hopping(10000, 4)
  

metropolis.ex01 <- function(n.iter=1e6, theta.init=0.1, sigma=0.2, plot.yn = T){
  # example from DBDA 2nd Ed. (Kruschke) ch. 07
  # metropolis algo. with 1 parameter 
  
  # "posterior of theta" function
  posterior.theta <- function(theta, N, z, a, b) {
    posterior = theta^z * (1-theta)^(N-z) * theta^(a-1) * (1 - theta)^(b-1) / beta(a,b)
  }
  
  # initialization
  theta.history = rep(0,nrow = n.iter,ncol = 1)
  theta.current = theta.init
  theta.history[1] = theta.current
  # values given in text
  mu = 0
  N = 20
  z = 14
  a = 1
  b = 1
  
  # main loop
  for (i_iter in 2:n.iter) {
    theta.proposed = theta.current + rnorm(1, mean=mu, sd=sigma)
    if (theta.proposed < 0 | theta.proposed > 1) {
      p.move = 0 
    } else {
      p.current = posterior.theta(theta.current, N, z, a, b)
      p.proposed = posterior.theta(theta.proposed, N, z, a, b)
      p.move = min(p.proposed/p.current, 1)
    }
    if (runif(1) < p.move) {
      theta.current = theta.proposed
    }
    theta.history[i_iter] = theta.current
  }
  
  # plotting results
  if (plot.yn == T) {
    par(mfrow = c(3, 1))
    hist(theta.history, nclass = 100, col = "orange", probability = T, xlab = "theta")
    den=density(theta.history)
    lines(den)
    plot(theta.history[(n.iter-100):n.iter], ((n.iter-100):n.iter), type = 'o', xlim = c(0,1), xlab="theta", ylab = "step in chain")
    plot(theta.history[1:100],1:100, type = 'o', xlim = c(0,1), xlab = "theta", ylab = "step in chain")
  }
  return(theta.history)
}

res = metropolis.ex01(10000, 0.1)

# initialization
n.iter = 10000; sigma = 0.1; counter = 0
z = c(6, 2); N = c(8, 7)
a = c(2, 2); b = c(2, 2); 
n.par = 2
th.hist = matrix(0, nrow = n.iter*n.par, ncol = n.par)
theta = runif(2)

# function to calc. prob. move
prob.gibbs <- function(theta, proposed, N, z, a, b, idx){
  p.th=dbeta(theta[idx], z[idx]+a[idx], N[idx]-z[idx]+b[idx])
  p.pro=dbeta(proposed, z[idx]+a[idx], N[idx]-z[idx]+b[idx])
  return(p.pro/p.th)
}

# main loop
for (i.iter in 1:n.iter){
  for (i.par in 1:n.par){
    proposed = theta[i.par] + rnorm(1, mean=0, sd=sigma)
    if (proposed > 1) {proposed = 1}
    if (proposed < 0) {proposed = 0}
    p.move= min(1, prob.gibbs(theta, proposed, N, z, a, b, i.par))
    if (runif(1) < p.move){
      theta[i.par] = proposed
    } 
    counter = counter + 1
    th.hist[counter, ] = theta
  }
}
par(mfrow=c(3,1))
HDI.plot(th.hist[,1])
HDI.plot(th.hist[,2])
plot(th.hist[,1],th.hist[,2], type='o',cex=0.1,xlim = c(0,1),ylim=c(0,1))
par(mfrow=c(1,1))

model.txt = "
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dbern( theta )
  }
  theta ~ dbeta( 1, 1 )
}"
writeLines(model.txt, "model.txt")

dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/z15N50.csv")
y=dat$y
Ntotal=length(dat$y)
datalist = list(y=y,Ntotal=Ntotal)

# jags
jagsModel = jags.model(file="model.txt",data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)

par(mfrow=c(2,2))
cols=c("orange", "skyblue","pink")

# chain
mcmc1<-as.mcmc(codaSamples[[1]])
mcmc2<-as.mcmc(codaSamples[[2]])
mcmc3<-as.mcmc(codaSamples[[3]])
plot(mcmcMat[1:5000],type='l')
lines(mcmcMat[5001:10000],col='red')
lines(mcmcMat[10001:15000],col='blue')

# autocorrelation
ac1=autocorr(mcmc1,lags=0:50)
ac2=autocorr(mcmc2,lags=0:50)
ac3=autocorr(mcmc3,lags=0:50)
plot(ac1, type='o', pch = 20, col = cols[1], ylab = "Autocorrelation", xlab = "Lag")
lines(ac2, type='o', pch = 20, col = cols[2])
lines(ac3, type='o', pch = 20, col = cols[3])

# shrink factor
resALL=mcmc.list(mcmc1,mcmc2,mcmc3)
gd1=gelman.plot(resALL, auto.layout = F)

# density
den1=density(mcmc1)
den2=density(mcmc2)
den3=density(mcmc3)
plot(den1, type='l', col = cols[1], main = " ", xlab = "param. value",lwd=2.5)
lines(den2, col = cols[2], lwd=2.5)
lines(den3, col = cols[3], lwd=2.5)

par(mfrow=c(1,1))

model.txt = "
model {
  for ( i_data in 1:Ntotal ) {
    y[ i_data ] ~ dbern( theta[s[i_data]] )
  }
  for ( i_s in 1:Nsubj) {
    theta[i_s] ~ dbeta( 2, 2 )
  }
}"
writeLines(model.txt, "model.txt")

dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/z6N8z2N7.csv")
y=dat$y
s=as.numeric(dat$s)
Ntotal=length(dat$y)
Nsubj=length(unique(s))
datalist = list(y=y,s=s,Ntotal=Ntotal,Nsubj=Nsubj)

jagsModel = jags.model(file="model.txt",data=datalist,n.chains=3,n.adapt=500)
update(jagsModel,n.iter=1000)
codaSamples=coda.samples(jagsModel,variable.names=c("theta"),n.iter=5000)
mcmcMat<-as.matrix(codaSamples)
par(mfrow=c(2,2))
cols=c("orange", "skyblue","pink")

# chain
mcmc1<-as.mcmc(codaSamples[[1]])
mcmc2<-as.mcmc(codaSamples[[2]])
mcmc3<-as.mcmc(codaSamples[[3]])
plot(temp1,type='l')
lines(temp2,col='red')
lines(temp3,col='blue')

# autocorrelation
ac1=autocorr(mcmc1,lags=0:50)
ac2=autocorr(mcmc2,lags=0:50)
ac3=autocorr(mcmc3,lags=0:50)
plot(ac1, type='o', pch = 20, col = cols[1], ylab = "Autocorrelation", xlab = "Lag")
lines(ac2, type='o', pch = 20, col = cols[2])
lines(ac3, type='o', pch = 20, col = cols[3])

# shrink factor
resALL=mcmc.list(mcmc1,mcmc2,mcmc3)
gd1=gelman.plot(resALL, auto.layout = F)

# density
den1=density(mcmc1)
den2=density(mcmc2)
den3=density(mcmc3)
plot(den1, type='l', col = cols[1], main = " ", xlab = "param. value",lwd=2.5)
lines(den2, col = cols[2], lwd=2.5)
lines(den3, col = cols[3], lwd=2.5)

par(mfrow=c(1,1))

### batting
dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/BattingAverage.csv")
N=dat$AtBats
z=dat$Hits
s=dat$PlayerNumber
Nsubj=length(unique(s))
c=as.numeric(dat$PriPos)
Ncat=length(unique(c))
dataList=list(z=z,N=N,c=c,Nsubj=Nsubj,Ncat=Ncat)
model.txt = "
model {
    for ( s in 1:Nsubj ) {
      z[s] ~ dbin( theta[s] , N[s] )
      theta[s] ~ dbeta( omega[c[s]]*(kappa[c[s]]-2)+1 , (1-omega[c[s]])*(kappa[c[s]]-2)+1 ) 
    }
    for (c in 1:Ncat) {
      omega[c] ~ dbeta(omegaO*(kappaO-2)+1 , (1-omegaO)*(kappaO-2)+1 )
      kappa[c] <- kappaMinusTwo[c] + 2
      kappaMinusTwo[c] ~ dgamma( 0.01 , 0.01 )  
    }
    omegaO ~ dbeta( 1 , 1 )
    kappaO <- kappaMinusTwoO + 2
    kappaMinusTwoO ~ dgamma( 0.01 , 0.01 )  
 }
"
writeLines(model.txt, "modelAVE.txt")

jagsModel=jags.model(file="ModelAVE.txt",data=dataList,n.chain=3,n.adapt=500)
update(jagsModel,n.iter=500)
codaSample=coda.samples(jagsModel,variable.names=c("theta","omega","kappa","omegaO","kappaO"),n.iter=10000)


model.txt = "
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dnorm( mu , 1/sigma^2  )
  }
  mu ~ dnorm( meanY , 1/(100*sdY)^2 )
  sigma ~ dunif( sdY/1000 , sdY*1000 )
}
"
writeLines(model.txt, "model.txt")

dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/TwoGroupIQ.csv")
y = dat$Score[dat$Group=="Smart Drug"]
Ntotal = length(y)

dataList = list(y = y, Ntotal = Ntotal, meanY = mean(y), sdY = sd(y))
jagsModel = jags.model("model.txt", data=dataList, n.chains=3, n.adapt=1000 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=c("mu","sigma"), n.iter=5000, thin=5)
mcmcMat<-as.matrix(codaSamples)

# calculating & plotting normality
par(mfrow=c(2,2))
HDI.plot(mcmcMat[,1])
hist(y,nclass=15,probability = T)
x.temp = seq(40,200,0.1)
n.plot = 100
randperm = sample(1:nrow(mcmcMat),n.plot)
for (i.plot in 1:n.plot){
  norm.temp=dnorm(x.temp,mean=mcmcMat[randperm[i.plot],1],sd=mcmcMat[randperm[i.plot],2])
  lines(x.temp,norm.temp,col='orange')
}
HDI.plot(mcmcMat[,2])

# calculating & plotting effect size
effect.size=(mcmcMat[,"mu"]-100)/mcmcMat[,"sigma"]
HDI.plot(effect.size)


#
dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/TwoGroupIQ.csv")
y = dat$Score[dat$Group=="Smart Drug"]
Ntotal = length(y)

model.txt="
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dt( mu , 1/sigma^2 , nu )
  }
  mu ~ dnorm( meanY , 1/(100*sdY)^2 )
  sigma ~ dunif( sdY/1000 , sdY*1000 )
  nu <- nuMinusOne+1
  nuMinusOne ~ dexp(1/30.0)
}"
writeLines(model.txt, "model.txt")

dataList = list(y = y, Ntotal = Ntotal, meanY = mean(y), sdY = sd(y))
jagsModel = jags.model("model.txt", data=dataList, n.chains=3, n.adapt=1000 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=c("mu","sigma","nu"), n.iter=5000, thin=5)
mcmcMat<-as.matrix(codaSamples)

# calculating & plotting normality
par(mfrow=c(3,2))
HDI.plot(mcmcMat[,1])
HDI.plot(mcmcMat[,3])
normality=log10(mcmcMat[,"nu"])
HDI.plot(normality)
effect.size=(mcmcMat[,"mu"]-100)/mcmcMat[,"sigma"]
HDI.plot(effect.size)

hist(y,nclass=20,probability = T)
n.plot = 100
randperm = sample(1:nrow(mcmcMat),n.plot)
for (i.plot in 1:n.plot){
  x.temp1 = seq(40,200,0.1)
  x.temp2 = (x.temp1 - mcmcMat[randperm[i.plot],1])/mcmcMat[randperm[i.plot],3]
  t.temp=dt(x.temp2,mcmcMat[randperm[i.plot],2])/mcmcMat[randperm[i.plot],3]
  lines(x.temp1,t.temp,col='orange')
}

par(mfrow=c(2,2))
plot(mcmcMat[,1],mcmcMat[,3],col='orange')
plot(mcmcMat[,1],log10(mcmcMat[,2]),col='orange')
plot(0,0,type='n')
plot(log10(mcmcMat[,2]),mcmcMat[,3],col='orange')


dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/TwoGroupIQ.csv")
y = dat$Score
group = as.numeric(dat$Group)
Ntotal = length(y)
Ngroup = length(unique(group))
model.txt="
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dt( mu[group[i]] , 1/sigma[group[i]]^2 , nu )
  }
  for (j in 1:Ngroup){
    mu[j] ~ dnorm( meanY , 1/(100*sdY)^2 )
    sigma[j] ~ dunif( sdY/1000 , sdY*1000 )
  }      
  nu <- nuMinusOne+1
  nuMinusOne ~ dexp(1/29)
}"
writeLines(model.txt, "model.txt")


dataList = list(y = y, Ntotal = Ntotal, meanY = mean(y), sdY = sd(y),Ngroup=Ngroup,group=group)
jagsModel = jags.model("model.txt", data=dataList, n.chains=3, n.adapt=5000 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=c("mu","sigma","nu"), n.iter=5000, thin=5)
mcmcMat<-as.matrix(codaSamples)

# plotting result
par(mfrow=c(5,2))
HDI.plot(mcmcMat[,1],xlabel="placebo Mean")

hist(y[dat$Group=="Placebo"],nclass=20,probability = T)
n.plot = 100
randperm = sample(1:nrow(mcmcMat),n.plot)
for (i.plot in 1:n.plot){
  x.temp1 = seq(40,200,0.1)
  x.temp2 = (x.temp1 - mcmcMat[randperm[i.plot],1])/mcmcMat[randperm[i.plot],4]
  t.temp=dt(x.temp2,mcmcMat[randperm[i.plot],3])/mcmcMat[randperm[i.plot],4]
  lines(x.temp1,t.temp,col='orange')
}

HDI.plot(mcmcMat[,2],xlabel="smart drug Mean")

hist(y[dat$Group=="Smart Drug"],nclass=20,probability = T)
n.plot = 100
randperm = sample(1:nrow(mcmcMat),n.plot)
for (i.plot in 1:n.plot){
  x.temp1 = seq(40,200,0.1)
  x.temp2 = (x.temp1 - mcmcMat[randperm[i.plot],2])/mcmcMat[randperm[i.plot],5]
  t.temp=dt(x.temp2,mcmcMat[randperm[i.plot],3])/mcmcMat[randperm[i.plot],5]
  lines(x.temp1,t.temp,col='orange')
}

HDI.plot(mcmcMat[,4],xlabel="placebo scale")

HDI.plot(mcmcMat[,2]-mcmcMat[,1],xlabel="Difference of Means")

HDI.plot(mcmcMat[,5],xlabel="smart drug scale")

HDI.plot(mcmcMat[,5]-mcmcMat[,4],xlabel="Difference of Scales")

HDI.plot(log10(mcmcMat[,3]),xlabel="Normality")

effect.size = (mcmcMat[,2]-mcmcMat[,1])/sqrt((mcmcMat[,5]^2+mcmcMat[,4]^2)/2)
HDI.plot(effect.size,xlabel="effect size")

Posted in UT

データ解析基礎論B 人工神経回路

library(nnet)
dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
set.seed(5)
x = dat[,1:3]
y = dat[,4]
dat.nnet = nnet(x,y, size = 150, linout= TRUE, maxit = 1000)
nnet.pred <- predict(dat.nnet,dat)
cor(dat.nnet$fitted.values,dat$sales)^2

n.data = nrow(dat);
n.sample = n.data*0.6; 
n.rep = 100
trainNN.cor = rep(0,n.rep); trainLM.cor = rep(0,n.rep)
testNN.cor = rep(0,n.rep); testLM.cor = rep(0,n.rep)
for (i.rep in 1:n.rep){  
  randperm = sample(1:n.data)  
  train.idx = randperm[1:n.sample]  
  test.idx = randperm[(n.sample+1):n.data]  
  dat.nnet <- nnet(sales~.,size = c(10), linout=T,decay= 0.01, maxit=1000, data = dat[train.idx,])  
  dat.lm <-lm(sales~.,data=dat[train.idx, ])  
  trainNN.cor[i.rep] <- cor(predict(dat.nnet,dat[train.idx, ]), dat[train.idx,]$sales)  
  trainLM.cor[i.rep] <- cor(predict(dat.lm,dat[train.idx, ]), dat[train.idx,]$sales)  
  testNN.cor[i.rep] <- cor(predict(dat.nnet,dat[test.idx, ]), dat[test.idx,]$sales)  
  testLM.cor[i.rep] <- cor(predict(dat.lm,dat[test.idx, ]), dat[test.idx,]$sales)
}
print(c(mean(trainNN.cor,na.rm=T),mean(testNN.cor,na.rm=T)))
print(c(mean(trainLM.cor,na.rm=T),mean(testLM.cor,na.rm=T)))
print(c(max(trainNN.cor,na.rm=T),max(testNN.cor,na.rm=T)))
print(c(max(trainLM.cor,na.rm=T),max(testLM.cor,na.rm=T)))
print(c(min(trainNN.cor,na.rm=T),min(testNN.cor,na.rm=T)))
print(c(min(trainLM.cor,na.rm=T),min(testLM.cor,na.rm=T)))

dat<-read.csv("http://matsuka.info/data_folder/tdkDA01.csv", header=T)
dat.nnet<-nnet(class~.,dat,size=30,maxit=1000,decay=0.05)
dat.pred<-predict(dat.nnet,dat,type="class")
table(dat.pred,dat$class)

dat.glm <- glm(class~., family="binomial",dat)
glm.pred <- predict(dat.glm, dat, type="response") < 0.5
table(glm.pred,dat$class)

dat<-read.csv("http://www.matsuka.info/data_folder/cda7-16.csv")
wts = rep(1,nrow(dat))
wts[which(dat$survival=="no")]=45
dat.nnet<-nnet(survival~., weights=wts, dat, size=100, maxit=1000, decay=0.1)
dat.pred<-predict(dat.nnet,dat,type="class")
table(dat.pred,dat$survival)

dat<-read.csv("http://matsuka.info/data_folder/tdkDA02.csv",header=T)
class.id<-class.ind(dat$class)
x = dat[,1:6]
dat.nnet<-nnet(x, class.id, size=30, maxit=1000, decay=0.01, softmax=TRUE)
max.id = apply(dat.nnet$fitted.values,1,which.max)
table(max.id,dat$class)

dat<-read.table("http://www.matsuka.info/data_folder/tdkPCA01.txt")
dat.nnet<-nnet(dat,dat,size=2, maxit=1000, decay=0.01, linout=TRUE)

院:認知情報解析演習

source("http://peach.l.chiba-u.ac.jp/course_folder/HDI_revised.txt")
dat <- read.csv("http://peach.l.chiba-u.ac.jp/course_folder/detectFakeSmile.csv")
dataList = list(y = dat$resp, x = dat$eye, s=dat$s, 
      Ntotal = nrow(dat), Nsubj=length(unique(dat$s)))

model.txt = "
  data {
    xm <- mean(x)
    xsd <- sd(x)
    for (i in 1:Ntotal){
      zx[i] = (x[i] - xm)/xsd
    }
  }
  model {
  for ( i in 1:Ntotal ) {
    y[ i ] ~ dbern(ilogit( zbeta0[s[i]] + zbeta1[s[i]] * zx[i]))
  }
  for (j in 1:Nsubj){
    zbeta0[j] ~ dnorm(0, 1/2^2)
    zbeta1[j] ~ dnorm(0, 1/2^2)
  }
  for (j in 1:Nsubj){
    beta1[j] <- zbeta1[j] / xsd
    beta0[j] <- zbeta0[j] - zbeta1[j] * xm / xsd
  }
}"
writeLines(model.txt, "model.txt")
parameters = c( "beta0" ,  "beta1")
jagsModel = jags.model( "model.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel , n.iter=1000)
codaSamples = coda.samples( jagsModel , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat<-as.matrix(codaSamples)
n.plot = 100
x= seq(30,200, length.out = 100)
eye = matrix(dat$eye, nrow=20)
resp = matrix(dat$resp, nrow=20) 
par(mfrow=c(5,4))
for (i.s in 1:n.s){
  randsamp = sample(1:nrow(mcmcMat),n.plot)
  plot(eye[i.s,],resp[i.s,],ylim=c(-0.1,1.1),xlim=c(30,200))
  for (i.plot in 1:n.plot){
    lines(x,1/(1+exp(-(mcmcMat[randsamp[i.plot],"beta0[20]"]+mcmcMat[randsamp[i.plot],"beta1[20]"]*x))),col='orange')
  }
}

model2.txt = "
  data {
    xm <- mean(x)
    xsd <- sd(x)
    for (i in 1:Ntotal){
      zx[i] = (x[i] - xm)/xsd
      }
    }
  model {
    for ( i in 1:Ntotal ) {
      y[ i ] ~ dbern(ilogit( zbeta0[s[i]] + zbeta1[s[i]] * zx[i]))
    }
    for (j in 1:Nsubj){
      zbeta0[j] ~ dnorm( zbeta0mu, 1/(zbeta0sigma)^2)
      zbeta1[j] ~ dnorm( zbeta1mu, 1/(zbeta1sigma)^2)
    }

    zbeta0mu ~ dnorm(0, 1/(10)^2)
    zbeta1mu ~ dnorm(0, 1/(10)^2)
    zbeta0sigma ~ dunif(1.0E-3, 1.0E+3)
    zbeta1sigma ~ dunif(1.0E-3, 1.0E+3)

    for (j in 1:Nsubj){
      beta1[j] <- zbeta1[j] / xsd
      beta0[j] <- zbeta0[j] - zbeta1[j] * xm / xsd
    }
    beta1mu <- zbeta1mu / xsd
    beta0mu <- zbeta0mu  -zbeta1mu * xm / xsd
}"
writeLines(model2.txt, "model2.txt")
parameters = c( "beta0" ,  "beta1", "beta0mu", "beta1mu")
library(rjags)
jagsModel2 = jags.model( "model2.txt", data=dataList, n.chains=3, n.adapt=500 )
update( jagsModel2 , n.iter=1000)
codaSamples2 = coda.samples( jagsModel2 , variable.names=parameters, n.iter=10000, thin=5)
mcmcMat2 <- as.matrix(codaSamples2)
n.plot = 100
x= seq(30,200, length.out = 100)
par(mfrow=c(5,4))
for (i.s in 1:n.s){
  randsamp = sample(1:nrow(mcmcMat),n.plot)
  plot(eye[i.s,],resp[i.s,],ylim=c(-0.1,1.1),xlim=c(30,200))
  for (i.plot in 1:n.plot){
    lines(x,1/(1+exp(-(mcmcMat2[randsamp[i.plot],"beta0[20]"]+mcmcMat2[randsamp[i.plot],"beta1[20]"]*x))),col='orange')
  }
}
par(mfrow=c(1,1))
HDI.plot(mcmcMat2[,"beta1mu"])

広域システム特別講義II W06

# 課題の解答例
ALC.init<-function(size) {
  # size[1] = input dimension
  # size[2] = number of exemplars
  # size[3] = number of categories
  alpha=matrix(1,nrow=1,ncol=size[1])/size[1]
  w=matrix(rnorm(size[3]*size[2])*0.1,ncol=size[2])
  c=2            # gain
  phi=3          # decision strength
  lrA=0.2        # learning rate for attentions
  lrW=0.1        # learning rate for associations
  return(list(alpha = alpha,
              w = w,
              lrA = lrA,
              lrW = lrW,
              c = c,
              phi = phi))
}

# alcove forward process
alcove.forward <- function(parmSet, inp, exemplar){
  diff= matrix(abs(matrix(1,n.exemp,1)%*%inp-exemplar),nrow=nrow(exemplar))
  exemp=exp(-parmSet$c*(diff%*%t(parmSet$alpha)))
  out=parmSet$w%*%exemp
  return(list(diff = diff, exemp = exemp, out = out))
}

# alcove backward process
alcove.backward <- function(res.forward,parmSet, target){
  err=target-res.forward$out
  dW=parmSet$lrW*res.forward$exemp%*%t(err)
  dA=-parmSet$lrA*(t(err)%*%parmSet$w*t(res.forward$exemp))%*%res.forward$diff*parmSet$c
  return(list(dW = dW, dA = dA))
}

# main function
alcove<-function(parmSet, inp, exemplar, targ,iter) {
  # ----ALCOVE for R ---
  # input arguments
  #  parmSet - parameter set
  # inp - input matrix (n_sample x n_dimension)
  # exemplar - exemplar matrix (n_sample x n_dimension)
  # targ - target matrix (n_sample x n_category)
  # iter - # of training epochs
  # ----ALCOVE for R ---
  # initialization
  inpDim = ncol(inp)
  n.inp = nrow(inp)
  n.exemp = nrow(exemplar)
  n.cat = ncol(targ)
  accu=rep(0,nrow=iter+1)
  attn=matrix(0,nrow=iter+1,ncol=inpDim)
  attn[1,]=alpha
  # end initialization
  
  # main loop
  for (i_iter in 1:iter) {
    ro=sample(1:n.inp, n.inp)
    prob_temp=0;
    for (i_tr in 1:n.inp) {
      res.forward <- alcove.forward(parmSet, inp[ro[i_tr],], exemplar)
      res.backward <- alcove.backward(res.forward, parmSet, targ[ro[i_tr],])
      parmSet$w = parmSet$w+t(res.backward$dW)
      parmSet$alpha = parmSet$alpha+res.backward$dA
      parmSet$alpha[which(parmSet$alpha<0)]=0
      pT=(exp(parmSet$phi*res.forward$out)/sum(exp(parmSet$phi*res.forward$out)))*targ[ro[i_tr],]
      prob_temp=prob_temp+pT[which(!pT==0)]
    }
    accu[i_iter+1]=prob_temp/n.inp
    attn[i_iter+1,]=alpha
  }
  attnN=attn/apply(attn,1, sum)
  out=matrix(0,nrow=n.cat,ncol=n.inp)
  for (i_tr in 1:n.inp) {
    res.forward<-alcove.forward(parmSet, inp[i_tr,], exemplar)
    out[,i_tr]=res.forward$out
  }
  return(list(attn=attn,attnN=attnN,accu=accu,out=out,ps=parmSet))
}


exemplar = matrix(c(0,0,0,
                    0,0,1,
                    0,1,0,
                    0,1,1,
                    1,0,0,
                    1,0,1,
                    1,1,0,
                    1,1,1),byrow=T,nrow=8)
inp = exemplar

target1 = matrix(c(0,0,0,0,1,1,1,1,
                   1,1,1,1,0,0,0,0),nrow=8)

target2 = matrix(c(1,1,0,0,0,0,1,1,
                   0,0,1,1,1,1,0,0),nrow=8)

target3 = matrix(c(1,1,1,0,0,1,0,0,
                   0,0,0,1,1,0,1,1),nrow=8)

target4 = matrix(c(1,1,1,0,1,0,0,0,
                   0,0,0,1,0,1,1,1),nrow=8)

target5 = matrix(c(1,1,1,0,0,0,0,1,
                   0,0,0,1,1,1,1,0),nrow=8)

target6 = matrix(c(1,0,0,1,0,1,1,0,
                   0,1,1,0,1,0,0,1),nrow=8)

seed.num = 1;n.train = 50
set.seed(seed.num)
parmSet<-ALC.init(c(3,8,2))
result1<-alcove(parmSet,inp,exemplar,target1,n.train)
plot(result1$accu,type='o',lwd=2,pch=20,col=1)
set.seed(seed.num)
parmSet<-ALC.init(c(3,8,2))
result2<-alcove(parmSet,inp,exemplar,target2,n.train)
lines(result2$accu,type='o',lwd=2,pch=20,col=2)
set.seed(seed.num)
parmSet<-ALC.init(c(3,8,2))
result3<-alcove(parmSet,inp,exemplar,target3,n.train)
lines(result3$accu,type='o',lwd=2,pch=20,col=3)

set.seed(seed.num)
parmSet<-ALC.init(c(3,8,2))
result4<-alcove(parmSet,inp,exemplar,target4,n.train)
lines(result4$accu,type='o',lwd=2,pch=20,col=4)
set.seed(seed.num)
parmSet<-ALC.init(c(3,8,2))
result5<-alcove(parmSet,inp,exemplar,target5,n.train)
lines(result5$accu,type='o',lwd=2,pch=20,col=5)
set.seed(seed.num)
parmSet<-ALC.init(c(3,8,2))
result6<-alcove(parmSet,inp,exemplar,target6,n.train)
lines(result6$accu,type='o',lwd=2,pch=20,col=6)
legend("bottomright",paste("type",1:6,sep=''),col=1:6,lwd=2,pch=20)
# word2vecなど
dat <- read.table("http://www.matsuka.info/data_folder/tdkPCA01.txt")
dat.pca <- princomp(dat)
dat.pca$score
dat.pca$loadings

library(nnet)
dat <- read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
set.seed(5)
x = dat[, 1:3]
y = dat[, 4]
dat.nnet = nnet(x, y, size = 150, linout = TRUE,maxit = 1000)
nnet.pred <- predict(dat.nnet, dat)
cor(dat.nnet$fitted.values,dat$sales)^2

dat.lm<-lm(sales~.,data=dat)
summary(dat.lm)

# pseudo-pca
dat<-read.table("http://www.matsuka.info/data_folder/tdkPCA01.txt")
dat.nnet<-nnet(dat,dat,size=2, maxit=1000, decay=0.01, linout=TRUE)
cor(dat.nnet$fitted.values,dat)

txt = "You say goodbye and I say hello.";txt = tolower(txt)
txt = gsub('[.]', ' .',txt)
words = unlist(strsplit(txt, " "))
uniq.words = unique(words)
corpus = rep(0,n.words)
corpus = match(words,uniq.words)

contxt <- function(corpus, uniq.words, words){
  cc = matrix(0, nrow=n.uniq, ncol=n.uniq)
  for (i.c in 1:n.uniq){
    loc = which(corpus==i.c)
    L = which(match(uniq.words, words[pmax(loc-1,1)])!='NA')
    R = which(match(uniq.words, words[pmin(loc+1,n.words)])!='NA')
    cc[i.c, c(L,R)]=cc[i.c, c(L,R)]+1
  }
  diag(cc)=0
  return(cc)
}

txt = "You say goodbye and I say hello.";
txt2corpus <- function(txt){
  txt = tolower(txt)
  txt = gsub('[.]', ' .',txt)
  words = unlist(strsplit(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)
}
corp = txt2corpus(txt)
corp2contxt <- function(corpus){
  len.corp = length(corpus)
  target = corpus[2:(len.corp-1)]
  col1 = corpus[1:(len.corp-2)]
  col2 = corpus[3:len.corp]
  context = cbind(col1,col2)
  return(list(context=context,target = target))
}

corp2contxt1S = function(corpus){
  len.corp = length(corpus)
  # creating target matrix
  idxT = cbind(1:(len.corp-2), corpus[2:(len.corp-1)])
  targ1S = matrix(0,nrow=len.corp-2,ncol=length(unique(corpus)))
  targ1S[idxT]=1
  
  # creating context matrices
  idxC1 = cbind(1:(len.corp-2),corpus[1:(len.corp-2)])
  idxC2 = cbind(1:(len.corp-2),corpus[3:len.corp])
  contxt1 = matrix(0,nrow=len.corp-2,ncol=length(unique(corpus)))
  contxt2 = matrix(0,nrow=len.corp-2,ncol=length(unique(corpus)))
  contxt1[idxC1]=1
  contxt2[idxC2]=1
  return(list(target=targ1S,contxt1=contxt1,contxt2=contxt2))
}


MatMult.forwd <- function(x, W){
  return(x%*%W)
}
init.W2V <- function(n.uniq,size.hidden){
  W.in = matrix(rnorm(n.uniq*size.hidden),nrow=n.uniq)*0.01
  W.out = matrix(rnorm(n.uniq*size.hidden),nrow=size.hidden)*0.01
  return(list(W.in = W.in, W.out= W.out))
}

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

MatMult.bckwd <- function(x, W, dout){
  dx = dout%*%t(W)
  dW = t(x)%*%dout
  return(list(dx = dx, dW = dW))
}

network<-init.W2V(7,5)
n.batch = 3;
n.data = nrow(dat$target)
n.iter = 2000;
lambda=0.1;
loss = rep(0, n.iter)
for (i.iter in 1:n.iter){
  samp.idx = sample(1:n.data,n.batch)
  batch.c1 = dat$contxt1[samp.idx,]
  batch.c2 = dat$contxt2[samp.idx,]
  batch.y = dat$target[samp.idx,]
  h1 <- MatMult.forwd(batch.c1, network$W.in)
  h2 <- MatMult.forwd(batch.c2, network$W.in)
  h = 0.5 * (h1 + h2)
  s = MatMult.forwd(h,network$W.out)
  z = softmax.forwd(s,batch.y)
  loss[i.iter] = z
  ds = softmax.bckwd(s, batch.y, 1)
  da = MatMult.bckwd(h,network$W.out,ds)
  da$dW = da$dW*0.5
  dIn1 = MatMult.bckwd(batch.c1,network$W.in,da$dx)
  dIn2 = MatMult.bckwd(batch.c2,network$W.in,da$dx)
  network$W.out = network$W.out - lambda*da$dW
  network$W.in = network$W.in - lambda*dIn1$dW
  network$W.in = network$W.in - lambda*dIn2$dW
}
plot(loss, type='l')
# more "effecient" way
embed.forwd <- function(W, idx){
  return(W[idx,])
}

embed.dot.forwd <- function(W, h, idx){
  return(O=rowSums(W[idx,]*h))
}

sigmoidWL.forwd <- function(O,target){
  y = 1/(1+exp(-O))
  loss=-sum(target*log(y)+(1-target)*log(1 - y))
  return(loss)
}

sigmoidWL.backwd <- function(O,target,dout=1){
  y = 1/(1+exp(-O)) 
  dW = (y - target)*dout
  return(dW)
}

embed.dot.backwd <- function(W,h,idx,dout){
  dh = dout * W[idx,]
  dW = dout * h
  return(list(dh = dh, dW = dW))
}

embed.backwd <- function(W, idx, dout){
  dW = matrix(0,nrow(W),ncol(W))
  for (i.idx in 1:length(idx)){
    dW[idx[i.idx], ] = dW[idx[i.idx], ] + dout[i.idx, ]
  }
  return(dW)
}

unigram.sampler <- function(corpus, target, power, sample.size){
  neg.corpus <- corpus[-which(corpus == target)]
  uniq.word = unique(neg.corpus)
  n.word = length(neg.corpus)
  prob = (as.vector(table(neg.corpus))/n.word)^power
  prob = prob/sum(prob)
  sample.idx = sample(uniq.word, sample.size, prob = prob)
  return(sort(sample.idx))
}

get.neg.sample <- function(corpus, target, power, sample.size){
  sample.id = matrix(0, nrow = length(target), ncol = sample.size)
  for (i.targ in 1:length(target)){
    sample.id[i.targ,] = unigram.sampler(corpus, target[i.targ], power, sample.size)
  }
  return(sample.id)
}

network<-init.W2V(8,5)
n.batch = 3;
txt = "You say goodbye and I say hello.."
corp = txt2corpus(txt)
corp = c(8,8,corp)
dat<-corp2contxt(corp)
n.data =  length(dat$target)

n.iter = 10000;
lambda=0.1;
loss = rep(0, n.iter)
sample.size = 2
for (i.iter in 1:n.iter){
  samp.idx = sample(1:n.data,n.batch)
  batch.c1 = dat$context[samp.idx,1]
  batch.c2 = dat$context[samp.idx,2]
  batch.y = dat$target[samp.idx]
  h1 <- embed.forwd(network$W.in,batch.c1)
  h2 <- embed.forwd(network$W.in,batch.c2)
  h = 0.5 * (h1 + h2)
  # positive
  s = embed.dot.forwd(network$W.out,h,batch.y)
  z = sigmoidWL.forwd(s,1)
  loss[i.iter] = z
  ds = sigmoidWL.backwd(s, 1, 1)
  dE = embed.dot.backwd(network$W.out,h, batch.y, ds)
  dh = dE$dh*0.5
  dIn1 = embed.backwd(network$W.in,dat$context[samp.idx,1], dh)
  dIn2 = embed.backwd(network$W.in,dat$context[samp.idx,2], dh)
  network$W.out[batch.y,] = network$W.out[batch.y,] - lambda*dE$dW
  network$W.in = network$W.in - lambda*dIn1
  network$W.in = network$W.in - lambda*dIn2
  # negative
  neg.set <- get.neg.sample(corp,batch.y, 0.75, sample.size)
  for (i.neg in 1:sample.size){
      s = embed.dot.forwd(network$W.out,h,neg.set[,i.neg])
      z = sigmoidWL.forwd(s,0)
      loss[i.iter] = loss[i.iter] + z
      ds = sigmoidWL.backwd(s, 0, dout=1)
      dE = embed.dot.backwd(network$W.out,h,neg.set[,i.neg], ds)
      dh = dE$dh*0.5
      dIn1 = embed.backwd(network$W.in, dat$context[samp.idx,1],dh)
      dIn2 = embed.backwd(network$W.in, dat$context[samp.idx,2],dh)
      network$W.out[neg.set[,i.neg],] = network$W.out[neg.set[,i.neg],] - lambda/sample.size*dE$dW
      network$W.in = network$W.in - lambda*dIn1
      network$W.in = network$W.in - lambda*dIn2
    }
}
par(mfrow=c(1,1))
plot(loss, type='l')

samp.idx = c(2:6,8,9,1)
batch.c1 = dat$context[samp.idx,1]
batch.c2 = dat$context[samp.idx,2]
batch.y = dat$target[samp.idx]
h1 <- embed.forwd(network$W.in,batch.c1)
h2 <- embed.forwd(network$W.in,batch.c2)
h = 0.5 * (h1 + h2)
s = MatMult.forwd(h,t(network$W.out))
z = sigmoid.forwd(s)
res=cbind(z,batch.y)

par(mfrow=c(4,1))
for (i in 1:4){
  col.spec = rep("black",8)
  col.spec[i]="orange"
  barplot(res[i, 1:8],col=col.spec)
}

par(mfrow=c(4,1))
for (i in 5:8){
  col.spec = rep("black",8)
  col.spec[i]="orange"
  barplot(res[i, 1:8],col=col.spec)
}
Posted in UT

データ解析基礎論B 判別分析

means = matrix(c(52.1,43.3,25.9,74.8,71,23.4),nrow=3,byrow=T)
var = array(c(355.7,203.8,203.8,252.5,180.4,-198.8,-198.8,359.4,435.1,212.6,212.6,166.4),dim=c(2,2,3))
library(plot3D)
library(matrixcalc)
x1 = seq(0,100,length.out=100)
M <- mesh(x1,x1)
inv.var1<-matrix.inverse(var[,,1])
inv.var2<-matrix.inverse(var[,,2])
inv.var3<-matrix.inverse(var[,,3])
dens1=M$x
dens2=M$x
dens3=M$x
for (i.r in 1:100){
  for (i.c in 1:100){
    dens1[i.r, i.c]=(c(M$x[i.r,i.c],M$y[i.r,i.c])- means[1,])%*%inv.var1%*%(rbind(M$x[i.r,i.c],M$y[i.r,i.c])- means[1,])
    dens2[i.r, i.c]=(c(M$x[i.r,i.c],M$y[i.r,i.c])- means[2,])%*%inv.var2%*%(rbind(M$x[i.r,i.c],M$y[i.r,i.c])- means[2,])
    dens3[i.r, i.c]=(c(M$x[i.r,i.c],M$y[i.r,i.c])- means[3,])%*%inv.var3%*%(rbind(M$x[i.r,i.c],M$y[i.r,i.c])- means[3,])
 }
}
contour(exp(-dens1),drawlabels = F, col='red')
contour(exp(-dens2),add=T,drawlabels =F,col='blue')
contour(exp(-dens3),add=T,drawlabels =F,col='green')
points(0.4,0.6,pch=18,cex=4)
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)
plot(dat$writing,dat$interview,xlab="writing",ylab="review",col=dat$cl,pch=20,cex=3,xlim=c(20,110),ylim=c(20,110))
abline(a=dat.lda$scaling[2],b=dat.lda$scaling[1])
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
a=intcpt/dat.lda$scaling[2]
b=-dat.lda$scaling[1]/dat.lda$scaling[2]
b2=dat.lda$scaling[2]/dat.lda$scaling[1]
x=0:120
lines(x,a+b*x,col='red')
lines(x,-60+b2*x,col='blue')

dat<-read.csv("http://matsuka.info/data_folder/tdkDA01.csv", header=T)
dat.lda<-lda(class~.,dat)
lda.pred<-predict(dat.lda,dat)
table(lda.pred$class, dat$class)

dat.lda<-lda(class~.,dat, CV=T)
cbind(dat.lda$posterior,dat$class)

dat<-read.csv("http://matsuka.info/data_folder/tdkDA02.csv",header=T)
dat.lda=lda(class~.,data=dat)
lda.pred<-predict(dat.lda,dat)
table(lda.pred$class, dat$class)
plot(dat.lda, dimen=2, col=as.numeric(lda.pred$class), cex=3) 
plot(dat.lda, dimen=3, col=as.numeric(lda.pred$class),cex=2)

dat.km<-kmeans(dat[,1:6],5)
table(lda.pred$class,dat.km$cluster)

dat<-read.csv("http://www.matsuka.info/data_folder/tdkCFA.csv")
dat1<-subset(dat,dat$popularity<5)
dat2<-subset(dat,dat$popularity>4 & dat$popularity<6)
dat3<-subset(dat,dat$popularity>6)
dat1$popularity="LP";dat2$popularity="MP";dat3$popularity="VP"
datT=rbind(dat1,dat2,dat3)
datT.lda<-lda(popularity~.,datT)

datT.pred<-predict(datT.lda,datT)
table(datT.pred$class,datT$popularity)   
plot(datT.lda,col=c(rep('red',20),rep('blue',28),rep('green',29)),cex=1.5)
plot(datT.lda,col=as.numeric(datT.pred$class),cex=1.5)
plot(datT.lda,col=as.numeric(datT.pred$class)+2,cex=1)
colors=rep("black",nrow(datT))
miss.idx = which(datT.pred$class!= datT$popularity)
colors[miss.idx]='red'
points(datT.pred$x,pch=20,col=colors)
legend("bottomright",c("correct pred.", "missed"), pch=20,col=(1:2))

dat<-read.csv("http://matsuka.info/data_folder/tdkDA01.csv", header=T)
dat.lda<-lda(class~.,dat)
dat.lda.pred<-predict(dat.lda)
dat.glm<-glm(class~.,family=binomial,data=dat)
dat.glm.pred<-predict(dat.glm,type="response")
head(cbind(dat.glm.pred,dat.lda.pred$posterior))
cor(dat.glm.pred,dat.lda.pred$posterior[,2])