Application Center - Maplesoft
  Submit your work
Application Center Applications Joint Cumulants of Polykays Preview

App Preview:

Joint Cumulants of Polykays

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application



 

Joint Cumulants of polykays

E. Di Nardo*
elvira.dinardo@unito.it
https://www.elviradinardo.it
Tel: +390116702862  Fax: +390116702878

 

G. Guarino**
giuseppe.guarino@rete.basilicata.it

  
* Mathematics Department “G. Peano”, University of Turin, Turin, Italy

**Local Health Authority of Potenza, Italy

Introduction

Abstract: Given a multivariate random sample, this set of functions returns the joint cumulants of multivariate polykays in terms of joint cumulants of the multivariate population underlying the sample.

Multivariate polykays are unbiased estimators with minimum variance of joint cumulant products. These estimators are usually indexed by a list of integer vectors: each vector corresponds to a joint cumulant.
When this list reduces to a list of single integers, the multivariate polykay  reduces to a univariate polykay. When this list reduces to a single integer vector, the multivariate polykay reduces to a multivariate k-statistic. When this list contains a single integer vector with only one element, the multivariate k-statistic reduces to a univariate k-statistic.

Application Areas/Subject: Combinatorics & algebraic methods in statistics.

Keyword: umbral calculus, symmetric polynomial, set partition, multiset, cumulant, k-statistic, polykay.

See Also:  Maple algorithm [1] and [2]

Remark: k-statistics, polykays and their multivariate generalization are commonly given in terms of power sums

 

S[j, k, () .. ()] = Sum(X[i]^j*Y[i]^k, i = i .. n) .. ()

with n the sample size.

Initialization

 

> 

restart

> 

with(combinat, partition, multinomial, stirling2)

[partition, multinomial, stirling2]

(2.1)

Background Functions

The polyk function handles four different algorithms. Depending on the input parameters, univariate or multivariate polykays as well as univariate or multivariate k-statistics are expressed in terms of power sums. Specifically, if the input is a list of integer vectors, univariate polykays are returned; if the input is a list of integer vectors, each containing only one integer, polykays are returned; if the input is a single vector of integers, multivariate k-statistics are returned; if the input is a single integer, univariate k-statistics are returned; if the input is a list with a one dimensional integer vector, univariate k-statistics are returned.  All these functions are in the MAPLE worksheet [1] and have been included here for user convenience. For the multivariate case, the algorithm listing all subdivisions of a multiset has also been included. This algorithm is given and discussed in the MAPLE worksheet [2]. For details on the procedures and mathematical backgroud see [3 - 6].

K-Statistics

The n-th k-statistic k[n] is the unique symmetric unbiased estimator of the n-th cumulant kappa[n] of a given statistical distribution, that is

                                                                  E[k[n]] =  kappa[n]

where E denotes the expectation. They are expressed in terms of power sums involving the random variables of a random sample.

 

> 

makeTab_ss := proc (N) [seq([factorial(N)*mul(k[i], i = y)/mul(factorial(x)^numboccur(y, x)*factorial(numboccur(y, x)), x = {op(y)}), mul(S[i], i = y)], y = partition(N))] end proc:

makeK_s := proc (N) [seq(k[i] = add(stirling2(i, j)*x^j*(-1)^(j-1)*factorial(j-1), j = 1 .. i), i = 1 .. N)] end proc:

fd := proc (j, k) expand(mul(n-i-k, i = 0 .. j-k)) end proc:

ks := proc (N) local u, v; v := expand(eval(makeTab_ss(N), makeK_s(N))); u := [seq(x^i = (-1)^(i-1)*factorial(i-1)*fd(N-1, i), i = 1 .. N)]; add(x[1]*x[2], x = expand(eval(v, u)))/mul(n-x, x = 0 .. N-1) end proc:

Exampe: The k-statistic of order 3 is

> 

ks(3)

(n^2*S[3]-3*n*S[1]*S[2]+2*S[1]^3)/(n*(n-1)*(n-2))

(3.1.1)

Polykays

The symmetric statistic k[r, s, `..`] is an unbiased estimator of the cumulant product kappa[r]*kappa[s]..., that is

E[k[r, s, `..`]]  = kappa[r]*kappa[s]...

where  E*denotes*the*expectation and kappa[r]  is the r-th cumulant. These statistics are called polykays or generalized k-statistics.

 

> 

makeCTR_s := proc (N) [seq(k[i] = add(mul(mu[k], k = v)*(-1)^(nops(v)-1)*factorial(nops(v)-1)*factorial(i)/mul(factorial(x)^numboccur(v, x)*factorial(numboccur(v, x)), x = {op(v)}), v = partition(i)), i = 1 .. N)] end proc:

makeMu := proc () local u, v, N, eu; N := add(i, i = args); eu := [seq(mu[i] = 1, i = 1 .. N)]; if nargs = 1 then u := [seq([[x]], x = partition(args[1]))] else u := [seq([seq(x, x = partition(i))], i = args)] end if; u := op(add(mul((-1)^(nops(i)-1)*factorial(nops(i)-1)*mul(mu[k], k = i)*countP(i), i = v)*fd(N-1, add(nops(k), k = v))/countP([seq(op(i), i = v)]), v = `if`(nargs = 1, comb(u, 1, []), [comb(u, 1, [])]))); seq(simplify(v/(eval(v, eu))) = eval(v, eu), v = u) end proc:

