広域システム特別講義II 10.11の講義

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

Posted in UT

広域システム特別講義II 課題1

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

z= \frac{1}{20}x^2+y^2
を最小化するxとyを慣性を含む勾配法を用いて求めてください。

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

広域システム 09.28

x=rnorm(n=1, mean=100, sd=15)
y=runif(n=3, min=1, max=10)
N = 10000

random.data = rnorm(N, mean=0, sd=1)
hist(random.data, nclass = 50, col = "navy", xlab = "Data",
     probability = T, main = "Histogram of Random Data")
dens = density(random.data)
lines(dens, col = "orange", lwd = 4)

sample(1:10,3)
sample(c(“gu”,“choki”,“pa”),1)
sample(1:10)
sample(0:1, 10, replace=T)

# FOR loop
for (i_loop in 1:5) {
  print(i_loop)
}

for (i_loop in 1:5) {
  print(c(i_loop, 2^i_loop))
}

counter <- 1
while(counter<=10){
  print(counter)
  counter<-counter+1
}

counter <- 1
while(counter^2 <= 10){
  print(c(counter, counter^2))
  counter<-counter+1
}

affil<-"cogsci"
if (affil=="cogsci") {
  print("you are wonderful")
}

affil<-"phil"
if (affil=="cogsci") {
  print("you are wonderful")
} else {
  print("still, you are wonderful")
}

counter=6
repeat{
  print(counter)
  counter = counter + 1
  if(counter>5){break}
}

counter=6
repeat{
  if(counter>5){break}
  print(counter)
  counter+counter+1
}

six.counter=0; N = 1000
for (i_loop in 1:N) {
  die<-sample(1:6,1)
  if (die==6) {six.counter=six.counter+1}
}
six.counter/N

N=1000; six.counter=rep(0,N);
for (i_loop in 1:N) {
  die<-sample(1:6,1)
  if (die==6) {six.counter[i_loop]=1}
}
plot(1:N,cumsum(six.counter)/(1:1000),type='l',ylim=c(0,1),
     lwd=2)
abline(h=1/6,lwd=2,col='red')

N = 1000
die.all <- sample(1:6,N,replace=T)
six.index <- die.all==6
par(mfrow = c(2,1))
par(oma=c(2,2,0,0),mar=c(4,4,1,1),mfrow=c(2,1))
plot(1:N, die.all, pch=20, col = 'red', ylim = c(0,7),
     ylab = "Result", xlab = "trial")
plot(1:N,cumsum(six.index)/(1:1000), type='l', ylim=c(0,1), lwd=2,
     ylab = "P(die = 6)", xlab = "trial")
abline(h=1/6,lwd=2,col='red')

# CLT w/ loop
N=10;nRep=10000;
means<-rep(0,nRep)
for (i_rep in 1:nRep) {
  dat<-runif(N)
  means[i_rep]=mean(dat)
}
hist(means,nclass=50,probability=T)
dens<-density(means)
lines(dens,col='skyblue',lwd=3)
xs=seq(-0,1,0.01)
theo.dens<-dnorm(xs,mean=0.5,sd=sqrt((1/12)/N))
lines(xs,theo.dens,col='orange',lwd=3,lty=2)

# CLT w/o loop
N=10
nRep=10000
dat<-matrix(runif(N*nRep),nrow=N)
means<-colMeans(dat)
hist(means,nclass=50,probability=T)
dens<-density(means); lines(dens,col='skyblue',lwd=3)
xs=seq(-0,1,0.01)
theo.dens<-dnorm(xs,mean=0.5,sd=sqrt((1/12)/N))
lines(xs,theo.dens,col='orange',lwd=3,lty=2)

r=-99;v1=1633;v2=355
while (r!=0){
  r=v1%%v2
  print(paste('v1 =',v1,', v2 = ',v2,',remainder = ',r))
  v1=v2
  v2=r
}

tol = 1e-7;grad = 1e10; lambda = 0.1;x = 10; x.hist = x
grad = 2*x+2
while (abs(grad)>tol){
  x = x - lambda*grad
  x.hist=c(x.hist,x)
  grad = 2*x+2
}
x.temp=seq(-10,10,0.1)
plot(x.temp, x.temp^2+2*x.temp,type='l',lwd=2)
lines(x.hist,x.hist^2+2*x.hist,type='o',pch=1,col='red',cex=1)
Posted in UT

UT 特別講義II week01

# SEC: random numbers
N = 10000
# N = 1000
random.data = rnorm(N, mean=0, sd=1)
hist(random.data, nclass = 50, col = "navy", xlab = "Data",
     probability = T, main = "Histogram of Random Data")
# density of generated data
dens = density(random.data)
lines(dens, col = "orange", lwd = 4)
# theoretical density
x = seq(-4,4,0.1)  
true.norm = dnorm(x, mean = 0, sd = 1)
lines(x,true.norm, col = "green", lty = 3, lwd = 4)
legend("topleft",c("empirical", "theoretical"), lty = c(1,3),
       col = c('orange','green'),lwd=4)

# SEC: law of large numbers
# simplest version
six.counter=0; N = 1000
for (i_loop in 1:N) {
  die<-sample(1:6,1)
  if (die==6) {six.counter=six.counter+1}
}
six.counter/N
# simpler version
N=1000; six.counter=rep(0,N);
for (i_loop in 1:N) {
  die<-sample(1:6,1)
  if (die==6) {six.counter[i_loop]=1}
}
plot(1:N,cumsum(six.counter)/(1:1000),type='l',ylim=c(0,1),lwd=2)
abline(h=1/6,lwd=2,col='red')
# simple version
N = 1000
die.all <- sample(1:6,N,replace=T)
six.index <- die.all==6
par(mfrow = c(2,1))
par(oma=c(2,2,0,0),mar=c(4,4,1,1),mfrow=c(2,1))
plot(1:N, die.all, pch=20, col = 'red', ylim = c(0,7), 
     ylab = "Result", xlab = "trial")
plot(1:N,cumsum(six.index)/(1:1000), type='l', ylim=c(0,1), lwd=2, 
     ylab = "P(die = 6)", xlab = "trial")
abline(h=1/6,lwd=2,col='red')

# CLT
par(mfrow=c(1,1))
N=10
nRep=10000
dat<-matrix(runif(N*nRep),nrow=N)
means<-colMeans(dat)
hist(means,nclass=50,probability=T)
dens<-density(means);lines(dens,col='skyblue',lwd=3)  
xs=seq(-0,1,0.01)
theo.dens<-dnorm(xs,mean=0.5,sd=sqrt((1/12)/N))
lines(xs,theo.dens,col='orange',lwd=3,lty=2)

# GCD
# script version
r=-99;v1=1633;v2=355
while (r!=0){
  r=v1%%v2
  print(paste('v1 =',v1,', v2 = ',v2,',remainder = ',r))
  v1=v2
  v2=r
}
# function version
GCD<-function(v1,v2){  
  real.v1=max(c(v1,v2))  
  real.v2=min(c(v1,v2))
  repeat{
    r=real.v1%%real.v2;real.v1=real.v2;real.v2=r    
    if (r==0){
      print(paste('GCD  is',real.v1)); 
      return(real.v1);
      break}
  }
}
GCD(1633,355)

# digit conversion
# binary to dec
bin2dec=function(bin) {  
  ones=which(rev(bin)==1)-1   
  dec=sum(2^ones)  
  return(dec)
} 
# dec 2 bin - script
num=150;bin=c()
while(num!=0) {
  rem=num%%2; 
  num=num%/%2
  bin=print(c(rem,bin))
}
# dec 2 bin - function
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)
			
Posted in UT