6
$\begingroup$

I've been trying to find information on the sampling distribution of the standard deviation for uniform distributions and have been having a heck of a time figuring out the expected value for the standard deviation of a sample. Hopefully someone can help out!

Basically, I'm trying to find the expected standard deviation given sample size n. Given a uniform, continuous distribution with population range (b-a), I know that the population standard deviation would be (b-a)/sqrt(12); however, what would be the expected standard deviation for a sample size n drawn from that continuous uniform distribution?

My hope is that given a Uniform Distribution with range W (or "PopRange" in the example below), I would like to calculate an expected standard deviation for a given sample size. Since the mean sample standard deviation changes as a function of the sample size, I'm hoping to quantify it. Here's some MATLAB code illustrating what I am hoping to achieve.

% Number of samples for the brute-force estimation
totalSamples = 10000;

% The population range (b-a)
PopRange = 30;

for n = 2:100
    % Draw samples of sample size n
    theSamples = rand(n,totalSamples)*PopRange;

    % Calculate the mean of the standard deviations and ranges.  This
    % should be  equivalent to the means of the sampling distributions
    % of the standard deviation and range
    uniformSampleSD(n) = mean(std(theSamples,[],1));
    uniformSampleRange(n) = mean(max(theSamples,[],1)-min(theSamples,[],1));

    % Calculate the expected standard deviation, note that this is just the
    % population standard deviation at the moment and does not incorporate
    % the sample size
    uniformSampleESD(n) = PopRange/sqrt(12);
    % SEE: https://en.wikipedia.org/wiki/Prediction_interval#Non-parametric_methods
    % Calculate the expected sample range
    uniformSampleERange(n) = (n-1)/(n+1)*PopRange;
end

% Note the red curve underestimating the population standard deviation for
% smaller sample sizes? This is fine, but I would like to find a numerical
% estimate for it
figure
plot(uniformSampleSD,'r')
hold on
plot(uniformSampleESD,'b')

Any help or pointers would be greatly appreciated!


Edit:

Perhaps I mis-worded what I am looking for? I am looking for the mean of the sampling distribution of the standard deviation for a continuous distribution.

If I draw 100000 sets from a uniform distribution [0,1] of sample size n, the distribution of standard deviations will not have a mean of sqrt((1-0)^2 / 12) = 0.289. Empirically, the mean of the sampling distribution of the sd will be:

n = 2, sd_mean = 0.235
n = 3, sd_mean = 0.263
n = 4, sd_mean = 0.273
n = 5, sd_mean = 0.278
n = 6, sd_mean = 0.280
n = 7, sd_mean = 0.282
n = 8, sd_mean = 0.283
n = 9, sd_mean = 0.284
n = 10, sd_mean = 0.285

Note that it asymptotically approaches sqrt((b - a)^2 / 12), but that this is inaccurate for small n. See these examples of the sampling distributions for those small n:

https://i.stack.imgur.com/AwnhE.jpg (note: the "population" and "observed" sd line labels were accidentally swapped)

These were created using this code: totalSamples = 100000;

for n = 2:10
    theSamples = rand(n,totalSamples);
    setSDs = std(theSamples,[],1);
    popSD = sqrt((1-0)^2/12);
    nelements = hist(setSDs,50);

    figure
    hist(setSDs,50)
    hold on
    plot([popSD popSD],[0,max(nelements)],'r')
    plot([mean(setSDs) mean(setSDs)],[0,max(nelements)],'g')
    xlabel('Count')
    ylabel('Sample Standard Deviation')
    title({'Sampling Distribution of the Standard Deviation with 100000 sets';sprintf('drawn from a uniform distribution [0,1] with sample size %u',n)})
    legend('Observations',sprintf('Pop SD = %f',popSD),sprintf('Observed Mean Sample SD = %f',mean(setSDs)))
end
$\endgroup$
5
  • 4
    $\begingroup$ An exact answer will likely be extremely complicated: see stats.stackexchange.com/questions/41467 for an analysis of the distribution of just the sum of uniforms. Perhaps an asymptotic approximation would be ok? $\endgroup$
    – whuber
    Mar 11, 2014 at 22:42
  • 2
    $\begingroup$ Any insight into the estimated mean for the sampling distribution for the standard deviation would be greatly appreciated, even if it is an asymptotic approximation. Looking at my bootstrapped estimates, the average standard deviation across 10000 sets of samples sized n quickly converges on (b-a)/sqrt(12) (by around n = 20), but it would be great to have a better understanding/characterization of the behavior when n is smaller. $\endgroup$
    – Steve
    Mar 11, 2014 at 22:50
  • $\begingroup$ As it is, your problem is not fully specified. The sample standard deviation is defined as the square root of the sample variance. However, the sample variance has two common alternative specifications, and the problem fails to mention which one you are using. $\endgroup$
    – wolfies
    Mar 12, 2014 at 7:05
  • 2
    $\begingroup$ Also, if $S^2$ denotes the sample variance (whichever defn you are using), and $S$ denotes the sample standard deviation, do you specifically REQUIRE $E[S]$, or would you be happy to have $E[S^2]$ ... the latter is easy to derive. $\endgroup$
    – wolfies
    Mar 12, 2014 at 7:10
  • 2
    $\begingroup$ I added an edit that includes (hopefully) more clarification about my question. It seems as though that I am trying to find an unbiased estimate of the standard deviation for a continuous uniform distribution. See Unbiased estimation of standard deviation $\endgroup$
    – Steve
    Mar 12, 2014 at 11:40

1 Answer 1

5
$\begingroup$

The integration is difficult even with as few as $3$ values. Why not estimate the bias in the sample SD by using a surrogate measure of spread? One set of choices is afforded by differences in the order statistics.

Consider, for instance, Tukey's H-spread. For a data set of $n$ values, let $m = \lfloor\frac{n+1}{2}\rfloor$ and set $h = \frac{m+1}{2}$. Let $n$ be such that $h$ is integral; values $n = 4i+1$ will work. In these cases the H-spread is the difference between the $h^\text{th}$ highest value $y$ and $h^\text{th}$ lowest value $x$. (For large $n$ it will be very close to the interquartile range.) The beauty of using the H-spread is that, being based on order statistics, its distribution can be obtained analytically, because the joint PDF of the $j,k$ order statistics $(x,y)$ is proportional to

$$x^{j-1}(1-y)^{n-k}(y-x)^{k-j-1},\ 0\le x\le y\le 1.$$

From this we can obtain the expectation of $y-x$ as

$$s(n; j,k) = \mathbb{E}(y-x) = \frac{k-j}{n+1}.$$

Set $j=h$ and $k=n+1-h$ for the H-spread itself. When $n=4i+1$, $j=i+1$ and $k=3i+1$, whence $s(4i+1; i+1, 3i+1)=\frac{2i}{4i+1}.$

At this point, consider regressing simulated (or even calculated values) of the expected SD ($sd(n)$) against the H-spreads $s(4i+1,i+1,3i+1) = s(n).$ We might expect to find an asymptotic series for $sd(n)/s(n)$ in negative powers of $n$:

$$sd(n)/s(n) = \alpha_0 + \alpha_1 n^{-1} + \alpha_2 n^{-2} + \cdots.$$

By spending two minutes to simulate values of $sd(n)$ and regressing them against computed values of $s(n)$ in the range $5\le n\le 401$ (at which point the bias becomes very small), I find that $\alpha_0 \approx 0.5774$ (which estimates $2\sqrt{1/12}\approx 0.57735$), $\alpha_1\approx 1.091,$ and $\alpha_2 \approx 1.$ The fit is excellent. For instance, basing the regression on the cases $n\ge 9$ and extrapolating down to $n=5$ is a pretty severe test and this fit passes with flying colors. I expect it to give four significant figures of accuracy for all $n\ge 5$.

Figure

#
# Expected spread of the j and kth order statistics (k > j) in n
# iid uniform values.
#
sd.r <- function(n,j,k) (k-j)/(n+1)
#
# Expected sd of n iid uniform values.
#
sim <- function(n, effort=10^6) {
  x <- matrix(runif(n * ceiling(effort/n)), ncol=n)
  y <- apply(x, 1, sd)
  mean(y)
}
#
# Study the relationship between sd.r and sim.
#
i <- c(1:7, 9, 15, 30, 300)
system.time({
  d <- replicate(9, t(sapply(i, function(i) c(4*i+1, sim(4*i+1), i))))
})
#
# Plot the results.
#
data <- as.data.frame(matrix(aperm(d, c(2,1,3)), ncol=3, byrow=TRUE))
colnames(data) <- c("n", "y", "i")
data$x <- with(data, sd.r(4*i+1,i+1,3*i+1))

plot(subset(data, select=c(x,y)), col="Gray", cex=1.2,
     xlab="Expected H-spread", ylab="Expected SD (via simulation)")

fit <- lm(y ~ x + I(x/n) + I(x/n^2) - 1, data=subset(data, n > 5))
j <- seq(1, 1000, by=1/4)
x <- sd.r(4*j+1, j+1, 3*j+1)
y <- cbind(x,  x/(4*j+1), x/(4*j+1)^2) %*% coef(fit)
lines(x[-(1:4)], y[-(1:4)], col="#606060", lwd=2, lty=2)
lines(x[(1:5)], y[(1:5)], col="#b0b0b0", lwd=2, lty=3)

points(subset(data, select=c(x,y)), col=rainbow(length(i)), pch=19)
#
# Report the fit.
#
summary(fit)
par(mfrow=c(2,2))
plot(fit)
par(mfrow=c(1,1))
#
# The fit based on all the data.
#
summary(fit <- lm(y ~ x + I(x/n) + I(x/n^2) - 1, data=data))
# 
# An alternative fit (fixing alpha_0).
#
summary(fit <- lm((y - sqrt(1/12))/x ~ I(1/n) + I(1/n^2) + I(1/n^3) - 1, data=data))
$\endgroup$
1
  • $\begingroup$ Thank you for the response and the interesting analysis! This was very helpful. $\endgroup$
    – Steve
    Apr 11, 2014 at 7:44

Your Answer

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

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