Michael Balch asserts this problem is solved, saying “[I]f X and Y are independent and have underlying marginal distributions conservatively represented by two Dempster-Shafer structures, then the Cartesian product formulation, for any operator g(X,Y), will yield a conservative Dempster-Shafer structure on Z = g(X,Y). See section 7.5.3 of my dissertation for proof. This theorem applies whether or not X or Y straddle zero. The question of zero-straddling does not come up in the proof. The proof does not even stipulate that X and Y need to be real-valued variables.”
Yager (1986) described a way to do arithmetic with Dempster-Shafer structures on the real line under the assumption of independence between the operands. Williamson and Downs (1990, pages 127 and 145) described essentially the same algorithm for horizontal discretizations of p-boxes. Here’s how it works. Consider these two p-boxes A and B:
These p-boxes can be sliced horizontally into pieces. In this case, we’ll slice them into three pieces each, but, in general, we could use any number of slices.
These slices represent Dempster-Shafer structures in that they are lists of intervals with corresponding probability masses that add to one. In his case, A = {([1,3], 1/3), ([2,4], 1/3), ([3,5], 1/3)}, and B = {([2,8], 1/3), ([6,10], 1/3), ([8,12], 1/3)}. The Yager algorithm is just to take the Cartesian product of such discretizations, where the intervals propagate according to the mathematical operation and, so long as A and B are independent, the probabilities multiply. If we want the sum A+B, we’d form the following Cartesian product:
The result is another list of intervals with associated probability masses which sum to one, thus, the result is a Dempster-Shafer structure which we can reassemble into a probability box like this:
Below are slow and fast implementations for R of this algorithm for use with p-boxes x and y. The slow version is easier to understand, but use the fast version if you actually want to compute, because it is much faster.
SLOWconv.pbox <- function(x,y) {
n = length(x@d)
c <- rep(0,n^2)
Zu <- rep(1,n)
Zd <- rep(1,n)
for (i in 1:n) for (j in 1:n) c[[(i-1)*n+j]] <- x@d[[i]] + y@d[[j]]
c <- sort(c)
for (k in 1:n) Zd[[k]] <- c[[(k-1)*n + n]]
for (i in 1:n) for (j in 1:n) c[[(i-1)*n+j]] <- x@u[[i]] + y@u[[j]]
c <- sort(c)
for (k in 1:n) Zu[[k]] <- c[[(k-1)*n + 1]]
pbox(u = Zu, d = Zd)
}
FASTconv.pbox <- function(x,y,op='+') {
if (op=='-') return(FASTconv.pbox(x,negate.pbox(y),'+'))
if (op=='/') return(FASTconv.pbox(x,reciprocate.pbox(y),'*'))
n = length(x@d)
Zu <- rep(1,n)
Zd <- rep(1,n)
Y <- rep(y@d, rep.int(n, n))
c <- sort(do.call(op,list(x@d,Y))) #c <- sort(x@d + Y)
k <- 1:n
Zd <- c[(k-1)*n + n]
Y <- rep(y@u, rep.int(n, n))
c <- sort(do.call(op,list(x@u,Y))) #c <- sort(x@u + Y)
Zu <- c[(k-1)*n + 1]
pbox(u = Zu, d = Zd)
}
# example with precise distributions
x = normal(5,1)
y = uniform(2,3)
z = FASTconv.pbox(x,y)
Williamson and Downs (1990, page 139, first paragraph of "Use of Mixtures for Nonmonotonic Operations") seem to specifically disclaim the use of the Cartesian product algorithm for computing the multiplicative or divisive sigma-convolution (i.e., assuming independence) whenever the operands straddle zero (i.e., are not sign-definite). We know the Williamson and Downs algorithms for Fréchet multiplication and division
don't work if either operand straddles zero, but we believe that our algorithms for both independent multiplication and independent division do work. We need to prove (or disprove) this belief.
The R code below gives some examples of precise distributions that straddle zero for which the Cartesian product algorithm seems to conform very well to the results of a Monte Carlo simulation for all four basic arithmetic operations. The p-boxes (shown in red and black) are well overwritten by the Monte Carlo distributions (shown in green). If it’s true for precise distributions, will it also be true for p-boxes?
source('ESTpbox.r')
pb = function(xs,x1,x2) {
if (xs=='u') return(U(x1,x2))
if (xs=='n') return(N(x1,x2))
if (xs=='l') return(logistic(x1,x2))
if (xs=='p') return(laplace(x1,x2))
if (xs=='g') return(gumbel(x1,x2))
if (xs=='c') return(cauchy(x1,x2))
}
pd = function(xs,x1,x2) {
if (xs=='u') return(runif(many,x1,x2))
if (xs=='n') return(rnorm(many,x1,x2))
if (xs=='l') return(rlogis(many,x1,x2))
if (xs=='p') return(qlaplace(runif(many),x1,x2))
if (xs=='g') return(qgumbel(runif(many),x1,x2))
if (xs=='c') return(qcauchy(runif(many),x1,x2))
}
shapes = c('u','n','l','p','g','c')
pbox.steps = 200
many = 10000
xs = shapes[order(runif(length(shapes)))][[1]]
x1 = runif(1,-10,2)
x2 = runif(1, 1,5)
ys = shapes[order(runif(length(shapes)))][[1]]
y1 = runif(1,-10,2)
y2 = runif(1, 1,5)
cat(xs,x1,x2,ys,y1,y2,'\n')
x = pb(xs,x1,x2)
y = pb(ys,y1,y2)
X = pd(xs,x1,x2)
Y = pd(ys,y1,y2)
old.par <- par(mfrow=c(2,2),mar=c(2,2,3,1))
z = x %|+|% y
plot(z)
Z = X + Y
edf(Z)
z = x %|-|% y
plot(z)
Z = X - Y
edf(Z)
z = x %|*|% y
plot(z)
Z = X * Y
edf(Z)
z = x %|/|% y
plot(z)
Z = X / Y
edf(Z)
par(old.par)
We need to prove (or disprove) that Yager’s Cartesian product works for computing product and quotient p-boxes from operands that straddle zero (despite the pessimism of Williamson and Downs).