---

We start with some R warm-up.

# Numerical value
v <- 23.5
v

## [1] 23.5

# Create a vector
apple <- c('red','green',"yellow")
apple

## [1] "red"    "green"  "yellow"

# Create a matrix
M <- matrix( c('a','a','b','c','b','a'), nrow = 2, ncol = 3, byrow = TRUE)
M

##      [,1] [,2] [,3]
## [1,] "a"  "a"  "b"
## [2,] "c"  "b"  "a"

# Removing a variable
rm(M)

Some basic statistics with R: normal and uniform distributions.

# Normal distribution

# Let's generate a vectors of random numbers from a normal distribution
x1 <- rnorm(10000, mean=0, sd=1)
# The breaks argument specifies how many bars are in the histogram
hist(x1, probability=TRUE, breaks = 100)

# pnorm returns the integral from −\inf to q of the pdf of the normal
# distribution where q is a Z-score
# a z-score is the number of standard deviations from the mean a data point is
pnorm(0)

## [1] 0.5

pnorm(2)

## [1] 0.9772499

# qnorm function is simply the inverse of the cdf, which you can also
# think of as the inverse of pnorm
qnorm(0.5)

## [1] 0

qnorm(0.98)

## [1] 2.053749

# Uniform distribution
punif(0.75, min = 0, max = 1)

## [1] 0.75

qunif(0.75, min = 0, max = 1)

## [1] 0.75

x2 <- runif(1000, min = 0, max = 1)
hist(x2, probability=FALSE, breaks = 100)

# Descriptive statistics
mean(x1)

## [1] -0.005080723

sd(x1)

## [1] 1.000038

var(x1)

## [1] 1.000076

min(x1)

## [1] -3.656309

max(x1)

## [1] 3.442262

median(x1)

## [1] -0.01164455

range(x1)

## [1] -3.656309  3.442262

quantile(x1)

##          0%         25%         50%         75%        100%
## -3.65630909 -0.67498496 -0.01164455  0.66671920  3.44226164

summary(x1)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
## -3.656309 -0.674985 -0.011645 -0.005081  0.666719  3.442262

---

The volesti package

volesti is a C++ package (with an R interface) for computing estimations of volume of polytopes given by a set of points or linear inequalities or Minkowski sum of segments (zonotopes). There are two algorithms for volume estimation and algorithms for sampling, rounding and rotating polytopes.

We can download the R package from https://CRAN.R-project.org/package=volesti

# first load the volesti library
#install.packages('volesti')
library(volesti)

## Loading required package: Rcpp

packageVersion("volesti")

## [1] '1.0.3'

You have access to the documentation of volesti functions like volume computation and sampling.

help("volume")
help("sample_points")

Full documentation here: https://cran.r-project.org/web/packages/volesti/volesti.pdf

Let’s try our first volesti command to estimate the volume of a 3-dimensional cube \(\{-1\leq x_i\leq 1, x_i\in\mathbb{R}\ |\ i=1,2,3\}\)

P <- GenCube(3,'H')
print(volume(P))

## [1] 8.023323

What is the exact volume of P? Did we obtain a good estimation?

Sampling

Sampling uniformly in the square.

library(ggplot2)

x1<-runif(1000, min = -1, max = 1)
x2<-runif(1000, min = -1, max = 1)

g<-ggplot(data.frame( x=x1, y=x2 )) + geom_point( aes(x=x, y=y))

g<-g+annotate("path",
   x=cos(seq(0,2*pi,length.out=100)),
   y=sin(seq(0,2*pi,length.out=100)),color="red")+coord_fixed()

plot(g)

Can we estimate the volume of the red ball via sampling? Solution: rejection sampling. The following computation illustrates that this will fail in (not so) high dimensions.

# run in around 1 min
for (d in 2:20) {
  num_of_points <- 100000
  count_inside <- 0

  points1 <- matrix(nrow=d, ncol=num_of_points)
  for (i in 1:d) {
    x <- runif(num_of_points, min = -1, max = 1)
    for (j in 1:num_of_points) {
      points1[i,j] <- x[j]
    }
  }

  for (i in 1:num_of_points) {
    if (norm(points1[,i], type="2") < 1) {
      count_inside <- count_inside + 1
    }
  }
  vol_estimation <- count_inside*2^d/num_of_points
  vol_exact <- pi^(d/2)/gamma(d/2+1)

  cat(d, vol_estimation, vol_exact, abs(vol_estimation- vol_exact)/
vol_exact, "\n")
}

## 2 3.14128 3.141593 9.952073e-05
## 3 4.1864 4.18879 0.0005706194
## 4 4.9336 4.934802 0.0002436168
## 5 5.23456 5.263789 0.005552847
## 6 5.09824 5.167713 0.01344362
## 7 4.6208 4.724766 0.02200447
## 8 4.04224 4.058712 0.004058461
## 9 3.14368 3.298509 0.04693906
## 10 2.42688 2.550164 0.04834357
## 11 1.86368 1.884104 0.0108401
## 12 1.6384 1.335263 0.2270244
## 13 0.73728 0.9106288 0.1903616
## 14 0.49152 0.5992645 0.1797946
## 15 0.32768 0.3814433 0.140947
## 16 0.65536 0.2353306 1.784848
## 17 0 0.1409811 1
## 18 0 0.08214589 1
## 19 0 0.0466216 1
## 20 0 0.02580689 1

Sampling via random walks

volesti supports 3 types of random walks

1. Ball walk

2. Random directions hit-and-run

3. Coordinate directions hit-and-run

There are two important parameters cost per step and mixing time that affects the accuracy and performance of the walks. Below we illustrate this by choosing different walk steps for each walk while sampling on the 100-dimensional cube.

#run in few secs
library(ggplot2)
library(volesti)
for (step in c(1,20,100,150)){
  for (walk in c("CDHR", "RDHR", "BW")){
    P <- GenCube(100, 'H')
    points1 <- sample_points(P, WalkType = walk, walk_step = step, N=1000)
    g<-plot(ggplot(data.frame( x=points1[1,], y=points1[2,] )) +
geom_point( aes(x=x, y=y, color=walk)) + coord_fixed(xlim = c(-1,1),
ylim = c(-1,1)) + ggtitle(sprintf("walk length=%s", step, walk)))
  }
}