# VolEsti (volume computation and sampling library) # Copyright (c) 2012-2020 Vissarion Fisikopoulos # Copyright (c) 2018-2020 Apostolos Chalkis # Copyright (c) 2020-2020 Marios Papachristou # Contributed and/or modified by Marios Papachristou, as part of Google Summer of Code 2020 program. # Licensed under GNU LGPL.3, see LICENCE file # Example script for using the logconcave sampling methods # Import required libraries library(ggplot2) library(volesti) # Sampling from logconcave density example # Helper function norm_vec <- function(x) sqrt(sum(x^2)) # Negative log-probability oracle f <- function(x) (0.5 * norm_vec(x)^2) # Negative log-probability gradient oracle grad_f <- function(x) (x) # Generate polytope P <- gen_cube(10, 'H') d <- dim(P\$A)[2] # Mode of logconcave density x_min <- 0 * c(1:d) # Smoothness and strong-convexity L <- 1 m <- 1 # Warm start point from truncated Gaussian warm_start <- sample_points(P, n = 1, random_walk = list("nburns" = 5000), distribution = list("density" = "gaussian", "variance" = 1/L, "mode" = x_min)) # Sample points n_samples <- 80000 n_burns <- n_samples / 2 pts <- sample_points(P, n = n_samples, random_walk = list("walk" = "HMC", "step_size" = 0.03, "nburns" = n_burns, "walk_length" = 50, "solver" = "leapfrog", "starting_point" = warm_start[,1]), distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f, "L_" = L, "m" = m)) # pts <- sample_points(P, n = n_samples, random_walk = list("walk" = "HMC", "step_size" = 0.3, "nburns" = n_burns, "walk_length" = 3, "solver" = "leapfrog", "starting_point" = warm_start[,1]), distribution = list("density" = "logconcave", "mode" = x_min, "variance" = 1)) # Plot histogram hist(pts, probability=TRUE, breaks = 100) cat("Sample mean is: ") sample_mean <- mean(pts) cat(sample_mean) cat("\n") cat("Sample variance is: ") sample_variance <- mean((pts - sample_mean)^2) cat(sample_variance)