We will try to get similar results as the previous report from 18 April, Probabilistic Simulation Plan with Bernoulli distribution but also others. This previous report suggest that we should draw consistent conclusions.
The classical distributions are: Bernoulli, Binomial, normal, exponential, Beta.
We discard the normal distribution as negative values are irrelevant in our case. All other distributions require to derive formulas. We will process them in the order of appearance above.
We assume that value M is taken with probability p (0 otherwise). The mean is mu=pM and the variance is p(1-p)M2. Hence, for a given coefficient of variation cov and a given mean, we have: p=1/(1+cov2) and M=mu/p.
rbern_cost <- function(n, mean, sigma) {
if (mean == 0 || sigma == 0)
return(rep(mean, n))
cov <- sigma/mean
p <- 1/(cov^2 + 1)
M <- mean/p
return(rbinom(n, 1, p) * M)
}
mean(rbern_cost(10000, 10, 1))
## [1] 10.01
sqrt(var(rbern_cost(10000, 10, 1)))
## [1] 1.005
There are n trials with probability of success p. The mean is npM and the variance is np(1-p)M2. It is preferable to fix the number of trials parameter and to control the mean with a scaling parameters M as this number is an integer. We get p=1/(1+cov2 x n).
rbino_cost <- function(n, mean, sigma) {
if (mean == 0 || sigma == 0)
return(rep(mean, n))
cov <- sigma/mean
size <- 100
p <- 1/(1 + cov^2 * size)
M <- mean/(size * p)
return(rbinom(n, size, p) * M)
}
mean(rbino_cost(10000, 10, 1))
## [1] 10.02
sqrt(var(rbino_cost(10000, 10, 1)))
## [1] 1.007
size <- 100
1/(1 + c(0.01, 0.1, 1, 10)^2 * size)
## [1] 9.901e-01 5.000e-01 9.901e-03 9.999e-05
This method may not be relevant when the coefficient of variation is greater than 1. Indeed, the number of successful trials is often zero due to the insufficient number of trials (and increasing n will be compensating by a decrease in p). But perhaps, XP will be OK.
We assume that there is a minimum value m and a rate lambda. The mean is mu=1/lambda+m and the variance is 1/lambda2. Hence, lambda=1/(mu x cov) and m=mu(1-cov).
rexp_cost <- function(n, mean, sigma) {
if (mean == 0 || sigma == 0)
return(rep(mean, n))
cov <- sigma/mean
lambda <- 1/(mean * cov)
m <- mean * (1 - cov)
return(rexp(n, lambda) + m)
}
mean(rexp_cost(10000, 10, 1))
## [1] 10
sqrt(var(rexp_cost(10000, 10, 1)))
## [1] 0.9586
rexp_cost(10, 10, 1)
## [1] 9.962 10.489 9.214 9.242 9.132 10.880 10.382 10.761 9.372 9.195
Here again, a coefficient of variation greater than 1 would lead to irrelevant results as negative values would be drawn.
We assume that there is a maximum value M. The mean is M x alpha/(alpha+beta) and the variance is M x alpha x beta/(alpha+beta)2/(alpha+beta+1). We want to fix the mean and alpha to a set of values. Then, we want to determine M and beta given the coefficient of variation. Maxima is needed to provide a correct derivation.
assume(alpha > 0, beta > 0, mu >0, cov > 0, M > 0);
declare(alpha, constant);
declare(mu, constant);
declare(cov, constant);
mean: mu = M*alpha/(alpha+beta);
variance: cov^2 = M^2*alpha*beta/(alpha+beta)^2/(alpha+beta+1)/mu^2;
solve([mean, variance], [beta, M]);
rbeta_cost <- function(n, mean, sigma, alpha) {
if (mean == 0 || sigma == 0)
return(rep(mean, n))
cov <- sigma/mean
M <- (cov^2 + 1) * mean/(1 - alpha * cov^2)
beta <- (M/mean - 1) * alpha
return(rbeta(n, alpha, beta) * M)
}
mean(rbeta_cost(10000, 10, 1, 50))
## [1] 10
sqrt(var(rbeta_cost(10000, 10, 1, 50)))
## [1] 0.9978
mean <- 10
cov <- 0.1
alpha <- 50
M <- (cov^2 + 1) * mean/(1 - alpha * cov^2)
beta <- (M/mean - 1) * alpha
curve(dbeta(x, 50, beta))
We will finally determine the set of values alpha such that relevant ratios of alpha/beta are tested (those between 0.01 and 1). First, we define the range of admissible values for alpha: alpha < 1/cov2.
cov <- c(0.01, 0.1, 1, 10)
1/cov^2
## [1] 1e+04 1e+02 1e+00 1e-02
We can use alpha=1000 up until cov=0.03, alpha=100 until cov=0.1, … and alpha=0.01 until the end (mention in the final article that all those choices cover a large sample of distributions and are likely to be representative). This represents 6 distributions (let's start with 3), plus the previous 3.
And of course, the previous distribution for checking the results.
rgamma_cost <- function(n, mean, sigma, alpha) {
if (mean == 0 || sigma == 0)
return(rep(mean, n))
shape <- mean^2/sigma^2
rate <- mean/sigma^2
return(rgamma(n, shape, rate))
}
mean(rgamma_cost(10000, 10, 1))
## [1] 9.978
sqrt(var(rgamma_cost(10000, 10, 1)))
## [1] 1.013
We take the previous simulation and plotting code:
get_all_makespans <- function(iterations, n, rdist, cov_min = 0.01, cov_max = 10) {
V_cov <- vector(length = iterations)
V_perf <- vector(length = 4 * iterations)
for (i in 1:iterations) {
set.seed(i)
cov <- exp(runif(1, log(cov_min), log(cov_max)))
V_cov[i] <- cov
rdist_cov <- function(n) return(rdist(n, 1, cov))
rdist_zero <- function(n) return(rep(0, n))
set.seed(i)
bino <- static_order_MC(tree_build(n, 1, 0), rdist_cov, rdist_zero,
1)
set.seed(i)
fibo <- static_order_MC(tree_build(n, 1, 1), rdist_cov, rdist_zero,
1)
set.seed(i)
dyn <- dyn_MC(n, rdist_cov, rdist_zero, 1)
set.seed(i)
dyn_com <- dyn_com_MC(n, rdist_cov, rdist_zero, 1)
V_perf[1:4 + (i - 1) * 4] <- c(bino, fibo, dyn, dyn_com)
}
V_cov <- rep(V_cov, each = 4)
V_method <- rep(factor(c("bino", "fibo", "dyn", "dyn_com")), times = iterations)
V <- data.frame(cov = V_cov, method = V_method, perf = V_perf)
return(V)
}
plot_all_ribbons <- function(V, limits = NULL) {
library(ggplot2)
library(plyr) # For rounding values
library(scales) # For comma formatter
# For grouping the statistic summary
V$pcov <- 10^round_any(log10(V$cov), 0.1)
tg <- ddply(V, c("pcov", "method"), summarise, pperf = median(perf), perf_inf = quantile(perf,
0.1), perf_sup = quantile(perf, 0.9))
tg$cov <- tg$pcov
tg$perf <- tg$pperf
p <- ggplot(data = tg, aes(x = cov, y = perf, group = method, shape = method,
fill = method))
p <- p + geom_ribbon(aes(ymin = perf_inf, ymax = perf_sup), alpha = I(1/5))
p <- p + geom_line(aes(color = method), alpha = I(1/2))
p <- p + geom_point(aes(color = method), alpha = I(1/2))
p <- p + geom_point(aes(color = method), color = alpha("black", 0.1))
p <- p + scale_x_log10(breaks = unique(10^round_any(log10(V$cov), 1)), labels = comma,
limits = limits)
p <- p + scale_y_log10(breaks = c(5, 10, 20, 50))
name <- "Method"
breaks <- c("bino", "fibo", "dyn", "dyn_com")
labels <- c("Binomial-stat", "Fibonacci-stat", "Tree-dyn", "Non-Commut-Tree-dyn")
p <- p + scale_colour_hue(name = name, breaks = breaks, labels = labels)
p <- p + scale_shape(name = name, breaks = breaks, labels = labels)
p <- p + scale_fill_hue(name = name, breaks = breaks, labels = labels)
p <- p + labs(x = "Coefficient of variation", y = "Schedule length", title = "Effect of the CV on the performances")
p <- p + coord_cartesian(ylim = exp(extendrange(log(tg$perf))))
p <- p + theme_bw()
p <- p + theme(legend.position = "bottom")
p <- p + annotation_logticks()
return(p)
}
Now, let's generate some figures:
n <- 64
iterations <- 1e+05 # 7*10 seconds for 1e2 -> 20 hours with 1e5
V_bern <- get_all_makespans(iterations, n, rbern_cost)
V_bino <- get_all_makespans(iterations, n, rbino_cost)
V_exp <- get_all_makespans(iterations, n, rexp_cost)
V_beta0_01 <- get_all_makespans(iterations, n, function(n, mean, sigma) return(rbeta_cost(n,
mean, sigma, 0.01)))
V_beta1 <- get_all_makespans(iterations, n, function(n, mean, sigma) return(rbeta_cost(n,
mean, sigma, 1)), cov_max = 1)
V_beta100 <- get_all_makespans(iterations, n, function(n, mean, sigma) return(rbeta_cost(n,
mean, sigma, 100)), cov_max = 0.1)
V_gamma <- get_all_makespans(iterations, n, rgamma_cost)
And let's plot them:
remove_min_plot <- function(V_dist, limits = NULL) {
V <- V_dist
# Zero value must be removed because of log scale (hopefully, they are
# rare enough to have no impact on median)
V[V$perf == min(V$perf), "perf"] <- unique(sort(V$perf))[2]
plot_all_ribbons(V, limits = limits)
}
remove_min_plot(V_bern)
remove_min_plot(V_bino)
remove_min_plot(V_exp)
remove_min_plot(V_beta0_01)
remove_min_plot(V_beta1, limits = exp(expand_range(log(c(0.01, 10)), 0.1)))
remove_min_plot(V_beta100, limits = exp(expand_range(log(c(0.01, 10)), 0.1)))
remove_min_plot(V_gamma)
Now, the same for the next figures:
get_all_makespans_over <- function(iterations, n, rdist, nb_pointx, nb_pointy = nb_pointx,
cov_min = 0.01, cov_max = 10) {
cd <- exp(seq(log(0.1), log(10), length.out = nb_pointy))
cov <- exp(seq(log(cov_min), log(cov_max), length.out = nb_pointx))
V_over <- data.frame()
for (i in 1:length(cd)) for (j in 1:length(cov)) {
rdistd <- function(n) return(rdist(n, 1, cov[j]))
rdistc <- function(n) return(rdist(n, cd[i], cd[i] * cov[j]))
set.seed(i + j * length(cov))
fibo <- static_order_MC(tree_build(n, 1, 1), rdistd, rdistc, iterations)
set.seed(i + j * length(cov))
dyn <- dyn_MC(n, rdistd, rdistc, iterations)
V_over <- rbind(V_over, data.frame(cd = cd[i], cov = cov[j], perf = mean(fibo)/mean(dyn),
perf_old = median(fibo/dyn)))
}
return(V_over)
}
plot_all_over <- function(V, limits = NULL) {
library(ggplot2)
library(plyr) # For rounding values
library(scales) # For darker lines and comma formatter
library(directlabels)
p <- ggplot(V, aes(cov, cd, fill = perf, z = perf))
p <- p + geom_tile()
p <- p + scale_x_log10(breaks = unique(10^round_any(log10(V$cov), 1)), labels = comma,
limits = limits)
p <- p + scale_y_log10(breaks = c(0.1, 0.2, 0.5, 1, 2, 5, 10), labels = percent)
p <- p + stat_contour(aes(color = ..level..), color = alpha("black", 1/2),
breaks = seq(0, 2, 0.1))
p <- p + labs(x = "Coefficient of variation", y = "Proportion of computation costs",
title = "Effect of the CV and computation costs\non the Fibonacci-stat performance")
p <- p + scale_fill_continuous("Performance ratio")
p <- p + theme_bw()
p <- p + theme(legend.position = "bottom")
p <- p + annotation_logticks()
p <- p + scale_color_continuous(low = "#000000", high = "#000000")
direct.label(p)
}
The computations:
n <- 64
iterations <- 1000 # 7*40 seconds with 1e1 & 11 points -> 100 hours with 1e3 & 33 points
nb_point <- 9 # odd because of symmetry and divisible by 3 for beta distribution
V_over_bern <- get_all_makespans_over(iterations, n, rbern_cost, nb_point)
V_over_bino <- get_all_makespans_over(iterations, n, rbino_cost, nb_point)
V_over_exp <- get_all_makespans_over(iterations, n, rexp_cost, nb_point)
V_over_beta0_01 <- get_all_makespans_over(iterations, n, function(n, mean, sigma) return(rbeta_cost(n,
mean, sigma, 0.01)), nb_point, cov_max = 9.99)
V_over_beta1 <- get_all_makespans_over(iterations, n, function(n, mean, sigma) return(rbeta_cost(n,
mean, sigma, 1)), nb_point * 2/3, nb_point, cov_max = 0.999)
V_over_beta100 <- get_all_makespans_over(iterations, n, function(n, mean, sigma) return(rbeta_cost(n,
mean, sigma, 100)), nb_point/3, nb_point, cov_max = 0.0999)
V_over_gamma <- get_all_makespans_over(iterations, n, rgamma_cost, nb_point)
and the plots:
plot_all_over(V_over_bern)
## Loading required package: grid
## Loading required package: quadprog
## Loading required package: proto
plot_all_over(V_over_bino)
plot_all_over(V_over_exp)
plot_all_over(V_over_beta0_01)
plot_all_over(V_over_beta1, limits = exp(expand_range(log(c(0.01, 10)), 0.1)))
plot_all_over(V_over_beta100, limits = exp(expand_range(log(c(0.01, 10)), 0.1)))
plot_all_over(V_over_gamma)
Now, the same for the next figures:
get_all_makespans_non_comm <- function(iterations, n, rdist, nb_pointx, nb_pointy = nb_pointx,
cov_min = 0.01, cov_max = 10) {
cd <- exp(seq(log(0.1), log(10), length.out = nb_pointy))
cov <- exp(seq(log(cov_min), log(cov_max), length.out = nb_pointx))
V_non_comm <- data.frame()
for (i in 1:length(cd)) for (j in 1:length(cov)) {
rdistd <- function(n) return(rdist(n, 1, cov[j]))
rdistc <- function(n) return(rdist(n, cd[i], cd[i] * cov[j]))
set.seed(i + j * length(cov))
bino <- static_order_MC(tree_build(n, 1, 0), rdistd, rdistc, iterations)
set.seed(i + j * length(cov))
fibo <- static_order_MC(tree_build(n, 1, 1), rdistd, rdistc, iterations)
set.seed(i + j * length(cov))
dyn_com <- dyn_com_MC(n, rdistd, rdistc, iterations)
# Comparing performance item-wise makes sense as the random costs are the
# same across all methods
b1 <- length(which(bino < fibo & bino < dyn_com))
b2 <- length(which(fibo < bino & fibo < dyn_com))
b3 <- length(which(dyn_com < bino & dyn_com < fibo))
best_old <- which.max(c(b1, b2, b3))
best_old <- factor(c("bino", "fibo", "dyn_com")[best_old])
# Ultimately, the mean could be sufficient
best <- which.min(c(mean(bino), mean(fibo), mean(dyn_com)))
best <- factor(c("bino", "fibo", "dyn_com")[best])
V_non_comm <- rbind(V_non_comm, data.frame(cd = cd[i], cov = cov[j],
best = best, best_old = best_old))
}
return(V_non_comm)
}
plot_all_non_comm <- function(V, limits = NULL) {
library(ggplot2)
library(plyr) # For rounding values
library(scales) # For comma formatter
p <- ggplot(V, aes(cov, cd, fill = best))
p <- p + geom_tile()
p <- p + scale_x_log10(breaks = unique(10^round_any(log10(V$cov), 1)), labels = comma,
limits = limits)
p <- p + scale_y_log10(breaks = c(0.1, 0.2, 0.5, 1, 2, 5, 10), labels = percent)
p <- p + labs(x = "Coefficient of variation", y = "Proportion of computation costs",
title = "Expected best method with\nvarying overlapping and variability")
breaks <- c("bino", "fibo", "dyn", "dyn_com")
labels <- c("Binomial-\nstat", "Fibonacci-\nstat", "Tree-dyn", "Non-Commut-\nTree-dyn")
p <- p + scale_fill_grey("Method", breaks = breaks, labels = labels)
p <- p + theme_bw()
p <- p + theme(legend.position = "bottom")
p <- p + annotation_logticks()
return(p)
}
The computations:
n <- 64
iterations <- 1000 # 7*80 seconds with 1e1 & 11 points -> 140 hours with 1e3 & 33 points
nb_point <- 9 # odd because of symmetry and divisible by 3 for beta distribution
V_non_comm_bern <- get_all_makespans_non_comm(iterations, n, rbern_cost, nb_point)
V_non_comm_bino <- get_all_makespans_non_comm(iterations, n, rbino_cost, nb_point)
V_non_comm_exp <- get_all_makespans_non_comm(iterations, n, rexp_cost, nb_point)
V_non_comm_beta0_01 <- get_all_makespans_non_comm(iterations, n, function(n,
mean, sigma) return(rbeta_cost(n, mean, sigma, 0.01)), nb_point, cov_max = 9.99)
V_non_comm_beta1 <- get_all_makespans_non_comm(iterations, n, function(n, mean,
sigma) return(rbeta_cost(n, mean, sigma, 1)), nb_point * 2/3, nb_point,
cov_max = 0.999)
V_non_comm_beta100 <- get_all_makespans_non_comm(iterations, n, function(n,
mean, sigma) return(rbeta_cost(n, mean, sigma, 100)), nb_point/3, nb_point,
cov_max = 0.0999)
V_non_comm_gamma <- get_all_makespans_non_comm(iterations, n, rgamma_cost, nb_point)
and the plots:
plot_all_non_comm(V_non_comm_bern)
plot_all_non_comm(V_non_comm_bino)
plot_all_non_comm(V_non_comm_exp)
plot_all_non_comm(V_non_comm_beta0_01)
plot_all_non_comm(V_non_comm_beta1, limits = exp(expand_range(log(c(0.01, 10)),
0.1)))
plot_all_non_comm(V_non_comm_beta100, limits = exp(expand_range(log(c(0.01,
10)), 0.1)))
plot_all_non_comm(V_non_comm_gamma)
This is positive as most of the figures are consistent. We just need to find a way to aggregate the results.