Application Center - Maplesoft
  Submit your work
Application Center Applications Spectral k-statistics Preview

App Preview:

Spectral k-statistics

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

Learn about Maple
Download Application



This Maple worksheet accompanies the papers:

 

(2013) Di Nardo E., McCullagh P., Senato D. Natural statistics for spectral samples. Annals of Statistics. 41(2), 982-1004.  http://arxiv.org/abs/1302.5892

  Spectral k-statistics

 

E. Di Nardo*

elvira.dinardo@unibas.it

http://www.unibas.it/utenti/dinardo/home.html;

Tel: +39 0971205890, Fax: +39 0971205896

 

G. Guarino**

giuseppe.guarino@rete.basilicata.it

 

* Dipartimento di Matematica e Informatica, Università degli Studi della Basilicata,Viale dell'Ateneo Lucano n.10, 85100 Potenza, Italy

**Medical Scool, Università del Sacro Cuore (Rome branch), Largo Agostino Gemelli n.8, 00168 Roma, Italy

 

 

Introduction

Abstract: The algorithm constructs natural statistics of a spectral sample, by using convolutions on the symmetric group and the Weingarten function. These statistics are unbiased estimators of cumulants of traces.


Application Areas/Subject: Computational statistics

Keyword: Random matrix, cumulant of traces,  polykays

 

Initialization

 

> 

restart

> 

with(combinat, Chi, partition, permute);

[Chi, partition, permute]

[invperm, mulperms]

[Flatten]

(2.1)

Background

 Consruction of Schur function

The procedure Sch takes in input an integer partition and returns the Schur polynomial in N indeterminates all evaluated in 1.

Sch := proc (lambda, N) local mu; mu := combinat[conjpart](lambda); Matrix(nops(mu), nops(mu), proc (i, j) options operator, arrow; `if`(mu[i]+i-j < 0, 0, binomial(N, mu[i]+i-j)) end proc); expand(linalg[det](%)) end proc:

 

 

Example: for the partition (1,2,3) of the integer 6

> 

Sch([1, 2, 3], N)

(1/45)*N^6-(1/9)*N^4+(4/45)*N^2

(3.1.1)

for the partition (1^2,2) of the integer 4

> 

Sch([1, 1, 2], N)

(1/8)*N^4-(1/4)*N^3-(1/8)*N^2+(1/4)*N

(3.1.2)

 

 Weingarten function``

The procedure Wg takes in input an integer partition and returns the Weingarten function as a rational function in N.

The algorithm makes use of Schur polynomials and the character of the symmetric group.

 

> 

Wg := proc (mu, N) local q, uno; q := add(x, x = mu); uno := [`$`(1, q)]; factor(add(Chi(lambda, uno)^2*Chi(lambda, mu)/Sch(lambda, N), lambda = partition(q))/factorial(q)^2) end proc:

 

Example: for the partition (1,2,3) of the integer 6

> 

Wg([1, 2, 3], N)

-(2*N^2+13)/(N*(N+5)*(N+4)*(N+2)*(N+1)^2*(N-1)^2*(N-2)*(N-4)*(N-5))

(3.2.1)

 

Example: for the partition (1^2,2) of the integer 4

> 

Wg([1, 1, 2], N)

-1/((N+3)*(N+1)*(N-1)*(N-3)*N)

(3.2.2)

Spectral k-statistics

The Maple routines

 Some details on secondary Maple routines

The procedure compldisjcyc takes as input a permutation and returns its decomposition in disjoint cycles.

In the output there are also the fixed points;

> 

compldisjcyc := proc (a) local v, S; v := convert(a, disjcyc); S := `minus`({op(a)}, {seq(op(c), c = v)}); [seq([i], i = S), op(v)] end proc:

 

Example: for the permutation  which fix 1 and 4 and switch 2 and 3  

> 

compldisjcyc([1, 3, 2, 4])

[[1], [4], [2, 3]]

(4.1.1.1)

Example: for the permutation  which sends 1 in 2, 2 in 3, 3 in 4 and 4 in 1

> 

compldisjcyc([2, 3, 4, 1])

[[1, 2, 3, 4]]

(4.1.1.2)

Example: for the identity permutations

> 

compldisjcyc([1, 2, 3, 4])

[[1], [2], [3], [4]]

(4.1.1.3)
> 

``

The procedure tipo takes as input a permutation in disjoint cycles and returns its cycle type, that is how many cycles of each length are present in the cycle decomposition of the permutation.

> 

tipo := proc () local n, v; if nargs = 1 then n := max(seq(op(x), x = args[1])) else n := args[2] end if; v := sort([seq(nops(c), c = args[1])]); [`$`(1, n-add(x, x = v)), op(v)] end proc:

 

Examples: for the permutation which fix 1 and 4 and switch 2 and 3

> 

tipo([[1, 4], [2], [3]])

[1, 1, 2]

(4.1.1.4)

Examples: for the permutation  which sends 1 in 2, 2 in 3, 3 in 4 and 4 in 1

> 

tipo([[1, 2, 3, 4]])

[4]

(4.1.1.5)

Examples: for the identity permutation

> 

