These are Matlab programs for the SJS paper "On the Expectations of Equivariant
Matrix-valued Functions of Wishart and Inverse Wishart Matrices" (with Grant
Hillier), Scandinavian Journal of Statistics (2024), DOI 10.1111/sjos.12707,
and for its extension to the beta-Wishart case (Hillier and Kan, 2026,
"Recurrence relations for the moments of beta-Wishart and inverse beta-Wishart
matrices").  A user manual (wishmom_manual_v3.pdf) is distributed with the
programs.  If you have questions, comments, or bug reports, please send them
to Raymond.Kan@rotman.utoronto.ca

Raymond Kan
Rotman School of Management
University of Toronto

Requirements: Matlab R2020b or later (the functions qk_coeff.m, qkn_coeff.m,
wish_ps.m and iwish_ps.m use pagemtimes).  The Symbolic Math Toolbox is needed
only for the *_sym.m routines that produce symbolic matrices (qk_coeff_sym.m,
qkn_coeff_sym.m, wish_ps_sym.m, iwish_ps_sym.m, dkimap.m without n1) and for
the example scripts example1.m-example6.m; everything else, including the
analytical-expression routines wishmom_sym.m and iwishmom_sym.m, runs in
plain double precision.

All programs take the parameter alpha = 2/beta describing the type of Wishart
distribution: alpha = 2 (real, default), alpha = 1 (complex), alpha = 1/2
(quaternion).  Other positive values of alpha are accepted; for those only
the scalar quantities E[p_kappa(W)] and E[p_kappa(W^{-1})] have a meaning
(see the manual, Section 2.4).  Throughout, n1 = n-m+1-alpha.

Version history
Version 1.0: 4/30/2022, initial release.
Version 1.1: 6/30/2024, update iwishmom.m to check non-existence of moments.
Version 1.2: 8/16/2024, introduce denpoly.m to compute the coefficients of
             the denominator polynomial of \tilde{H}_k and \tilde{C}_k.
Version 2.0: 8/18/2024, introduce dkimapr.m, iwish_psr.m, and qkn_coeffr.m.
             These programs allow us to compute the numerator and denominator
             polynomials of \tilde{D}_k, \tilde{H}_k, and \tilde{C}_k without
             using symbolic toolbox. example4.m to example6.m are updated to
             use the new programs.
Version 2.1: 8/28/2024, rewrite dkimapr.m to make use denpoly.m.  The resulting
             expressions for the numerator and denominator polynomials of
             \tilde{D}_k are shortened as a result.
Version 2.2: 9/13/2024, fixed an error in the algorithm of dkimapr.m
Version 2.3: 9/14/2024, improve the algorithm of dkimapr.m
Version 3.0: 9/9/2026
   * iwishmom.m: corrected the existence test for the inverse moments.  The
     moment E[prod_j tr(W^{-j})^{f_j} W^{-iw}] of total degree
     k = sum_j j*f_j + iw exists if and only if (n-m+1)/alpha > k, i.e.
     n1 > (k-1)*alpha (real case: n-m+1 > 2k, complex case: n-m+1 > k).
     Versions 1.1-2.3 used the stricter test n1 > 2k and returned NaN for
     some moments that exist (e.g. E[W^{-1}] for a real Wishart matrix with
     n = m+1.5).  New helper iwish_exist.m.
   * wishmom.m, iwishmom.m: the degenerate case k = 0 now returns m (or I_m).
   * wishmom_rec.m: now correct for every alpha (the recursion used the real
     and complex coefficients n+r and n only; the general coefficient is
     n+(alpha-1)r).
   * qk_coeff.m, qk_coeff_sym.m, qkn_coeff.m, qkn_coeff_sym.m: alpha now
     defaults to 2 when omitted (as in all other routines).
   * dkimapr.m, qkn_coeffr.m, iwish_psr.m: the numerator/denominator
     polynomials are now verified numerically at a test value of n1; a
     warning is issued if the check fails (loss of precision for large k).
   * New: wishmom_sym.m and iwishmom_sym.m produce the analytical expressions
     (coefficients and LaTeX strings) of E[prod tr(W^j)^{f_j} W^{iw}] and of
     E[prod tr(W^{-j})^{f_j} W^{-iw}] without the symbolic toolbox (the same
     functionality as in the R package wishmom).  New helpers wm_index.m and
     wm_polystr.m.
   * New: wishmom_test.m, a script that checks the programs against each
     other, against the identities of the paper, and against Monte Carlo.
   * Documentation: comment headers corrected (several typos, iwish_psn.m
     returns \tilde{H}_k^{-1}); LaTeX manual wishmom_manual_v3.pdf.
   * Scratch files test.m, test1.m, test2.m and unidl.m removed; unid.m kept.
Version 3.0a: 9/10/2026 (documentation and tests only; no change to any
             computed result)
   * New: wm_weights.m, the reversibility weights 1/(alpha^l(lambda)*z_lambda)
     of the mapping matrix D_k.  With M = dkmap(k,alpha) one has
     diag(w)*M = M'*diag(w), so M is similar to a symmetric matrix and is
     therefore diagonalizable with real eigenvalues, all at most k*alpha;
     this is what makes D_{-k} = n1*I-M invertible for n1 > k*alpha.
   * wishmom_test.m: new check 5a (reversibility of D_k-n*I, and agreement of
     its spectrum with the root list of unid.m).
   * unid.m, dkimapr.m: corrected comment headers.  The root list is indexed
     by the cells (i,j) with (i+1)*(j+1) <= k+1 EXCLUDING (i,j) = (0,0): the
     empty cell is addable only to the empty diagram.  The code was already
     correct; only the comments described the set wrongly.  Note that
     denpoly.m uses a different list (all cells of the partitions of k, which
     does include the content 0 once, from Ct_1 = 1/n1).
   * verification/dkalgebra.py added.

Moments of Wishart
dkmap.m: creates D_k, but without the n on the diagonal
qk_coeff.m: program to compute coefficients for \mathcal{C}_k.
qk_coeffn.m: program to compute coefficients for \mathcal{C}_k using
             double precision when n is numeric.
qk_coeff_sym.m: program to compute coefficients for \mathcal{C}_k
                symbolically.
wish_ps.m: program to compute coefficients of \mathcal{H}_k
wish_psn.m: program to compute coefficients of \mathcal{H}_k using
            double precision when n is numeric.
wish_ps_sym.m: program to compute coefficients of \mathcal{H}_k
               symbolically.
wishmom.m: Computes E[\prod_{j=1}^r tr(W^j)^{f_j}*W^{iw}] using an
           analytical formula.
wishmom_sym.m: Analytical expression (polynomials in n) of
           E[\prod_{j=1}^r tr(W^j)^{f_j}*W^{iw}], with a LaTeX string.
wishmom_rec.m: Computes E[\prod_{j=1}^r tr(W^j)^{f_j}*W^{iw}] using a
               recursive algorithm.  Note that this is for demonstration
               purpose, and it is very slow when \sum_{j=1}^r j*f_j+iw
               is large.  Users are advised to use wishmom.m instead.
example1.m: a program that shows how to compute D_k using dkmap.m
example2.m: a program that shows how to compute \mathcal{H}_k using wish_ps.m
example3.m: a program that shows how to compute \mathcal{C}_k using qk_coeff.m

Moments of inverse Wishart
denpoly.m: compute the denominator polynomial for \tilde{H}_k and \tilde{C}_k
dkimap.m: creates \tilde{D}_k
dkimapr.m: creates coefficients of the numerator and denominator polynomials
           of \tilde{D}_k without using symbolic toolbox
qkn_coeff.m: program to compute coefficients for \tilde{C}_k^{-1}.
qkn_coeffn.m: program to compute coefficients for \tilde{C}_k^{-1}
              using double precision when n1 is numeric.
qkn_coeffr.m: program to compute coefficients for the numerator and
              denominator polynomials of \tilde{C}_k without using
              symbolic toolbox.
qkn_coeff_sym.m: program to compute coefficients for \tilde{C}_k^{-1}
                 symbolically.
iwish_ps.m: program to compute coefficients of \tilde{H}_k^{-1}
iwish_psn.m: program to compute coefficients of \tilde{H}_k^{-1} using
             double precision when n1 is numeric.
iwish_psr.m: program to compute coefficients for the numerator and
             denominator polynomials of \tilde{H}_k without using
             symbolic toolbox.
iwish_ps_sym.m: program to compute coefficients of \tilde{H}_k^{-1}
                symbolically.
iwishmom.m: Computes E[\prod_{j=1}^r tr(W^{-j})^{f_j}*W^{-iw}] using an
            analytical formula (NaN if the moment does not exist).
iwishmom_sym.m: Analytical expression (rational functions of n1) of
            E[\prod_{j=1}^r tr(W^{-j})^{f_j}*W^{-iw}], with a LaTeX string.
iwish_exist.m: existence test for the inverse moments of degree k.
example4.m: a program that shows how to compute \tilde{D}_k using dkimapr.m
example5.m: a program that shows how to compute \tilde{H}_k^{-1} using iwish_psr.m
example6.m: a program that shows how to compute \tilde{C}_k^{-1} using qkn_coeffr.m

Auxiliary programs:
ip_desc.m: Create integer partitions for an integer k, sorted in reverse
           lexicographical order
wm_index.m: position of a generator W^{iw} prod tr(W^j)^{f_j} in Q_k or q_k
wm_polystr.m: string representation of a polynomial
unid.m: the distinct eigenvalues of D_{-k} (minus n1), used by dkimapr.m
wm_weights.m: reversibility weights for D_k (see Version 3.0a above)
wishmom_test.m: consistency checks (run it after installation)

verification/: Python scripts (sympy) used to verify the programs and the
           results of the companion note: an independent implementation of
           the recurrences (wishrec.py), symbolic checks against the paper
           (chk_paper.py), against the Jack-polynomial moment formulae for
           general alpha (jackchk*.py), and Monte Carlo checks for the
           real, complex and quaternion cases (mc.py).
