How to derive the combined frequency distribution from two independent normally distributed population?

Question

We have two different independent internal function calls for an application module, max of the time taken by two function calls individually determines how much the overall module takes. Individually, both the internal function calls are normally distributed. So, it means for a single execution, time taken by internal_function_call_1 can be at P90 datapoint while for internal_function_call_2 it could be P10 datapoint for their respective normal distributions. It means we cannot just take max of P10 of both calls and say what’s the P10 time taken by overall module. But I cannot take max of P10 and P90 of two functions to say what’s P10 or P90th time taken by P90 datapoint of overall module execution either. Need help to find what statistical methodology needs to be followed here to find the frequency distribution of overall module execution. Hope what I wrote makes sense, looking for help. Thanks!

Answer 1

If the normal distributions are identical you could use Ordered Statistics theory(see casella and berger section 5.4). But I think a better(correct?) approach here is from Multiple Random Variables(See 4.2.10). In your case assuming mean and the standard deviation of the functions is different and not known so the formulation is the following:

Let $X_3 = max(X_1, X_2)$. Here the $P(X_3 < x) = P(max(X_1, X_2) < x) = P(X_1 < x) * P(X_2 < x) = F_{X_1}(x)*F_{X_2}(x) $ which comes from the observation that the set $\{max(X_1, X_2) < x\} <=> \{X_1<x, X_2<x\}$ and so if you want to calculate $P(\{X_1<x, X_2<x\})$ it becomes the above form. You can computationally use this as follows with hardcoded values for the mean and variance of the two functions($\mu = 0, \sigma= 0.5$ and $\mu = 1, \sigma= 0.5$):

def f1_cdf(x):
  return stats.norm.cdf(x, 0, 0.5)
def f2_cdf(x):
  return stats.norm.cdf(x, 1, 0.5)
def analytic_nth_percentile(x, p=90):
  return f1_cdf(x)*f2_cdf(x) - p/100.0
Closed form percentiles using CDF function we formulated above
analytic_percentiles = [root_scalar(analytic_nth_percentile, args=p, bracket=[-10, 10], method='brentq').root for p in range(1, 100)]
Lets grab a large sample and compute the sample percentiles to compare
def random_sample(size=100):
  return np.stack((stats.norm.rvs(0, 0.5, size=size), stats.norm.rvs(1, 0.5, size=size)), axis=1)
s = random_sample(size=100000)
maxs = np.max(s, axis = 1)
sample_percentiles = [np.percentile(maxs, p) for p in range(1, 100)]
comparison_percentiles = np.stack((sample_percentiles, analytic_percentiles), axis=-1)
Lets compare the two
comparison_percentiles = np.stack((sample_percentiles, analytic_percentiles), axis=-1)
np.mean((comparison_percentiles[:, 0] - comparison_percentiles[:, 1])**2)
Out[197]: 8.115858961563016e-06

You can check this process is empirically correct by using KS test or a QQ plot which was like: