4
$\begingroup$

I sample $n$ points from standard normal and need the mean and variance of 3rd and 4th sample cumulants.

@JimB suggested that variance of 4th sample cumulant is given by the expression below. This seems like a very algebraically nice expression, is there way to derive/verify it with the help of Mathematica?

$$ \frac {24(n-1)((n-6)n+24)}{n^4} $$

Edit, Sep 5

The following gives a way to verify the formula for a fixed $n$, but not helpful for deriving the formula in terms of $n$

n = 5;
Clear[v];
rvec = Array[v, n];
sampleDist = ProductDistribution @@ Table[NormalDistribution[], n];
cumulantDist = 
  TransformedDistribution[Cumulant[rvec, 4], 
   rvec \[Distributed] sampleDist];
Variance[cumulantDist] == (24 (n - 1) ((n - 6) n + 24))/n^4
Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ Have a look at the mathstatica package, they also have a free book mathstatica.com/book/bookcontents.html $\endgroup$
    – yarchik
    Sep 6, 2020 at 21:03
  • $\begingroup$ @yarchik thanks for the pointer, it seems this issue is discussed in Section 7.3 "moments of moments". I actually bought 1.0 version of their package a while back, good to see they are still maintaining it $\endgroup$
    – Yaroslav Bulatov
    Sep 7, 2020 at 4:32

5 Answers 5

Reset to default
3
$\begingroup$

Here is a solution which is completely general, and works for any random variable whose moments exist. The standard Normal is then just a special and simple case.

The OP is interested in the $3^\text{rd}$ and $4^\text{th}$ sample cumulants (a type of sample moment, similar to sample central moments). The lingua franca for working with sample moments, and moments of sample moments, is to structure everything in terms of power sum notation, $s_r = \sum_{i=1}^n X_i^r$.

The mathStatica package for Mathematica has a function to express sample cumulants into power sum notation. Here is the $4^\text{th}$ sample cumulant, expressed in terms of power sums:

enter image description here

This statistic is a random variable, and the OP wishes to find some moments of it.

The mean of k4 is simply the $1^\text{st}$ Raw Moment of k4:

enter image description here

The variance of k4 is the $2^\text{nd}$ Central Moment of k4:

enter image description here

... where the solution is expressed in terms of the central moments $\mu_i$ of the population of the parent random variable $X$.

Note that this solution is completely general, and holds for any random variable whose moments exist.

The standard Normal case

In the simple standard Normal case $X \sim N(0,1)$, the first 8 central moments $\mu_r$ of the population are:

enter image description here

Substituting these values into the solutions above yields the mean of k4 as:

enter image description here

and the variance of k4 as:

enter image description here

All done.

Improve this answer
$\endgroup$
2
$\begingroup$

This is similar approach using FindSequenceFunction using the sample cumulant definition that Mathematica uses. (That sample cumulant definition is not always unbiased for the population cumulant.)

First define a function that finds the variance for the k-th cumulant for a sample size of $n$ from a Normal[0,1] distribution.

variance[k_, n_] := Module[{e1, e2, expectations},
  (* Finds variance of k-th cumulant for a sample of size n from a Normal[0,1] distribution *)
  (* Raw moments of a Normal[0,1] *)
  expectations = Table[x[i_]^i -> Expectation[z^i, z \[Distributed] NormalDistribution[0, 1]], 
    {i, 2 k, 1, -1}];
  (* Expectation of sample cumulant *)
  e1 = (Cumulant[Table[x[i], {i, n}], k] // Expand) /. expectations;
  (* Expecation of square of sample cumlant *)
  e2 = (Cumulant[Table[x[i], {i, n}], k]^2 // Expand) /. expectations;
  e2 - e1^2]

Now find the variance for the 3rd and 4th cumulants:

v3 = Table[{n, variance[3, n]}, {n, 2, 10}];
FindSequenceFunction[v3, n] // FullSimplify

$$\frac{6 (n-2) (n-1)}{n^3}$$

v4 = Table[{n, variance[4, n]}, {n, 2, 10}];
FindSequenceFunction[v4, n] // FullSimplify

$$\frac{24 (n-1) ((n-6) n+24)}{n^4}$$

Improve this answer
$\endgroup$
2
$\begingroup$

mathStatica (for me) has the more straightforward way to getting what you asked (as @wolfies has demonstrated). Here is how to accomplish that in a somewhat similar fashion in Mathematica.

(* Sample statistic *)
(k4 = MomentConvert[Cumulant[4], "SampleEstimator"] /. PowerSymmetricPolynomial[0] -> n) // TraditionalForm

Sample estimator for the 4th cumulant

(* Mean of sample statistic *)
(meank4 = MomentConvert[k4, "CentralMoment"] /. PowerSymmetricPolynomial[0] -> n // FullSimplify) // TraditionalForm

Mean of sample statistic

(* Variance of sample statistic *)
(vark4 = MomentConvert[k4^2 - meank4^2, "CentralMoment"] /. PowerSymmetricPolynomial[0] -> n // Simplify) // TraditionalForm

Variance of sample 4th cumulant

These can be simplified under the assumption of random samples from a Normal(0,1) distribution:

meank4 /. CentralMoment[i_] -> CentralMoment[NormalDistribution[0, 1], i] // FullSimplify

Mean under N(0,1)

vark4 /. CentralMoment[i_] -> CentralMoment[NormalDistribution[0, 1], i] // FullSimplify

Variance under N(0,1)

Improve this answer
$\endgroup$
1
$\begingroup$

An idea from JimB was to use FindSequenceFunction

getExpr[n_] := Module[{v},
  rvec = Array[v, n];
  sampleDist = ProductDistribution @@ Table[NormalDistribution[], n]; 
  cumulantDist = 
   TransformedDistribution[Cumulant[rvec, 4], 
    rvec \[Distributed] sampleDist]; Variance[cumulantDist]
  ];
values = Table[{n, getExpr[n]}, {n, 1, 9}];
FindSequenceFunction[values] (* (24 (-24 + 30 #1 - 7 #1^2 + #1^3))/#1^4 & *)
Improve this answer
$\endgroup$
1
$\begingroup$

An alternative way to use MomentConvert + MomentEvaluate:

$\kappa _3$

Mean

μκ3 = FullSimplify @ 
   MomentConvert[Cumulant[3], {"SampleEstimator", n}, {"CentralMoment", n}];

TraditionalForm @ μκ3

enter image description here

MomentEvaluate[μκ3, NormalDistribution[]]

0

Variance

σκ3 = FullSimplify[
   MomentConvert[Cumulant[3]^2, {"SampleEstimator", n}, {"CentralMoment", n}] - μκ3^2];

TraditionalForm @ σκ3

enter image description here

MomentEvaluate[σκ3, NormalDistribution[]]

enter image description here

$\kappa _4$

Mean

μκ4 = FullSimplify @ 
   MomentConvert[Cumulant[4], {"SampleEstimator", n}, {"CentralMoment", n}];

TraditionalForm @ μκ4

enter image description here

MomentEvaluate[μκ4, NormalDistribution[]]

enter image description here

Variance

σκ4 = FullSimplify[
   MomentConvert[Cumulant[4]^2, {"SampleEstimator", n}, {"CentralMoment", n}] - μκ4^2];

TraditionalForm @ σκ4

enter image description here

MomentEvaluate[σκ4, NormalDistribution[]]

enter image description here

Improve this answer
$\endgroup$

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.