computeBound: Solves programs (5) and (EC.19) over 2-point masses...

In cmottet/RobustTail: RobustTail

Description

This function solves programs (5) and (EC.19), in the case when their respective feasible region is restricted to distribution functions with at most two point supports.

Usage

computeBound(H, mu, sigma, lambda, nu = 1, direction = c("max", "min"))

Arguments

H	The function defined in the objective value of program (5) and (EC.19)
mu	Either an ordered vector containing bounds of the first moment in program (EC.19), or a scalar value as in program (5)
sigma	Either an ordered vector containing bounds of the second moment in program (EC.19), or a scalar value as in program (5)
lambda	A real number giving the limit value of the ratio H(x)/x^2 when x goes to infinity
nu	A real number as defined in program (5) and (EC.19) (actually denoted nubar in the latter case)
direction	A string either min or max identifying the type of wether program (5) should be a min or a max program. Default is max.

Value

bound	the optimal objective value
P	the optimal distribution function with point masses p = (p1,p2) and point supports x = (x1,x2)

Examples

##############################################################################
#### Solving a program alike  to (5)
##############################################################################
####
#### max P(X > c)
#### s.t. sum(px) = 1
####      sum(px^2) = 2
####      (p,x) is a two point mass distribution function where x => 0
####
#### where c is some positive number. We point out that the solution to this
#### problem is given in Theorem 3.3 of
#### "Optimal Inequalities in Probability Theory: a Convex Optimization Approach"
#### by D. Bertsimas and I. Popescu.
####
################################################################################

c <- qexp(0.9)
H <- function(x) as.numeric(c <= x)
mu <- 1
sigma <- 2
lambda <- 0

output <- computeBound(H,mu, sigma, lambda)

# Check that optimal upper bound is equal to  the analytical solution
CMsquare <- (sigma- mu^2)/mu^2
delta <-  c/mu-1
data.frame(Algorithm = output$bound, Analytical = CMsquare/(CMsquare + delta^2))

# Check that the output is feasible
with(output$P, data.frame(moments = c(sum(p),sum(p*x), sum(p*x^2)), truth = c(1,mu,sigma) ))

##############################################################################
#### Solving a program alike  to (EC.19)
##############################################################################
####
#### max P(X > c)
#### s.t. sum(p)  = 1
####      1 <= sum(px) <= 1
####      2 <= sum(px^2) <= 2
####      (p,x) is a two point mass distribution function where x => 0
####
#### where c is some positive number. This program is the same as the first one
#### but formulated with inequalities rather than equalities.
####
################################################################################
mu <- c(1,1)
sigma <- c(2,2)
lambda <- 0

output <- computeBound(H,mu, sigma, lambda)

# Check that optimal upper bound is equal to  the analytical solution
data.frame(Algorithm = output$bound, Analytical = CMsquare/(CMsquare + delta^2))

# Check that the output is feasible
with(output$P, data.frame(lowerMomentBound = c(1, mu[1], sigma[1]), moments = c(sum(p),sum(p*x), sum(p*x^2)), upperMomentBound = c(1, mu[2], sigma[2]) ))

##############################################################################
#### Solving a program alike  to (1)
##############################################################################
####
#### max P(x >b)
#### s.t. f(a) = eta
####      1-F(a) = beta
####      f'(a) => -nu
####      f(x) is convex for all x => a
####      f(x) is non-negative for all x => a
####
#### where b is some real number larger than a. This problem is of the form of
#### program (1). By Theorem 4, it is equivalent to solve
####
#### max sum(p H(x))
#### s.t sum(px) = mu
####     sum(px^2) = sigma
####     (p,x) is a two point mass distribution function where x => 0
####
#### where mu = eta/nu, sigma = beta/nu, H(x) = 1/2(x + a -b)I(x =a =>b)
####
################################################################################

# Assume the true distribution function is a standard exponential
a <- qexp(0.7)
eta <- dexp(a)
nu <- dexp(a)
beta <- 1-pexp(a)
mu <- eta/nu
sigma <- 2*beta/nu
lambda <- 1/2

# Compute the optimal upper bound for various valus of b
b <- qexp(seq(0.7,0.99, by = 0.01))

runFunc <- function(b){
H <- function(x) 1/2*(x + a - b)^2*( x +a >= b)
bound <- RobustTail::computeBound(H,mu,sigma,lambda,nu)$bound
output <- data.frame(b = b, bound = bound)
}

dataPlot <- plyr::ldply(lapply(X = b, FUN = runFunc))

library(ggplot2)
ggplot(dataPlot, aes(x = b, y = bound)) +
geom_line() +
labs(y = "Optimal Upper Bound") +
ylim(c(0,0.3))

cmottet/RobustTail documentation built on May 13, 2019, 8:44 p.m.