Disclaimer

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

データ解析基礎論B review02

n.subj = 5
set.seed(10)
groupA1 = rnorm(n.subj, mean = 100, sd =15)
groupA2 = rnorm(n.subj, mean = 100, sd =15)
groupID = c(rep("A1",n.subj),rep("A2",n.subj))
dat = data.frame(score = c(groupA1,groupA2),ID = groupID)
grand.mean = mean(dat$score)
mean.A1 = mean(groupA1)
mean.A2 = mean(groupA2)
plot(dat$ID, dat$score)
points(as.numeric(dat$ID), dat$score, col=c(rep('red',n.subj),rep("blue",n.subj)),pch=20)

idx = 1:(n.subj+n.subj)
plot(idx,dat$score,col=c(rep('red',n.subj),rep("blue",n.subj)),
     pch=20,ylim=c(65,145),ylab ="score",xlab="subjID")
abline(h=grand.mean,lwd=2)
abline(h=mean.A1,lty=2,col='red',lwd=2)
abline(h=mean.A2,lty=2,col='blue',lwd=2)
legend("topright",c("overall mean", "G1 mean","G2 mean"),
       lty=c(1,2,2),col=c('black','red','blue'),lwd=2)

ss.total = sum((dat$score-grand.mean)^2)
ss.group1 = n.subj*(mean.A1-grand.mean)^2
ss.group2 = n.subj*(mean.A2-grand.mean)^2
ss.group = ss.group1+ss.group2
ss.error1 = sum((groupA1-mean.A1)^2)
ss.error2 = sum((groupA2-mean.A2)^2)
ss.error = ss.error1+ss.error2
df.group = 2 - 1
df.error = (n.subj-1)+(n.subj-1)
ms.group = ss.group/df.group
ms.error = ss.error/df.error
f.value = ms.group/ms.error
p.value = 1 - pf(f.value, df.group, df.error)
summary(aov(score~ID, data=dat))
n.rep = 1e5
n.subj = 5
f.hist = rep(0,n.rep)
for (i.rep in 1:n.rep){
  groupA1 = rnorm(n.subj, mean = 100, sd =15)
  groupA2 = rnorm(n.subj, mean = 100, sd =15)
  groupID = c(rep("A1",n.subj),rep("A2",n.subj))
  dat = data.frame(score = c(groupA1,groupA2),ID = groupID)
  grand.mean = mean(dat$score)
  mean.A1 = mean(groupA1)
  mean.A2 = mean(groupA2)
  ss.total = sum((dat$score-grand.mean)^2)
  ss.group = n.subj*(mean.A1-grand.mean)^2+n.subj*(mean.A2-grand.mean)^2
  ss.error = sum((groupA1-mean.A1)^2)+sum((groupA2-mean.A2)^2)
  df.group = 2 - 1
  df.error = (n.subj-1)+(n.subj-1)
  ms.group = ss.group/df.group
  ms.error = ss.error/df.error
  f.hist[i.rep] = ms.group/ms.error
}

hist(f.hist,1000,xlim=c(0,15),probability = T)
f.temp = seq(0,15,0.05)
f.true = df(f.temp,df.group,df.error)
lines(f.temp,f.true,col='red',lty=2,lwd=2)
thres.f = qf(0.95,df.group,df.error)
abline(v = thres.f, col='green')
prop.FP = sum(f.hist > thres.f)/n.rep

データ解析基礎論B Reivew01

# central limit theorem
n.sample = 1e6
n.data = 10
dat = matrix(runif(n.sample*n.data), nrow = n.data)
means = colMeans(dat)
sds = apply(dat,2,sd)

