-
Notifications
You must be signed in to change notification settings - Fork 0
/
plotSurface.R
157 lines (135 loc) · 4.81 KB
/
plotSurface.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
if(!require(lattice)) { install.packages("lattice"); library(lattice) }
if(!require(mvtnorm)) { install.packages("mvtnorm"); library(mvtnorm) }
if(!require(numDeriv)) { install.packages("numDeriv"); library(numDeriv) }
if(!require(ellipse)) { install.packages("ellipse"); library(ellipse) }
if(!require(grDevices)) { install.packages("grDevices"); library(grDevices) }
if(!require(wesanderson)) { install.packages("wesanderson"); library(wesanderson) }
plotTrajectory <- function (start, G, scale= 2) {
current_gen = start
es = eigen(cov2cor(G))$values
v1 = sqrt(es[1])/1.2 * eigen(cov2cor(G))$vectors[,1]
v2 = sqrt(es[2])/1.2 * eigen(cov2cor(G))$vectors[,2]
polygon(ellipse(G[1,2], centre = current_gen, level = 0.3), col = wes_palette("Royal1")[1])
segments(current_gen[1] - v1[1], current_gen[2] - v1[2],
current_gen[1] + v1[1], current_gen[2] + v1[2], lwd = 2)
segments(current_gen[1] - v2[1], current_gen[2] - v2[2],
current_gen[1] + v2[1], current_gen[2] + v2[2], lwd = 2)
points(current_gen[1], current_gen[2], pch = 19)
net_beta = c(0, 0)
for(i in 1:gen){
beta = grad(W_bar, t(current_gen))
net_beta = net_beta + beta
next_gen = current_gen + (G/scale)%*%beta
# arrows(current_gen[1], current_gen[2],
# next_gen[1], next_gen[2], pch = 18, length = 0.14, lwd = 2.5, col = "black")
# arrows(current_gen[1], current_gen[2],
# current_gen[1] + beta[1]/5, current_gen[2] + beta[2]/5,
# pch = 18, length = 0.14, lwd = 3, col = wes_palette("Rushmore")[3])
current_gen = next_gen
}
net_delta = (G/scale) %*% net_beta
arrows(start[1], start[2],
start[1] + net_delta[1], start[2] + net_delta[2],
pch = 18, length = 0.14, lwd = 2.5, col = 'black')
arrows(start[1], start[2],
start[1] + net_beta[1]/5, start[2] + net_beta[2]/5,
pch = 18, length = 0.14, lwd = 2.5, col = wes_palette("FantasticFox")[5])
}
w_cov = matrix(c(1.0, 0.7,
0.7, 1.0), ncol = 2)
G = matrix(c(1.0, 0.8,
0.8, 1.0), ncol = 2)
gen = 15
step = 0.1
mypalette = colorRampPalette(c(wes_palette(10, name = "Zissou", type = "continuous"), "darkred"))
#################################
# Multiple Peaks
#################################
W_bar = function(x) {
log(
dmvnorm(x, mean = c(4, 5), sigma = w_cov) +
1.1*dmvnorm(x, mean = c(1, 5), sigma = w_cov) +
dmvnorm(x, mean = c(2, 2), sigma = w_cov) +
dmvnorm(x, mean = c(7, 5), sigma = w_cov) +
dmvnorm(x, mean = c(5, 2), sigma = w_cov))
}
x <- seq(-2, 9.6, step) ## valores para mu
y <- seq(-2, 9.6, step)
X <- as.matrix(expand.grid(x, y))
Z <- vector()
for(i in 1:nrow(X)){
Z[i] <- W_bar(c(X[i,1], X[i,2]))
}
Z = exp(Z - log(sum(exp(Z))))
b <- matrix(Z, length(x))
mypalette = colorRampPalette(c("white", wes_palette(10, name = "Zissou", type = "continuous"), "darkred"))
png("multipeaklandscape.png", width = 1080, height = 900)
filled.contour(x, y, z = b, color.palette = mypalette,
plot.axes = {
axis(1);
axis(2);
plotTrajectory(c( 1.8,4.03), G)
plotTrajectory(c(0.25,0.25), G)
plotTrajectory(c(5 ,4 ), G)
})
dev.off(dev.cur())
#################################
# Single Peak
#################################
W_bar = function(x) {
log(dmvnorm(x, mean = c(3, 3), sigma = w_cov))
}
x <- seq(0, 6.5, step) ## valores para mu
y <- seq(-0.5, 6, step)
X <- as.matrix(expand.grid(x, y))
Z <- vector()
for(i in 1:nrow(X)){
Z[i] <- W_bar(c(X[i,1], X[i,2]))
}
Z = exp(Z - log(sum(exp(Z))))
b <- matrix(Z, length(x))
png("singlepeaklandscape.png", width = 1080, height = 900)
filled.contour(x, y, z = b, color.palette = mypalette,
plot.axes = {
axis(1);
axis(2);
plotTrajectory(c(2,4), G)
plotTrajectory(c(1,1), G)
plotTrajectory(c(5.5,3), G)
}
)
dev.off(dev.cur())
#################################
# Gradient
#################################
G_1 = matrix(c(1.0, -0.8,
-0.8, 1.0), ncol = 2)
G_2 = matrix(c(1.0, 0,
0, 1.0), ncol = 2)
G_3 = matrix(c(1.0, 0.8,
0.8, 1.0), ncol = 2)
gen = 30
step = 0.1
W_bar = function(x) {
log(x[2])
}
x <- seq(1, 9, step) ## valores para mu
y <- seq(1, 9, step)
X <- as.matrix(expand.grid(x, y))
Z <- vector()
for(i in 1:nrow(X)){
Z[i] <- W_bar(c(X[i,1], X[i,2]))
}
Z = exp(Z - log(sum(exp(Z))))
b <- matrix(Z, length(x))
png("linearlandscape.png", width = 1080, height = 900)
filled.contour(x, y, z = b, color.palette = mypalette,
plot.axes = {
axis(1);
axis(2)
plotTrajectory(c(3,4), G_1, 3)
plotTrajectory(c(5,4), G_2, 3)
plotTrajectory(c(7,4), G_3, 3)
}
)
dev.off(dev.cur())