# Computing statistics of sample average over $n$ IID variables

I need to plot $$E[\log|1-a X|]=0$$ as a function of $$a$$ where $$X$$ is sample average of $$n$$ IID draws from standard $$\chi$$^2. The following is a symbolic approach (borrowed from here) which works for $$n=1$$ and $$n=2$$ but breaks down for higher values.

What is the nicest way to get this plot for higher values of $$n$$?

The obstacle is that PDF of chi-squared sample mean doesn't evaluate for $$n>2$$, hence need an automatic way of constructing the appropriate integral, tips?

(* distribution of average of n chi-squared variables *)
sampleAverage[n_] := (
xvec = Array[x, n];
varspec =
Table[x[i] \[Distributed] ChiSquareDistribution[1], {i, 1, n}];
TransformedDistribution[Mean[xvec], varspec]
);
n = 2;
SF = StringForm;
pdf = PDF[sampleAverage[n], x];
expectation =
Assuming[a > 0,
Integrate[Log[Abs[1 - a*x]]*pdf, {x, 0, Infinity}]];
Plot[expectation, {a, 0, 8},
PlotLabel ->
SF["zero crossing for n= at a=", n,
a /. First@FindRoot[expectation, {a, 25/10}]]]


• "The obstacle is that PDF of chi-squared sample mean doesn't evaluate for n>2." Not true : t = Table[{n, FullSimplify[ PDF[TransformedDistribution[Sum[x[i], {i, n}]/n, Table[x[i] \[Distributed] ChiSquareDistribution[1], {i, n}]], z], Assumptions -> z > 0]}, {n, 1, 8}]. You'll find that the general pdf can be defined as pdf[n_, z_] := Exp[-z n/2] z^(n/2 - 1) n^(n/2)/(2^(n/2) Gamma[n/2]).
– JimB
Jan 10 at 4:27

## 1 Answer

One can find an expression for $$E[\log|1-a X|]$$ for $$n=2,3,4$$ using Integration but NIntegration appears necessary for $$n>4$$.

The pdf for the mean of $$n$$ $$\chi^2_1$$ independent random variables is given by

pdf[n_?IntegerQ, z_] := Exp[-z  n/2] z^(n/2 - 1) n^(n/2)/(2^(n/2) Gamma[n/2])


The means of $$\log|1-a X|$$ for $$n=2,3,4$$ are found as follows:

mean[n_?IntegerQ, a_] := FullSimplify[
Integrate[Log[1 - a  z] pdf[n, z], {z, 0, 1/a}, Assumptions -> a > 0] +
Integrate[Log[a  z - 1] pdf[n, z], {z, 1/a, ∞}, Assumptions -> a > 0], Assumptions -> a > 0]

m2 = mean[2, a]
(* -E^(-1/a) ExpIntegralEi[1/a] *)

m3 = mean[3, a] // Simplify
(* 2 - EulerGamma - (3 HypergeometricPFQ[{1, 1}, {1/2, 2}, -(3/(2 a))])/a + Log[a/6] *)

m4 = mean[4, a] // Simplify
(* 1 - ((2 + a) E^(-2/a) ExpIntegralEi[2/a])/a *)


To solve for values of $$a$$ that make the means equal to zero:

Solve[m2 == 0, a][[1]]
(* {a -> 2.68451} *)
FindRoot[m3 == 0, {{a, 2}}]
(* {a -> 2.69888} *)
FindRoot[m4 == 0, {{a, 2}}]
(* {a -> 2.663} *)


For larger values of $$n$$ numerical integration seems necessary:

n0 = 500;
FindRoot[(NIntegrate[Log[1 - a z] pdf[n0, z], {z, 0, 1/a}, WorkingPrecision -> 30] +
NIntegrate[Log[a z - 1] pdf[n0, z], {z, 1/a, ∞}] // Quiet) == 0, {{a, 2.6}}, WorkingPrecision -> 30] // Quiet
(* {a -> 2.00808} *)


It looks like $$a \rightarrow 2$$ as $$n \rightarrow \infty$$.

• Thanks! How did you get the PDF? Jan 10 at 17:07
• From the table created in my comment in your post, I saw that each pdf was a multiple of $z^{\frac{n}{2}-1} \exp \left(-z n/2\right)$. Then the coefficient to make the pdf integrate to 1 was found with 1/Integrate[Exp[-z n/2] z^(n/2 - 1), {z, 0, \[Infinity]}, Assumptions -> n > 0].
– JimB
Jan 10 at 17:20
• OK, ChatGPT was useful here, sum of chi-squared is another chi-squared, so alternative way is to do use PDF of TransformedDistribution[x/n, x \[Distributed] ChiSquareDistribution[n]] Jan 10 at 17:23
• That's even better. (I should have done it that way.)
– JimB
Jan 10 at 17:26