hist(means,100,probability = T)
x.temp  = seq(0,1, 0.01)
true.sd = sqrt(1/12)
true.dens = dnorm(x.temp,mean=0.5, sd = true.sd/sqrt(n.data))
lines(x.temp, true.dens,col='red',lwd=2)
ci.low = qnorm(0.025,mean=0.5, sd = true.sd/sqrt(n.data))
ci.high = qnorm(0.975,mean=0.5, sd = true.sd/sqrt(n.data))
abline(v=ci.low,col='blue',lwd=2)
abline(v=ci.high,col='blue',lwd=2)

# standardizing data
z = (means-0.5)/(true.sd/sqrt(n.data))
hist(z,100,probability = T)
x.tempZ = seq(-4,4,0.01)
z.dens = dnorm(x.tempZ,mean=0, sd = 1)
lines(x.tempZ, z.dens,col='red',lwd=2)
ci.lowZ = qnorm(0.025,mean=0, sd = 1)
ci.highZ = qnorm(0.975,mean=0, sd = 1)
abline(v=ci.lowZ,col='blue',lwd=2)
abline(v=ci.highZ,col='blue',lwd=2)

# t distribution
t.temp = seq(-4,4,0.01)
t.dens = dt(t.temp, df=(n.data-1))
lines(t.temp, t.dens,col='green',lwd=2)
ci.lowT = qt(0.025,df=(n.data-1))
ci.highT = qt(0.975,df=(n.data-1))
abline(v=ci.lowT,col='cyan',lwd=2)
abline(v=ci.highT,col='cyan',lwd=2)


### cell automata
nCell=201
nGen=100
res=matrix(0,nrow=nGen,ncol=nCell)
res[1,100]=1
for (i_gen in 2:nGen) {
  for (i_cell in 2:(nCell-1)) {
    res[i_gen,i_cell]=sum(res[i_gen-1,(i_cell-1):(i_cell+1)])%%2;
  } 
 # 左端
  res[i_gen,1]=(res[i_gen-1,nCell]+res[i_gen-1,1]+res[i_gen-1,2])%%2 
 # 右端
  res[i_gen,nCell]=(res[i_gen-1,(nCell-1)]+res[i_gen-1,nCell]+res[i_gen-1,1])%%2;
} 
image(res)


# general version
dec2bin<-function(num, digits=8) {
  bin=c()
  if (num==0){
    bin=0
  } else {
    while(num!=0){
      rem=num%%2
     num=num%/%2
      bin=c(rem,bin)
    }
  }
  if (length(bin) < digits){
    res=matrix(0,nrow=1,ncol=digits)
    res[(digits-length(bin)+1):digits]=bin
  } else {res=bin}
  return(res)
}

transFUN<-function(st,ruleID){
  output=dec2bin(ruleID,8);
  a=matrix(c(1,1,1,1,1,0,1,0,1,1,0,0,0,1,1,0,1,0,0,0,1,0,0,0),nrow=8,byrow=T)
  newSt=output[which(apply(a,1,function(x) {all(x==st)}))]
  return(newSt)
}
ECA<-function(nCell, nGen,ruleID){
  res=matrix(0,nrow=nGen,ncol=nCell)
  res[1,ceiling(nCell/2)]=1;
  for (i_gen in 2:nGen) {
    for (i_cell in 2:(nCell-1)) {
      res[i_gen,i_cell]=transFUN(res[i_gen-1,(i_cell-1):(i_cell+1)],ruleID)
     }
   res[i_gen,1]=transFUN(c(res[i_gen-1,nCell],res[i_gen-1,1],res[i_gen-1,2]),ruleID)
   res[i_gen,nCell]=transFUN(c(res[i_gen-1,(nCell-1)],res[i_gen-1,nCell],res[i_gen-1,1]),ruleID)
  }
  return(res)
}
res<-ECA(200,100,99)
image(res) 

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

### first bayesian analysis
theta = seq(0,1,0.1)
prior  = c(0:5,4:0)
prior = prior/sum(prior)
likelihood = theta
posterior = likelihood*prior
posterior = posterior/sum(posterior)
par(mfrow = c(3, 1))
barplot(prior)
barplot(likelihood)
barplot(posterior)