tipo([], 4)

[1, 1, 1, 1]

(4.1.1.6)

 

 The master function

 

The procedure CXX takes as input a permutation and returns the formula (5.6) in Theorem 5.2, see [1]. This formula corresponds to the convolution between products of traces of a spectral sample X and the inverse of a function giving the spectral sample size powered by the number of disjoint cycles.

> 

CX := proc () local b, n, binv; b := args[1]; binv := invperm(b); if nargs = 1 then n := max(seq(op(x), x = b)) else n := args[2] end if; add(Wg(tipo(mulperms(binv, convert(a, disjcyc)), n), N)*E(mul(Tr(mul(X[i], i = c)), c = compldisjcyc(a))), a = permute(n)); expand(%) end proc:

 

 

Examples: for the permutation  which fix 1 and 4 and switch 2 and 3  

> 

CXX([[1, 4], [2], [3]])

E(Tr(X)^2*Tr(X^2))/((N+3)*(N+2)*(N+1)*(N-1)*(N-2)*(N-3))-10*E(Tr(X^2)^2)/((N+3)*(N+2)*(N+1)*(N-1)*(N-2)*(N-3)*N)-20*E(Tr(X)*Tr(X^3))/((N+3)*(N+2)*(N+1)*(N-1)*(N-2)*(N-3)*N)-E(Tr(X^2)^2)/((N+3)*(N+1)*(N-1)*(N-3)*N)-4*E(Tr(X)*Tr(X^3))/((N+3)*(N+1)*(N-1)*(N-3)*N)-E(Tr(X)^4)/((N+3)*(N+1)*(N-1)*(N-3)*N)+10*E(Tr(X^4))/((N+3)*(N+2)*(N+1)*(N-1)*(N-2)*(N-3))+N^2*E(Tr(X)^2*Tr(X^2))/((N+3)*(N+2)*(N+1)*(N-1)*(N-2)*(N-3))

(4.1.2.1)

 

Examples: for the permutation  which sends 1 in 2, 2 in 3, 3 in 4 and 4 in 1  

> 

CXX([[2, 3, 4, 1]])

10*E(Tr(X)^2*Tr(X^2))/((N+3)*(N+2)*(N+1)*(N-1)*(N-2)*(N-3))-5*E(Tr(X^2)^2)/((N+3)*(N+2)*(N+1)*(N-1)*(N-2)*(N-3)*N)-20*E(Tr(X)*Tr(X^3))/((N+3)*(N+2)*(N+1)*(N-1)*(N-2)*(N-3)*N)-2*E(Tr(X^2)^2)/((N+3)*(N+1)*(N-1)*(N-3)*N)-4*E(Tr(X)*Tr(X^3))/((N+3)*(N+1)*(N-1)*(N-3)*N)-5*E(Tr(X)^4)/((N+3)*(N+2)*(N+1)*(N-1)*(N-2)*(N-3)*N)+E(Tr(X^4))/((N+3)*(N+2)*(N+1)*(N-1)*(N-2)*(N-3))+N^2*E(Tr(X^4))/((N+3)*(N+2)*(N+1)*(N-1)*(N-2)*(N-3))

(4.1.2.2)

 

The procedure eTr takes as input the output of the procedure CXX and replaces traces with power sums indexed by their powers.

> 

eTr := proc (y) options operator, arrow; S[degree(y)] end proc:eE := proc (y) options operator, arrow; y end proc:
eCXX := proc () local Ris, var; Ris := eval(CXX(args), [Tr = eTr, E = eE]); Ris := simplify(Ris); var := `minus`(indets(Ris), {N}); mul(factorial(nops(x)-1), x = compldisjcyc(op(args)))*collect(numer(Ris), var)/factor(denom(Ris)) end proc:

 

 

Example: for the permutation which sends 4 in 1, with fixed 2 and 3

> 

eCXX([[4, 2, 3, 1]])

(-5*S[1]^4+10*S[1]^2*S[2]*N+(-4*N^2-4)*S[3]*S[1]+(3-2*N^2)*S[2]^2+(N^3+N)*S[4])/((N+3)*(N+2)*(N+1)*(N-1)*(N-2)*(N-3)*N)

(4.1.2.3)

Example: for the permutation which sends 1 in 2, 2 in 3, 3 in 4 and 4 in 1  

> 

eCXX([[2, 3, 4, 1]])

6*(-5*S[1]^4+10*S[1]^2*S[2]*N+(-4*N^2-4)*S[3]*S[1]+(3-2*N^2)*S[2]^2+(N^3+N)*S[4])/((N+3)*(N+2)*(N+1)*(N-1)*(N-2)*(N-3)*N)

(4.1.2.4)
> 

NULL

Conclusions

This algorithm extends the symmetric functions k-statistics and polykays to spectral sampling. Spectral samples are eigenvalues of freely randomized classical sample.  The notion of freely randomized classical sample has been introduced for the first time in [1].

References

1] Di Nardo E., McCullagh P., Senato D. (2013) Natural statistics for spectral samples. Annals of Statistics. 41(2), 982-1004.  http://arxiv.org/abs/1302.5892

 

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

 

> 

NULL