# PDF of sum of squared UniformDistribution (TransformedDistribution)

When drawing a number of samples from a NormalDistribution, the distribution of the sample variance is related to the ChiSquareDistribution.

I'm trying to find the pdf of the distribution for the sample variance when samples are drawn from a UniformDistribution.

In the simplest case

n=2;
dvar = TransformedDistribution[(Total[Table[x[i]^2, {i, 1, n}]] - Total[Table[x[i], {i, 1, n}]]^2/n)/(n - 1), Table[x[i] \[Distributed] UniformDistribution[{-a, a}], {i, 1, n}]]
PDF[dvar,x]


it works. However for, n=3 or higher, I cannot evaluate the PDF explicitly -- didn't count on it; no 'UniformSumSquareDistribution' function -- or numerically (e.g. PDF[dvar /. a->1, 0.1]).

I can calculate the CentralMoments fine. It doesn't help with fixed limits in the UniformDistribution, e.g. [{0,1}].

What else to try?

Edit: Background

In order to avoid misunderstandings, I'll just add a little background to my query.

The symmetric uniform distribution has

d = UniformDistribution[{-a, a}]; {Mean[d], Variance[d]}
(* {0, a^2/3} *)


But if you don't know the distribution, and can only draw a limited amount of samples, say 3, you can only get an estimate of the variance by calculating the sample variance. I.e.

n = 3; Variance[RandomVariate[d /. a -> 1, n]]
(* 0.0616734 *)


Repeat it, and you get something else. In any case, not the "expected" 1/3.

You can see how the distribution of the sample variances look, e.g. by

Histogram[Variance /@
Partition[RandomVariate[UniformDistribution[{-1, 1}], 3000000],3],
100, "PDF"] As it can be seen, it is going to be hard to find the "true" variance of the underlying (in principle unknown) distribution. The reliability of one of the calculated sample variances can be inferred from the properties of the distribution shown.

The basic properties can be found, e.g. (as above)

n=3;
dvar = TransformedDistribution[(Total[Table[x[i]^2, {i, 1, n}]] - Total[Table[x[i], {i, 1, n}]]^2/n)/(n - 1),
Table[x[i] \[Distributed] UniformDistribution[{-a, a}], {i, 1, n}]];
{Mean[dvar], Variance[dvar]}
(* {a^2/3, a^4/15} *)


So if you could do the experiment (draw three samples and calculate the sample variance) an infinite number of times, the mean would give the correct value, but if you do it only once, you will get a value which you must assign an uncertainty. E.g. the 95 % confidence interval - this is where you would use 2*Sqrt[Variance] if it was a normal distribution - but it requires that you can perform calculations on the distribution dvar. Optimally, the PDF or CDF or the ability to calculate Quantile.

If you want to increase the reliability of your estimate, you could take more samples (increase n above). That all good, you still find the Mean[dvar(n)] = a^2/3 and you can easily infer that Variance[dvar(n)] = 2 (2 n + 3) a^4/(45 n (n-1)). But what is e.g. the 95 % confidence interval?

Getting Mathematica to evaluate PDF[dvar, 0.1] or more relevant, CDF[dvar, 0.025] would be a first step.

All this is easy when the underlying distribution was the normal distribution; that's given by the properties of the ChiSquareDistribution.

• If you haven't already done so, you should become familiar with the Irwin-Hall distribution (sum of uniform random variables: en.wikipedia.org/wiki/Irwin%E2%80%93Hall_distribution) to see an approach for determining the distribution of the sample variance. – JimB Nov 22 '16 at 16:19
• Thanks, but isn't Irwin-Hall just a UniformSumDistribution[]? Nothing squared there, as far as I can see. I would expect to be more in the direction of BetaDistribution, since Uniform^2 ~ BetaDistribution. – HJensen Nov 22 '16 at 16:25
• Yes, it is just about the sum of uniform random variables. But how that distribution is derived would be useful for what you want to do with the sample variance. – JimB Nov 22 '16 at 16:56
• For those new to this topic (the sampling distribution of the sample variance, given a non-Normal parent), there is a nice video on youtube: youtube.com/watch?v=V4Rm4UQHij0 – wolfies Dec 3 '16 at 16:16

Given: $$X \sim \text{Uniform}(-a,a)$$ with pdf $$f(x)$$:

Let $$(X_1, \dots, X_n)$$ denote a random sample of size $$n$$ on parent $$X$$. The OP then defines the following statistic, namely the sample variance:

$$h_2 = \frac{1}{n-1}\sum _{i=1}^n \left(X_i-\bar{X}\right){}^2$$

where $$\bar{X}$$ denotes the sample mean. Clearly, $$h_2$$ is itself a random variable, and the OP is interested in the sampling distribution and variance of $$h_2$$. The statistic $$h_2$$ is also known as the $$2^\text{nd}$$ h-stastistic. As background, h-statistics have the elegant property that they are unbiased estimators of population central moments, so $$E[h_2] = \mu_2 = Var(X)$$.

The Variance of $$h_2$$

This brings us into the realm of finding moments of moments ... for which the natural language is to use power sum notation, where the $$r^\text{th}$$ power sum is defined as:

$$s_r=\sum _{i=1}^n X_i^r$$

We can express $$h_2$$ in terms of power sums by any number of ways, but perhaps the simplest is to use the HStatistic function in the mathStatica add-on to Mathematica (see Rose and Smith, Mathematical Statistics with Mathematica at Chapter 7 ... a free download of the chapter is available here).

The variance of $$h_2$$ is just the $$2^\text{nd}$$ CentralMoment of $$h_2$$, and so can be found with the mathStatica package function (expressed here in terms of central moments of $$X$$):

where $$\mu_2$$ and $$\mu_4$$ denote the $$2^\text{nd}$$ and $$4^\text{th}$$ central moments of the parent random variable $$X$$. Note that this result holds for ANY distribution whose moments exist ... not just for the OP's symmetric Uniform. In the case of the OP's symmetric Uniform, calculating $$\mu_2$$ and $$\mu_4$$ yields the exact variance of $$h_2$$ as:

... using here the Expect function from the mathStatica add-on.

The sampling distribution of $$h_2$$

Calculating the exact distribution of $$h_2$$ sounds like it will be a messy and complicated process, if indeed a closed form is attainable. But, if the sample size $$n$$ is not too small, an elegant and good approximation is likely be attained, whether by appealing to the Central Limit Theorem or otherwise, by:

$$h_2 \overset{a}{\sim } N\big( E[h_2], Var[h_2] \big)$$

where:

• $$E[h_2] = \mu_2 = \text{Var}[X]$$:
• $$\text{Var}(h_2)$$ was derived above

All done.

So how does it work?

The following diagram compares:

• our suggested Normal approximation (dashed red curve) to a ...
• Monte Carlo simulation of the true pdf of $$h_2$$ (blue squiggly curve)

... at 4 different sample sizes $$n = 10,30,50$$ and $$100$$. When $$n = 10$$, there is a small noticeable skewness offset, but for larger values of $$n$$, it appears to work very neatly. In all plots, parameter $$a= 1$$.

• Thank, @wolfies, it seems I need to check out the mathStatica package. I'll give you answer-credits for the hint. It's nice to see that the calculation of the variance of h2 gives the result I inferred. Can mathStatica also calculate things like confidence intervals? Concerning the sampling examples, the interesting part is n < 10, e.g. n=3, far from the CLT-normal distributions. It's quite relevant for, say, reliability of results from expensive testing like fire-testing of products. – HJensen Nov 26 '16 at 21:10
• Seems the closest I'll get to e.g. the n=3 distribution is something like PearsonDistribution[1, 1, 0.0411311, -0.681234, 0.82862, 0.00239931] – HJensen Nov 27 '16 at 9:17
dist[a_] = UniformDistribution[{-a, a}];


Since your UniformDistribution has zero Mean, the distribution of its Variance is

dVar1[a_] = TransformedDistribution[x^2,
x \[Distributed] dist[a]];

PDF[dVar1[a], x] The distribution of the variance can then also be written

dVar2[a_] = ProbabilityDistribution[
1/(2 a Sqrt[x]), {x, 0, a^2}];

PDF[dVar2[a], x] Assuming[{0 < x < a^2},
PDF[dVar1[a], x] == PDF[dVar2[a], x] //
Simplify]

(*  True  *)

Plot[
Evaluate[
Tooltip[PDF[dVar2[#], x], #] & /@
Range],
{x, 0, 10},
PlotRange -> {-0.05, 1},
Frame -> True,
FrameLabel -> (Style[#, 14, Bold] & /@
{"x", "PDF[dVar2[a], x]"}),
Axes -> False,
PlotLegends -> Automatic] #[dVar1[a]] & /@ {Mean, Variance}

(*  {a^2/3, (4*a^4)/45}  *)

#[dVar2[a]] & /@ {Mean, Variance}

(*  {a^2/3, (4*a^4)/45}  *)


Either form of the distribution can be used with RandomVariate

Show[
Histogram[RandomVariate[dVar1, 10000], Automatic, "PDF"],
Plot[PDF[dVar1, x], {x, 0, 1},
PlotRange -> {0, 5},
PlotStyle -> Thick]] Show[
Histogram[RandomVariate[dVar2, 10000], Automatic, "PDF"],
Plot[PDF[dVar2, x], {x, 0, 1},
PlotRange -> {0, 5},
PlotStyle -> Thick]] • But the question (I think) is about the distribution of the sample variance for sample sizes larger than 2. That statistic won't be the sum of the squares of independent uniform random variables because the square of the sample mean is also involved. But as the sample size gets larger, then the distribution of the sample variance should be very close to the sum of squares of independent uniform random variables (when the true mean is zero). – JimB Nov 22 '16 at 20:23
• Thanks for the effort, @Bob Hanlon, but as Jim mentions, I am interested in the PDF of the sample variance for n > 2. As I mention in my post, n = 2 can be calculated directly in Mathematica. I can also easily get a visual of the shape, e.g. from: Histogram[ Table[Variance /@ Partition[RandomVariate[UniformDistribution[{0, 1}], n 1000000],n], {n, 2, 8}], 100, "PDF"]. You are calculating something different. – HJensen Nov 23 '16 at 8:58
• Since your UniformDistribution has zero Mean, the distribution of its Variance is ... Population variance is not a random variable (it is a constant), and thus has no distribution. Sample variance is something quite different. The OP has not, in fact, defined which estimator of sample variance he is using, which would be a good starting point. – wolfies Nov 23 '16 at 12:59
• @wolfies The standard unbiased estimator V(X)=1/(n-1) Sum (Xi - E(X))^2. Same as implemented in Mathematica Variance[]. The aim is to look at the reliability of the calculated variance (or standard deviation) of a few samples drawn from - in this case - a uniform distribution. – HJensen Nov 23 '16 at 18:33