### plotting beta dist.
par(mfrow = c(3, 1))
theta = seq(0, 1, 0.001)
prior = dbeta(theta, 0.1, 1)
plot(prior,type ='l',lwd=3)
prior = dbeta(theta, 2, 2)
plot(prior,type ='l',lwd=3)
prior = dbeta(theta, 4, 1)
plot(prior,type ='l',lwd=3)

### jisshu 1
posterior_beta <- function(a, b, z, N){
  par(mfrow = c(3, 1))
  theta = seq(0, 1, 0.001)
  prior = dbeta(theta, a, b)
  likelihood = theta^z * (1-theta)^(N - z)
  posterior = theta^(z+a-1) * (1-theta)^(N-z+b-1) / beta(z+a, N-z+b)
  # normalizing posterior
  posterior = posterior/sum(posterior)

  # plotting prior
  plot(theta,prior,type='l',lwd=2)
  polygon(theta,prior,col='orange')

  # plotting likelihood
  plot(theta,likelihood,type='l',lwd=2)
  polygon(theta,likelihood,col='orange')

  # plotting posterior
  plot(theta,posterior,type='l',lwd=2)
  polygon(theta,posterior,col='orange')
  par(mfrow=c(1,1))
}

###   jisshu 2 - island hopping
# initialization
n.iter = 10000
location = 4
location.history = rep(0, n.iter)
location.history[1] = location
location.counter = rep(0, 7)
location.counter[location] = location.counter[location] + 1
# end initialization

# main loop
for (i.iter in 2:n.iter){
  proposed = location + sample(c(-1,1),1)
  if (proposed == 8){ proposed = 7}
  # if (proposed = 0){ proposed = 1}
  p.move = min(1, proposed/location)
  if (runif(1) < p.move){
    location = proposed
  } # else {location = location}
  location.history[i.iter] = location
  location.counter[location] = location.counter[location] + 1
}
par(mfrow = c(2,1))
plot(location.history, 1:n.iter, type = 'o', log = 'y')
barplot(location.counter)


