# Distribution of the average of a sample of standard deviations

I have a sample of standard deviations. How is its average distributed?

The sample is not big so normality is hardly possible.

• How was this data set constructed? Do you know if the data comes from any particular distribution? – Blue Marker Oct 12 '14 at 6:54
• The data are weight of coins from the trials of the pyx collected over several years. I do assume normality of initial distribution. I can find averages and standard deviations for each year. Taking the sequence of standard deviations as my new data I found its average and want to find a confidence interval for it. For the latter I need to find its distribution. – Ari Belenkiy Oct 12 '14 at 17:09

So you're looking for the distribution of the average of standard deviations of sets of Normal draws. Let's say that there are $M$ sets, and for each of these $M$ sets, there are $N_i$ draws from an i.i.d. $\mathcal{N}(\mu_i,\sigma_i^2)$ for $i=1,\dots,M$.

Take any draw $i \in \{1,\dots,M\}$. The sample standard deviation is $$s_i = \sqrt{\frac{1}{N_i} \sum_{j=1}^{N_i} (x_{ji}-\bar{x}_i)}$$ where $x_{ji}$ is the $j^\text{th}$ observation of set $i$ and $\bar{x}_i$ is the average of observations in set $i$. Then a result of Kenney and Keeping (1951) (referenced here http://mathworld.wolfram.com/StandardDeviationDistribution.html) is that $s_i$ will have the distribution $$f_{N_i}(s_i) = 2 \frac{\left(\frac{N_i}{2 \hat\sigma_i^2}\right)^{(N_i-1)/2}}{\Gamma\left(\frac{1}{2}(N_i-1)\right)} \exp \left( -N_i s_i^2 / (2 \hat\sigma_i^2) \right) s_i^{N_i-2}$$ where $\Gamma(\cdot)$ is the Gamma function, and $\hat\sigma_i^2 = N_is_i^2/(N_i-1)$. Let $V_i$ denote the random variable with the density given by the equation above. Then the distribution that you are looking for is the distribution of $$W = \frac{V_1 + \dots + V_{M}}{M} .$$ So how do you compute the density of $W$?

If $X$ and $Y$ are independent random variables with density functions $f_{X}(x)$ and $f_{Y}(y)$ respectively, then $Z =X + Y$ has density equal to its convolution, $$f_{Z=X+Y}(z) = (f_{X} * f_{Y})(z) = \int_{-\infty}^\infty f_{X}(z-y) f_{Y}(y) \,dy .$$

Apply this fact recursively to find the distribution of $W$, and after a ton of incredibly ugly math, you'll get the density of $W$.

• Something appears wrong with this density -it simplifies to the point of equaling $f(s) = A\cdot(1/s)$, with $A$ a constant. – Alecos Papadopoulos Oct 12 '14 at 22:08
• While yes, s must follow a kind of chi-2 distribution, the authors of the formula seems to miss adding <s> (mean) instead of s proper. Besides, they seem to draw a different distribution for different Ni - different sizes of initial samples. But the final answers seems to boil down to a convolution which is indeed an ugly stuff. Is there a good approximation for it? – Ari Belenkiy Oct 13 '14 at 2:30

You have a sample of standard deviations? i.e. more than one standard deviation sampled from the same population. Do you also have a sample of means? If yes, maybe plot the means on a graph and see what the mean distribution looks like.

Maybe if you mention what you are trying to find out from the data set, someone might be able to help you find a way around.

This seems like an ideal situation for bootstrapping.

Suppose you have $N$ standard deviations stored in some vector devs. Select a sample of size $N$ with replacement from devs, using the sample command.

First, I'm going to create some dummy data.

# Create dummy data drawn from a normal distribution
# with mean 0 and standard deviation 1

data <- rnorm(10000, 0, 1) # population


Next, generate 100 standard deviations, each based on 50 observations from the population we created above.

devs <- replicate(100,
sd(sample(data, 50, replace=TRUE)),
simplify='vector')


Vector devs is our vector of standard deviations. We're now going to get a sample of size $N$ from devs and take the mean. It is very important to sample with replacement: otherwise, you will simply get the same sample over and over.

means <- replicate(length(devs),
mean(sample(devs, length(devs), replace=TRUE)),
simplify='vector')


You can do all the stuff you normally do with a vector, including plotting the distribution and finding a confidence interval.

quantile(means, c(0.025, 0.975)) # 95% CI of sampling distribution

# Plot to see distribution
hist(means, xlim = c(0.9, 1.1), col='blue')


The output of the function is below. It will vary somewhat, since I didn't set a seed. Notice, however, how nicely it estimates our population standard deviation!