What is the correct way to average random variables and get correct quantiles

Question

Say I have two random variables $A$ and $B$ which may or may not be independent. I also have their $0.95$ quantiles $Q95_A$ and $Q95_B$. What is a valid way to average these densities and obtain valid quantiles?

Example Use Case

I am able to get draws from $A$ and $B$. In reality these random variables are probabilistic forecasting methods. I would like to obtain a combination forecast by averaging these densities.

More Clairity - Is this a question about mixture densities or about convolution operations?

From the wikipedia page on mixture distributions https://en.wikipedia.org/wiki/Mixture_distribution there is a distinction between when one would desire a mixture distribution vs a convolution operation.

Mixture densities are used when the the random variable's density is the sum of components

Convolution operators are used when the random variable's value is the sum of the underlying random variable values.

An example is given on the wikipedia page

As an example, the sum of two jointly normally distributed random variables, each with different means, will still have a normal distribution. On the other hand, a mixture density created as a mixture of two normal distributions with different means will have two peaks provided that the two means are far enough apart, showing that this distribution is radically different from a normal distribution.

The goal here is convolution not a mixture density.

Answer 1

The most flexible method I am aware of is to use the bootstrap approach. I have seen some recommend concatenation as a way to average densities see this link, but from my simulations I do not see that method working.

I will compare the concatenation (concat) approach and the bootstrap (boot_avg) approach here via simulation. My hope is that if there are other methods available they can also be run in a similar simulation.

Some math

To have a flexible method in hand which could work in general it is nice to start with an analytical solution and see which methods approach the analytical truth.

If random variable A is:

$$ f_A(x) = \int_{y=\infty}^{x} N(y\mid 10, 5) \, dy $$

and random variable B is:

$$ f_B(x) = \int_{y=\infty}^{x} N(y\mid 20, 5) \, dy $$

where

$$ N(x\mid\mu,\sigma) = \frac{e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2}}{\sigma \sqrt{2\pi}} $$

Then the analytical solution would be:

$$ f_Z(x) = \int_{y=\infty}^{x} N(y\mid 15,\frac{5}{\sqrt{2}}) \, dy $$

Conclusion, with $Z$ in hand, the true density of the average of two normals, we can now see which methods approach this density.

In the simulations below I show a bootstrap approach which seems to work well to average the densities and the quantiles generated from this method give valid quantiles which match the true $Z$.

Simulation

from numpy.lib.function_base import meshgrid
import numpy as np
from itertools import product, repeat
from numpy.random import triangular, normal
from plotnine import ggplot, aes, facet_wrap, geom_histogram, stat_ecdf, labs
from plotnine.themes import theme_538
from plotnine.scales import scale_fill_manual, scale_color_manual
import pandas as pd

rng = np.random.default_rng(42)

red = "#DD0031"
white = "#FFFFFF"
burgandy = "#AF272F"
dark_blue = "#004F71"
light_blue = "#3EB1C8"
blackish = "#202336"

n, u1, u2, sd = 500, 10, 20, 5
u_avg = np.mean([u1,u2])
a = rng.normal(u1, sd, size=n)
b = rng.normal(u2, sd, size=n)
concat = rng.choice(np.concatenate([a,b]), size = n)
z = rng.normal(u_avg, sd/np.sqrt(2), size=n)
def boot(tuple_a_and_b):
  return np.mean([rng.choice(tuple_a_and_b[0]), rng.choice(tuple_a_and_b[1])])
tuple_a_and_b = (a,b)
boot_replications = repeat(tuple_a_and_b, n) # n replications
boot_avg = np.array(list(map(boot, boot_replications)))

df_distributions = pd.DataFrame({'A': a, 'B': b, 'boot_avg': boot_avg, 'concat': concat, 'true': z}).unstack().reset_index().rename(columns={"level_0": "distribution", 0: "value"}).drop(columns = "level_1")
df_distributions

(ggplot(df_distributions)
 + geom_histogram(aes(x = 'value', fill = 'distribution'), color = 'white')
 +facet_wrap('~distribution',ncol = 1)
 +theme_538()
 +scale_fill_manual(values = [red, burgandy, light_blue, dark_blue, blackish])
  +labs(title = "The Bootstrap Average Approaches The True Average")
)

(ggplot(df_distributions)
 + stat_ecdf(aes(x = 'value', color = 'distribution'))
 +theme_538()
 +scale_color_manual(values = [red, burgandy, light_blue, dark_blue, blackish])
 +labs(title = "The Bootstrap Average Approaches The True Average")
)

Now lets see the $0.95$ quantile from the different distributions.

data = {
  'calculation':['quantile(a, 0.95)','quantile(b, 0.95)','quantile(concat, 0.95)','quantile(boot_avg, 0.95)','mean(q95(a), q95(b))','true'],
  'q95': [  
  np.round(np.quantile(a, .95), 2),
  np.round(np.quantile(b, .95), 2),
  np.round(np.quantile(concat, .95), 2),
  np.round(np.quantile(boot_avg, .95), 2),
  np.round(np.mean([np.quantile(a, .95),np.quantile(b, .95)]), 2),
  np.round(np.quantile(true, .95), 2)
  ]
}
pd.DataFrame(data).to_markdown(index=False)

calculation	q95
quantile(a, 0.95)	18.02
quantile(b, 0.95)	28.09
quantile(concat, 0.95)	26.68
quantile(boot_avg, 0.95)	20.9
mean(q95(a), q95(b))	23.06
true	20.58

The the plots and table show that the concatenation method does not approach the true average density while the bootstrap approach does.