### jisshu 3 metropolis method
metropolis.ex01 <- function(n.iter=1e6, theta.init=0.1, sigma=0.2, plot.yn = T){

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


### exact beta posterior 2 variable
Exact_beta_posterior2V <- function(a1,b1,a2,b2,z1,N1,z2,N2){
  # DBDA ch. 7 exact posterior distribution
  # WARNING: posterior is not nomarlized
  try(library(plot3D))
  try(library(MASS))
  temp.thetas=seq(0,1,0.01)
  M=mesh(temp.thetas,temp.thetas)
  theta1=M$x
  theta2=M$y
  prior=dbeta(theta1,a1,b1)*dbeta(theta2,a2,b2)
  likelihood=theta1^z1 * (1-theta1)^(N1-z1)*theta2^z2 * (1-theta2)^(N2-z2)
  posterior=prior*likelihood

  par(mfrow=c(3,2))
  par(mar=c(4,1,1,1))
  persp(prior,theta =30,phi=20,lwd=0.1,col="orange")
  contour(prior,labels="",lwd=1.5,col="navy")

  persp(likelihood,theta =30,phi=20,lwd=0.1,col="orange")
  contour(likelihood,labels="",lwd=1.5,col="navy")

  persp(posterior,theta =30,phi=20,lwd=0.1,col="orange")
  contour(posterior,labels="",lwd=1.5,col="navy")

  par(mfrow=c(1,1))
}

### metropolis 2 var
metropolis.ex02 <- function(n.iter=1e6, theta.init=c(0.5,0.5), sigma=0.2){
  # example from DBDA 2nd Ed. (Kruschke) ch. 07
  # metropolis algo. with 2 parameters

  # "posterior of theta" function
  posterior.theta2<-function(theta, N1, z1, a1, b1, N2, z2, a2, b2) {
    posterior = theta[1]^z1 * (1-theta[1])^(N1-z1) *
      theta[2]^z2*(1-theta[2])^(N2-z2) *
      theta[1]^(a1-1)*(1 - theta[1])^(b1-1) *
      theta[2]^(a2-1)*(1 - theta[2])^(b2-1) /
      (beta(a1,b1)*beta(a2,b2))
  }

  # initialization
  try(library(MASS))
  theta.history=matrix(0,nrow=n.iter,ncol=2)
  theta.current=theta.init
  theta.history[1,]=theta.current
  # values given in text
  mu=0
  N1=8
  z1=6
  a1=2
  b1=2
  N2=7
  z2=2
  a2=2
  b2=2

  # main loop
  for (i_iter in 2:n.iter) {
    theta.proposed=theta.current+rnorm(2,mean=mu,sd=sigma)
    if (any(theta.proposed<0) | any(theta.proposed>1)) {
      p.move=0
    } else {
      p.current=posterior.theta2(theta.current, N1, z1, a1, b1, N2, z2, a2, b2)
      p.proposed=posterior.theta2(theta.proposed, N1, z1, a1, b1, N2, z2, a2, b2)
      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
  par(mfrow=c(2,2))
  hist(theta.history[,1],nclass=100, probability = T, col="orange")
  den=density(theta.history[,1])
  lines(den)
  den3d <- kde2d(theta.history[,1],theta.history[,2], n = 100)
  par(mar=c(1,1,1,1))
  persp(den3d, box=FALSE, col = "orange", lwd=0.1, theta=10, phi=20)
  plot(theta.history[,1],theta.history[,2],type='o', ylim=c(0,1), xlim=c(0,1), col="orange")
  hist(theta.history[,2], nclass=100, probability = T, col="orange")
  den=density(theta.history[,2])
  lines(den)

  return(theta.history)
}

res=metropolis.ex02(10000,c(0.1,0.1),0.2)

### jisshu 4 Gibbs sampling
gibbs.ex01 <- function(n.iter=1e6, theta.init=c(0.5,0.5), sigma = 0.2) {

  # posterior of theta function
  posterior.thetaG<-function(theta, N, z, a, b, idx) {
    p = dbeta(theta[idx],z[idx]+a[idx],N[idx]-z[idx]+b[idx])
  }

  # initialization
  try(library(MASS))
  n.theta = 2
  theta.history = matrix(0, nrow = n.iter, ncol = n.theta)
  theta.current = theta.init
  theta.proposed = rep(0, n.theta)
  theta.history[1,] = theta.current
  # values given in text
  mu = 0
  N = c(8,7)
  z = c(6,2)
  a = c(2,2)
  b = c(2,2)

  # main loop
  for (i_iter in 2:n.iter) {
    for (i_theta in 1:n.theta) {
      theta.proposed[i_theta] = theta.current[i_theta] + rnorm(1, mean=mu, sd=sigma)
      if (theta.proposed[i_theta]<0 | theta.proposed[i_theta]>1) {
        p.move=0
      } else {
        p.current = posterior.thetaG(theta.current, N, z, a, b, i_theta)
        p.proposed = posterior.thetaG(theta.proposed,N, z, a, b, i_theta)
        p.move = min(p.proposed/p.current, 1)
      }
      if (runif(1) < p.move) {
        theta.current[i_theta] = theta.proposed[i_theta]
      }
      theta.history[i_iter, i_theta] = theta.current[i_theta]
    }
  }
  # plotting results
  par(mfrow=c(2,2))
  hist(theta.history[,1], nclass = 100, probability = T, col="orange")
  den=density(theta.history[,1])
  lines(den)
  den3d <- kde2d(theta.history[,1],theta.history[,2], n = 100)
  par(mar=c(1,1,1,1))
  persp(den3d, box=FALSE, col = "orange", lwd=0.1, theta=10, phi=20)
  plot(theta.history[,1],theta.history[,2],type='o', ylim=c(0,1), xlim=c(0,1), col="orange")
  hist(theta.history[,2], nclass=100, probability = T, col="orange")
  den=density(theta.history[,2])
  lines(den)

  return(theta.history)
}

res=gibbs.ex01(10000,c(0.1,0.1),0.2)
Posted in UT

広域システム特別講義II

# ryukou
for (t in 1:max_time) {
  x[(t+1)] = x[t] - alpha*x[t]*y[t]*dt
  y[(t+1)] = y[t] + (alpha*x[t]*y[t] - beta*y[t])*dt
  z[(t+1)] = z[t] + beta*y[t]*dt
}

# renai
initX=c(rep(-40,5),rep(-20,5),rep(0,5),rep(20,5),rep(40,5))
initY=c(rep(c(-40,-20,0,20,40),5))

# cake game
n.cake = 10
pay.mat = matrix(0,11,11)
for (i.cake in 1:n.cake){
  pay.mat[(i.cake+1),1:(n.cake-i.cake+1)] =i.cake
}

p.cake = runif(n.cake+1)
p.cake = p.cake/sum(p.cake)
max.time = 50
dt = 0.01
t = seq(0,max.time,dt)
n.iter = length(t)
p.hist = matrix(0,nrow = n.iter, ncol = (n.cake+1))
p.hist[1,] = p.cake
for (i.time in 2:n.iter){
  W = colSums(p.cake*t(pay.mat))
  W.ave = sum(outer(p.cake,p.cake)*pay.mat)
  p.cake = p.cake + p.cake*(W - W.ave)/W.ave * dt
  p.hist[i.time,] = p.cake
}


plot(p.hist[,1],ylim=c(0,1),type='l',lwd=2,ylab = 'Proportion',xlab="time")  
for (i.strat in 2:(n.cake+1)){
  lines(p.hist[,i.strat],col=i.strat,lwd=2)
}
legend("topleft",paste("request = ",0:10),col=1:(n.cake+1),lwd =2,cex=0.75)

Posted in UT

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

実習課題の解答例

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)

Qlearn.ex<-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.ex(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) 
# MC value evaluation
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=rep(0.25,4);V = rep(0,14)max.iter = 10000;
state.count=rep(0,15)
for (i.iter in 1:max.iter){
  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]  
  }
  uniq.seq = unique(state.seq)
  for (i.uniq in 1:(length(uniq.seq)-1)){
    first.visit = which(state.seq == uniq.seq[i.uniq])[1]
    V[uniq.seq[i.uniq]] = V[uniq.seq[i.uniq]] + r*(length(state.seq)-first.visit-1)
  }
  state.count[uniq.seq] = state.count[uniq.seq] + 1
}
V = matrix(c(0,V/state.count[1:14],0),nrow=4)



