PDF of sum of squared UniformDistribution (TransformedDistribution)

Question

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.

Answer 1

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

_{(source: tri.org.au)}

_{(source: tri.org.au)}

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).

_{(source: tri.org.au)}

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$):

_{(source: tri.org.au)}

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:

_{(source: tri.org.au)}

... 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]$:

_{(source: tri.org.au)}

• $\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)

_{(source: tri.org.au)}

... 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$.

Answer 2

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[4]],
 {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[1], 10000], Automatic, "PDF"],
 Plot[PDF[dVar1[1], x], {x, 0, 1},
  PlotRange -> {0, 5},
  PlotStyle -> Thick]]

Show[
 Histogram[RandomVariate[dVar2[1], 10000], Automatic, "PDF"],
 Plot[PDF[dVar2[1], x], {x, 0, 1},
  PlotRange -> {0, 5},
  PlotStyle -> Thick]]