In [1]:

    

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


    

    




Loading required package: lattice
Loading required package: mvtnorm
Loading required package: numDeriv
Loading required package: ellipse
Loading required package: wesanderson

In [2]:

    

plotTrajectory <- function (start) {
  current_gen = start
  polygon(ellipse(0.8, 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%*%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
  }
  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])
  #  net_delta = G %*% 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')
}


    

In [4]:

    

#w_cov = matrix(c(1, 0.7, 0.7, 1), ncol = 2)
w_cov = 1.3
G = matrix(c(1, 0.8, 0.8, 1)/2, ncol = 2)

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]

gen = 15

W_bar = function(x) {
  a = 0.5
  b = 5
  # a * x -1 * y + b = 0
  dist = (a * x[1] - x[2] + b)/sqrt(a*a + 1)
  log(dnorm(dist, mean = 0, sd = w_cov))
}
step = 0.05
x <- seq(-1.5, 8.5, step) ## valores para mu
y <- seq(-1.5, 8.5, step)
X <- as.matrix( expand.grid(x, y))
colnames(X) <- c("mu","var")
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(wes_palette(10, name = "Zissou", type = "continuous"), "darkred"))

filled.contour(x, y, z = b, color.palette = mypalette,
               plot.axes = { 
                 axis(1); 
                 axis(2);
                 plotTrajectory(c(2,0))
               }
)