# black jack problem 01
# function to calc. value of cards
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)}
}}}
 
# 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) {
    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)
    }
    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)
}
 
# function 2 plot results
plot.BJ_MC<-function(V){
 par(mfrow=c(1,2))
 image(V[,,1],main="with usable Ace",xaxt='n',yaxt='n',
   xlab="Dealer showing",ylab="Player sum")
 axis(1,at=seq(0,1,length.out=10),label=c("A",paste(2:10)))
 axis(2,at=seq(0,1,length.out=10),label=12:21)
 image(V[,,2],main="without usable Ace",xaxt='n',yaxt='n',
   xlab="Dealer showing",ylab="Player sum")
 axis(1,at=seq(0,1,length.out=10),label=c("A",paste(2:10)))
 axis(2,at=seq(0,1,length.out=10),label=12:21)
}
 
V=BJ_MC(17)
plot.BJ_MC(V)

# black jack problem 2
# function to calc. value of cards
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)}
}}}
 
# main function
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) {
    # initial draw
    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)
    }
    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)
    }
    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))
}
 
 
# function 2 plot results
plot.BJ_MC2<-function(V){
 par(mfrow=c(2,2))
 image(V[[2]][,,1],main="Utility with usable Ace",xaxt='n',yaxt='n',
  xlab="Dealer showing",ylab="Player sum")
 axis(1,at=seq(0,1,length.out=10),label=c("A",paste(2:10)))
 axis(2,at=seq(0,1,length.out=10),label=12:21)
 image(V[[2]][,,2],main="Utility without usable Ace",xaxt='n',yaxt='n',
  xlab="Dealer showing",ylab="Player sum")
 axis(1,at=seq(0,1,length.out=10),label=c("A",paste(2:10)))
 axis(2,at=seq(0,1,length.out=10),label=12:21)
 image(V[[1]][,,1],main="Policy with usable Ace",xaxt='n',yaxt='n',
  xlab="Dealer showing",ylab="Player sum")
 axis(1,at=seq(0,1,length.out=10),label=c("A",paste(2:10)))
 axis(2,at=seq(0,1,length.out=10),label=12:21)
 image(V[[1]][,,2],main="Policy without usable Ace",xaxt='n',yaxt='n',
  xlab="Dealer showing",ylab="Player sum")
 axis(1,at=seq(0,1,length.out=10),label=c("A",paste(2:10)))
 axis(2,at=seq(0,1,length.out=10),label=12:21)
}
 
