  A pair of lower and upper cumulative distribution functions, also called
probability box or p-box, is among the most popular models used in imprecise
probability theory. They arise naturally in expert elicitation, for instance in
cases where bounds are specified on the quantiles of a random variable, or when
quantiles are specified only at a finite number of points. Many practical and
formal results concerning p-boxes already exist in the literature. In this
paper, we provide new efficient tools to construct multivariate p-boxes and
develop algorithms to draw inferences from them. For this purpose, we formalise
and extend the theory of p-boxes using Walley's behavioural theory of imprecise
probabilities, and heavily rely on its notion of natural extension and existing
results about independence modeling. In particular, we allow p-boxes to be
defined on arbitrary totally preordered spaces, hence thereby also admitting
multivariate p-boxes via probability bounds over any collection of nested sets.
We focus on the cases of independence (using the factorization property), and
of unknown dependence (using the Fr\'echet bounds), and we show that our
approach extends the probabilistic arithmetic of Williamson and Downs. Two
design problems---a damped oscillator, and a river dike---demonstrate the
practical feasibility of our results.
