10
$\begingroup$

Forgive me if I've missed something rather obvious.

I'm a physicist with what is essentially a (histogram) distribution centered about a mean value that approximates to a Normal distribution. The important value to me is the standard deviation of this Gaussian random variable. How would I go about trying to find the error on the sample standard deviation? I have the feeling its something to do with the error on each bin in the original histogram.

$\endgroup$
  • $\begingroup$ A hint is provided at stats.stackexchange.com/questions/26924. In general, the sampling error of a variance can be computed in terms of the first four moments of the distribution and therefore the sampling error of the SD can at least be estimated from those moments. $\endgroup$ – whuber May 7 '12 at 21:33
11
$\begingroup$

It sounds like you're asking for a calculation of the standard deviation of the sample standard deviation. That is, you're asking for ${\rm SD}(s) = \sqrt{ {\rm var}(s) }$, where

$$ s = \sqrt{ \frac{1}{n-1} \sum_{i=1}^{n} (X_i - \overline{X}) }, $$

$X_1, ..., X_n \sim N(\mu, \sigma^2)$ and $\overline{X}$ is the sample mean.

First, we know from the basic properties of variance that

$$ {\rm var}(s) = E(s^2) - E(s)^2 $$

Since the sample variance is unbiased, we know $E(s^2) = \sigma^2$. In Why is sample standard deviation a biased estimator of $\sigma$?, $E(s)$ is calculated, from which we can infer

$$ E(s)^2 = \frac{2 \sigma^2 }{n-1} \cdot \left( \frac{ \Gamma(n/2) }{ \Gamma( \frac{n-1}{2} ) } \right)^2 $$

therefore

$$ {\rm SD}(s) = \sqrt{ E(s^2) - E(s)^2 } = \sigma \sqrt{ 1 - \frac{2}{n-1} \cdot \left( \frac{ \Gamma(n/2) }{ \Gamma( \frac{n-1}{2} ) } \right)^2 } $$

$\endgroup$
  • $\begingroup$ Good point. I got an estimate of the variance of s^2. Taking the square root gives an estimate of the standard deviation of s^2. But you answered the actual question which was to get the standard deviation of s. I would assume that for practical reasons you too woulde replace σ with s to get an estimate using the formula. $\endgroup$ – Michael Chernick May 8 '12 at 22:32
  • $\begingroup$ Yes, that's right, you can replace $\sigma$ with $s$ and this approximation performs well even for modest sample sizes - I did some testing with $n=20$. $\endgroup$ – Macro May 8 '12 at 23:24
10
$\begingroup$

The quantity $X=(n-1) s^2/\sigma^2$ has a chi-squared distribution with $n-1$ degrees of freedom when the samples are independent and distributed with the same normal distribution This quantity can be used to get confidence intervals for the variance of the normal and its standard deviation. If you have the raw values and not just the central value of the bins you can calculate $s^2$.

It is known that if $X$ has a chi-squared distribution with $n-1$ degrees of freedom its variance is $2(n-1)$. Knowing this and the fact the $\mathrm{Var}(cX) = c^2 \mathrm{Var}(X)$ we get that $s^2$ has a variance equal to $$\frac{2(n-1)\sigma^4}{(n-1)^2} =\frac{2\sigma^4}{n-1} \>.$$ Although $\sigma^4$ is unknown you can approximate it by $s^4$ and you have a rough idea of what the variance of $s^2$ is.

$\endgroup$
  • $\begingroup$ I was going to post this at the beggining, but the problem as I see it here is that $\sigma^2$ is unknown. Given that fact, I don't know if it is valid to approximate $s^4\approx \sigma^4$ if we don't even know the sample size. I recall that one can show that the fourth moment can have serious problems with outliers. $\endgroup$ – Néstor May 8 '12 at 0:40
  • $\begingroup$ $s^4$ is a consistent estimator of $\sigma^4$ (provided $\sigma^4$ exists), right @Nesp? I think this is usually what is meant when people said "approximate" or "rough idea". $\endgroup$ – Macro May 8 '12 at 0:44
  • 2
    $\begingroup$ Maybe is the lack of sleep, but, isn't that like circular reasoning? $\endgroup$ – Néstor May 8 '12 at 0:47
  • $\begingroup$ We assumed from the onset that the data came from a normal distribution so there is no outlier issue. I meant rough in the way Macro suggests. I agree that the sample size affects how close s^4 is to σ^4. But the worry about outliers is offbase Nesp. If you downvoted me for that I think it is very unfair. What I presented was the standard way of estimating the standard deviation for s^2 when data are NORMALLY DISTRIBUTED. $\endgroup$ – Michael Chernick May 8 '12 at 0:55
  • $\begingroup$ @Nesp, Michael has given a consistent estimator of the variance of the sample standard deviation from a normally distributed sample - for large samples it will do well - simulate it and find out. I'm not sure why you think this is circular reasoning. $\endgroup$ – Macro May 8 '12 at 1:01
7
$\begingroup$

There are several ways of quantifying the error of the standard deviation in the normal case. I am going to present the profile likelihood of $\sigma$ which can be used for approximating confidence intervals.

Let $x=(x_1,...,x_n)$ be a sample from a Normal$(\mu,\sigma)$. The corresponding likelihood function is given by

$${\mathcal L}(\mu,\sigma) \propto \dfrac{1}{\sigma^n}\exp\left(-\dfrac{1}{2\sigma^2}\sum_{j=1}^n(x_j-\mu)^2\right)$$

Then, the Maximum Likelihood Estimators are given by $(\hat\mu,\hat\sigma)=(\bar x,s)$, where $s=\sqrt{\dfrac{1}{n}\sum_{j=1}^n(x_j-\bar x)^2}$. Given that you are interested on quantifying the error on $\sigma$, you can then calculate the normalised profile likelihood of this parameter as follows.

$$R_p(\sigma)=\dfrac{\sup_{\mu}{\mathcal L}(\mu,\sigma)}{{\mathcal L}(\hat\mu,\hat\sigma)} = \left(\dfrac{\hat\sigma}{\sigma}\right)^n\exp\left[\dfrac{n}{2}\left(1-\left(\dfrac{\hat\sigma}{\sigma}\right)^2\right)\right]$$

Note that $R_p:{\mathbb R}_+\rightarrow (0,1]$. An interval of level $0.147$ has an approximate confidence of $0.95$. Next I attach an $R$ code that can be used for calculating these intervals. You can modify it accordingly in your context (or if you post the data I can include these changes).

data = rnorm(30)
n = length(data)
sg = sqrt(mean((data-mean(data))^2))
# Profile likelihood
rp = function(sigma) return( (sg/sigma)^n*exp(0.5*n*(1-(sg/sigma)^2))  )
vec = rvec = seq(0.5,1.5,0.01)
for(i in 1:length(rvec)) rvec[i] = rp(vec[i])
plot(vec,rvec,type="l")
rpc = function(sigma) return(rp(sigma)-0.147)
# Approximate 95% confidence interval
c(uniroot(rpc,c(0.7,0.8))$root,uniroot(rpc,c(1.1,1.3))$root)

An advantage of this sort of intervals is that they are invariant under transformations. In this case if you calculate an interval for $\sigma$, $I=(L,U)$, then the corresponding interval for $\sigma^2$ is simply $I^{\prime}=(L^2,U^2)$.

$\endgroup$
  • $\begingroup$ I think he really just wanted the standard deviation of s. $\endgroup$ – Michael Chernick May 8 '12 at 0:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.