otherTests/multivariateMixture.R

In lucas-castillo/SampleR: Tools To Compare Human Performance To Sampling Algorithms

require(distr6)
require(ggplot2)
require(SampleR)


set.seed(1)

M <- matrix(0, nrow = 15, ncol = 2)
for (i in 1:15){
  M[i,] <- runif(2) * 18 - 9
}
# m = pracma::rand(15,2) * 18-9
listDistr <- list()
for (i in 1:15){
  listDistr[[i]] <- MultivariateNormal$new(mean = c(m[i,1], m[i,2]))
}

weights <- rep(1/length(listDistr), length(listDistr))

seq_weights <- vector()

pdfs <- list()
for (i in 1:length(listDistr)){
  dist <- listDistr[[i]]
  pdfs[[i]] <- dist$pdf
}

for (i in 1:length(weights)){
  if (i == 1){
    seq_weights[i] <- weights[i]
  } else {
    seq_weights[i] <- weights[i] + seq_weights[i-1]
  }
}

myPDF <- function(x){
  densities <- vector()
  for (i in 1:length(pdfs)){
    densities[i] <- pdfs[[i]](data = matrix(x,nrow=1))
  }
  return(sum(densities * weights))
}
myPDF2 <- function(x, log = FALSE){
  densities <- vector()
  for (i in 1:length(weights)){
    densities[i] <- mvtnorm::dmvnorm(x, mean = c(m[i,1], m[i,2]), log = log)
  }
  return(sum(densities * weights))
}

system.time(myPDF(c(5,-1)))
system.time(myPDF2(c(5,-1)))

dims <- 2
myRand <- function(x){
  randNums <- matrix(ncol = dims, nrow = x)
  for (j in 1:x){
    randomN <- stats::runif(1)
    finish = F
    i = 0
    while (!finish){
      i = i + 1
      if (randomN <= seq_weights[i]){
        dist <- listDistr[[i]]
        randNums[j,] <- unlist(unname(dist$rand(1)))
        finish = T
      }
    }
  }
  return(randNums)
}

mapDensity1 <- function(pdf, start, size, cellsPerRow = 50){
  # start is a vector <- c(x, y)
  #  is a number n so that the map ranges from x, y to x + n, y + n
  xRange <- seq(from = start[1], to = start[1] + size, length.out = cellsPerRow)
  xxRange <- rep(xRange, cellsPerRow)

  yRange <- seq(from = start[2], to = start[2] + size, length.out = cellsPerRow)
  for (i in 1:cellsPerRow){
    if (i == 1){
      yyRange <- rep(yRange[i], cellsPerRow)
    } else {
      yyRange <- c(yyRange, rep(yRange[i], cellsPerRow))
    }
  }

  density <- vector()


  for (i in 1:length(yyRange)){
    density[i] <- pdf(c(xxRange[i],yyRange[i]))
    }

  df <- data.frame(x = xxRange, y = yyRange, density = density)
  return(df)
}


mapDensity2 <- function(start, size, cellsPerRow = 50){
  # start is a vector <- c(x, y)
  #  is a number n so that the map ranges from x, y to x + n, y + n
  xRange <- seq(from = start[1], to = start[1] + size, length.out = cellsPerRow)
  xxRange <- rep(xRange, cellsPerRow)

  yRange <- seq(from = start[2], to = start[2] + size, length.out = cellsPerRow)
  for (i in 1:cellsPerRow){
    if (i == 1){
      yyRange <- rep(yRange[i], cellsPerRow)
    } else {
      yyRange <- c(yyRange, rep(yRange[i], cellsPerRow))
    }
  }

  densities <- matrix(ncol = length(listDistr), nrow = cellsPerRow ** 2)

  for (i in 1:length(pdfs)){
    densities[,i] <- pdfs[[i]](xxRange, yyRange) * weights[i]
  }
  density <- rowSums(densities)
  df <- data.frame(x = xxRange, y = yyRange, density = density)
  return(df)
#
#   for (i in 1:length(yyRange)){
#     density[i] <- pdf(c(xxRange[i],yyRange[i]))
#   }
#
#   df <- data.frame(x = xxRange, y = yyRange, density = density)
#   return(df)
}


# hills_df <- mapDensity1(myPDF2, c(-10,-10), 20, 150)
# write.csv(hills_df, "hillsDF.csv")

hills_df <- read.csv("otherTests/hillsDF.csv")
hill_map <- ggplot(hills_df) + geom_raster(mapping = aes(x = x, y = y, fill = density)) + scale_fill_viridis_c()

DS <- myRand(1024)
DS_df <- data.frame(x = DS[,1], y = DS[,2])
DS_path <- hill_map + geom_path(DS_df, mapping = aes(x,y))

MCMC <- mcmc(myPDF2, start = c(0,0), sigma_prop = diag(2) / 8)
MCMC_df <- data.frame(x = MCMC[[1]][,1], y = MCMC[[1]][,2])
MCMC_path <- hill_map + geom_path(MCMC_df, mapping = aes(x,y)) + theme_minimal()
# ggplot2::ggsave(filename = "docs/figures/MCMC_path.jpg", plot = MCMC_path)

MC3 <- mc3(myPDF2, start = c(0,0), sigma_prop = diag(2) / 8, swap_all = FALSE)
MC3_df <- data.frame(x = MC3[[1]][,1,1], y = MC3[[1]][,2,1])
MC3_path <- hill_map + geom_path(MC3_df, mapping = aes(x, y)) + theme_minimal()
# ggplot2::ggsave(filename = "docs/figures/MC3_path.jpg", plot = MCMC_path)

HMC <- hmc(myPDF2, epsilon = 4, L = 500, start = c(-5,0), iterations = 10)
HMC_df <- data.frame(x = HMC[[1]][,1], y = HMC[[1]][,2])
HMC_path <- hill_map + geom_path(HMC_df, mapping = aes(x = x, y = y))


set.seed(1)
NUTS <- hmc_nuts(myPDF2, start = c(0,0), epsilon = .5)
NUTS_df <- data.frame(x = NUTS[,1], y = NUTS[,2])
NUTS_path <- hill_map + geom_path(NUTS_df, mapping = aes(x, y))