Figure illustrating the distributions (April 30, 2014)

We plot on one figure all the distributions that are used with a CV of 0.5 (or 0.4).

library(ggplot2)

dgamma_cost <- function(x, mean, sd) {
    shape <- mean^2/sd^2
    rate <- mean/sd^2
    return(dgamma(x, shape, rate))
}

dbern_cost <- function(x, mean, sd) {
    cov <- sd/mean
    p <- 1/(cov^2 + 1)
    M <- mean/p
    return(dbinom(x/M, 1, p)/M)
}

dbino_cost <- function(x, mean, sd) {
    cov <- sd/mean
    size <- 100
    p <- 1/(1 + cov^2 * size)
    M <- mean/(size * p)
    return(dbinom(x/M, size, p)/M)
}

dexp_cost <- function(x, mean, sd) {
    cov <- sd/mean
    lambda <- 1/(mean * cov)
    m <- mean * (1 - cov)
    return(dexp(x - m, lambda))
}

dpois_cost <- function(x, mean, sd) {
    lambda <- sd^2
    m <- mean - lambda
    return(dpois(x - m, lambda))
}

dunif_cost <- function(x, mean, sd) {
    min <- mean - sd * sqrt(12)/2
    max <- mean + sd * sqrt(12)/2
    return(dunif(x, min, max))
}

library(triangle)
dtriangle_cost <- function(x, mean, sd) {
    a <- mean - sqrt(6) * sd
    b <- mean + sqrt(6) * sd
    return(dtriangle(x, a, b))
}

dbeta_cost <- function(x, mean, sd, alpha) {
    cov <- sd/mean
    M <- (cov^2 + 1) * mean/(1 - alpha * cov^2)
    beta <- (M/mean - 1) * alpha
    return(dbeta(x/M, alpha, beta)/M)
}

mean <- 1
sd <- 0.5

x <- seq(-1, 5, 0.001)
suppressWarnings(dists <- data.frame(x = x, bino = dbino_cost(x, mean, sd), 
    gamma = dgamma_cost(x, mean, sd), exp = dexp_cost(x, mean, sd), unif = dunif_cost(x, 
        mean, sd), triangle = dtriangle_cost(x, mean, sd), pois = dpois_cost(x, 
        mean, sd), beta1 = dbeta_cost(x, mean, sd, 1), beta0_01 = dbeta_cost(x, 
        mean, sd, 0.01), bern = dbern_cost(x, mean, sd)))

library(reshape)

dists_data <- melt(dists, id = "x")
dists_data$continuous <- dists_data$variable %in% c("gamma", "exp", "unif", 
    "beta0_01", "beta1")
p <- ggplot(dists_data, aes(x = x, y = value, color = variable, size = variable, 
    linetype = variable))
p <- p + facet_wrap(~continuous)
p <- p + geom_line(alpha = I(1/2))
p <- p + scale_size_discrete(range = c(0.5, 1))
p <- p + coord_cartesian(xlim = c(-0.5, 3), ylim = c(0, 2))
p

plot_dist <- function(p, ymax) {
    p <- p + geom_line(size = 5)
    p <- p + scale_size_discrete(range = c(0.5, 1))
    p <- p + coord_cartesian(xlim = c(-0.5, 3), ylim = c(0, ymax))
    p <- p + theme(panel.background = element_blank(), panel.grid.minor = element_blank(), 
        panel.grid.major = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(), 
        axis.title.x = element_blank(), axis.title.y = element_blank())
    return(p)
}
plot_dist(ggplot(dists, aes(x = x, y = bern)), ymax = max(dists$bern))

plot_dist(ggplot(dists, aes(x = x, y = beta1)), ymax = min(3, max(dists$beta1)))

plot_dist(ggplot(dists, aes(x = x, y = beta0_01)), ymax = min(3, max(dists$beta0_01)))

plot_dist(ggplot(dists, aes(x = x, y = bino)), ymax = max(dists$bino))

plot_dist(ggplot(dists, aes(x = x, y = exp)), ymax = max(dists$exp))

plot_dist(ggplot(dists, aes(x = x, y = gamma)), ymax = max(dists$gamma))

plot_dist(ggplot(dists, aes(x = x, y = pois)), ymax = max(dists$pois))

plot_dist(ggplot(dists, aes(x = x, y = triangle)), ymax = max(dists$triangle))

plot_dist(ggplot(dists, aes(x = x, y = unif)), ymax = max(dists$unif))