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

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.

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.

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), rng.choice(tuple_a_and_b)])
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.

• Nice code, but I'm puzzled: if you have say 100 samples from $A$ and 100 from $B$, then you have 200 samples of $A \cup B$ together. Can't you just sort the 200, and read off whatever quantiles you want ? Could you describe the real problem ? Mar 25 at 15:15
• I think AUB would be a 50:50 mixture not a convolution?
– Alex
Mar 25 at 17:33
• why one or the other -- don't you just have 200 data points ? Mar 25 at 18:49
• @Alex do you have a reference for the bootstrap method? I'm interested on what is the procedure if the sample sizes are not the same, or if there are more than two variables. Apr 5 at 23:33
• @J.R.C. I do not. I looked in a number of my stats books and could not find much. I think this is because the problem may be more probability than stats. In my two probability texts the most I was able to get was the difference between mixtures and convolutions, but not a lot of code or simulations on implementing arithmetic in both cases. I learned that mixtures are the easier of the two because you just decide what proportion you want from each distribution. Convolutions are something that seem a little more complicated so I just tested the bootstrap here.
– Alex
Apr 12 at 13:34