Calculating probability distribution that describes the mean and variance of a normal distribution given N samples

Question

Given N samples from a normal distribution with unknown mean and variance, how to estimate the probability distribution of the mean and variance (separately) of that normal distribution?

To clarify, the question is about calculating the probability distribution that describes the mean, and the probability distribution that describes the variance of the normal distribution.

Answer 1

For a random sample of size $n$ from $\mathsf{Norm}(\mu,\sigma),$ $\bar X = \frac 1n \sum_{i=1}^n X_i$ estimates $\mu.$ Also, $S^2 = \frac{1}{n-1}\sum_{i=1}^n (X_i - \bar X)^2$ estimates $\sigma^2.$ Both estimates are unbiased: $E(\bar X) = \mu,$ $E(S^2) = \sigma^2.$

The quantity $T = \frac{\bar X = \mu}{S/\sqrt{n}} \sim \mathsf{T}(\nu=n-1),$ Student's t distribution with $\nu = n-1$ degrees of freedom. Therefore, if numbers $-t^*$ and $t^*$ cut probability $0.025$ from the lower and upper tails, respectively, of this distribution, then a 95% confidence interval for $\mu$ is of the form $\bar X \pm t^*S/\sqrt{n}.$

Also, the quantity $\frac{(n-1)S^2}{\sigma^2} \sim \mathsf{Chisq}(\nu = n-1).$ Therefore, if numbers $L$ and $U$ cut probability $0.025$ from the lower and upper tails, respectively, of this distribution, then a 95% confidence interval for $\sigma^2$ is of the form $\left(\frac{(n-1)S^2}{U},\, \frac{(n-1)S^2}{L}\right).$ Also, a 95% confidence interval for $\sigma$ is found by taking square roots of the endpoint of the CI for $\sigma^2.$

To illustrate computation of these confidence intervals in R, we begin with a random sample of size $n = 100$ from $\mathsf{Norm}(\mu, \sigma),$ and use this random sample to find CIs for $\mu$ and $\sigma.$

n = 100;  mu = 50;  sg = 7
set.seed(2021)
x = rnorm(n, mu, sg)
summary(x);  length(x);  sd(x)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  34.21   43.10   49.06   48.78   53.42   64.84 
[1] 100        # sample size
[1] 7.200381   # sample standard deviation

Then a 95% CI for $\mu$ is $(47.35,\, 40.21),$ which happens to include the value of $\mu = 50$ in the simulation. Such a confidence interval will contain the population mean for about 95% of samples.

ci.mu = mean(x) + qt(c(.025,.975), n-1)*sd(x)/sqrt(n)
ci.mu
[1] 47.34990 50.20732

In R, the one-sample t test procedure also makes a confidence interval which may be printed using $-notation without showing the entire printout. This is the same CI shown just above.

t.test(x)$conf.int
[1] 47.34990 50.20732
attr(,"conf.level")
[1] 0.95

A 95% CI for $\sigma$ is $(6.32,\, 8.36).$ Again here our CI happens to contains the estimated population parameter $\sigma = 7.$

ci.sg = sqrt( (n-1)*var(x)/qchisq(c(.975,.025), n-1) ) 
ci.sg
[1] 6.321984 8.364505

Notes:

(1) The CI for $\mu$ is centered at the point estimate $\bar X.$ However, CIs for $\sigma^2$ and for $\sigma$ are not centered at $S^2$ and $S$ because the chi-squared distribution is not symmetrical.

(2) For normal data, $\bar X \sim\mathsf{Norm}(\mu,\sigma/\sqrt{n}).$ Then, if $\sigma$ happens to be known, a 95% CI for $\mu$ is of the form $\bar X \pm 1.96\frac{\sigma}{\sqrt{n}}.$

(3) If $\mu$ happens to be known, then $\sigma^2$ is estimated by $\widehat{\sigma^2} = \frac 1n\sum_{i=1}^n(X_i - \mu)^2$ and $\frac{n\,\widehat{\sigma^2}}{\sigma^2} \sim\mathsf{Chisq}(\nu = n).$ This leads to the 95% CI $\left(\frac{n\,\widehat{\sigma^2}}{U}, \frac{n\,\widehat{\sigma^2}}{L}\right),$ where $L$ and $U$ cut probability $0.025$ from the lower and upper tails, respectively, of the chi-squared distribution. Again here, a CI for $\sigma$ can be obtained by taking square roots of endpoints of the CI for $\sigma^2.$

(4) Even though $E(S^2) = \sigma^2,$ the sample standard deviation $S$ is not exactly an unbiased estimator for $\sigma,$ (Unbiasedness does not 'survive' the nonlinear square root operation.) However, for large $n$ the bias is small so $E(S) \approx \sigma.$ For sample sizes as small as $n = 15,$ the bias may be noticeable: Then $E(S) \approx 0.982\sigma.$ So if $\sigma = 7,$ then $E(S) \approx 6.87 \ne 7.$

set.seed(213)
s = replicate(10^6, sd(rnorm(15, 50, 7)))
mean(s)
[1] 6.874778

s = replicate(10^6, sd(rnorm(150, 50, 7)))
mean(s)
[1] 6.988528