comb := proc (V, ptr, Y) if ptr = nops(V)+1 then return Y end if; seq(comb(V, ptr+1, [op(Y), L]), L = V[ptr]) end proc:

countP := proc (u) factorial(add(x, x = u))/mul(x^numboccur(u, x)*factorial(numboccur(u, x)), x = {op(u)}) end proc:

ps := proc () local u, v, N; N := add(x, x = args); u := expand(eval(makeTab_ss(N), makeCTR_s(N))); v := expand(eval(u, [makeMu(args), mu = 0])); add(x[1]*x[2], x = v)/mul(n-x, x = 0 .. N-1) end proc:

Example; The*polykay*of*order*(2, 1)*is

> 

ps(2, 1)

(-S[1]^3+(1+n)*S[1]*S[2]-n*S[3])/(n*(n-1)*(n-2))

(3.2.1)

Multiset subdivision

The makeTab function lists all subdivisions of a multiset (for a thorough discussion see [2]). With multivariate k-statistics and/or multivariate polykays, this function speeds up the overall procedure (see [7]).  For details on the procedure and mathematical background see also [3] and [4].

 

> 

nRep := proc (u) mul(factorial(x[2]), x = convert(u, multiset)) end proc:

makeTab := proc () local U; if add(x, x = args) = 0 then return 0 end if; U := [seq(`if`(args[i] = 0, NULL, [seq([[seq(P || i^z, z = y)], multinomial(args[i], seq(r, r = y))], y = partition(args[i]))]), i = 1 .. nargs)]; if nops(U) = 1 then [seq([x[1], x[2]/nRep(x[1])], x = op(U))] else [seq(URmV(op(x)), x = [comb(U, 1, [])])] end if end proc:
> 

NULL

Multivariate k-statistics

The multivariate k-statistic is the unique symmetric unbiased estimator of the n-th joint cumulant kappa[i, j, () .. ()] of a given multivariate statistical distribution, that is
  

E[k[i, j, () .. ()]] = kappa[i, j, () .. ()]

where E denotes the expectation. They are expressed in terms of multivariate power sums involving the random vectors of a random sample.

 

> 

makeTab_sm := proc () local vP, U; U := makeTab(args); vP := sort([seq(P || i, i = 1 .. nargs)]); [seq([mul(k[add(degree(x, vP[i]), i = 1 .. nops(vP))], x = y[1])*y[2], mul(S[seq(degree(x, vP[i]), i = 1 .. nops(vP))], x = y[1])], y = U)] end proc:

km := proc () local u, v, N; N := add(x, x = args); v := expand(eval(makeTab_sm(args), makeK_s(N))); u := [seq(x^i = (-1)^(i-1)*factorial(i-1)*fd(N-1, i), i = 1 .. N)]; add(x[1]*x[2], x = expand(eval(v, u)))/mul(n-x, x = 0 .. N-1) end proc:

Example; The*multivariate*k-statistics*of*order*(2, 1)is

> 

km(2, 1)

(n^2*S[2, 1]-n*S[0, 1]*S[2, 0]-2*n*S[1, 0]*S[1, 1]+2*S[0, 1]*S[1, 0]^2)/(n*(n-1)*(n-2))

(3.4.1)

Note that the power sums in (3.4.1) are:

                S[2, 1] = Sum(X[i]^2*Y[i], i = i .. n)    S[0, 1] = Sum(Y[i], i = i .. n)    S[2, 0] = Sum(X[i]^2, i = i .. n)    S[1, 0] = Sum(X[i], i = i .. n)   
S[1, 1] = Sum(X[i]*Y[i], i = i .. n)

where (X[1], () .. (), X[n]) refers to the population X and  (Y[1], () .. (), Y[n]) refers to the population "Y. "

 

Multivariate Polykays

The multivariate polykay "k[r,s,...;i,j,...;...]" is an unbiased estimator of the cumulant product kappa[r, s, () .. ()]*kappa[i, j, () .. ()]..., that is

E["k[r,s,...;i,j,...;...]"]  = kappa[r, s, () .. ()]*kappa[i, j, () .. ()]...

where  E denotes the expectation and  kappa[r, s, () .. ()]*kappa[i, j, () .. ()]  are joint cumulants of order (r,s,...) and (i,j,..) respectively.

 

> 

makeTab_mm := proc () local vP, U; U := makeTab(args); vP := sort([op(indets(U))]); [seq([mul(k[seq(degree(x, vP[i]), i = 1 .. nops(vP))], x = y[1])*y[2], mul(S[seq(degree(x, vP[i]), i = 1 .. nops(vP))], x = y[1])], y = U)] end proc:

ctr := proc () local vP, U; U := makeTab(args); vP := [seq(P || i, i = 1 .. nargs)]; add(v[2]*(-1)^(nops(v[1])-1)*factorial(nops(v[1])-1)*mul(mu[seq(degree(x, vP[i]), i = 1 .. nops(vP))], x = v[1]), v = makeTab(args)) end proc:

makeCTR_m := proc () [seq(k[op(i)] = ctr(op(i)), i = comb([seq([seq(x, x = 0 .. y)], y = args)], 1, []))] end proc:

ricVtab := proc (v, V, N) local u, vv; vv := sort(v[1]); for u in V do if sort(u[1]) = vv then return v[2]*fd(N-1, nops(v[1]))/u[2] end if end do; return 0 end proc:

Example; The*multivariate*polykay*of*order*([1, 1], [1, 0])is

> 

pm([1, 1], [1, 0])

(n*S[1, 0]*S[1, 1]-S[0, 1]*S[1, 0]^2-n*S[2, 1]+S[0, 1]*S[2, 0])/(n*(n-1)*(n-2))

(3.5.1)

Note that the vectors in the input list must be the same length: if any vector is "smaller," the function adds 0 until it reaches the correct length.

Master function "polyk" handling all cases

This function allows us to recall all functions for generate k-statistics, polykays and their multivariate generalizzations.
The input is the following:

-  for generate  k-statistics k[r]  the parameter is:  [ r ]

-  for generate  polykays k[`r;s`] the parameter is:  [ r ], [ s ]

-  for generate multivariate k-statistics k[r, s]  the parameter is:  [ r, s ]

-  for generate multivariate polykays k[`r,s; u,v`]  the parameter is:  [ r, s ],  [ u, v]

> 

polyk := proc () if nargs = 1 then if nops(args[1]) = 1 then ks(op(args[1])) else km(op(args[1])) end if elif add(`if`(nops(x) = 1, 0, 1), x = args) = 0 then ps(seq(op(x), x = args)) else pm(args) end if end proc:

 

Esample: The k-statistic of order 3 is

> 

polyk([3])

(n^2*S[3]-3*n*S[1]*S[2]+2*S[1]^3)/(n*(n-1)*(n-2))

(3.6.1)

Example; The*polykay*of*order*(2, 1)*is

> 

polyk([2], [1])

(-S[1]^3+(1+n)*S[1]*S[2]-n*S[3])/(n*(n-1)*(n-2))

(3.6.2)

Example; The*multivariate*k-statistics*of*order*(2, 1)is

> 

polyk([2, 1])

(n^2*S[2, 1]-n*S[0, 1]*S[2, 0]-2*n*S[1, 0]*S[1, 1]+2*S[0, 1]*S[1, 0]^2)/(n*(n-1)*(n-2))

(3.6.3)

Example; The*multivariate*polykay*of*order*([1, 1], [1, 0])is

> 

polyk([1, 1], [1])

(n*S[1, 0]*S[1, 1]-S[0, 1]*S[1, 0]^2-n*S[2, 1]+S[0, 1]*S[2, 0])/(n*(n-1)*(n-2))

(3.6.4)

Replacing symbols with numerical data

The following functions allow the random variables of a random sample  to be replaced with corresponding numerical observations.  


The powS function  returns the sum of the r-th powers of the observations.

> 

powS := proc () if nargs = 1 then Sum('X'[i]^args[1], i = 1 .. 'n') else Sum(mul('X'[i, j]^args[j], j = 1 .. nargs), i = 1 .. 'n') end if end proc:

 

Example

> 

powS(5, 3, 1)

Sum(X[i, 1]^5*X[i, 2]^3*X[i, 3], i = 1 .. n)

(3.6.1.1)
> 

powS(9);

Sum(X[i]^9, i = 1 .. n)

(3.6.1.2)

 

The npolyk function NULLcomputes a numerical value of a k-statistic or a polykay (including the multivariate cases),  replacing the random variables of the random sample with the corresponding numerical observations. The parameters are the following:

-  for generate  k-statistics k[r]  the parameters are:  [ r ],  [ [ n1, n2, ...] ]  

-  for generate  polykays k[`r;s`] the parameters are:  [ [ r ], [ s ] ],  [ [ n1, n2, ...] ]  

-  for generate multivariate k-statistics k[r, s]  the parameters are: [ [ r , s ] ],  [ [ n1a, n2a], [ n1b, n2b] , ... ]    

-  for generate multivariate polykays k[`r,s; u,v`]  the parameters are:  [ [ r , s ],  [ u , v] ],  [ [ n1a, n2a], [ n1b, n2b] , ... ]    

> 

npolyk := proc (V, data) local res, ind, N, vE, Si; Si := false; if nops(V) = 1 then if nops(op(V)) = 1 then Si := true; N := nops(op(data)); res := ks(op(op(V))) else N := nops(data); res := km(op(op(V))) end if elif add(`if`(nops(x) = 1, 0, 1), x = V) = 0 then Si := true; N := nops(op(data)); res := ps(seq(op(x), x = V)) else N := nops(data); res := pm(op(V)) end if; ind := `minus`(indets(res), {n}); vE := seq(y[1] = Sum(mul('X'[`if`(Si, op([j, i]), op([i, j]))]^y[2, j], j = 1 .. nops(y[2])), i = 1 .. N), y = seq([x, [op(x)]], x = ind)); eval(res, [eval(evalf(vE), [X = data]), n = N]) end proc:

Examples: suppose to have the following (univariate) numerical sample

> 

data := [[16.34, 10.76, 11.84, 13.55, 15.85, 18.20, 7.51, 10.22, 12.52, 14.68, 16.08, 19.43, 8.12, 11.20, 12.95, 14.77, 16.83, 19.80, 8.55, 11.58, 12.10, 15.02, 16.83, 16.98, 19.92, 9.47, 11.68, 13.41, 15.35, 19.11]]

[[16.34, 10.76, 11.84, 13.55, 15.85, 18.20, 7.51, 10.22, 12.52, 14.68, 16.08, 19.43, 8.12, 11.20, 12.95, 14.77, 16.83, 19.80, 8.55, 11.58, 12.10, 15.02, 16.83, 16.98, 19.92, 9.47, 11.68, 13.41, 15.35, 19.11]]

(3.6.1.3)

 

An estimation of the mean (k-statistic of first order) is

> 

npolyk([[1]], data)

14.02166667

(3.6.1.4)

An estimation of the variance (k-statistic of second order) is

> 

npolyk([[2]], data)

12.65006954

(3.6.1.5)

An estimation of the skewness is  k[3] / sqrt(k[2]^3)

> 

npolyk([[3]], data)/sqrt(npolyk([[2]], data)^3)

-0.3216232203e-1

(3.6.1.6)

An estimation of the kurtosis is  k[4] / k[2]^2

> 

npolyk([[4]], data)/npolyk([[2]], data)^2

-.8852926052

(3.6.1.7)

The estimation of the product kappa[3]*kappa[2] is    

> 

npolyk([[3], [2]], data)

-15.56102316

(3.6.1.8)

 

Examples: suppose to have the following (bivariate) numerical sample

> 

data := [[5.31, 11.16], [3.26, 3.26], [2.35, 2.35], [8.32, 14.34], [13.48, 49.45], [6.25, 15.05], [7.01, 7.01], [8.52, 8.52], [.45, .45], [12.08, 12.08], [19.39, 10.42]]

[[5.31, 11.16], [3.26, 3.26], [2.35, 2.35], [8.32, 14.34], [13.48, 49.45], [6.25, 15.05], [7.01, 7.01], [8.52, 8.52], [.45, .45], [12.08, 12.08], [19.39, 10.42]]

(3.6.1.9)

 

An estimation of the joint cumulant kappa[2, 1] is  

> 

npolyk([[2, 1]], data)

-23.73790404

(3.6.1.10)

An estimation of the product of the joint cumulant  kappa[2, 1] kappa[1, 0]is

> 

npolyk([[1, 1], [1, 0]], data)

294.2657621

(3.6.1.11)

An estimation of the product of the joint cumulant kappa[2, 1] kappa[2, 0]kappa[1, 0]is

> 

npolyk([[2, 1], [2, 0], [1, 0]], data)

12369.47450

(3.6.1.12)

 

Joint Cumulant of polykays: auxiliary functions

The expression produced by the following set of functions corresponds to formula (4.10) pp. 96 in [8]. This formula expresses the cumulants of multivariate polykays in terms of joint cumulants of the population underlying the  multivariate sample. This choice is motivated by the reduction of the overall final result complexity.
Referring to the simplest case of k-statistics, the procedure is as follows:
1) the cumulants of the k-statistics are expressed in terms of their moments using the function ctr ;
2) the moments of the k-statistics are expressed in terms of the moments of the population underlying the sample using the functions pSS and E;
3) the moments of the population underlying the sample are expressed in terms of its cumulants using the function rct.

The same procedure works for multivariate polykays.

In the following details are given for each step.  

 

Cumulants and Moments

The ctr function (see [1]) expresses cumulants of a random vector in terms of joint moments mu[i, j, () .. ()]. Details on the procedure and mathematical background are given in [1] and [3]

Example: The joint cumulant kappa[2, 1]in terms of joint moments is

> 

ctr(2, 1)

2*mu[0, 1]*mu[1, 0]^2-mu[0, 1]*mu[2, 0]-2*mu[1, 0]*mu[1, 1]+mu[2, 1]

(4.1.1)

The rct function expresses joint moments of a random vector in terms of joint  cumulants kappa[i, j, () .. ()]. Details on the procedure and mathematical background are given in [3].

> 

rtc := proc () local vP, U; U := makeTab(args); vP := [seq(P || i, i = 1 .. nargs)]; add(v[2]*mul(kappa[seq(degree(x, vP[i]), i = 1 .. nops(vP))], x = v[1]), v = makeTab(args)) end proc:

Example: The joint moment mu[2, 1]in terms of joint cumulants is

> 

rtc(2, 1)

kappa[0, 1]*kappa[1, 0]^2+kappa[0, 1]*kappa[2, 0]+2*kappa[1, 0]*kappa[1, 1]+kappa[2, 1]

(4.1.2)

Product of Augmented Symmetric Functions  

Products of multivariate polykays return products of power sums. For example the expectation of k[2] k[1] is:
E[k[2] k[1]] = E[(n*S[2]-S[1]^2)/(n*(n-1))*(S[1]/n)] = (n*E[S[2]*S[1]]-E[S[1]^3])/(n^2*(n-1))
where Typesetting[delayDotProduct](E*denotes*the*expectation.To, express, true)*E[S[2]*S[1]]and E[S[1]^3] in terms of moments, a fundamental result of estimation might be used if these products are transformed in augmented symmetric functions, see [9]. For univariate sample, an augmented symmetric function involves products of variables with different indexes, that is
A[[i, j, k, () .. ()]]"=(∑) X[a]^(i)X[b]^(j)"X[c]^k...
The pSS function carries out this task.

> 

uRv := proc (u, v) local w; w := [sort([op(u[1]), v*'X']), u[2]], seq(`if`(has(u[1, i], 'X'), NULL, [sort([op(u[1, 1 .. i-1]), u[1, i]+v*'X', op(u[1, i+1 .. nops(u[1])])]), u[2]]), i = 1 .. nops(u[1])); w end proc:
NULL
pSS := proc () local u, v, U, V; V := sort([args], proc (x, y) evalb(nops(y) < nops(x)) end proc); U := pS(S[V[1]], S[V[2]]); for v in V[3 .. nops(V)] do U := add(pS(u, S[v]), u = [op(U)]) end do; U end proc:

Example: product of power sums S[1]*S[2]^3 in terms of  augmented symmetric functions:

> 

pSS([1], [2], [2], [2])

A[[1, 2, 2, 2]]+3*A[[2, 2, 3]]+3*A[[1, 2, 4]]+3*A[[3, 4]]+3*A[[2, 5]]+A[[1, 6]]+A[[7]]

(4.2.1)

Example: product of power sums S[2, 1]*S[1, 1]*S[2, 3] in terms of  multivariate augmented symmetric functions:

> 

pSS([[2, 1]], [[1, 1]], [[2, 3]])

A[[[1, 1], [2, 1], [2, 3]]]+A[[[2, 3], [3, 2]]]+A[[[1, 1], [4, 4]]]+A[[[2, 1], [3, 4]]]+A[[[5, 5]]]

(4.2.2)

The following is a test explaing how to transform the power sum product S[1, 1]*S[1, 2] with in augmented symmetric functions with n=3:

> 

F1 := (Sum(x[i]*y[i], i = 1 .. 3))*(Sum(x[i]*y[i]^2, i = 1 .. 3)):F1 := sort(expand(evalf(F1)))

x[1]^2*y[1]^3+x[1]*x[2]*y[1]^2*y[2]+x[1]*x[2]*y[1]*y[2]^2+x[1]*x[3]*y[1]^2*y[3]+x[1]*x[3]*y[1]*y[3]^2+x[2]^2*y[2]^3+x[2]*x[3]*y[2]^2*y[3]+x[2]*x[3]*y[2]*y[3]^2+x[3]^2*y[3]^3

(4.2.3)
> 

pSS([[1, 1]], [[1, 2]])

A[[[1, 1], [1, 2]]]+A[[[2, 3]]]

(4.2.4)
> 

F2 := Sum(Sum(x[i]*y[i]*x[j]*y[j]^2, i = 1 .. 3), j = 1 .. 3)+Sum(x[i]^2*y[i]^3, i = 1 .. 3):F2 := sort(add(add(`if`(i = j, 0, x[i]*y[i]*x[j]*y[j]^2), i = 1 .. 3), j = 1 .. 3)+add(x[i]^2*y[i]^3, i = 1 .. 3))

x[1]^2*y[1]^3+x[1]*x[2]*y[1]^2*y[2]+x[1]*x[2]*y[1]*y[2]^2+x[1]*x[3]*y[1]^2*y[3]+x[1]*x[3]*y[1]*y[3]^2+x[2]^2*y[2]^3+x[2]*x[3]*y[2]^2*y[3]+x[2]*x[3]*y[2]*y[3]^2+x[3]^2*y[3]^3

(4.2.5)
> 

evalb(F2 = F1);

true

(4.2.6)

Expected values of Augmented Simmetric Functions

The fundamental result of estimation we refer to allows to express expectation of augmented symmetric functions in terms of moments [10], that is
E[A[[i, j, k, () .. ()]]] = n(n-1)...(n-s+1)mu[i]*mu[j]*mu[k] .. ()
where s denotes the number of elements in the list "[i,j,k,...]. "The product n(n-1)...(n-s+1) is the lower factorial n[s]. The df  function computes the lower factorial. This result can be suitably generalized to the multivariate case.

> 

df := proc (m) options operator, arrow; mul(n-i, i = 0 .. m-1) end proc

proc (m) options operator, arrow; mul(n-i, i = 0 .. m-1) end proc

(4.3.1)

Example: the decreasing factorial n[3] 

> 

df(3)

n*(n-1)*(n-2)

(4.3.2)

The E function computes the expectation of augmented symmetric functions.

> 

E := proc () local v, vArgs; if nargs = 1 then return n*mu[op(args)] end if; vArgs := seq([u], u = [args]); v := seq([(1/2)*op(u)[1], op(op(u)[2])], u = [op(2*pSS(vArgs))]); add(df(nops(u[2]))*u[1]*mul(mu[op(r)], r = u[2]), u = v) end proc:

 

Example: the expectation of A[[[1, 1], [2, 1]]]+A[[[3, 2]]] is

> 

E([1, 1], [2, 1])+E([3, 2])

n*(n-1)*mu[1, 1]*mu[2, 1]+2*n*mu[3, 2]

(4.3.3)


The fE function repeatedly executes the E function on an expression.

> 

fE := proc (s) local sn, sd, vI, vE; sd := denom(s); sn := expand(numer(s)); vI := `minus`(indets(s), {n}); vE := [seq(u = 1, u = vI)]; if nops(vI) = 1 or type(sn, `*`) then sn := [seq([[seq([[op(w)], degree(u, w)], w = vI)], eval(u, vE)], u = [sn])] else sn := [seq([[seq([[op(w)], degree(u, w)], w = vI)], eval(u, vE)], u = sn)] end if; sn := [seq([[seq(`if`(w[2] = 0, NULL, w), w = u[1])], u[2]], u = sn)]; expand(add(E(seq(`$`(w[1], w[2]), w = u[1]))*u[2], u = sn))/sd end proc:
> 

``

Example: the expectation of ((n*S[2]-S[1]^2)/(n*(n-1))*(S[1]/n))*is

> 

fE(polyk([2])*polyk([1]));

(-n^3*mu[1]^3+n^3*mu[1]*mu[2]+3*n^2*mu[1]^3-4*n^2*mu[1]*mu[2]-2*n*mu[1]^3+n^2*mu[3]+3*n*mu[1]*mu[2]-n*mu[3])/(n^2*(n-1))

(4.3.4)

Example: the expectation of  S[2]*S[3]^2+is*S[4]

> 

fE(S[2]*S[3]^2+S[4])

n^3*mu[2]*mu[3]^2-3*n^2*mu[2]*mu[3]^2+n^2*mu[2]*mu[6]+2*n^2*mu[3]*mu[5]+2*n*mu[2]*mu[3]^2-n*mu[2]*mu[6]-2*n*mu[3]*mu[5]+n*mu[4]+n*mu[8]

(4.3.5)

Joint Cumulant of polykays: the main function

The jcks function returns the joint cumulants of multivariate polykays according to the procedure described in the previous section.

> 

jcks := proc (vK, vC) local s, vE, nv, dv, vS, evK; evK := [seq(polyk(vK[i]), i = 1 .. nops(vK))]; s := ctr(op(vC)); vE := [seq(u = fE(mul(evK[i]^[op(u)][i], i = 1 .. nops([op(u)]))), u = indets(s))]; s := simplify(eval(s, vE)); vE := [seq(u = rtc(op(u)), u = [seq(`if`(op(0, u) = mu, u, NULL), u = indets(s))])]; s := simplify(eval(s, vE)); nv := numer(s); dv := denom(s); vS := [seq(v = 1, v = `minus`(indets(s), {n}))]; if nops(vS) = 1 then return s end if; vS := [seq([u[1]/u[2], u[2]], u = [seq([v, eval(v, vS)], v = nv)])]; vS := {seq([v[1], add(`if`(u[1] = v[1], u[2], 0), u = vS)], v = vS)}; add(v[1]*simplify(v[2]/dv), v = vS) end proc:

Examples

 

kappa[1](k[2])is the first cumulant of k[2], that*can*be*computed*as

> 

jcks([[2]], [1])

kappa[2]

(5.1.1)

kappa[2](k[2])is the second cumulant of k[2], that*can*be*computed*as

> 

jcks([[2]], [2])

kappa[4]/n+2*kappa[2]^2/(n-1)

(5.1.2)

kappa[3](k[2])is the third cumulant of k[2], that*can*be*computed*as

> 

jcks([[2]], [3])

kappa[6]/n^2+8*kappa[2]^3/(n-1)^2+4*kappa[3]^2*(n-2)/(n*(n-1)^2)+12*kappa[2]*kappa[4]/(n*(n-1))

(5.1.3)


In the following the variances of the first few k-statistics are computed, see [11]. Note that "Var(X)=Cov(X,X)."

> 

var(k[1]) = jcks([[1], [1]], [1, 1])

var(k[1]) = kappa[2]/n

(5.1.4)
> 

var(k[2]) = jcks([[2], [2]], [1, 1])

var(k[2]) = kappa[4]/n+2*kappa[2]^2/(n-1)

(5.1.5)
> 

var(k[3]) = jcks([[3], [3]], [1, 1])

var(k[3]) = kappa[6]/n+6*kappa[2]^3*n/((n-1)*(n-2))+9*kappa[3]^2/(n-1)+9*kappa[2]*kappa[4]/(n-1)

(5.1.6)
> 

var(k[4]) = jcks([[4], [4]], [1, 1])

var(k[4]) = kappa[8]/n+24*kappa[2]^4*n*(n+1)/((n-1)*(n-2)*(n-3))+34*kappa[4]^2/(n-1)+144*kappa[2]*kappa[3]^2*n/((n-1)*(n-2))+16*kappa[2]*kappa[6]/(n-1)+72*kappa[2]^2*kappa[4]*n/((n-1)*(n-2))+48*kappa[3]*kappa[5]/(n-1)

(5.1.7)

 

 

More examples are given in [10], pag. 265.
The joint cumulant "kappa(2^2 1)" of order (2, 1) of (k[2], k[1]) is

> 

jcks([[2], [1]], [2, 1])

kappa[5]/n^2+4*kappa[2]*kappa[3]/(n*(n-1))

(5.1.8)

The joint cumulant kappa(2^2*1^2) of order (2, 2) of (k[2], k[1])*is

> 

jcks([[2], [1]], [2, 2])

kappa[6]/n^3+4*kappa[3]^2/(n^2*(n-1))+4*kappa[2]*kappa[4]/(n^2*(n-1))

(5.1.9)

 

 

The joint cumulant kappa(3*2) of order (1, 1) of Typesetting[delayDotProduct]((k[3], k[2])*is, see[11]*pag.271, true)

> 

jcks([[3], [2]], [1, 1])

kappa[5]/n+6*kappa[2]*kappa[3]/(n-1)

(5.1.10)

NULL

The joint cumulant kappa(4*2^2) of order (1,2) of (k[4], k[2])*is see[11]*pag.272

> 

jcks([[4], [2]], [1, 2])

kappa[8]/n^2+4*kappa[4]^2*(7*n-10)/(n*(n-1)^2)+120*kappa[2]*kappa[3]^2/(n-1)^2+20*kappa[2]*kappa[6]/(n*(n-1))+80*kappa[2]^2*kappa[4]/(n-1)^2+8*kappa[3]*kappa[5]*(5*n-7)/(n*(n-1)^2)

(5.1.11)

 

The joint cumulant of order (1,1) of "(k^(r,s),k^(t,u)) is "(see[8] pag.94)

> 

jcks([[1, 1, 0, 0], [0, 0, 1, 1]], [1, 1])

kappa[1, 1, 1, 1]/n+kappa[1, 0, 0, 1]*kappa[0, 1, 1, 0]/(n-1)+kappa[1, 0, 1, 0]*kappa[0, 1, 0, 1]/(n-1)

(5.1.12)


The joint cumulant of order (1,1,1) of "kappa[k](r,s|t,u|v,w) is "(see[8] pag.94)

> 

jcks([[1, 1, 0, 0, 0, 0], [0, 0, 1, 1, 0, 0], [0, 0, 0, 0, 1, 1]], [1, 1, 1]);

kappa[1, 0, 0, 0, 1, 0]*kappa[0, 1, 1, 1, 0, 1]/(n*(n-1))+kappa[1, 0, 0, 1, 0, 0]*kappa[0, 1, 1, 0, 1, 1]/(n*(n-1))+kappa[1, 0, 0, 1, 1, 1]*kappa[0, 1, 1, 0, 0, 0]/(n*(n-1))+kappa[1, 0, 1, 0, 0, 0]*kappa[0, 1, 0, 1, 1, 1]/(n*(n-1))+kappa[1, 0, 1, 0, 1, 1]*kappa[0, 1, 0, 1, 0, 0]/(n*(n-1))+kappa[1, 0, 1, 1, 0, 1]*kappa[0, 1, 0, 0, 1, 0]/(n*(n-1))+kappa[1, 0, 1, 1, 1, 0]*kappa[0, 1, 0, 0, 0, 1]/(n*(n-1))+kappa[1, 1, 0, 1, 0, 1]*kappa[0, 0, 1, 0, 1, 0]/(n*(n-1))+kappa[1, 1, 0, 1, 1, 0]*kappa[0, 0, 1, 0, 0, 1]/(n*(n-1))+kappa[1, 1, 1, 0, 0, 1]*kappa[0, 0, 0, 1, 1, 0]/(n*(n-1))+kappa[1, 1, 1, 0, 1, 0]*kappa[0, 0, 0, 1, 0, 1]/(n*(n-1))+kappa[1, 0, 0, 0, 0, 1]*kappa[0, 1, 0, 1, 0, 0]*kappa[0, 0, 1, 0, 1, 0]/(n-1)^2+kappa[1, 0, 0, 0, 0, 1]*kappa[0, 1, 1, 0, 0, 0]*kappa[0, 0, 0, 1, 1, 0]/(n-1)^2+kappa[1, 0, 0, 0, 1, 0]*kappa[0, 1, 0, 1, 0, 0]*kappa[0, 0, 1, 0, 0, 1]/(n-1)^2+kappa[1, 0, 0, 0, 1, 0]*kappa[0, 1, 1, 0, 0, 0]*kappa[0, 0, 0, 1, 0, 1]/(n-1)^2+kappa[1, 0, 0, 1, 0, 0]*kappa[0, 1, 0, 0, 0, 1]*kappa[0, 0, 1, 0, 1, 0]/(n-1)^2+kappa[1, 0, 0, 1, 0, 0]*kappa[0, 1, 0, 0, 1, 0]*kappa[0, 0, 1, 0, 0, 1]/(n-1)^2+kappa[1, 0, 1, 0, 0, 0]*kappa[0, 1, 0, 0, 0, 1]*kappa[0, 0, 0, 1, 1, 0]/(n-1)^2+kappa[1, 0, 1, 0, 0, 0]*kappa[0, 1, 0, 0, 1, 0]*kappa[0, 0, 0, 1, 0, 1]/(n-1)^2+kappa[1, 0, 0, 0, 0, 1]*kappa[0, 1, 1, 1, 1, 0]/(n*(n-1))+kappa[1, 0, 0, 1, 0, 1]*kappa[0, 1, 1, 0, 1, 0]*(n-2)/(n*(n-1)^2)+kappa[1, 0, 0, 1, 1, 0]*kappa[0, 1, 1, 0, 0, 1]*(n-2)/(n*(n-1)^2)+kappa[1, 0, 1, 0, 0, 1]*kappa[0, 1, 0, 1, 1, 0]*(n-2)/(n*(n-1)^2)+kappa[1, 0, 1, 0, 1, 0]*kappa[0, 1, 0, 1, 0, 1]*(n-2)/(n*(n-1)^2)+kappa[1, 1, 1, 1, 1, 1]/n^2

(5.1.13)

Estimation of joint cumulants of polykays

Given a random numerical sample, the njcks function allows estimating the joint cumulants of multivariate polykays. The function invokes the nployk function in the background section. The first two parameters are the same as in the njcks function, while the data vector must be entered as third parameter. The procedure consists in expressing the joint cumulant of a multivariate polykay in terms of joint cumulants of the underlying population and then replacing the occurrences of these joint cumulants and/or their products with the corresponding multivariate numerical polykays.

 

