Say I have draw N values from a normal distribution [$\mu_1$, $\sigma_1$]. Below are 10 sampled points compared to the normal distribution they're sampled from

Step 1: Sampling N points from a gaussian

I then create a normal mixture of normals [$\mu_2$, $\sigma_2$], where $\mu_2$ are the values sampled in step 1, and $\sigma_2$ is constant. Below are the 10 individual gaussians, and the normal mixture.

Step 2: Creating a function of N gaussians with means from step 1

Then, I calculate the 'percentile range', only I do so by integrating up to N percent of the total integral of the function. Below is the function with vertical lines defining the first 25% and last 25% of the integral.

Step 3: Calculating the positions where 25% and 75% of the integrated function are contained

Here's a psuedocode example:

def gaussian_sum(x, M, STD):
    norms = [norm.pdf(x, M[i], STD[i]) for i in range(len(M))]
    return sum(norms)/len(M)

testX = np.linspace(-4000, 4000, 4000)

means = np.random.normal(loc=0, scale = M*np.sqrt(940), size=100)
stds = [C]*len(means)

# Step 1
plt.plot(testX, Gauss(testX, 1, 0, M*np.sqrt(940)), "--", color='red')
for mean in means:
    plt.axvline(mean, ymin = 0, ymax = 0.2)

# Step 2
plt.plot(testX, [gaussian_sum(t, means, stds) for t in testX])
for i in range(len(means)):
    plt.plot(testX, norm.pdf(testX, means[i], stds[i]))
# Step 3
plt.plot(testX, [gaussian_sum(t, means, stds) for t in testX])
aftlow = scipy.optimize.bisect(lambda x: scipy.integrate.quad(gaussian_sum, -4000, x, args=(means, stds))[0] - 0.25, -4000, 4000)
afthigh = scipy.optimize.bisect(lambda x: scipy.integrate.quad(gaussian_sum, -4000, x, args=(means, stds))[0] - 0.75, -4000, 4000)

# Final result
print("Width:", afthigh-aftlow)

How would I derive the distribution of widths for a given value of N? I've been reading about order statistics and saw this excellent example here and I've seen examples using the percentile range, but this doesn't take into account the fact that this would be a continuous distribution.

  • 1
    $\begingroup$ Welcome to CV. Start by simplifying your problem. Based on your code, I believe you are just generating N iid values from a Normal distribution with parameters $\mu_1+\mu_2, \sqrt{\sigma_1^2+\sigma_2^2}$ (in an unnecessarily complicated way). Thus, you seek the distribution of the difference of two order statistics of a sample of N from that distribution. Apply the formula as indicated in your linked reference. $\endgroup$
    – whuber
    Feb 17, 2023 at 17:27
  • $\begingroup$ @whuber I don't think that's right, since the mean of the distribution should still be at 0. I'm specifically interested in how the functions generated don't resemble a normal distribution for low values of N. I've added some diagrams to help with this, in case it wasn't clear. $\endgroup$
    – Hunty2312
    Feb 18, 2023 at 20:37
  • 1
    $\begingroup$ It is very unclear. The pseudocode looks like you are generating Gaussian values conditional on Gaussian means. The graphics look more like mixture distributions, but you haven't shown how you create the graphics. $\endgroup$
    – whuber
    Feb 18, 2023 at 21:49
  • 2
    $\begingroup$ Again, this reads like an overly complicated method of doing something simpler -- but the text description is unclear. As far as I can tell, you are not summing the results: you are summing their density functions. As such, you appear to be creating a Normal mixture of Normals. The ambiguity lies in what it means to "create a function of a sum of ... Gaussians." $\endgroup$
    – whuber
    Feb 19, 2023 at 14:20
  • 1
    $\begingroup$ You're right, this is generating a normal mixture of normals - I just didn't know such a term existed. I've edited the question to reflect this $\endgroup$
    – Hunty2312
    Feb 19, 2023 at 14:44


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.