V=BJ_MC(5e6)
plot.BJ_MC2(V)


# Sarsa learning
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) {
      max.id = which(Qs[state[1],state[2],]==max(Qs[state[1],state[2],]))
      if (length(max.id)>1){
        move = sample(max.id,1)
      } else {move=max.id}
    } 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) {
        max.id = which(Qs[state[1],state[2],]==max(Qs[state[1],state[2],]))
        if (length(max.id)>1){
          move = sample(max.id,1)
        } else {move=max.id}
      } 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)
}  

Qs=sarsa.ex6.5(10000,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)
Posted in UT

データ解析基礎論B W10

install.packages("nnet")
library(nnet)

# regression
dat<-read.csv("http://www.matsuka.info/data_folder/tdkReg01.csv")
for (i.iter in 1:10){
set.seed(5)
print(i.iter)
x = dat[,1:3]
y = dat[,4]
dat.nnet = nnet(x,y, size = 150, linout= TRUE,maxit = 1000)
nnet.pred <-predict(dat.nnet,dat)
print(c(cor(nnet.pred,dat$sales),cor(nnet.pred,dat$sales)^2))
plot(dat.nnet$fitted.values, dat$sales,pch=20,col='black')
points(predict(dat.lm), dat$sales,pch=20,col='red')
abline(h=max(dat$sales),col='green')
abline(v=max(dat$sales),col='green')
abline(h=min(dat$sales),col='green')
abline(v=min(dat$sales),col='green')
abline(h=mean(dat$sales),col='green')
abline(v=mean(dat$sales),col='green')


n.data = nrow(dat);n.sample = n.data*0.6; n.rep = 1000 
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 = 10, linout=T,decay= 0.1, 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)))

# logistic regression
dat<-read.csv("http://www.matsuka.info/data_folder/tdkDA01.csv")
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.glm<-glm(class~social*coop.*dilig.*enterp., family="binomial",dat)


dat<-read.csv("http://www.matsuka.info/data_folder/cda7-16.csv")
dat.nnet<-nnet(survival~., dat, size=30, maxit=1000, decay=0.01)
dat.pred<-predict(dat.nnet,dat,type="class")
table(dat.pred,dat$survival)

Ns = summary(dat$survival)
(Ns[1]/Ns[2])^-1

wts = rep(1,nrow(dat))
wts[which(dat$survival=="no")]=45
dat.nnet<-nnet(survival~., weights=wts, dat, size=30, maxit=1000, decay=0.01)
dat.pred<-predict(dat.nnet,dat,type="class")
table(dat.pred,dat$survival)

