---
address:
- 'The Alan Turing Institute, 96 Euston Road, London NW1 2DB'
- 'London School of Economics and Political Science, Houghton Street, London, WC2A 2AE'
author:
- Nikki Sonenberg
- Edward Wheatcroft
- Henry Wynn
bibliography:
- 'References.bib'
title: Majorisation as a theory for uncertainty
---

Introduction {#sec:intro}
============

Majorisation, also called rearrangement inequalities, yields a type of stochastic ordering in which two or more distributions are rearranged in decreasing order of their probability mass (discrete case) or probability density (continuous case) and then compared. When methods of majorisation are applied to probabilities and probability distributions, they provide a concept of the peakedness. This is independent of the ‘location’ of the probabilities, i.e., of the support of the distribution. This geometry-free property makes majorisation a good candidate as a foundation for the idea of uncertainty which is the focus of this paper. Majorisation is a partial, not total, ordering and implies that one or more of the class of order-preserving functions with respect to the ordering might be used as an entropy or ‘uncertainty metric’. Many well-known metrics fall into this category, one of which is Shannon entropy, widely used in information theory.

Consider the question, ‘is uncertainty geometric?’. If we think our friend is in London, Birmingham or Edinburgh with probabilities $p_1,p_2, p_3$ respectively (where $p_1+p_2+p_3=1$), does it make any difference to our uncertainty, however we measure it, if the locations are changed to Reading, Manchester and Glasgow with the same probabilities, respectively? In fact, could we just permute the order of the first three cities? If our answer is no, i.e., that there is no difference, then we are in the realm of entropy and information. In the above cases, the Shannon entropy is $-\{p_1 \log(p_1) + p_2 \log(p_2) + p_3 \log(p_3)\}$ and we see that the subscript simply serves as a way to collect the probabilities, not to locate them in the geography of the UK.

Another element of the majorisation approach is that it is, in a well-defined sense, dimension-free. In this paper, we show how this approach enables us to create, for a multivariate distribution, a univariate decreasing rearrangement (DR) by considering a decreasing threshold and ‘squashing’ all of the multivariate mass for which the density is above the threshold to a univariate mass adjacent to the origin. This creates the possibility of comparing multivariate distributions with different numbers of dimensions.

We introduce a set of operations that can be applied to study uncertainty in a range of settings and illustrate these with examples. We see this work as a merging of methods used in applied mathematics and statistics with general methodology for the study of uncertainty. The methods discussed provide a foundation for the extension to Bayesian probabilities, a topic for further work.

There is a large literature on majorisation. The classical results of Hardy, Littlewood and Polya [@Hardy1988] led to developments in a wide variety of fields. Marshall and Olkin’s [@Marshall2011] key volume on majorisation (later extended in [@Marshall2011]) built on these results. Applications in mathematical economics can be found in portfolio theory and income distributions built on classical work by Lorenz [@Lorenz1905] and Gini [@Gini1914] (see recent work by Arnold and Sarabia [@Arnold2018]). Majorisation has also been used in chemistry for mixing liquids and powders [@Klein1997] and in quantum information [@Partovi2011]. Statistical applications include experimental design [@Giovagnoli1987; @Pukelsheim1987], and in application to testing [@eaton1974monotonicity; @tong1988some]. Majorisation has been employed in the area of proper scoring rules by considering the partial and total ordering in the class of well-calibrated experts [@Degroot1986; @Degroot1985; @Degroot1988]. The common feature of these studies is the need to compare and quantify the degree of variation between distributions. We note that the theory of the decreasing rearrangement of functions, which we have used to a limited extent for probability densities, can be considered an area of functional analysis particularly in the area of inequalities of the kind which say that a rearrangement of a function increases or decreases some special functional [@lieb2001graduate].\
This paper is organised as follows. In the remainder of this section we introduce the concept of majorisation of probabilities and present related concepts and previous work. In Section \[sec:cont\_major\], we present results for the continuous case. In Section \[sec:multivariate\], we present the key idea of reducing multivariate distributions into a one dimensional decreasing rearrangement and illustrate this with examples. In Section \[sec:operations\], we collect together operations for the study of uncertainty and, in Section \[sec:algebra\], a lattice and an algebra for uncertainty. In Section \[sec:empirical\], we discuss empirical applications. Concluding remarks are given in Section \[sec:conclusion\].

Discrete majorisation and related work {#sec:discrete}
--------------------------------------

We introduce majorisation for discrete distributions following Marshall *et al.* [@Marshall2011].

Consider two discrete distributions with $n$-vectors of probabilities $$p_1=(p^{(1)}_1,\ldots, p_n^{(1)}) \quad \text{ and } \quad p_2=(p^{(2)}_1,\ldots, p_n^{(2)}),$$ where $\sum_i p_i^{(1)}=\sum_i p_i^{(2)}=1$. Placing the probabilities in decreasing order: $$\tilde{p}^{(1)}_1 \geq \ldots \geq \tilde{p}_n^{(1)}\quad \text{ and } \quad \tilde{p}^{(2)}_1 \geq \ldots \geq \tilde{p}_n^{(2)},$$ it is then said that $p_2$ majorises $p_1$, written $p_1 \preceq p_2$ when, for all $n$, $$\sum_{i=1}^n \tilde{p}_i^{(1)} \leq \sum_{i=1}^n \tilde{p}_i^{(2)}.$$

This definition of majorisation is a partial ordering, that is, not all pairs of vectors are comparable. As argued by Partovi [@Partovi2011], this is not a shortcoming of majorisation, rather a consequence of its rigorous protocol for ordering uncertainty. Marshall *et al* [@Marshall2011] provide several equivalent conditions to $p_1\preceq p_2$. We consider three (A1-A3) of the best known in detail below.\
(A1) There is a doubly stochastic $n\times n$ matrix $P$, such that $$\begin{aligned}
\label{equiv:doubly}
p_1 = P p_2.\end{aligned}$$ This is a well known result by Hardy, Littlewood and Pólya [@Hardy1988]. The intuition of this result is that a probability vector which is a mixture of the permutations of another is more disordered. The relationship between a stochastic matrix $P$ and the stochastic transformation function in the refinement concept was presented by DeGroot [@Degroot1988].\
(A2) Schur [@Schur1923] demonstrated that, if (A1) holds for some stochastic matrix $P$, this leads to the following equivalent condition. For all continuous convex functions $h( \cdot )$, $$\begin{aligned}
\label{condition3}
\sum_{i=1}^n h(\tilde{p}_i^{(1)}) \leq  \sum_{i=1}^n h(\tilde{p}_i^{(2)}),\end{aligned}$$ for all $n$.

The sums in (\[condition3\]) are special cases of the more general Schur-convex functions on probability vectors. Details on the characteristics and properties of Schur-convex functions are provided by Marshall *et al.* [@Marshall2011]. In particular, entropy functions such as Shannon information, for which $h(y)=y\log(y)$, are Schur-convex. We also highlight a special case of the Tsallis information for which $$h(y)=\frac{y^{\gamma}-1}{\gamma}, \quad\gamma>0,$$ where, in the limit $\gamma\rightarrow 0$, Shannon information is obtained. The condition (A2) is equivalent to the majorisation ordering for distributions, and we consider it as a continuous extension to Equation (\[condition3\]) (see Section \[sec:cont\_major\] for details). The condition (A2) indicates that the ordering imposed by majorisation is stronger than the ordering by any single entropic measure and, in a sense, is equivalent to all such (entropic) measures taken collectively [@Partovi2011].\
(A3) Let $\pi(p) = (p_{\pi(1)}, \ldots, p_{\pi(n)})$ be the vector whose entries are a permutation $\pi$ of the entries of a probability vector $p$, with symmetric group $S$, then $$\begin{aligned}
p_1 \in \mbox{conv}_{\pi \in S} (\{\pi(p_2)\}).\end{aligned}$$ That is to say, $p_1$ is in the convex hull of all permutations of entries of $p_2$. Majorisation is a special case of group-majorisation (G-majorisation) for the symmetric (permutation) group [@Giovagnoli1985]. The general type of groups for which the theory really works are reflection (Coxeter) groups, e.g. permutation and sign changes in $n$ dimensions (called the $B_n$ series). The main work on G-majorisation was by Eaton and Perlmann [@Eaton1977].

Rearrangements can be viewed as special instances of transportation plans, which move a given mass distribution to another distribution of the same total mass (see Buchard [@Burchard2009]). The recent use of transport mapping in UQ is to improve Markov chain Monte Carlo (MCMC) and sequential Monte Carlo (SMC) samplers [@Marzouk2016; @Parno2018].

Continuous majorisation {#sec:cont_major}
=======================

Following Hardy *et al.* [@Hardy1988], in this section we describe continuous majorisation, extending on the discrete case presented in Section \[sec:discrete\]. Ryff [@Ryff1965] provided continuous analogues to discrete majorisation by replacing vectors and matrices with integrable functions and linear operators.

\[drdefn\] Let $f(x)$ be a (univariate) pdf and define $m(y)=\mu\{z: f(z) \geq y\}$. The decreasing rearrangement of $f(x)$ is then $$\begin{aligned}
\tilde{f}(z)=\mbox{sup}\{t: m(t) >z\},\; z >0.\end{aligned}$$

Let $\tilde{f}_1(z)$ and $\tilde{f}_2(z)$ be the DR of two pdfs $f_1(x)$ and $f_2(x)$, respectively and $\tilde{F}_1(z)$ and $\tilde{F}_2(z)$ their corresponding cdfs. We say that $f_2(x)$ majorises $f_1(x)$, written $f_1 \preceq f_2$, if and only if $$\tilde{F}_1(z) \leq  \tilde{F}_2(z),\;\; \mbox{for all} \quad  z > 0.$$

Similarly to the discrete case, we give three equivalent conditions for majorisation, $\preceq$:\
(B1) For some non-negative doubly stochastic kernel $P(x,t)$, $$\begin{aligned}
f_1(x) = \int P(x,t) f_2(t) dt.\end{aligned}$$ (B2) For all continuous convex functions $h(\cdot)$, $$\begin{aligned}
\int h(f_1(z)) dz \leq \int h(f_2(z))dz.\end{aligned}$$ (B3) Slice condition: $$\begin{aligned}
\int(f_1(x)-c)_+ dx \leq \int(f_2(x)-c)_+dx, \quad c>0. \label{eq:slice}\end{aligned}$$

\[example|\_beta\]

Consider the Beta$(a,b)$ distribution with $(a, b)=(3, 2)$ and pdf $p(z)=12(1-z)z^2$. We look for $z_1$ and $z_2$ (where $z_1<z_2$) such that $p(z_1)=p(z_2)=c$, that is the points at which $p(z)$ intersects the line $y=c$. The pdf and the horizontal line $y=c$ are both shown in Figure \[fig:DensityPlot\]. Setting $z=z_2-z_1$, this gives us the system of equations: $$\begin{aligned}
\label{example_beta_system}
\begin{cases}
p(z_1)=12(1-z_1)z_1^2=y ,\\
p(z_2)=12(1-z_2)z_2^2=y,\\
z_2-z_1=z,\\
0\leq z\leq 1.
\end{cases}\end{aligned}$$

![The identification of $z_{1}$ and $z_{2}$ for the probability density function of Beta(3,2).[]{data-label="fig:DensityPlot"}](DensityFunctionBeta.png){width="34.00000%"}

The DR can be obtained from Equation (\[example\_beta\_system\]) by eliminating $z_1$ and $z_2$ and setting $\tilde{f} = y$. The elimination variety is $ 48z^6 - 96z^4 + 9y^2 + 48z^2 - 16y$. We obtain the explicit solution, $$\begin{aligned}
\tilde{f}(z) = \left\{ 
 \begin{array}{l}  
\frac{8}{9} 
+ \frac{4}{9} (-27z^6 + 54z^4 - 27z^2 + 4)^{\frac{1}{2}}, \quad 0 \leq z \leq \frac{1}{\sqrt{3}}  \\                                              \frac{8}{9} -  \frac{4}{9} (-27z^6 + 54z^4 - 27z^2 + 4)^{\frac{1}{2}}, \quad \frac{1}{\sqrt{3}} \leq z \leq 1.
\end{array}
\right.\end{aligned}$$

This variety is shown in Figure \[fig:beta\] and the DR $\tilde{f}(z)$ is the section of the curve over $[0,1]$ decreasing from $\left(0, \frac{16}{9}\right)$ to $(1,0)$. The cdf of the DR of Beta$(2,3)$ is $\tilde{F}(z)=-4z^3 + 6z^2$. If $\tilde{F}(z)$ is the cdf corresponding to $\tilde{f}$ obtained by adjoining the equations $Y = F(z_2)-F(z_1)$, we obtain the second variety $4z^4 + 3Y^2 - 8Yz = 0,$ illustrated in Figure \[fig:beta\], and $\tilde{F}(z)$ is the upper portion of the curve from $(0,0)$ to $(1,1)$, namely $$\tilde{F}(z) = 2\left(\frac{2}{3} + \frac{\sqrt{-3z^2 + 4}}{3}\right)z.$$

It is hard to derive when $f_1(x) \leq f_2(x)$ for the general case in which $(a_1,b_1)$ and $(a_2,b_2)$ and $f_i(x) \sim \mbox{Beta}(a_i,b_i),\; i=1,2$. However, we can prove the following.

Assume $a_1 , b_1 , a_2 , b_2 > 1$. If pdfs $f_1 (x)\sim \text{Beta}(a_1 , b_1 )$ and $f_2 (x)\sim\text{Beta}(a_2 , b_2 )$, have the same mode, then $\max_x  f_1(x) \leq  \max_x f_2(x)$ if and only if $X_1 \preceq X_2.$

We first prove that, under the same mode condition, $f_1(x)$ and $f_2(x)$ intersect at two distinct $x$-values at which the values of $f_1(x)$ and $f_2(x)$ are the same. Setting the modes equal, that is $$\frac{a_1-1}{a_1+b_1-2} = \frac{a_2-1}{a_2+b_2-2},$$

We find that, and assuming without loss of generality, that $a_2  > a_1$, we obtain $$\frac{f_1(x)}{f_2(x)} = \left\{x (1-x)^u \right\}^v C,$$

where $u= a_2-a_1,v = \frac{b_1-1}{b_2-1}$ and $C$ is a constant. If we set this equal to 1, we obtain two solutions given by $$x (1-x)^u = C^{-\frac{1}{v}}.$$ It is then straightforward to verify that the common value of $f_1(x)$ and $f_2(x)$ is the same at the two solutions. The proof of the theorem is completed by using the slice condition in Equation .

Multivariate case: matching of uncertainty {#sec:multivariate}
==========================================

The following construction shows how to induce a one dimensional DR from a multidimensional distribution. This is in fact the continuous multidimensional version of the following discrete version. Thus, suppose we have a 2-dimensional table with probabilities $p_{i,j}, i=1, \ldots, I; j= 1, \ldots, J$. Line up the probabilities in decreasing (non-increasing) order from $1$ to $n =IJ$ with support on $1, \ldots, n$, respectively. If we now count the number of $\{i,j\}$ such that the probability is greater than or equal to a constant $c$, it is the same as the original table as for the DR.

\[multitouni\] A univariate decreasing rearrangement $\tilde{f}(z)$, compatible with $f(x)$, is, for all constants $c\geq 0$, $$\begin{aligned}
\label{decreaseRearrangemult}
\int_{\{x:f(x)\geq c \}}f(x)dx=\int_{\{z:\tilde{f}(z)\geq c\}}\tilde{f}(z)dz.\end{aligned}$$

Following [@Burchard2009]: as $$\begin{aligned}
     {\{x:f(x)\geq c \}} = {\{z:\tilde{f}(z)\geq c\}},\end{aligned}$$ then the volume of these sets are consistent.

The notation on the lemma will be important to us. We shall use $f(x)$ for the multivariate pdf for a random variable $X = (X_1, \ldots, X_n)$ of dimension $n$, and $\tilde{f}(z)$ will be its one-dimensional pdf and we reserve $\tilde{F}(z)$ for the corresponding cdf, which will be the basis for the majorisation.

The following lemma shows that the information/entropy for $X \sim f(x)$ and $Z \sim \tilde{f}(z)$ are the same. This is a crucial result and gives us confidence in the matching.

Let $f(x)$ be a multidimensional pdf and $\tilde{f}(z)$ on $[0, \infty]$ its decreasing rearrangement. Then, given a convex function $\varphi(x)$, we have $$\int_S \varphi(f(x))  dx = \int_0^{\infty} \varphi(\tilde{f}(z)) dz,$$ where $S$ is the support of $f(x)$.

(sketch) The proof consists of matching volume to length elements in $S$ and $[0,1)$. For $c>0$ and small $\delta c  > 0$ we have $$\int_{x: f(x) \geq c, x \in S} f(x)dx -  \int_{x: f(x) \geq c + \delta c, x \in S} f(x) dx =  \int_{z: \tilde{f}(z) \geq c, z \in [0, \infty)} \tilde{f}(z) dz  -\int_{z: \tilde{f}(z) \geq c +\delta c, z \in [0, \infty)} \tilde{f}(z) dz.$$ We can then write, approximately, $$u(c)A(c, \delta c)=u(c)L(c, \delta c),$$ where $A(c, \delta c)$ and $L(c, \delta c)$ are the corresponding increments in volume and length, respectively, as corresponding to the interval $[c,c+\delta c)$, that is $ f^{(-1)}([c,c+ \delta))$ and $ \tilde{f}^{(-1)}([c,c+ \delta))$, respectively. Cancelling $c$, we can equate $A(c, \delta c)$ and $L(c, \delta c)$, and this allows us to recapture and equate the integrals of any measurable function $u(\cdot)$: $$u(c)A(c,\delta c)  =  u(c) L(c, \delta c).$$ In particular, we can write $u(c) = \varphi(f(c)).$

In Examples \[multi\_norm\_ex\] and \[indep\_exp\_ex\], we demonstrate how to obtain a DR for the standard multivariate Normal distribution and the $n$-fold independent standard exponential distribution. In the examples below, we use the following idea to carry out computations. There may be cases in which, for a given $c$, the inverse set $f^{(-1)}(c)$ is described by some useful quantity $\delta$. Moreover $\delta$, expressed as a function of $x$, then becomes a random variable with a known (univariate) distribution. In that case, we can write Definition \[multitouni\] as $$\tilde{F}\big(\tilde{f}^{-1}(c) \big) =F_{\delta}\big(f_{{X}}^{-1}(c) \big)\nonumber \label{Fstuff},$$ then $$\begin{aligned}
 \label{eq:valid_DR}
\tilde{f}(r) & =  f_{\delta} \left( f_X^{(-1)}(\tilde{f}(r)) \right) \frac{\partial}{\partial r}\left(f_X^{(-1)}(\tilde{f}(r))\right).\end{aligned}$$

\[multi\_norm\_ex\] We provide a representation of Definition \[multitouni\] for a two-dimensional multivariate normal in Figure \[fig:DRM2\].

![*Left panel:* Density plot of a two-dimensional standard multivariate normal. The dashed line and blue shaded region correspond to $f(x)=c$ and $\int_{\{x:f(x)\geq c \}}f(x)dx$ respectively. *Central panel:* A plot to demonstrate the connection between ${X}=(X_1, X_2)^T$ and $R^2$. The radius $r$ of a blue circle corresponds to $x=f^{-1}(c)$. *Right panel:* A DR $\tilde{f}(z)$ obtained for a multivariate normal. The blue shaded region corresponds to $\int_{\{z: \tilde{f}(z)\geq c\}}\tilde{f}(z) dz$.[]{data-label="fig:DRM2"}](MultivariateNorm1.png "fig:"){width="30.00000%"} ![*Left panel:* Density plot of a two-dimensional standard multivariate normal. The dashed line and blue shaded region correspond to $f(x)=c$ and $\int_{\{x:f(x)\geq c \}}f(x)dx$ respectively. *Central panel:* A plot to demonstrate the connection between ${X}=(X_1, X_2)^T$ and $R^2$. The radius $r$ of a blue circle corresponds to $x=f^{-1}(c)$. *Right panel:* A DR $\tilde{f}(z)$ obtained for a multivariate normal. The blue shaded region corresponds to $\int_{\{z: \tilde{f}(z)\geq c\}}\tilde{f}(z) dz$.[]{data-label="fig:DRM2"}](MultivariateNorm2.png "fig:"){width="30.00000%"} ![*Left panel:* Density plot of a two-dimensional standard multivariate normal. The dashed line and blue shaded region correspond to $f(x)=c$ and $\int_{\{x:f(x)\geq c \}}f(x)dx$ respectively. *Central panel:* A plot to demonstrate the connection between ${X}=(X_1, X_2)^T$ and $R^2$. The radius $r$ of a blue circle corresponds to $x=f^{-1}(c)$. *Right panel:* A DR $\tilde{f}(z)$ obtained for a multivariate normal. The blue shaded region corresponds to $\int_{\{z: \tilde{f}(z)\geq c\}}\tilde{f}(z) dz$.[]{data-label="fig:DRM2"}](MultiVariateNorm3.png "fig:"){width="30.00000%"}

Let a real random vector ${X}=(X_1, \dots, X_n)^T$ be an $n$-variate standard normal distribution with density $$f_{{X}}(x_1, \dots, x_n)=\frac{1}{(2\pi)^{\frac{n}{2}}}\exp\Bigg\{-\frac{\sum_{i=1}^n x_i}{2} \Bigg\}.$$ We refer to a real random vector ${X}$ as a spherical Gaussian random vector with ${X}\sim\text{N}_n({0}, I_n)$, where ${0}$ is an $n$-vector of zeros and $I_n$ is an $n\times n$ identity matrix. Define the square of the radius of a spherical Gaussian random vector, that is, $$R^2 = \sum_{i=1}^nX_i^2.$$ To construct a DR, we slice the pdf of a multivariate normal at $f_{{X}}(x_1, \dots, x_n)=c$, then $$\begin{aligned}
    \label{eq:DR_derivation1}
    &c = \frac{1}{(2\pi)^{n/2}}\exp\Big\{-\frac{1}{2}\sum_{i=1}^n x_i^2 \Big\}.
    \end{aligned}$$ Given the relationship between the vector ${X}$, and a random variable $R^2$, define $r^2 =\sum_{i=1}^n x_i^2$, then $$\begin{aligned}
        r=\Big(-2\log\big((2\pi)^{n/2} c\big) \Big)^{1/2},\end{aligned}$$ where the volume of the $n$-dimensional Euclidean ball of radius $r$ is $$\label{eq:ball_volume}
    V_n(r)=\frac{\pi^{n/2}}{\Gamma\Big(\frac{n}{2}+1 \Big)}r^n.$$ Substituting Equation (\[eq:DR\_derivation1\]) into Equation (\[eq:ball\_volume\]), we obtain $$\begin{aligned}
c=\frac{1}{(2\pi)^{n/2}}\exp\bigg\{-\frac{1}{2}\bigg(\frac{V_n(r)\Gamma(n/2+1)}{\pi^{n/2}} \bigg)^{2/n} \bigg\},\end{aligned}$$ noting the values of $c$ and $V_n(r)$ are dependent on each other. To generalise the expression above, we replace $c$ and $V_n(r)$ with $\tilde{f}(z)$ and $z$, respectively. The final form of the DR is $$\label{eq:DRM_mult_normal}
\tilde{f}(z) = \frac{1}{(2\pi)^{n/2}}\exp\bigg\{-\frac{1}{2}\bigg(\frac{z}{V_n} \bigg)^{2/n} \bigg\},$$ where $V_n$ is the volume of the unit sphere in $\mathbb{R}^n$.

We validate the form of the DR in Equation (\[eq:DRM\_mult\_normal\]) by employing the construction from Equation (\[eq:valid\_DR\]). Here, $R^2=\sum_{i=1}^n X_i^2$ follows a Chi-squared distribution with $n$ degrees of freedom, $R^2\sim \chi_n^2$, with pdf and cdf, $$\begin{aligned}
f_{R^2}(y)=\frac{1}{2^{n/2}\Gamma(n/2)}y^{n/2-1}\exp\left\{-\frac{y}{2} \right\}, \qquad
F_{R^2}(y)=\frac{1}{\Gamma(n/2)}\gamma\left(\frac{n}{2}, \frac{y}{2} \right),\end{aligned}$$ respectively. The result for the DR in Equation (\[eq:DRM\_mult\_normal\]) can be verified by employing Equation (\[eq:valid\_DR\]) with $\delta = R^2$ and directly substituting the required functions.

\[indep\_exp\_ex\] Let the real random vector ${X}=(X_1, \dots, X_n)^T$ be an $n$-fold independent standard exponential distribution with density $$f_X(x_1, \dots, x_n)=\exp\left\{-\sum_{i=1}^n x_i \right\}.$$ As $f_X(x_1, \dots, x_n) = f_1(x_1) f_2(x_2) \cdots f_n(x_n)$, slicing the pdf at $c=f_{{X}}(x_1, \dots, x_n)$ yields $$\begin{aligned}
 \label{eq:DR_derivation2}
    &-\log(c)=\sum_{i=1}^n x_i.\end{aligned}$$ The volume of an $n$-dimensional simplex in which all $n$ variables are greater than $0$ but with sum less than $R$ is $$\label{eq:simplex_volume}
V_n = \frac{R^n}{n!}.$$ Substituting Equation (\[eq:DR\_derivation2\]) into Equation (\[eq:simplex\_volume\]), we obtain $$\begin{aligned}
c=\exp\left\{-(n!V_n)^{1/n} \right\},\end{aligned}$$ where the values of $c$ and $V_n$ depend on each other. To generalise this expression, we replace $c$ and $V_n$ with $\tilde{f}(z)$ and $z$, respectively. The DR can then be written as $$\label{eq:DRM_mult_exp}
\tilde{f}(z)=\exp\left\{-(n!z)^{1/n} \right\}.$$ To verify the form of the DR in Equation (\[eq:DRM\_mult\_exp\]), we employ the construction in Equation (\[eq:valid\_DR\]). Define the random variable $R=\sum_{i=1}^n X_i$, such that $R\sim\text{Gamma}(n, 1)$, with pdf and cdf, $$\begin{aligned}
    f_R(y)=\frac{1}{\Gamma(n)}y^{n-1}e^{-y}, \qquad F_R(y)=\frac{1}{\Gamma(n)}\gamma(n, y).\end{aligned}$$ Similarly to Example \[multi\_norm\_ex\], we derive the DR, $\tilde{f}(z)$, for a Gamma-distributed random variable $R$, that is, $$\label{eq:DR_mult_exp}
\tilde{f}(r)=f_R\big(f_{{X}}^{-1}(\tilde{f}(r))\big)\frac{\partial}{\partial r}(f_{{X}}^{-1}(\tilde{f}(r))).$$

The results for the DR in Equation (\[eq:DRM\_mult\_exp\]) can be verified by employing Equation (\[eq:valid\_DR\]) with the required functions.

Some operations with $\preceq$ {#sec:operations}
==============================

Inverse Mixing {#subsec:InverseMix}
--------------

We present *inverse mixing* as a method for combining uncertainty given two distributions over two different populations.

\[def:inversemixing\] Define the inverse mixture $$\tilde{f}_1\; [+] \;\tilde{f}_2=
\left(\tilde{f}_1^{(-1)}(z)+\tilde{f}_2^{(-1)}(z) \right)^{(-1)},$$

and the $\alpha$-inverse mixture $$\tilde{f}_1\; [+]_{\alpha}\; \tilde{f}_2=\left(\tilde{f}_1\left(\frac{z}{1-\alpha}\right)^{(-1)} + \tilde{f}_2\left(\frac{z}{\alpha}\right)^{(-1)} \right)^{(-1)},$$ where $0 < \alpha < 1$ is the mixing parameter.

For the case in which $\alpha= 1/2$, we claim that, in the discrete case, the *inverse mixture* corresponds to a combination of all the probabilities in both populations scaled by a factor $\alpha=\frac{1}{2}$, i.e., $\frac{1}{2}p_i\cup\frac{1}{2}q_i, i=1, \dots, 5$ and sorting them in a decreasing order. In the following examples, we demonstrate inverse mixing for the discrete and continuous cases.

\[exampleworkplace\] Consider two distinct groups of people in a work place. Denote the probability of the $i$-th member of group one and the $i$-th member of group two obtaining a promotion to be $p_i$ and $q_i$, respectively. Let the probabilities $p_{1},...,p_{5}$ and $q_{1},...,q_{5}$ be those denoted in Table \[tab:Probabilities\_MF\], noting that $   p_1\geq p_2\geq\cdots\geq p_5$, $q_1\geq q_2\geq \cdots q_5$ and $\sum_{i}p_i=1$, $\sum_{i}q_i=1$,

  ----------------- ------- ------- ------- ------- -------
         $i$           1       2       3       4       5
   \[0.5ex\] $p_i$   0.577   0.192   0.128   0.064   0.038
        $q_i$        0.730   0.219   0.036   0.007   0.007
  ----------------- ------- ------- ------- ------- -------

  : Ordered probabilities for Example \[exampleworkplace\].[]{data-label="tab:Probabilities_MF"}

To perform inverse mixing, we take the inverse of each pmf, combine them and sorting them into ascending order (demonstrated in panel (a) of Figure \[fig:inversetilt\]). The inverse is taken to obtain a pmf (panel (b), noting the change of scale on the $y$ axis). The result of direct mixing, obtained by summing ordered probabilities of two populations and scaling by a factor $\alpha$, i.e., $\frac{1}{2}(p_i+q_i), i=1, \dots, 5$, is shown in panel (c). Both mixtures provide information about the joint population, but the inverse mixture also preserves information about the individual subpopulations.

![(a) the addition of two inverse pmfs, (b) inverse mixture distribution with $\alpha=\frac{1}{2}$, (c) direct mixture distribution with $\alpha=\frac{1}{2}$.[]{data-label="fig:inversetilt"}](combine_inv.png){width="100.00000%"}

We now consider inverse mixing for the continuous case. In the continuous case, we need to pay attention to the maximum values (modes) of the probability distributions. We demonstrate the importance of this condition in the following two examples.

Given univariate and bivariate exponential distributions with the following form of decreasing rearrangements, $$\tilde{f}_1(z)=\exp\{-z\}, \quad \tilde{f}_2(z)=\exp\{-(2z)^{1/2}\},$$ we observe that $0<\tilde{f}_1(z), \tilde{f}_2(z)\leq 1$ and obtain the functional inverses, $$\tilde{f}_1^{(-1)}(z) = -\log(z), \quad \tilde{f}_2^{(-1)}(z)=\frac{1}{2}(\log(z))^2, \quad z\in(0, 1].$$ The left and central panels in Figure \[fig:inverse\_mixing\_con1\] show $\tilde{f}_1(z)$, $\tilde{f}_2(z)$, $\tilde{f}_1^{(-1)}(z)$ and $\tilde{f}_2^{(-1)}(z)$. The maximum value of these two probability functions occurs at the same point, $\tilde{f}_1(0)=\tilde{f}_2(0)=1$, so there is no kink in the inverse mixing of these two distributions. The inverse mixture of the two distributions is then $$f^{(1)}(z) =  \left\{ \tilde{f}_1^{(-1)}\left(\frac{z}{1-\alpha}\right) +  \tilde{f}_2^{(-1)}\left(\frac{z}{\alpha}\right) \right\}^{(-1)} = \left\{-\log\Big(\frac{z}{1-\alpha}\Big) +\frac{1}{2} \Big(-\log\Big(\frac{z}{\alpha}\Big) \Big)^2 \right\}^{(-1)} ,$$ for $0 \leq \alpha \leq 1$. The direct averaging of $f_1(x)$ and $f_2(x)$ gives: $$f^{(2)}(z) =   \left\{ (1-\alpha)\tilde{f}_1^{(-1)}(z)+  \alpha \tilde{f}_2^{(-1)}(z)\right\} ^{(-1)}= \left\{(1-\alpha)(-\log(z))+\frac{\alpha}{2}(-\log(z))^2 \right\}^{(-1)}.$$

We specify $\alpha=1/2$ to obtain the following expression: $$f^{(1)}(z)=\Big\{-\log(2z)+\frac{1}{2}\big[\log(2z)\big]^2 \Big\}^{(-1)}.$$ Since the expression above is a quadratic in $\log(2z)$, we obtain the two solutions $f^{(1)}(z)=\frac{1}{2}\exp\{1+\sqrt{1+2z}\}$ and $f^{(1)}(z)=\frac{1}{2}\exp\{1-\sqrt{1+2z}\}$. Proceeding with the second solution, as the first solution does not integrate to 1, we obtain the mean, variance and Shannon entropy: $\frac{7}{2}, \frac{99}{4}$ and $\frac{3}{2}+\log(2)$, respectively.

We perform similar calculations for direct mixing with $\alpha=1/2$. The pdf has the form, $$f^{(2)}(z)=\exp\left\{1-\sqrt{1+4z}\right\},$$ and, in this case, the values of the mean, variance and Shannon entropy are $\frac{7}{4}, \frac{99}{16}$ and $\frac{3}{2}$. Based on the pdfs, we have the following relationship for inverse mixing and direct averaging when $\alpha=1/2$, $$\label{eq:inverse_direct}
f^{(2)}(z)=2f^{(1)}(2z),$$ illustrated in the right panel of Figure \[fig:inverse\_mixing\_con1\]. Here, $f^{(1)}(z)$ (red line) stretches the support of the distributions, and lowers the overall maximum, whereas $f^{(2)}(z)$ (blue line) preserves the maximum and shrinks the support, confirmed by Equation (\[eq:inverse\_direct\]).

![*Left panel:* plot of DR functions $\tilde{f}_1(z)$ (solid line) and $\tilde{f}_2(z)$ (dashed line). *Central panel:* plot of functional inverses of the DR functions, i.e. $\tilde{f}_1^{(-1)}(z)$ (solid line) and $\tilde{f}_2^{(-1)}(z)$ (dashed line). *Right panel:* plot of inverse mixing and direct averaging, $f^{(1)}(z)$ (red line) and $f^{(2)}(z)$ (blue line).[]{data-label="fig:inverse_mixing_con1"}](inverse_mixing_con1){width="100.00000%"}

\[example:6\] We consider exponential distributions with means 1 and 2 and note that they are already DRs: $$\tilde{f}_1(z)=\exp\{-z\}, \quad\tilde{f}_2(z)=\frac{1}{2}\exp\{-z/2\}.$$ The left panel in Figure \[fig:inverse\_mixing\_con2final\] shows $\tilde{f}_1(z)$ and $\tilde{f}_2(z)$, depicted by solid and dotted lines respectively. We note that $0<\tilde{f}_1(z)\leq 1$ and $0<\tilde{f}_2(z)\leq\frac{1}{2}$, which indicates that there is different support for the functional inverses, i.e. $$\begin{aligned}
&\tilde{f}_1^{(-1)}(z)=-\log(z), \quad z\in (0, 1],\\
&\tilde{f}_2^{(-1)}(z)=-2\log(2z), \quad z\in (0, 1/2].\end{aligned}$$ The inverse mixing of these two distributions with $\alpha=\frac{1}{2}$ is defined as $$f^{(1)}(z)=\tilde{f}_1\; [+]_{\frac{1}{2}}\; \tilde{f}_2=\left\{-\log(2z)-2\log(4z) \right\}^{(-1)}.$$ To avoid negative values of the expression inside the functional inverse, we propose the following modification: $$\label{eq:arg_inverse}
\tilde{f}_1^{(-1)}(2z)+\tilde{f}_2^{(-1)}(2z) =\max\{0, -\log(2z)\}+\max\{0, -2\log(4z)\}.$$

In the central plot of Figure \[fig:inverse\_mixing\_con2final\], the dotted line corresponds to the function in Equation (\[eq:arg\_inverse\]). We note that the introduced modification results in a kink in $\tilde{f}_1^{(-1)}(2z)+\tilde{f}_2^{(-1)}(2z)$ at $z=0.25$. To obtain the inverse mixture, we are required to take another functional inverse by swapping the abscissa and ordinate. Therefore we observe a kink in $f^{(1)}(z)$ at $z=\log(2)$ in the right panel of Figure \[fig:inverse\_mixing\_con2final\] (blue line). We obtain the final form of the inverse mixing: $$\begin{aligned}
f^{(1)}(z)&=\begin{cases}
\frac{1}{2}\exp\{-z\}, &\mbox{if } 0<z<\log(2) \\
\frac{1}{2}\exp\{ \frac{-2\log(2)-z}{3} \}, &\mbox{if } z \geq \log(2),
\end{cases}\end{aligned}$$ and we obtain the values of the mean, variance and Shannon entropy: $2.85$, $8.91$ and $1 + \frac{3\log(2)}{2}$.

Similarly to example 5, we consider the direct averaging of these distributions with $\alpha=\frac{1}{2}$, i.e. $$f^{(2)}(z)=\big\{\frac{1}{2}\tilde{f}_1^{(-1)}(z)+\frac{1}{2}\tilde{f}_2^{(-1)}(z) \big\}^{(-1)}=
\left\{-\frac{1}{2}\log(z)-\log(2z) \right\}^{(-1)}.$$ As with inverse mixing, to avoid negative values, we modify the argument of the functional inverse: $$\label{eq:arg_direct}
\frac{1}{2}\tilde{f}_1^{(-1)}(z)+\frac{1}{2}\tilde{f}_2^{(-1)}(z)=\max\left\{ 0,-\frac{1}{2}\log(z)\right\}+\max\left\{0, -\log(2z) \right\}.$$

The solid line in the central plot of Figure \[fig:inverse\_mixing\_con2final\] corresponds to the function in Equation (\[eq:arg\_direct\]). We observe a kink in the function at $z=\frac{1}{2}$. As a result, we obtain a kink in $f^{(2)}(z)$ at $z=-\frac{1}{2}\log{\frac{1}{2}}$ in the right panel of Figure \[fig:inverse\_mixing\_con2final\] (red line). The final form of the direct averaging is $$\begin{aligned}
f^{(2)}(z)
&=\begin{cases}
\exp\{-2z\}, &\mbox{if } 0<z<-\frac{1}{2}\log(\frac{1}{2}), \\
\exp\Big(\frac{-2z-2\log(2)}{3}\Big), &\mbox{if } z\geq -\frac{1}{2} \log(\frac{1}{2}),
\end{cases}\end{aligned}$$ where the values of the mean, variance and Shannon entropy are $1.42$, $2.23$ and $1+\frac{\log(2)}{2}$.

![*Left panel:* plot showing DR functions $\tilde{f}_1(z)$ and $\tilde{f}_2(z)$. *Central panel:* plot showing the function $\tilde{f}_{12}^{(-1)}(z)$ which is employed in direct averaging and inverse mixing. *Right panel:* plot showing pdfs obtained from inverse mixing and direct averaging, $f^{(1)}(z)$ and $f^{(2)}(z)$.[]{data-label="fig:inverse_mixing_con2final"}](inverse_mixing_con2final){width="100.00000%"}

From the representation of inverse mixing and direct averaging in Figure \[fig:inverse\_mixing\_con2final\], we can see that $f^{(1)}(z)$ stretches the support of the distribution, whilst $f^{(2)}(z)$ shrinks it. The maximum (mode) from the direct averaging is double the maximum of the inverse mixing.

This example shows how distributions can be approximated using DR from higher dimensional distributions. Consider the pdf $$f(x)=3(1-x)^2,\quad x\in[0, 1],$$ which has the functional inverse $f^{(-1)}(x)=1-\sqrt{ x/3}$. We need to perform an expansion in $t=\log(x)$, substituting $x=\exp\{t\}$ in the test function, defined as $$h(t)=f^{(-1)}(\exp\{t\})=1-\sqrt{\exp\{t\}/3}.$$ We provide the first three terms of the Taylor series, $$h(t)\approx \Big(1-\frac{\sqrt{3}}{3}\Big)-\frac{\sqrt{3}}{6}t-\frac{\sqrt{3}}{24}t^2,$$ and obtain the approximation $$\hat{f}^{(-1)}(x) \approx  \left(1-\frac{\sqrt{3}}{3}\right)-\frac{\sqrt{3}}{6}\log(x)-\frac{\sqrt{3}}{24}(\log(x))^2,$$ from which we obtain the functional inverse, $$\hat{f}(x)=\exp\left( -2+2\sqrt{2\sqrt{3}-1-2\sqrt{3}x}\right).$$ From the left panel in Figure \[fig:series\_plot\], we observe that the functional inverse (solid line) and its approximation (dotted line) are close to each other and, from the right panel, that $\hat{f}(x)$ (red line) is a good approximation to $f(x)$ (blue line).

![*Left panel:* plot of the functional inverse, $f^{(-1)}(x)$(solid line), and its approximation, $\hat{f}^{(-1)}(x)$ (dotted line). *Right panel:* plot of $f(x)$ (blue line) and its approximation $\hat{f}(x)$ (red line).[]{data-label="fig:series_plot"}](series_plot){width=".7\textwidth"}

Independence, conditional distributions
---------------------------------------

It is well known and axiomatic that Shannon information $S$ and entropy $(-S)$ are additive under independence: if $ X$ and $Y$ are independent random variables of the same dimension then, $$S(X+Y) = S(X) + S(Y).$$ It is a natural conjecture that entropy is a maximum in some sense when random variables are independent. This result holds for Shannon entropy.

\[lem:lem4\_2\] Within the class of bivariate random variables $(X_1,X_2)$ with given marginal distributions $X_1 \sim f_1(x) $ and $X_2 \sim f_2(x)$, the maximum Shannon entropy is uniquely achieved when $X_1$ and $X_2$ are independent.

Let $H( \cdot )$ be Shannon information. For random variables $X,Y$, we have the well-known expansion for the joint information: $$H(X,Y) = H(X) + \mbox{E}_X H(Y |X).$$ Also well known is the inequality, which follows from Jensen’s inequality, for the second term $$\mbox{E}_X H(Y |X) \geq H(Y),$$ with equality, and uniquely for Shannon information. The resulting additivity, $H(X,Y) = H(X) + H(Y)$, characterises Shannon information.

We can argue by Lemma \[lem:lem4\_2\] that, within the class with fixed marginals, the independence case cannot be uniformly dominated within the ordering $\preceq$.

This example shows that if we change the type of entropy, in this case to Tsallis, then the independent case may no longer achieve the minimum. Take $X_1,X_2 = 0,$ with probabilities $$p_{00} = \alpha \beta,\; p_{10} = (1-\alpha)\beta,\; p_{01} = \alpha (1-\beta),\; p_{11} = (1-\alpha)(1-\beta)$$ We can generate all distributions with the same margins with a perturbation $\epsilon$: $$p_{00} = \alpha \beta + \epsilon,\; p_{10} = (1-\alpha)\beta-\epsilon,\; p_{01} = \alpha (1-\beta)-\epsilon,\; p_{11} = (1-\alpha)(1-\beta)+\epsilon,$$ with the restriction that $|\epsilon| < \min(p_{00}, p_{10}, p_{01}, p_{11}).$ Taking the Tsallis entropy with $\gamma = 1$, we find the minimum when $$\epsilon = - \frac{(2\beta-1)(2\alpha-1)}{4},$$ which is zero if at least one of $\alpha$ or $\beta$ is $\frac{1}{2}$, which, interestingly, is a little less restricted than the case in which the distribution must be uniform i.e., all $p_{ij} = \frac{1}{4}$.

Note that independence and inverse mixing are closely related. If $X_1 =\{0,1\}$ is a binary random variable with $\mbox{prob}\{X_1=1\} = \alpha$ and $X_2$ has $J$ levels, the one-dimensional DR of $(X_1,X_2)$ is the inverse mixture with mixing parameter $\alpha$. Thinking in terms of a $2 \times J$ table, we combine the top row of probabilities, the distribution of $X_2$ multiplied by $\alpha$, with the bottom row in which they are multiplied by $1-\alpha$. In the continuous case, if the DR pdfs of $X_1$ and $X_2$ are $\tilde{f}_1(z)$ and $\tilde{f}_2(z)$, respectively then, for the joint distribution, we can think of $\tilde{f}_1$ weighting $\tilde{f}_2$ (or vice versa) in a similar way. This leads to the formula which we write informally: $$\tilde{f}_{12} (z) = \int \tilde{f}_1 \left( \frac{z}{\tilde{f}_2^{-1}(x)} \right) dx = \int \tilde{f}_2 \left( \frac{z}{\tilde{f}_1^{-1}(x)} \right) dx.$$ Iterating this operation, we can cover independent random variables in several dimensions and recapture results, such as those in Section \[sec:cont\_major\]. As with Example \[example:6\], care has to be taken in handling supports and limits of integration. The following results enable us to propagate the ordering $\preceq$ via conditional distributions.

Consider two pairs of joint distributed random variable $(X_1, X_3)$ and $(X_2,X_3)$ and suppose that the conditional distributions satisfy $X_1\; \vline \; X_3 \preceq X_2 \; \vline \; X_3$, for all values of the conditioning random variable, $X_3$. Then, for the joint distributions, $(X_1,X_3) \preceq (X_2,X_3)$.

This follows since, with simple notation, the joint distributions are $f(x_1,x_3) = f(x_1| x_3) f(x_3)$ and $f(x_2,x_3) = f(x_2| x_3) f(x_3)$. By assumption, we have $\tilde{F}(x_1|x_3) \leq \tilde{F}(x_2|x_3)$, for all $x_3$. An inverse mixing with respect to $f(x_3)$, as mentioned above, completes the proof.

This result can be understood easily in the discrete case using tables. It says that if there are two tables of probabilities with the same row margins then if every row of one table dominates the corresponding row of the other table then the whole table dominates the other whole table.

Volume-contractive mappings
---------------------------

We have seen that the area of the support is a key component of studying $\preceq$. For example, in the discrete case, if $n=0$ and $p= (p_1, p_2, p_3,0)$ are our probabilities with $p_1+p_2+p_3 = 1$, we can say we have support size 3. If we then split $p_3$ to form $q = (p_1, p_2, \frac{p_3}{2}, \frac{p_3}{2})$ then $q \preceq p$. In the continuous case, we refer to such an operation as dilation: locally, we have the same amount of density but stretch the support. This is a dilation in the continuous case via a special kind of transformation of the random variable, whose inverse we can call contractive. The result implies that the volume contractible mappings decrease uncertainty.

A volume differential invertible mapping $ h:  \mathbb R \rightarrow \mathbb R$, $y = h(x)$ will be called a volume-contractive mapping if the determinant of its Jacobian: $J = \vline \max{det} \left\{ \frac{\partial y_i}{\partial x_j}\right\} \vline$ satisfies $ 0 < J \leq 1$ for $x \in \mathbb R$.

If $h( \cdot) $ is a volume contractible inverse mapping $\mathbb R \rightarrow \mathbb R$, then, for any random variable $X \sim f_X(x)$, it holds that: $$X \preceq Y = h(X).$$

We give a proof for the one dimensional case and, in addition, assume $f_X(x)$ and $f_Y(y)$ are invertible. Using the slice condition, we want to show that $$\mbox{prob} \{ f_X(X) \geq c\} \geq  \mbox{prob} \{ f_Y(Y) \geq c\}.$$ Developing the left hand side, we see that $$\begin{aligned}
\{ f_X (X) \geq C \} & \Leftrightarrow & \{X \geq f_X^{-1}(c) \} \\
                              & \Leftrightarrow & \{Y \geq h(f_X^{-1}(c)) \} \\
                              & \Leftrightarrow & \{X \geq f_X^{-1}(c) \}.\end{aligned}$$ Computing the density of $Y$ as $$f_Y(y) = |J^{-1}| f_X(h^{-1} (y)),$$ gives $$f_Y^{-1}(c) =   h\left(f_X^{-1} ( J c)\right).$$ We thus need to establish whether $h(f_X^{-1} (c)) \geq h( f_X^{-1} ( J c)).$ We see the statement reduces to $c \geq J^{-1} c,$ which holds by assumption.

Contractive flows in sensitivity analysis
-----------------------------------------

A motivation for the previous subsection was to provide a method of analysis for systems within the general area of uncertainty quantification, one aspect of which is sensitivity analysis: the study of the propagation of variability through systems from input to output. We have seen above how a volume contractive mapping $Y=h(U)$ can decrease uncertainty. This subsection covers a closely related idea: how to show that, for two different inputs, the outputs may have more or less uncertainty. There are several areas of study in which it is hoped to decrease the variability of $Y$ via different types of intervention on the input $U$. Examples include the theory of antithetic variables in classical simulation, stochastic control, portfolio optimisation and robust design [@bates2002optimisation]. More generically, trying different random inputs $U$ is central to Monte Carlo simulation.

We represent the intervention on $u$ as a deterministic transformation and we seek a way to reduce the variability of $Y$ by shifting $U$ in some way. Consider the composition: $$y_0  = G(u_0), \quad 
u_1  =  h(u_0), \quad
y_1  =   G(u_1).$$ For some function intervention on $u$ given by $h(\cdot)$, we express schematically $$\begin{array}{ccc}
y_0 & \leftarrow & u_0  \\
 \rotatebox{-90}{$\dashrightarrow$} &       & \rotatebox{-90}{$\rightarrow$} \\
y_1 & \leftarrow & u_1  
\end{array}$$ Here, we are interested in what kind of interventions will result in volume contractive mappings from $y_0$ to $y_1$, induced by the intervention $h(u)$ in the dashed arrow in the diagram. When this holds, we can say that there is less uncertainty about the stochastic output $Y_1$ than about $Y_0$, for any input $U_0$.

A note of caution is that the induced function $y_0 \rightarrow y_1$ needs to be properly defined, in which case we can say that $h(u)$ is compatible. It is somewhat easier in explanation when the $u$-space and $y$-space have the same dimension. In addition, as we now see, local developments are easier. It is convenient to express $h(u)$ as a [*flow*]{} of the form: $$\begin{aligned}
h(u,\epsilon) = h(u_0) + \epsilon \xi(u_0) +\mbox{O}(\epsilon^2).\end{aligned}$$ In one dimension, $$\begin{aligned}
y_1 = G(u_0 + \epsilon \xi(u_0) ) + \mbox{O}(\epsilon^2),\end{aligned}$$ then $$\begin{aligned}
\frac{d y_1}{d y_0} & =  G'(u_0 + \epsilon \xi (u_0)) \frac{d u_0}{dy_0} + \mbox{O}(\epsilon^2),\\
& =  1 +  \epsilon\frac{\xi(u_0) G'' (u_0) + G'(u_0)\xi'(u_0)}{G'(u_0)} + \mbox{O}(\epsilon^2).\end{aligned}$$ Thus, if $\xi(u) = c  > 0 $, that is a constant, then $\frac{d y_1}{d y_0} < 1$ (locally) if and only if $$\begin{aligned}
\frac{G''(u_0)}{G'(u_0)} = \frac{d}{du_0} \log G'(u_0) < 1.\end{aligned}$$

For the multivariate linear case, we assume that the dimension of the $u$- and $y$- space are the same, namely $n$, for $n \times n$ matrices $A,B$ write: $$\begin{aligned}
y_0  = A u_0, \quad 
u_1  =  u_0 + \epsilon B u_0  +  \mbox{O}(\epsilon^2),\quad
y_1  =   Au_1,\end{aligned}$$ so that $$\begin{aligned}
y_1 = y_0 + \epsilon A B A^{-1}y_0 + \mbox{O}(\epsilon^2).\end{aligned}$$ Then locally we want $$\begin{aligned}
|\mbox{det} (I + \epsilon A B A^{-1})| < 1,\end{aligned}$$ for small $\epsilon$. If $\{\lambda_i\}$ are the eigenvalues of $ABA^{-1}$ then the condition is $$\begin{aligned}
| \prod_{i=1}^n (1+ \epsilon \lambda_i)| < 1.\end{aligned}$$

There is a particular problem when the input space and output space have different dimensions. Thus let the $y$-space above be one dimensional and the $u$-space $n$-dimensional. We can write for an $n$-vector $a$ $$y_0  = a^Tu_0, \quad 
u_1  =  u_0 + \epsilon B u_0  + \mbox{O}(\epsilon^2)\quad
y_1  =   a^Tu_1,$$ so that $$y_1 = (a^T + \epsilon a^T B)u_0 +  \mbox{O}(\epsilon^2).$$ Then, when we consider the random input version, and with the additional condition that the components of the random $U_0$ are independent, we have that $Y_0 \prec Y_1 $, when $B^T a \leq 0$ (componentwise). Geometrically, this says that the vector $a$ must lie in the dual cone generated by the columns of $B$.

Consider the problem when the input space and output space have different dimensions. We study the quadratic case, whose importance stems from the fact that, under suitable smoothness conditions, all function are locally quadratic. Let $A$ be an $n \times n$ symmetric positive definite matrix and consider the quadratic form $$y = u^TAu.$$ The function $y$ has a minimum at the origin. Starting at the point $u_0$, the natural shift which would contract the output is the direction of steepest descent which is given by the negative of the gradient at $u_0$, namely $$\frac{\partial y}{\partial u}\; \vline_{_{\;u=u_0}} =  2Au.$$ The flow in that direction would be expressed by $$u_1 = u_0 - \epsilon Au_0 +\mbox{O}(\epsilon^2),\quad y_0=u_0^TAu_0,$$ which gives $$\begin{aligned}
y_1 & = & u_1^T A u_1 + O(\epsilon^2),\\
       &  = & y_0  - 2\epsilon u_0^T A^2 u_0 + \mbox{O} (\epsilon^2), \\
       & =  & u_0^T (A - 2\epsilon A^2) u_0 + \mbox{O}(\epsilon^2).
 \end{aligned}$$ Taking the spectral decomposition of $A: A= \sum \lambda_i z_i z_i^T$, where the $\{\lambda_i\}$ are non-negative eigenvalues and the $\{z_i\}$ are standard unit length eigenvectors of $n\times 1$ dimension, we have $$y_1 = \sum_{i=1}^n (\lambda_i - 2\epsilon \lambda_i^2) (u_0^Tz_i)^2 + O(\epsilon^2).$$ If we take the input $U_0$ to be a vector of iid standard $N(0,1)$ random variables, then the variables $\{Z_i = (U_0^Tz_i)^{2}\}$ are independent $\chi^2$ random variables with 1 degree of freedom. Then it is straightforward to show that $$\sum_{i=1}^n \lambda_i Z_i \prec  \sum_{i=1}^n (\lambda_i - 2\epsilon \lambda_i^2)\;  Z_i,$$ so that approximately $Y_0 \prec Y_1$ where $Y_0 = U_0^T A U_0$, $Y_1 = U_1^T A U_1$ and $U_1 = U_0 - \epsilon AU_0 +\mbox{O}(\epsilon^2)$.

Algebra for uncertainty {#sec:algebra}
=======================

Let us recall some of the notation we have used. The majorisation ordering $\preceq$ is defined by: $$X_1 \preceq X_2  \Leftrightarrow \tilde{F}_1(z)  \leq \tilde{F}_2(z)\; \mbox{for all}\; x>0,$$ where $\tilde{F}$ means the (one dimensional) cdf for the decreasing rearrangement density $\tilde{f}$. We may write, equivalently, $\tilde{F}_1(z) \preceq \tilde{F}_2(z).$ Recall also that we can compare distributions in different dimensions via the construction given in Section \[sec:multivariate\].

We have the general concept of inverse mixing $ \tilde{h}_1\; [+] \; \tilde{h}_2$, as in Definition \[def:inversemixing\], applied to invertible functions $ h_1(x)$ and $h_2(x)$, whether or not they are densities. Denote by $\tilde{F}_1 \otimes \tilde{F}_2$ the associated cdf of the density, $$\begin{aligned}
 \label{eq:x}
 \frac{1}{2} ( \tilde{f}_1(z) \; [+] \; \tilde{f}_2(z)) =  \left(\tilde{h}_1(2z)^{(-1)} + \tilde{h}_2(2z)^{(-1)} \right)^{(-1)}.
 \end{aligned}$$ The first type of algebra is based on the following definition.

\[def:lattice\] For any two DR cdfs $\tilde{F}_1$ and $\tilde{F}_2 (z)$, we define $$\tilde{F}_1(z) \vee \tilde{F}_2 (z)  = \max(\tilde{F}_1(z), \tilde{F}_2(z)),$$ $$\tilde{F}_1(z) \wedge \tilde{F}_2 (z)  = \min(\tilde{F}_1(z), \tilde{F}_2(z)),$$ which themselves are cdfs. The partial ordering $\preceq$ under the meet and join $\vee$ and $\wedge$ defines a lattice which we refer to as the [*uncertainty lattice*]{}. It is satisfying that the ‘meet’ and ‘join’ which are defined once $\preceq$ is established can be manifested by the max and min of Definition \[def:lattice\].

We stress that, because we can embed a multidimensional distribution with density $f(x)$ into the one dimensional DR with pdf $\tilde{f}(z)$ and cdf $\tilde{F}(z)$, we can claim that the lattice is universal.

The inverse mixing cdf $\otimes$ can be combined with $\vee$ (or $\wedge$). To make the notation more appropriate, we will replace $\vee$ by $\oplus$ and then $\oplus$ and $\otimes$ yield a so-called max-plus (also called tropical) algebra [@maclagan2015introduction]. For this to be valid, we need distributivity.

\[dist\] We have $$\tilde{F}_3 \otimes ( \tilde{F}_1 \oplus \tilde{F}_2) =  (\tilde{F}_3 \otimes \tilde{F}_1) \oplus  (\tilde{F}_3 \otimes \tilde{F}_2),$$ where $\tilde{F}_1,\tilde{F}_2$ and $\tilde{F}_3$ are all cdfs for DR where $\tilde{F}_1(z) \oplus \tilde{F}_2 (z) = \max(\tilde{F}_1(z),\tilde{F}_2(z))$ and $\otimes$ is given in Equation (\[eq:x\]).

This follows from the useful expression of inverse mixing for cdfs, namely: $$\tilde{F}_1(z) \otimes  \tilde{F}_2(z)  = \left( \tilde{F}_1(2z)^{(-1)} +  \tilde{F}_2(2z)^{(-1)} \right)^{(-1)}.$$ Swapping min for max, we obtain $$\min\left(\tilde{F}_3(2z)^{(-1)} + \tilde{F}_1(2z)^{(-1)}, \tilde{F}_3(2z)^{(-1)} + \tilde{F}_2(2z)^{(-1)} \right) =$$ $$\tilde{F}_3(2z)^{(-1)} + \min \left(\tilde{F}_1(2z)^{(-1)}, \tilde{F}_2(2z)^{(-1)}\right),$$ which is the key step and common to maxplus (tropical) proofs.

We can now define a (semi) ring which we shall call the [*uncertainty ring*]{}. We first note that the $\otimes$ unit element will be $0$ and the $\oplus$ unit element will be $-\infty$. We immediately have the issue that, to define the ring, we need to remove the fact that $\tilde{f}$ is a density and $\tilde{F}$ is a non-negative function with $\tilde{F}(z) \rightarrow 1$ as $z \rightarrow \infty$. The polynomials which comprise the ring require powers and monomials.\
Consider the pdf arising from $\tilde{F}_1 \otimes \tilde{F}_1$. This is: $$\begin{aligned}
\tilde{f}(z) & = & \left( \tilde{f}(2z)^{(-1)} + \tilde{f}(2z)^{(-1)} \right)^{(-1)}, \\
                & = & \frac{1}{2} \tilde{f}\left(\frac{z}{2} \right), \end{aligned}$$ which is the pdf for the scaled random variable $Y= 2 X$ where $X \sim \tilde{f}(z)$. The pdf for $Y$ is $\tilde{F}\left(\frac{z}{2}\right)$. In general the $k$-th $\otimes$ power is $$\otimes^n \tilde{F} (z) = \tilde{F}\left( \frac{z}{n} \right).$$ The intuition of this expression is that increasing powers represent increasing dilation and form a decreasing chain with respect to our $\preceq$ ordering. A monomial with respect to $\otimes$ takes the form $$\prod_{i=1}^{m} \otimes^{\alpha_i} \tilde{F}_i(z).$$ Adjoining the base field $\mathbb R$, and appealing to Lemma \[dist\], we can define a ring of tropical polynomials [@glicksberg1959convolution].

The uncertainty ring is the toric semi ring of non-decreasing, twice-differentiable functions on $[0, \infty)$ on $\otimes$ and $\oplus$ and with $\oplus$ identity as $-\infty$.

To obtain proper decreasing densities we need to impose the additional condition that $\tilde{f}(z)$ is decreasing, $\tilde{F}(0) = 0$ and $\tilde{F}(z) \rightarrow 1$ as $z \rightarrow \infty$. Assuming that we have proper pdfs, we summarise the operations that we have:

1.  Scalar multiplication $\tilde{F} \rightarrow \beta\tilde{F}$, $ \beta \in \mathbb R$.

2.  Inverse mixing $\tilde{F}_1 \otimes \tilde{F}_2$.

3.  Maximum and minimum of $\tilde{F}_1, \tilde{F}_2$, denoted $\vee$ and $\wedge$, respectively. Noting that $\vee$ is written $\oplus$, when discussing the ring. We can also define a min-plus algebra and may use $\oplus$.

4.  Powers and monomials.

5.  Convolution $\tilde{F}_1 * \tilde{F}_2$. This refers to the DR pdf of the sum of independent random variables $X_1 \sim F_1$ and $X_2 \sim F_2$.

Further natural developments using ring concepts such as ideals are the subject of further work. In fact, convolutions themselves form a semi-group [@feller2008introduction], but we do not delve into the relationship between our ring and that semi-group. It is instructive to work over the binary field so that we do not have to use full scales from $\mathbb R$, but only $\{0,1\}$. This also has the advantage that, in every polynomial, we have proper pdfs and cdfs. An analogy is Boolean algebra. In Example \[sec6:example\_exp\] we illustrate these operations and the complexity obtained from a single distribution.

\[sec6:example\_exp\] Exponential with unit mean. Let $X_1 \sim \exp\{-x\}, $ with $x>0$. We have, $$\begin{aligned}
\tilde{F}_1 & = &   1-e^{-z}, \\
\tilde{F}_2  =  \tilde{F}_1 *  \tilde{F}_1 & = &   1 - (1+\sqrt{2z}) e^{-\sqrt{2z}}, \\
\tilde{F}_3  =  \otimes^2 \tilde{F}_1 & = & \frac{z}{2} \left( 1- e^{-z} \right), \\
\tilde{F}_4 =  \otimes^2 \tilde{F}_2 & = & \frac{1}{2} \left(  z- e^{-\sqrt{2z}}   \right)  \\
\tilde{F}_5 =  \tilde{F}_1 \otimes \tilde{F}_2 & = & \frac{z}{2} - \left( \frac{z}{2} + \frac{z^{\frac{3}{2}}}{\sqrt{2}} \right) e^{-\sqrt{2z}}.\end{aligned}$$ For the uncertainty lattice in Definition \[def:lattice\], we have $$\tilde{F}_4 \preceq \tilde{F}_5 \preceq \tilde{F}_2 \preceq \tilde{F}_1, \quad \text{and} \quad \tilde{F}_5 \preceq \tilde{F}_3 \preceq \tilde{F}_1,$$ and neither $\tilde{F}_2$ or $\tilde{F}_3$ dominate the other. Under $\vee$ and $\wedge$ we can include $  \tilde{F}_3 \vee \tilde{F}_2$ and $\tilde{F}_3 \wedge \tilde{F}_2,$ and show that $$\tilde{F}_4 \preceq \tilde{F}_3 \wedge \tilde{F}_2    \preceq  \tilde{F}_3 \vee \tilde{F}_2   \preceq \tilde{F}_1.$$

Algebra and entropies
---------------------

Given the equivalence of the majorisation in Section \[sec:cont\_major\], we study how the above structures affect the manipulation of uncertainty measured by a single metric, $$H(f(x)) = \int h(f(x)) dx,$$ for some convex function $h( \cdot)$. Since $$H(f(x)) = \int_{0}^{\infty} h(\tilde{f}(z)) dz,$$ we can, without loss of generality, consider $H$ as a functional of the DR with the advantage that the operations $\otimes,\oplus$ can be applied. The following, then, is a collection of operations and results, which we can claim as a toolbox for handling uncertainty, which can be applied to any uncertainty ($-H$), in this class. We omit the proofs.

1.  $\tilde{f}_1 \preceq \tilde{f_2} \Rightarrow H(\tilde{f}_1) \leq H(\tilde{f}_2).$

2.  $H\left(\tilde{f}_1 \otimes \tilde{f}_2 \right) = H\left(\tilde{\frac{1}{2}f}_1\right) + H\left(\frac{1}{2}\tilde{f}_2\right).$

3.  $H\left(\tilde{f}_1 \oplus \tilde{f}_2 \right) \geq \max \left\{H\left(\tilde{f}_1\right), H\left(\tilde{f}_2\right) \right\}.$

4.  $H\left(\tilde{f}_1 \otimes (\tilde{f}_2 + \tilde{f}_3) \right) \geq H\left(\frac{1}{2}\tilde{f}_1\right) + \max\left\{H\left(\frac{1}{2}\tilde{f}_2\right), H\left(\frac{1}{2}\tilde{f}_3\right)\right\}.$

Uncertainty toolbox: future scenarios
-------------------------------------

We present two practical situations in which the results in Section \[sec:algebra\] are employed to handle the uncertainty in practical situations, by combining them in two different ways.

Suppose there is initially to be $i=1,2$ races with a different number of horses $n_i$ ranked in order of their probability of winning (in the mind of a punter or bookmaker): $$p_{(i,1)} \geq \cdots \geq p_{(i,n_i)},$$ such that $\sum_{j=1}^{n_i}p_{i,j}=1, i=1,2$. If the two sets of horses should be combined into a single race, the issue then is how to combine the probabilities. The simplest approach is to divide each probability by two and combine all the probabilities, ranking them in the process, i.e, $\tilde{F}_1 \otimes \tilde{F}_2$, that is inverse mixing with $\alpha = \frac{1}{2}$. One can imagine some effect which may lead to having $\alpha \neq \frac{1}{2}$ such as the track being wet and not suited to one set of horses.

Consider the case in which we have a single race (say race 1) with two punters who rank the horses in the same order, but with different probabilities, i.e., $p_{(1,1)}\geq \ldots \geq p_{(1,n_1)}$ and $q_{(1,1)} \geq \ldots \geq \ldots q_{(1,n_1)}$, respectively. To define a joint betting strategy, a set of odds combining each of their own odds could take different approaches: an optimistic, more certain, approach with $\max(\tilde{F}_1(z),\tilde{F}_2(z))$ or a more pessimistic, uncertain, approach with $\min(\tilde{F}_1(z),\tilde{F}_2(z))$.

We note that the argument in the last paragraph is predicated on the two punters’ initial rank order being the same. If not, there is a danger that the same horse may appear twice in the min or max ordering. The min or max may then refer to a kind of hypothetical race. Nonetheless, we suggest that they are useful notionally. The same issue arises if one considers an average of the two actual probabilities, direct mixing, $\frac{1}{2}(p_{(1,1)} +q_{(1,1)} )\geq \ldots \geq\frac{1}{2}(p_{(1,n_1)}+q_{(1, n_1)}),$ which corresponds to taking $\frac{1}{2}(\tilde{F}_1(z) + \tilde{F}_2(z))$.

Now suppose there are two sets of horses of size $n_1$ and $n_2$ and two punters. We see then that the toolbox provides various ways to combine uncertainties within rows and columns. With obvious notation, we have

                punter 1            punter 2
  -------- ------------------- ------------------
   race 1   $\tilde{F}_{11}$    $\tilde{F}_{21}$
   race 2   $ \tilde{F}_{12}$   $\tilde{F}_{22}$

\[commiteeexample\] Consider a scenario with two committees, each with the same number of members. Each committee covers a different area of oversight, for example different technologies to solve a particular problem, say ‘electricity’ and ‘gas’. We assume that, within each committee, there is a common agreement about a range of possible future scenarios in the committee’s area. We define two future scenario types:

1.  an active situation: to choose between the two technologies on some grounds, in which the two committees’ assessments would not be combined.

2.  a passive situation: where only the probability of a particular technology being used can be assessed, which may be a consequence of unforeseen events.

In the latter case, the inverse mixing is a way to combine the uncertainties into an overall assessment. The actual future would be one gas-based alternative or one electricity alternative, the “horse” that won the combined race. The difficulties that may arise, as pointed out already, are that the private initial order of members of the [*same*]{} committee may not be the same as one familiar in subjects such as choice theory and rater assessment.

Note that independence between committees tends to increase uncertainty, and, in our development, the unidimensional DR is equivalent to inverse mixing. This shows, heuristically, that, when two committees act independently, there is more, rather than less, uncertainty.

Empirical decreasing rearrangements {#sec:empirical}
===================================

We extend the applications of DR to analyse a data set collected from an experiment in a large number of trials or from a product of computer simulations. We present two approaches for deriving the empirical DR and its associated cdf from a data set. In Section \[subsec:Discrete\_Climate\], we assess the uncertainties associated with climate projections, in which we transform continuous variables into discrete ones by grouping values, which is easy to perform in two dimensions. In Section \[subsec:Cont\_Heat\], algorithms are presented to obtain approximations for $\tilde{f}(z)$ and $\tilde{F}(z)$. This allows us to perform majorisation over data sets in higher dimensions, which we demonstrate analysing the risk associated with individual infrastructure decisions in an energy system.

Discrete Majorisation - climate projections {#subsec:Discrete_Climate}
-------------------------------------------

The UKCP18 climate projections [@UKCP18] consider four different scenarios of greenhouse gas concentrations in the atmosphere, called Representative Concentration Pathways (RCP). We consider variables over the 12km gridboxes that cover the UK. UKCP18 uses ensemble methods in which the model is run multiple times with slightly differing initial conditions and parameter values to account for observational and parametric uncertainty. We consider two variables:

1.  Increase in mean air temperature at a height of 1.5 metres,

2.  Percentage increase in precipitation,

where each variable is relative to the baseline period of 1981-2010. The projections correspond to mean daily values over the period from 2050 to 2079. The data, illustrated in Figure \[fig:climate\_scatter\], is discretised by dividing each variable into ranges and counting the number of ensemble members that fall into each category in the two dimensions. The temperature anomaly is divided into five categories, whilst the increase in precipitation is divided into four categories, therefore an ensemble member falls into one of 20 categories.

![Scatterplot showing how ensemble members are categorised according to the value each variable. Each point represents an ensemble member and each colour represents a different RCP.[]{data-label="fig:climate_scatter"}](climate_scatter.png)

We present the joint distribution for two contrasting scenarios for RCP2.6 and RCP8.5, of temperature anomaly and percentage change in precipitation is shown in Table \[table:clim\_dist\]. The ordered probabilities for each RCP are shown in Figure \[fig:climate\_example\_pdf\_cdf1\]. We use this to obtain empirical cdfs of DR together with the maximum and minimum of the cdfs depicted in the left and central panel plots in Figure \[fig:climate\_example\_pdf\_cdf2\]. We observe that the empirical cdf for RCP2.6 lies above that of RCP8.5. The uncertainty is therefore related to the level of assumed greenhouse gas concentration in the atmosphere. In particular, RCP2.6 carries the lowest level of uncertainty among the considered scenarios since its cdf corresponds to $\tilde{F}_1(z) \lor\tilde{F}_2(z)\lor\tilde{F}_3(z)\lor\tilde{F}_4(z)$, where the subscript indicates the scenario. In contrast, RCP8.5 carries the most uncertainty, since its cdf corresponds to $\tilde{F}_1(z)  \land\tilde{F}_2(z)\land\tilde{F}_3(z)\land\tilde{F}_4(z)$.

[|ll|&gt;cc&gt;cc&gt;cc&gt;cc&gt;cc|]{} & &\
& & & & & &\
& [**$<$0%**]{} & 0.01 & 0.00 & 0.18 & 0.02 & 0.17 & 0.09 & 0.01 & 0.10 & 0.00 & 0.04\
& [**0-5%**]{} & 0.01 & 0.00 & 0.22 & 0.03 & 0.28 & 0.16 & 0.04 & 0.20 & 0.00 & 0.08\
& [**5-10%**]{} & 0.00 & 0.00 & 0.02 & 0.01 & 0.05 & 0.07 & 0.01 & 0.11 & 0.00 & 0.06\
& [**$>$10%**]{} & 0.00 & 0.00 & 0.00 & 0.00 & 0.00 & 0.01 & 0.00 & 0.02 & 0.00 & 0.02\

![Ordered probabilities under each RCP.[]{data-label="fig:climate_example_pdf_cdf1"}](climate_example_pdf_cdf1.pdf)

![*Left panel*: Empirical cdf for the DR for each scenario. *Central panel*: the representation of $\max(\tilde{F}_1(z), \tilde{F}_2(z), \tilde{F}_3(z), \tilde{F}_4(z))$ and $\min(\tilde{F}_1(z), \tilde{F}_2(z), \tilde{F}_3(z), \tilde{F}_4(z))$. *Right panel*: cdf of inverse mixing with equal weights specified for each scenario.[]{data-label="fig:climate_example_pdf_cdf2"}](climate_example_pdf_cdf2.pdf){width="100.00000%"}

We can apply inverse mixing, described in Section \[subsec:InverseMix\], to combine the uncertainties in the four scenarios. For the discrete case, we specify an equal weighting for each pdf, and we illustrate cdf of the inverse mixing in the right panel in Figure \[fig:climate\_example\_pdf\_cdf2\].

Majorisation in higher dimensions {#subsec:Cont_Heat}
---------------------------------

We denote a data set $x_{ij}, i=1, \dots, m$ and $j=1, \dots, n$, where $m$ and $n$ correspond to the total number of data points and dimensions respectively. To obtain the DR, we require the density function values to construct the distribution function $m(y)$ (see Section \[sec:cont\_major\]). We assume that the observed data are a sample from a population with unknown pdf $f_X(x_1, \dots, x_n)$, from which we estimate the pdf $\hat{f}_X(x_1, \dots, x_n)$. In our examples, we employ kernel density estimation to obtain $\hat{f}_X(x_1, \dots, x_n)$ using the `ks` package in `R`, which automatically selects the bandwidth parameters [@ks2020]. Algorithm 1 describes the approach adopted in this paper to obtain empirical DRs $\tilde{f}_{\hat{f}}(z)$. We note from Definition \[drdefn\] that obtaining the DR is a two-stage process. At the first stage, we are required to obtain a measure (distribution) function $m(y)$. For instance, in 1D, a distribution function $m(y)$ returns the length of the intervals on the $x$-axis at which $f(x)\geq y$. At the second stage, we employ the measure function to derive the DR.

Based on data $x_{ij}\in R, i=1, \dots, m$ and $j=1, \dots, n$, fit a pdf $\hat{f}_X(x_1, \dots, x_n)$ using kernel density estimation Produce a uniform and/or space-filling set $S$ of size $N$ across the input space $R$, with “sites” $s\in S$ Plot the estimated measure function values, $m_{\hat{f}}(y)$ against $y$ Swap the abscissa with ordinate, so that $\tilde{f}_{\hat{f}}(z)$ and $z$ correspond to $y$ and $m_{\hat{f}}(y)$.

As part of this algorithm, we employ Monte Carlo integration [@Fok1989] to estimate the volume of domain $S_y$ to derive the measure function $m_{\hat{f}}(y)$. In particular, [@Fok1989] proposed specifying another domain $R$ (a hypercube or a hyperplane) of known volume $\text{Vol}(R)$, such that $
S_y\in R.
$ The ratio of two volumes, $p=\text{Vol}(S_y)/\text{Vol}(R)$, and the volume $\text{Vol}(S_y)$ are estimated by $$\hat{p}=N_y/N\quad \text{and}\quad \hat{\text{Vol}}(S_y)=\hat{p}\text{Vol}(R).$$ We employ Algorithm 1 on a bivariate data set. We start by generating a random sample of size $m=200$, i.e. $x_{ij}, i=1, \dots, 200$ and $j=1, 2$, from a standard bivariate normal distribution, since we have the closed form expression for DR of a standard 2-dimensional normally distributed random variable $X$: $$\label{eq:DR2D}
    \tilde{f}(z)=\frac{1}{2\pi}\exp\Big\{-\frac{1}{2}\frac{z}{\pi} \Big\}.$$ To perform the algorithm, we produce a uniform sample of points of size $N=2500$ across the domain $R=[-5, 5]\times [-5, 5]$ of $\text{Vol}(R)=10^2$.

Figure \[fig:EmpDRFinal1\] demonstrates the implementation of Algorithm 1 as well as comparing $\tilde{f}_{\hat{f}(z)}$ to $\tilde{f}(z)$. In the left panel, we depict the estimated values of the distribution function $m_{\hat{f}}(y)$ against $y$. Note that the smoothness of the estimated distribution function depends on $M$, i.e. we expect to obtain a smooth representation of the distribution function with a large value of $M$ (cutoffs in density). We observe that the empirical DR (red dashed line) is overlapping with the DR in equation (\[eq:DR2D\]) (blue solid line) in the right panel.

![*Left panel*: A plot of the estimated measure function $m_{\hat{f}}(y)$ against $y$. *Right panel*: A DR plot where the blue solid line and the red dashed line correspond to $\tilde{f}(z)$ and $\tilde{f}_{\hat{f}}(z)$.[]{data-label="fig:EmpDRFinal1"}](EmpDRFinal.pdf){width=".8\textwidth"}

We present Algorithm 2 for obtaining an empirical cdf of the DR, an approximation to $\tilde{F}(z)$, denoted as $\tilde{F}_{\hat{f}}(z)$.

Specify an equally spaced vector $\boldsymbol{z}^*=(z_1^*, z_2^*, \dots, z_l^*)$ Fit a linear interpolator (spline) through $\{z_i, \tilde{f}_{\hat{f}}(z_i)\}_{i=1}^M$ (these values were derived in Algorithm 1) to obtain values of $\tilde{f}_{\hat{f}}(z_i^*), i=1, \dots, l$ Plot $\tilde{F}_{\hat{f}}(z^*)$ against $z^*$.

The weighting of computed probabilities by $$\frac{1}{\sum_{k=1}^{n-1}P(z_k<z<z_{k+1})},$$ in Algorithm 2 comes from the assumption that $z$ is upper bounded and we can only compute probabilities at the values specified in $\boldsymbol{z}^*$. Therefore, we expect $\sum_{k=1}^{n-1}P(z_k<z<z_{k+1})=1$. However, we tend to observe this sum to be slightly less than one due to errors introduced by numerical integration. Our proposed weighting is similar to normalisation performed as part of the construction of a histogram.

We proceed to demonstrate the implementation of Algorithm 2 in Figure \[fig:CDFDR2D\] on the bivariate data set used previously. We have a closed form expression for $\tilde{F}(z)$ of a standard 2-dimensional normally distributed random variable $X$: $$\label{eq:CDF_DR2D}
    \tilde{F}(z)=1-\exp\Big\{-\frac{z}{2\pi} \Big\}.$$

![*Left panel*: binned probability representation of the empirical DR, $\tilde{f}_{\hat{f}}(z)$, obtained as part of Algorithm 2. *Right panel*: empirical cdf of the DR, $\tilde{F}_{\hat{f}}(z)$ obtained from Algorithm 2. The blue solid line and red dashed line correspond to $\tilde{F}(z)$ and $\tilde{F}_{\hat{f}}(z)$ respectively.[]{data-label="fig:CDFDR2D"}](CDFModifiedScheme2D.pdf){width=".75\textwidth"}

Based on the right panel in Figure \[fig:CDFDR2D\], we conclude that by employing Algorithm 2, the empirical cdf $\tilde{F}_{\hat{f}}(z)$ is an accurate representation of $\tilde{F}(z)$.

![Empirical cdf of DRs obtained for random samples from the standard normal distribution with $n=1, \dots, 4$ (dimension).[]{data-label="fig:EmpiricalCDF"}](EmpiricalCDF.png){width="50.00000%"}

We used Algorithms 1 and 2 to construct $\tilde{F}_{\hat{f}}(z)$ in Figure \[fig:EmpiricalCDF\] for a random sample from the multivariate standard normal with varying dimensions. We observe that the empirical cdf of the DR based on a sample from the univariate normal distribution majorises the remaining empirical cdfs. This is an expected result and confirms our expectation that uncertainty increases as we increase the number of dimensions. The accuracy of the representation of the DR is sensitive to the number of data points in the sample, $m$, while the smoothness of DR depends on $M$.

### District heating example {#subsubsec:DHE}

We compare the uncertainty associated with three potential design options for supplying heat in a model based on a system in Brunswick, Germany [@REUSEHEAT]. District heating networks allow heat from a centralised source to be distributed to buildings through a network of insulated pipes [@werner2013district]. This allows for a range of potential sources of heat to be connected and, in recent years, the idea of ‘reusing’ excess heat produced by nearby sources such as factories has gained traction as part of efforts to decarbonise the energy sector. This is done by using excess heat to heat water which is then distributed through the network. Traditionally, high temperature sources such as from heavy industry have been used, there has been increased interest in recent years in low temperature sources such as data centres, metro systems and sewage which require electric heat pumps to ‘upgrade’ the temperature before being suitable for use in the system.

The Brunswick case is a demonstrator on the EU funded ReUseHeat project [@REUSEHEAT] that aims to demonstrate the use of low temperature sources of heat for use in district heating networks. The city’s existing district heating network is powered by a Combined Heat and Power (CHP) plant, which uses natural gas as a fuel and outputs both heat for use in the network and electricity. The network in the newly constructed area of interest will be connected to the CHP and, in addition, there is an option to use excess heat from a data centre to provide at least some of the heat to the district. There are therefore three potential design options to be considered and these are shown in Table \[tab:Design\]. The question of interest is then which design option should be chosen, taking into consideration both cost and carbon emissions.

Although not used in actual decision making for the city or in the ReUseHeat project itself, a simple model was described in [@Volodina2020] which outputs both Net Present Cost (NPC) in € and $\text{CO}_2$-equivalent emissions (in metric tonnes). We use these simulations to demonstrate majorisation as a tool for decision making. Local and global sensitivity analysis was performed for each design option by varying a number of inputs to the model, resulting in a total of 81 simulations. In addition, three scenarios were considered and these are shown in Table \[tab:Scenarios\], further details of which are given in [@Volodina2020].

In Figure \[fig:ScatterPlot\], we illustrate the outputs of the model under sensitivity analysis by plotting emissions against NPC for each of the three design options and under each scenario. The highest level of emissions comes from design option 1 due to the use of natural gas whilst design option 3 is shown to be the most environmentally friendly. However, there is an inverse relationship between carbon emissions and cost with design option 3 being the most expensive.

![Net Present Costs against carbon emissions for all three design options under the three scenarios.[]{data-label="fig:ScatterPlot"}](ScatterPlot.png){width="70.00000%"}

We demonstrate majorisation and DR in the context of the model outputs under the three scenarios. In particular, we obtain $\tilde{f}_{\hat{f}}(z)$ and $\tilde{F}_{\hat{f}}(z)$ based on the distribution of points under the sensitivity analysis for each design option and under each scenario by employing Algorithms 1 and 2. To apply equal importance to both outputs, we scale the data on $[0, 1]$ and generate a uniform set $S$ across $[0, 1]\times [0, 1]$ of size $N=2500$. To produce a smooth representation of the DR and its cdf, we use a value of $M=5000$, which corresponds to the number of cutoffs in the density values.

  ----------------- -------------------------------
   **Design type**  **Description**
   Design option 1  Combined Heat and Power (CHP)
   Design option 2  CHP and Heat Pump
   Design option 3  Heat Pump
  ----------------- -------------------------------

  : Description of design options in the District Heating study [@Volodina2020].

\[tab:Design\]

  -------------- ---------------------- --------------------- ------------------------------------------
   **Scenario**  **Emission Penalty**   **Consumer demand**   **Commodity prices**
      Green      100€/metric tonne      -1% annual change     $\uparrow$ gas, $\downarrow$ electricity
     Neutral     40€/metric tonne       small fluctuations    small fluctuations
      Market     no penalty             +1% annual change     $\downarrow$ gas, $\uparrow$ electricity
  -------------- ---------------------- --------------------- ------------------------------------------

  : Description of scenarios in the District Heating study [@Volodina2020].

\[tab:Scenarios\]

![Empirical cdf for decreasing rearrangements $\tilde{F}_{\hat{f}}(z)$ for (i) all three design options plotted together for each individual scenario (*first row*), (ii) all three scenarios plotted together for each individual design option (*second row*).[]{data-label="fig:CDF_HeatExample"}](ScenarioPlotHeat.pdf "fig:"){width="100.00000%"} ![Empirical cdf for decreasing rearrangements $\tilde{F}_{\hat{f}}(z)$ for (i) all three design options plotted together for each individual scenario (*first row*), (ii) all three scenarios plotted together for each individual design option (*second row*).[]{data-label="fig:CDF_HeatExample"}](DesignPlotHeat.pdf "fig:"){width="100.00000%"}

Plots of the empirical cdfs of the DR $\tilde{F}_{\hat{f}}(z)$ are shown in Figure \[fig:CDF\_HeatExample\]. In the first row, the cdfs obtained for all three design options are considered together under individual scenarios. A feature here is that, under the green and neutral scenarios, the empirical cdf for design option 3 lies above that for design option 2, which lies above that for design option 1 whilst, under the market scenario, the ordering of the empirical cdfs changes and the empirical cdf for design option 1 lies above that for design option 2. We conclude that under all three scenarios, the (unknown) distribution function associated with design option 3 majorises the cdfs for both design options 1 and 2. We therefore consider that, for the outputs of interest, design option 3 is less uncertain than the alternatives.

The second row of plots shows empirical cdfs for each scenario plotted together under individual design options. This gives a slightly different view, allowing us to compare how the uncertainty under each design option varies under the different scenarios. For example, under design option 1, there is clear ordering of the cdfs which implies that the Market scenario is less uncertain than the Neutral scenario which, in turn, is less uncertain than the Green scenario.

We now demonstrate the uncertainty tools from Section \[sec:algebra\] in order to combine the uncertainty under different scenarios and produce orderings of design options. In particular, under each design option, we find the minimum of the empirical cdfs associated with individual scenarios to obtain an approximation to $\tilde{F}_1(z)\lor\tilde{F}_2(z)\lor\tilde{F}_3(z)$. This is shown in the left panel of Figure \[fig:maxminDR\] and can be considered to represent an ‘optimistic’ approach. We find that design option 3 majorises the other design options.

![*Left panel*: representation of $\max(\tilde{F}_1(z), \tilde{F}_2(z), \tilde{F}_3(z))$. *Right panel*: representation of $\min(\tilde{F}_1(z), \tilde{F}_2(z), \tilde{F}_3(z))$.[]{data-label="fig:maxminDR"}](maxminDR.pdf){width=".7\textwidth"}

We also produce an approximation to $\tilde{F}_1(z)\land\tilde{F}_2(z)\land\tilde{F}_3(z)$, which corresponds to a ‘pessimistic’ approach. The results are shown in the right panel of Figure \[fig:maxminDR\] in which we obtain the maximum of the empirical cdfs associated with individual scenarios. In this case, we observe a clear ordering between design options: design option 3 majorises design option 2, which majorises design option 1. Under both the pessimistic and optimistic outlooks, we therefore conclude that design option 3 is less uncertain than the two alternatives.

![cdfs from inverse mixing with different weightings on each scenario: *Left panel*: equal weights on each scenario. *Central panel*: $\alpha_{G}=0.7$, $\alpha_{N}=0.2$ and $\alpha_{M}=0.1$. *Right panel*: $\alpha_{G}=0.05$, $\alpha_{N}=0.05$ and $\alpha_{M}=0.9$[]{data-label="fig:HeatInverseMixing"}](HeatInverseMixing.pdf){width="100.00000%"}

We employ inverse mixing to combine the uncertainty associated with all three scenarios under individual design options. This is done by estimating probabilities from the empirical DR, by Algorithm 2 and ordering the probabilities. We consider the defined weights in inverse mixing as probabilities of the occurrence of each scenario. Let $\alpha_{G}, \alpha_{N}$ and $\alpha_{M}$ be the weights applied to the Green, Neutral and Market scenarios, respectively. We consider three different cases for the weights: (i) equal weights, (ii) $\alpha_G=0.7$, $\alpha_M=0.15$ and $\alpha_N=0.15$ and (iii) $\alpha_G=0.05$, $\alpha_M=0.9$ and $\alpha_N=0.05$. The results of the inverse mixing in the three cases are shown in the left, middle, and right panels of Figure \[fig:HeatInverseMixing\], respectively. In cases (i) and (ii), there are clear orderings in which the cdf of design option 3 lies above the cdf of design option 2 which lies above that of design option 1. In case (iii), however, there is no ordering between the empirical cdfs. In particular, the cdfs for design options 1 and 2 cross. However, the pdf associated with design option 3 majorises the pdfs for both design options 1 and 2 and we conclude that design option 3 is the least risky option in all three cases. It is important to note that, whilst the above results provide useful guidance for comparing uncertainty, the uncertainty is only one aspect of such decisions and one would want to take into account the actual costs and carbon emissions (rather than just their variability) in each case. However, here we have demonstrated majorisation to be an intuitive approach to comparing uncertainty and ultimately aiding informed decisions in such settings.

Concluding remarks {#sec:conclusion}
==================

The concept of uncertainty is the subject of much discussion, particularly at the technical interface between scientific modelling and statistics. It is our contention that uncertainty is close in spirit to entropy, but that restriction to a limited definition of entropy can be lifted. The fact that there are different types points to the existence of a wider framework which may then widen the scope of uncertainty. The idea presented is that a candidate for a wider framework is a stochastic ordering under for which most, if not all, types of entropy are order preserving. Here, we suggest majorisation, which only compares the rank order of probability mass, continuous or discrete. We have shown that any two distributions can be compared, and consider this to be a principal contribution of the paper. We demonstrated this approach to assess the uncertainty in applications of climate projections and energy systems.
