Summation of uncertainty intervals

Question

I have the following data based on the 6 independent regions that make up the world. I want to compile them to find a global estimate. I can add best estimates from each region to give the best estimate for the the world. I am not sure how to deal with the uncertainty intervals as I feel adding the low and high estimates does not seem correct

Region	Best Estimate	Low	High
African Region	379	321	441
Region of the Americas	19	18	20
Eastern Mediterranean Region	80	68	92
European Region	21	20	22
South-East Asia Region	698	650	747
Western Pacific Region	87	80	93

Answer 1

I agree with Allan that having the data that makes up the estimates would be ideal. However, if you don't have the data I'm not sure you'd want to make the assumption of these all being Uniform distributions since this gives equal weight of all values between the lower/upper limits given. Since the column name is Best Estimate then I would think that you may not agree with that assumption.

Making a different assumption that these could each be Normally distributed and are independent of each other would allow you to take advantage of the sum of normally distributed random variables.

$$ \begin{aligned} X \sim \mathcal N(\mu_X, \sigma^2_X) \\ Y \sim \mathcal N(\mu_Y, \sigma^2_Y) \\ Z = X + Y \end{aligned} $$

Then $$Z \sim \mathcal N(\mu_X + \mu_Y, \sigma^2_X + \sigma^2_Y)$$

What you don't have is the variance for each but you could make another rough estimate using the range rule. For this you simply take the range for each and divide it by $4$ to have a rough estimate of your standard deviation.

Putting this all together in Julia you would have:

using Distributions

means = [379, 19, 80, 21, 698, 87]
lows = [321, 18, 68, 20, 650, 80]
highs = [441, 20, 92, 22, 747, 93]

# using the range rule
stds = (highs - lows) / 4
vars = stds.^2

# using the assumption of independent Normal R.V.s
world = Normal(sum(means), sqrt(sum(vars)))

# what's our 95% interval for the grand total?
quantile(world, [0.025, 0.5, 0.975])

Results in:

 1207.2072954499056
 1284.0
 1360.7927045500944

Answer 2

While you cannot blindly add any intervals, assuming that those are lower and higher bounds, it's easy. If you only care about adding the intervals and learning what would be the interval for the sum, this can be solved with interval arithmetic. The sum of the intervals is

$$ [x_1, x_2] + [y_1, y_2] = [x_1+y_1, x_2+y_2] $$

so you only need to sum lower bounds together and upper bounds together for the lower and upper bound of the sum. This gives you a conservative estimate (worst case, best case).

For data represented as the lower end, higher end, and the mode, there is a very simple probability distribution to represent it, the triangular distribution. It is popular because of its simplicity. You can use the distribution to conduct a Monte Carlo simulation to answer the question about the possible range of the values given the information provided. For the simulation, we will sample from the triangular distributions parametrized by lower bound, higher bound, and the mode and sum the results, to learn the distribution of the sums.

Let me also use Julia like @joshualeond.

using Random
using Distributions
Random.seed!(42)

dist = TriangularDist.(lows, highs, means)
world = [sum(rand.(dist)) for i in 1:500_000]

quantile(world, [0.025, 0.5, 0.975])
## 3-element Vector{Float64}:
##  1223.6470001955206
##  1285.2884768397992
##  1347.739670270986

As you can see, the average ends up very close to the result obtained by @joshualeond (+1) using the normal distribution, but the spread of the results is lower because the normal distribution is unbounded, while with triangular distribution we define the hard bounds that limit the variability.

The triangular distribution is hardly the most realistic choice, but it can be useful as a very simple model for simulation given the three-point summaries like the bounds and mode.

Answer 3

Ideally you'd have the raw data concatenate everything and with that you can make the uncertainty intervals. But if you only have this data, you could make assumption about the shape of distributions and simulate the result. For example assuming everthing is uniform:

import numpy as np
import pandas as pd

sim = []
for n in range(int(1e4)):
    africa = np.random.uniform(321, 421)
    americas = np.random.uniform(18, 20)
    medi = np.random.uniform(68, 92)
    europe = np.random.uniform(20, 22)
    sasia = np.random.uniform(650, 747)
    wpacific = np.random.uniform(80, 93)
    estimate = np.mean([africa, americas, medi, europe, sasia,wpacific])
    print(estimate)
    sim.append(estimate)

And with that is possible to estimate the intervals, and you could make this better weighting by the populations.