Distribution of sample sets drawn from a normal distribution

Question

I'm not very well versed in statistics, but I have a data analysis project that tests a certain type of data fitting algorithm with simulated data. I want to validate my initialization logic of the algorithm against theory, but I can't wrap my head around it. I'm drawing sets of samples from a normal distribution $f(x,\mu,\sigma)$ with, let's say, $n$ samples each, which are then binned to a histogram as my simulated data.

My problem boils down to this: I need to evaluate how the number of samples in a certain histogram bin fluctuates across the sample sets. I can calculate the mean number of samples in the histogram bin $[x_1,x_2]$ with $n\cdot(\Phi(x_2)-\Phi(x_1))$. But how do I get the distribution of the number of samples in the bin around this mean value?

My intuition says that the number of samples in a certain histogram bin in the sample sets are still normally distributed around the mean value calculated above. Still, I don't know how to calculate the standard deviation of the distribution.

Edit: I may have explained things poorly, so I try to elaborate in more detail with an example.

If I draw $n = 1000$ random samples form a normal distribution $f(x,\mu = 500,\sigma = 5)$ and bin the samples to a histogram with integer bins, I get a following histogram with 77 samples in the bin containing the $\mu$. The mean number of samples in the bin [499.5,500.5] calculated using CDFs is $1000\cdot(\Phi(500.5)-\Phi(499.5))\approx 79.7$.

And sure enough, if draw $m=10^5$ sets of $n=1000$ random samples, bin them to histograms, and note the number of samples in the "center bin", I get the following normalized distribution: And if I fit a gaussian function to the distribution, I get $\mu_{fit} = 79.1$ and $\sigma_{fit} = 8.52$ that kinda confirms the mean value calculated with the CDF. But how do I get, or is it possible to get, the $\sigma$ theoretically?

The idea is that if I can deduce the "center bin", I could use the number of samples in that bin ($c$) to limit my search space for $n$: Let us say, for all intents and purposes, the value of $c$ is at most $6\sigma_{fit}$ away from the mean $c$, and the probable range for $n$ can be calculated from that, which can be used as a search space for the fitting algorithm.

And yes, I know there are several valid ways to fit Gaussian functions to data. The key idea is to test out a different way of fitting.

Answer 1

To answer your question directly, I'll start by considering a specific example. Suppose $X \sim \mathsf{Norm}(\mu=100,\sigma=15)$ and the sample size is $n=100.$ Then it is reasonable that your histogram will have a bin $(80,90).$

The number $Y$ of observations in this bin is binomial: $Y\sim\mathsf{Binom}(n=100, p),$ where $p = P(80 < X < 90) = 0.1613,$ as computed in R below.

Also $E(Y) = np = 16.13,$ as you say, and $Var(Y) = np(1-p) = 13.53.$ For $n$ as large as $n = 100,$ you are right to say that $Y \stackrel{aprx}{\sim} \mathsf{Norm}(16.13,\, 3.678).$

p = diff(pnorm(c(80,90), 100, 15));  p
[1] 0.1612813
100*p;  100*p*(1-p), sqrt(100*p*(1-p))
[1] 16.12813
[1] 13.52697
[1] 3.677902

A simulation of a million such realizations of $Y$ gives the following result. Centers of small red circles show the PDF of $\mathsf{Binom}(n,p).$

set.seed(1124)
m = 10^6;  y = numeric(m)
for(i in 1:m) {
  x = rnorm(100, 100, 15)
  y[i] = sum((x>80) & (x<90))  }
mean(y);  sd(y)
[1] 16.12329
[1] 3.676105
mn = min(y)-.5; mx = max(y)
hist(y, prob=T, br = (0:40)+.5, col="skyblue2")
 points(0:100, dbinom(0:100, 100,p), col="red")

If, as suggested by @Dave, your ultimate goal is to check whether your simulated sample is normal, then there are some goodness-of-fit tests that you can use to test that. I will illustrate two of many.

In R, the Shapiro-Wilk test is implemented in the procedure shapiro.test. It tests whether your simulated data have a normal distribution.

set.seed(2020)
x = rnorm(100, 100, 15)
shapiro.test(x)

        Shapiro-Wilk normality test

data:  x
W = 0.98906, p-value = 0.5895

The P-value above 0.05 indicates that the simulated data are consistent with some normal distribution. By contrast, a Kolmogorov-Smirnov test, implemented in R as ks.test confirms that my data are specifically consistent with sampling from a population distributed $\mathsf{Norm}(\mu=100,\sigma=15):$

ks.test(x, pnorm, 100, 15)

        One-sample Kolmogorov-Smirnov test

data:  x
D = 0.088848, p-value = 0.4088
alternative hypothesis: two-sided

Note: If you are simulating your data via a well-vetted pseudorandom generator such as the Mersenne Twister using a method such as rnorm in R, then you can have high confidence that the sample is consistent with normal as simulated. If you are using some untested algorithm to simulate your data, then it is prudent to test several samples to see if you are getting samples as intended.