> 

njcks := proc (vK, vC, vD) local s, u, sn, sd, vI, vE, N; s := jcks(vK, vC); if nops(vK[1]) = 1 then N := nops(op(vD)) else N := nops(vD) end if; sd := denom(s); sn := expand(numer(s)); vI := `minus`(indets(s), {n}); vE := [seq(u = 1, u = vI)]; if nops(vI) = 1 or type(sn, `*`) then sn := [seq([[seq([[op(w)], degree(u, w)], w = vI)], eval(u, vE)], u = [sn])] else sn := [seq([[seq([[op(w)], degree(u, w)], w = vI)], eval(u, vE)], u = sn)] end if; sn := [seq([[seq(`if`(w[2] = 0, NULL, w), w = u[1])], u[2]], u = sn)]; eval(expand(add(npolyk([seq(`$`(w[1], w[2]), w = u[1])], vD)*u[2], u = sn))/sd, ['n' = N]) end proc:

 

 

Example: suppose to have the following (univariate) numerical sample

> 

data1 := [[16.34, 10.76, 11.84, 13.55, 15.85, 18.20, 7.51, 10.22, 12.52, 14.68, 16.08, 19.43, 8.12, 11.20, 12.95, 14.77, 16.83, 19.80, 8.55, 11.58, 12.10, 15.02, 16.83, 16.98, 19.92, 9.47, 11.68, 13.41, 15.35, 19.11]];

[[16.34, 10.76, 11.84, 13.55, 15.85, 18.20, 7.51, 10.22, 12.52, 14.68, 16.08, 19.43, 8.12, 11.20, 12.95, 14.77, 16.83, 19.80, 8.55, 11.58, 12.10, 15.02, 16.83, 16.98, 19.92, 9.47, 11.68, 13.41, 15.35, 19.11]]

(5.2.1)

 

The estimation of the joint cumulant of order (2,2) of (k[2], k[1])*isNULL

> 

njcks([[2], [1]], [2, 2], data1)

0.1918448804e-1

(5.2.2)

Example: suppose to have the following (bivariate) numerical sample

> 

data2 := [[5.31, 11.16], [3.26, 3.26], [2.35, 2.35], [8.32, 14.34], [13.48, 49.45], [6.25, 15.05], [7.01, 7.01], [8.52, 8.52], [.45, .45], [12.08, 12.08], [19.39, 10.42]];

[[5.31, 11.16], [3.26, 3.26], [2.35, 2.35], [8.32, 14.34], [13.48, 49.45], [6.25, 15.05], [7.01, 7.01], [8.52, 8.52], [.45, .45], [12.08, 12.08], [19.39, 10.42]]

(5.2.3)

 

An estimation of the joint cumulant of order (1,1) of k[1, 1] k[1, 2]is

> 

njcks([[1, 1], [1, 2]], [1, 1], data2)

18250.79716

(5.2.4)
> 

``

Conclusions

The jcks function computes the cumulants of multivariate polykays. These routines are particularly useful for assessing the goodness of the estimate obtained through these estimators. For example, the variance of first- and fourth-order k-statistics provides insights into the goodness-of-fit estimation of the mean, variance, skeweness, and kurtosis, respectively. Further examples are provided in the last section of this worksheet.  

References

[1] (2009) Di Nardo E., Guarino G., Senato D. Fast algorithms for k-statistics, polykays and their multivariate generalizations. (Worksheet Maple Software: available at https://www.maplesoft.com/Applications/Detail.aspx?id=33041)

 

[2] (2009) Di Nardo E., Guarino G., Senato D. Multiset subdivisions. (Worksheet Maple Software: : available at https://www.maplesoft.com/Applications/Detail.aspx?id=33039).

 

[3] (2015) Di Nardo E. Symbolic calculus in mathematical statistics: a review. Seminaire Lotharingien de Combinatoire Vol. 67 (B67a), pp. 72

 

[4] (2009) Di Nardo E., Guarino G., Senato D. A new method for fast computing unbiased estimators of cumulants. Statistics and Computing 19, 155--165.

 

[5] (2008) Di Nardo E., Guarino G., Senato D. A unifying framework for k-statistics, polykays and their multivariate generalizations. Bernoulli 14, 440--468.

 

[6] (2006) Di Nardo E., Senato D. A symbolic method for k-statistics. Applied Mathematics Letters 19, 968--975.

 

[7] (2011) Di Nardo E., Guarino G., Senato D. A new algorithm for computing the multivariate Faa di Bruno's formula. Applied Mathematics and Computation 217, 6286--6295


[8] P. McCullagh, Tensor Methods in Statistics, Monographs on Statistics and Applied Probability, Taylor & Francis Group, 22 December 2017

 

[9] (2008) Di Nardo E., Guarino G., Senato D. Symbolic computation of moments of sampling distributions. Computational Statistics and Data Analysis 52, 4909--4922.

 

[10] Stuart, A., Ord, J.K., 1987. Kendall’s Advanced Theory of Statistics 1. Charles Griffin and Company Limited, London.

 

[11] Weisstein, Eric W. "k-Statistic." From MathWorld--A Wolfram Web Resource. https://mathworld.wolfram.com/k-Statistic.html

 


Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities

> 

``