0
$\begingroup$

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}]]]

enter image description here

$\endgroup$
1
  • $\begingroup$ "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]). $\endgroup$
    – JimB
    Commented Jan 10 at 4:27

1 Answer 1

3
$\begingroup$

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

$\endgroup$
4
  • $\begingroup$ Thanks! How did you get the PDF? $\endgroup$ Commented Jan 10 at 17:07
  • $\begingroup$ 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]. $\endgroup$
    – JimB
    Commented Jan 10 at 17:20
  • $\begingroup$ 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]] $\endgroup$ Commented Jan 10 at 17:23
  • $\begingroup$ That's even better. (I should have done it that way.) $\endgroup$
    – JimB
    Commented Jan 10 at 17:26

Your Answer

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

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