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

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

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


install.packages("plot3D")
library(plot3D)
x = seq(-10,10,0.2)
y = seq(-5,5,0.2)
M = mesh(x,y)
Z = as.vector(1/20*M$x^2)+as.vector(M$y^2)
Z.mesh = matrix(Z,nrow(M$x)) contour(x, y, Z.mesh, drawlabels = F, nlevels=40)  その後linesで更新履歴を追加してみると良いかもしれません。 # gdは勾配法でのx、yの更新履歴 lines(gd, type='o', col = 'green', pch=20) # gdmは慣性ありの勾配法でのx、yの更新履歴 lines(gdm, type='o', col = 'blue', pch=20)  実装例 fun01 = function(x,y){ z = 1/20*x^2 + y^2 return(z) } fun01.grad <- function(x,y){ grad.x = 1/10*x grad.y = 2*y return(list(gx=grad.x, gy = grad.y)) } # script n.loop = 100; lambda = 0.9 x = y = rep(0,n.loop) x[1] = -7; y[1] = 2 for (i.loop in 2:n.loop){ g = fun01.grad(x[i.loop - 1],y[i.loop - 1]) x[i.loop] = x[i.loop - 1] - lambda*g$gx
y[i.loop] = y[i.loop - 1] - lambda*g$gy } # function min.fun01 <- function(init.x, init.y, n.loop, labmda){ x = y = rep(0,n.loop) x[1] = init.x; y[1] = init.y for (i.loop in 2:n.loop){ g = fun01.grad(x[i.loop - 1],y[i.loop - 1]) x[i.loop] = x[i.loop - 1] - lambda*g$gx
y[i.loop] = y[i.loop - 1] - lambda*g$gy } return(list(x = x, y = y)) } res <- min.fun01(-7, 2, 100, 0.9) temp.x = seq(-10,10,0.2);temp.y = seq(-5,5,0.2) M = mesh(temp.x,temp.y) Z = as.vector(1/20*M$x^2)+as.vector(M$y^2) Z.mesh = matrix(Z,nrow(M$x))
contour(temp.x, temp.y, Z.mesh, drawlabels = F, nlevels=40)
lines(res$x, res$y, type='o', col="green")

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

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

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

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

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


wikipediaの記事

1つ目の方策はgreedy法で、リワードの推定値が最も高いスロットマシーンを選択するといったものです。
2つ目の方策は基本的にはgreedy法で、ある確率epsilonでランダムにスロットマシーンを選択するものです。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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