# discriminant analysis
dat<-read.csv("http://www.matsuka.info/data_folder/tdkDA02.csv")
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)

基礎実習MB02

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
 

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

# 5x5 cellular workd
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,]]))
  }
}

#### 4 x 4 celluler world
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;


# value iteration
V=c(rep(0,100),1);V.hist=c()
p=c(0.4,0.6);
gamma=1;delta=1; tol=1e-20
max.a=rep(0,101)
while (delta>tol) {
  delta=0;
  for (i_state in 1:99) {
    v=V[i_state+1]
    temp=matrix(0,nrow=1,ncol=i_state)
    for (i_action in 1:i_state) {
       temp[i_action]=sum(p*(gamma*c(V[(min(i_state+i_action,100)+1)],
         V[(max(i_state-i_action,0)+1)])))
     }
    V[i_state+1]=max(temp)
    max.a[i_state+1]=which.max(round(temp,8))
    delta=max(delta,abs(v-V[i_state+1]))
  }
  V.hist=rbind(V.hist,V)
}
# plotting results
par(mfrow=c(1,2))
plot(V.hist[1,],type='l',lwd=2,xlab="Capital",ylab="Value Estimates",col='red')
lines(V.hist[2,],lwd=2,col='blue')
lines(V.hist[3,],lwd=2,col='green')
lines(V.hist[32,],lwd=2,col='black')
legend("topleft",c("sweep 1","sweep 2","sweep 3", "sweep 32"),
 col=c("red","blue","green","black"),lwd=2)
barplot(max.a,xlab="Capital",ylab="Final Policy",col="white")
Posted in UT

データ解析基礎論B W08

dat<-data.frame(x1=c(50,69,93,76,88,43,56,38,21,25),
                x2=c(15.5,18.4,26.4,22.9,18.6,16.9,21.6,12.2,16.0,10.5), 
                cl=c(rep("h",5),rep("d",5)))
library(MASS)
dat.lda<-lda(cl~.,data=dat)
 2, col='skyblue')

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


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

plot(dat.lda, dimen=2, col=as.numeric(dat$class), cex=3)
plot(dat.lda, dimen=5, col=as.numeric(lda.pred$class),cex=2)


dat.km<-kmeans(dat[,1:6],5)

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)
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=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.glm<-glm(class~.,family=binomial,data=dat)

#MDS
dat<-data.frame(p1=c(4,1,5,1,5),p2=c(1,5,4,3,1))
rownames(dat)<-c('a','b','c','d','e')
dat.mds<-cmdscale(dist(dat),2)
plot(dat.mds[,1],dat.mds[,2], type='n')
text(dat.mds[,1],dat.mds[,2],labels=row.names(dat))
dat<-read.csv("http://matsuka.info/data_folder/tdkMDS02.csv", row.name=1)
dat.mds<-cmdscale(dist(t(dat)),2,eig=T)

plot(dat.mds$points[,1],dat.mds$points[,2], type='n')
text(dat.mds$points[,1],dat.mds$points[,2],labels=colnames(dat), cex=2)

dat<-read.csv("http://www.matsuka.info/data_folder/tdkCFA.csv")

データ解析基礎論B WAB08

データ解析基礎論B WAB08 IRTなど
提出期限 2017.12.12 講義開始前まで

WAB08.1
このデータは、ある試験のデータ(計30問)です。まず、ltmパッケージのdescriptなどを使用して、各問題の記述統計値やクロンバックのアルファなどを求めてください。
次に、1パラメターのIRTモデル(rasch)及び2パラメターのIRTモデル(ltm)を用いて分析してください。そして、モデルを比較しより適切なモデルはどちらか説明してください。最後に、最適なモデルの問題の難しさ(Dffclt)の推定値と、各問題の統計量(平均値)の相関を求めてください。

WAB08.2
このデータを用いて、WAB08.1と同じ分析を行ってください。