Using the wishmom Package

Introduction

The wishmom package provides functions to compute the expectation of matrix-valued functions of β-Wishart and inverse β-Wishart distributions (β = 1: Real Wishart, β = 2: Complex Wishart, β = 4: Quaternion Wishart, β = 8: Octonion Wishart). The main functions in this package are wishmom() and iwishmom(), which handle the numerical computation of moments of β-Wishart and inverse β-Wishart distributions, respectively. The corresponding functions for generating analytical expressions of moments of β-Wishart and inverse β-Wishart distributions are wishmom_sym() and iwishmom_sym(). These programs are developed based on the results in Letac and Massam (2004) and Hillier and Kan (2024).

You can install the package and load it using:

install.packages("wishmom")
library("wishmom")


Mathematical Background

β-Wishart Distribution

The β-Wishart distribution is a fundamental distribution in multivariate statistics. When β = 1, 2, 4, 8, it is the real Wishart, complex Wishart, quaternion Wishart, and octonion Wishart, respectively. The density function of W ∼ Wmβ(n, Σ), i,e., a β-Wishart distribution with n degrees of freedom and an m × m covariance matrix Σ, is given by (when n > m − 1) (see Díaz-García and Gutiérrez-Jáimez (2011, Corollary 1))

$$ f(W) = \frac{\left(\frac{\beta}{2}\right)^\frac{mn\beta}{2}}{\Gamma_m^{(\beta)}\left(\frac{n \beta}{2}\right) |\Sigma|^\frac{n \beta}{2}}|W|^{\frac{(n-m+1)\beta}{2}-1} \mbox{etr}\left(-\frac{\beta \Sigma^{-1}W}{2}\right), $$

where

$$ \Gamma_m^{(\beta)}(a) = \pi^\frac{m(m-1)\beta}{4}\prod_{i=1}^m \Gamma\left(a-\frac{(i-1)\beta}{2}\right). $$

Note that we do not require n to be an integer but the definition of the density of W requires β = 1, 2, 4 or 8. However, if our interest is only on the functions of eigenvalues of W, we can generalize this to any real β > 0. Therefore, for any symmetric functions (say power-sum) of the eigenvalues of W, they can be well defined even when β is not equal to 1, 2, 4 or 8. For W ∼ Wmβ(n, Σ), the joint density of its eigenvalues λ1 ≥ ⋯ ≥ λm is given by (see Dresnky, Edelman, Genoar, Kan, and Koev (2021))

$$ f(\lambda_1,\ldots,\lambda_m) = \frac{| \Sigma|^{-\frac{n\beta}{2}}}{\mathcal{K}_m^ {(\beta)}\left(\frac{n\beta}{2}\right)} |\Lambda|^{\frac{(n-m+1)\beta}{2}-1} {}_0^{}F_0^{(\beta)}\left(-\frac{\beta}{2}\Lambda,\Sigma^{-1}\right)\prod_{1 \leq i < j \leq m}(\lambda_i-\lambda_j)^{\beta}, $$ where Λ = Diag(λ1, …, λm),

$$ \mathcal{K}_m^{(\beta)}(a) = \frac{\left(\frac{2}{\beta}\right)^{ma}} {\pi^{\frac{m(m-1)\beta}{2}}} \frac{\Gamma_m^{(\beta)}\left(\frac{m\beta}{2}\right)\Gamma_m^{(\beta)}(a)}{\left[\Gamma\left(\frac{\beta}{2}\right)\right]^m}, $$

$$ {}_0^{}F_0^{(\beta)}(A,B) = \sum_{k=0}^\infty \sum_{\kappa \vdash k} \frac{C_\kappa^{(\beta)}(A)C_\kappa^{(\beta)}(B)} {k!C_\kappa^{(\beta)}(I_m)}, $$

and Cκ(β)(X) is the Jack function of the eigenvalues of X.

Instead of using β, we use α = 2/β in our programs to describe the type of Wishart distribution. Therefore, α = 2 is for real Wishart, α = 1 is for complex Wishart, and α = 1/2 is for quaternion Wishart.

Moments of Matrix-valued Functions of β-Wishart and Inverse β-Wishart Distributions

Let λ = (λ1, …, λk) be an integer partition of a positive integer k, where |λ| = λ1 + … + λk = k, with λ1 ≥ λ2 ≥ ⋯ ≥ λk ≥ 0. The power-sum symmetric function pλ(W) of W corresponding to a partition λ is defined as

$$ p_{\lambda}(W) = \prod_{i=1}^{\ell(\lambda)}p_{\lambda_i}(W), $$

where ℓ(λ) is the number of non-zero parts of λ, and pi(W) = tr(Wi). We are interested in computing

𝔼[Wrpλ(W)]    and    𝔼[W−rpλ(W−1)],

where W ∼ Wmβ(n, Σ).

The method that we use is based on a generalization of the recurrence relations given in Hillier and Kan (2024) for which the cases of β = 1 and 2 were developed. Specifically, we have the following recurrence relations:

where ñ = n − m + 1 − (2/β) and λ(i) is λ with its i-th element removed. Together with the boundary conditions 𝔼[W] = nΣ and 𝔼[W−1] = Σ−1/ñ, we can obtain 𝔼[Erpλ(W)] and 𝔼[W−rpλ(W−1)]. Note that 𝔼[W−rpλ(W−1)] exists if and only if ñ > 2(r + |λ|).

Let k = r + |λ|. Hillier and Kan (2024) show that where cλ, ρ and c̃λ, ρ are constants that depend on n and ñ, respectively, but they do not depend on Σ. In addition, we have where hκ and h̃κ are constants that depend on n and ñ, repsectively, but they do not depend on Σ.

Using the recurrence relations, Hillier and Kan (2024) develop efficient algorithms for computing the constants cλ, ρ, c̃λ, ρ, hκ and h̃κ.


Main Functions in the Package

There are four main functions in this package: wishmom(), iwishmom(), wishmom_sym(), and iwishmom_sym(). The first two are used to numerically compute 𝔼[Wrpλ(W)] and 𝔼[W−rpλ(W−1)], respectively. The last two are used to generate an analytical expression of 𝔼[Wrpλ(W)] and 𝔼[W−rpλ(W−1)], respectively.

Moments of β-Wishart:

wishmom()

The function wishmom() computes $\mathbb{E}\left[\prod_{j=1}^r \mbox{tr}(W^j)^{f_j}W^{iw}\right]$ numerically, where W ∼ Wmβ(n, Σ). When iw = 0, it computes $\mathbb{E}[\prod_{j=1}^r \mbox{tr}(W^j)^{f_j}]$.

Arguments

  • n: degrees of freedom of the β-Wishart distribution
  • S: covariance matrix of the β-Wishart distribution
  • f: a vector of nonnegative integers fj that represents the power for tr(Wj), j = 1, …, r
  • iw: Power of W
  • alpha: The type of Wishart distribution (α = 2/β):
    • 1/2: Quaternion Wishart
    • 1: Complex Wishart
    • 2: Real Wishart (default)

Output

When iw = 0, it returns $\mathbb{E}[\prod_{j=1}^r \mbox{tr}(W^j)^{f_j}]$. When iw ≠ 0, it returns $\mathbb{E}[\prod_{j=1}^r \mbox{tr}(W^j)^{f_j}W^{iw}]$.

Examples

# Example 1: For E[tr(W)^4] with W ~ W_m^1(n,S), where n and S are defined below:
n <- 20
S <- matrix(c(25, 49,
              49, 109), nrow=2, ncol=2)
wishmom(n, S, 4) # iw = 0, for real Wishart distribution
#> [1] 8.705084e+13

# Example 2: For E[tr(W)^2*tr(W^3)*W^2] with W ~ W_m^1(n,S), where n and S, are defined below:
n <- 20
S <- matrix(c(25, 49,
              49, 109), nrow=2, ncol=2)
wishmom(n, S, c(2, 0, 1), 2, 2) # iw = 2, for real Wishart distribution
#>              [,1]         [,2]
#> [1,] 9.039462e+23 1.956948e+24
#> [2,] 1.956948e+24 4.258714e+24

# Example 3: For E[tr(W)^2*tr(W^3)] with W ~ W_m^2(n,S), where n and S are defined below:
n <- 20
S <- matrix(c(25, 49 + 2i,
              49 - 2i, 109), nrow=2, ncol=2)
wishmom(n, S, c(2, 0, 1), 0, 1) # iw = 0, for complex Wishart distribution
#> [1] 2.078126e+17

# Example 4: For E[tr(W)*tr(W^2)^2*tr(W^3)^2*W] with W ~ W_m^2(n,S), where n, S, are defined below:
n <- 20
S <- matrix(c(25, 49 + 2i,
              49 - 2i, 109), nrow=2, ncol=2)
wishmom(n, S, c(1, 2, 2), 1, 1) # iw = 1, for complex Wishart distribution
#>                            [,1]                       [,2]
#> [1,] 3.418999e+41+5.014362e+20i 6.943130e+41-2.833930e+40i
#> [2,] 6.943130e+41+2.833930e+40i 1.532151e+42-2.882805e+22i

wishmom_sym()

The function wishmom_sym() generates an analytical expression of $\mathbb{E}\left[\prod_{j=1}^r \mbox{tr}(W^j)^{f_j}W^{iw}\right]$, where W ∼ Wmβ(n, Σ). When iw = 0, it generates an analytical expression of $\mathbb{E}[\prod_{j=1}^r \mbox{tr}(W^j)^{f_j}]$.

Arguments

  • f: a vector of nonnegative integers fj that represents the power for tr(Wj), j = 1, …, r
  • iw: Power of W
  • alpha: The type of Wishart distribution (α = 2/β):
    • 1/2: Quaternion Wishart
    • 1: Complex Wishart
    • 2: Real Wishart (default)
  • latex: The type of output
    • TRUE: LaTeX expression
    • FALSE: Dataframe (default)

Output

When iw = 0, the output is an analytical expression of $\mathbb{E}[\prod_{j=1}^r \mbox{tr}(W^j)^{f_j}]$. When iw ≠ 0, the output is an analytical expression of $\mathbb{E}[\prod_{j=1}^r \mbox{tr}(W^j)^{f_j}W^{iw}]$. If latex = FALSE (default), the output is a data frame that stores the coefficients for the analytical expression. If latex = TRUE, the output is a $\LaTeX$ formatted string of the result in terms of n and Σ.

Examples

# Example 1: For E[tr(W)^4] with W ~ W_m^1(n,Sigma), represented as a dataframe:
wishmom_sym(4) # iw = 0, for real Wishart distribution
#>     kappa h_kappa
#> 1       4      48
#> 2     3,1     32n
#> 3     2,2     12n
#> 4   2,1,1   12n^2
#> 5 1,1,1,1     n^3

# Example 2: For E[tr(W)*tr(W^2)W] with W ~ W_m^1(n,Sigma), represented as a dataframe:
wishmom_sym(c(1,1), 1) # iw = 1, for real Wishart distribution
#>   i   rho          c
#> 1 1     3    4n^2+4n
#> 2 1   2,1 n^3+n^2+4n
#> 3 1 1,1,1        n^2
#> 4 2     2  2n^2+2n+8
#> 5 2   1,1         6n
#> 6 3     1 4n^2+4n+16
#> 7 4  <NA>     24n+24

# Example 3: For E[tr(W)^4] with W ~ W_m^2(n,Sigma), represented as a LaTeX string:
writeLines(wishmom_sym(4, 0, 1, latex=TRUE)) # iw = 0, for complex Wishart distribution
#> (6)p_{(4)}(\Sigma)
#> +(8n)p_{(3,1)}(\Sigma)
#> +(3n)p_{(2,2)}(\Sigma)
#> +(6n^2)p_{(2,1,1)}(\Sigma)
#> +(n^3)p_{(1,1,1,1)}(\Sigma)

# Example 4: For E[tr(W)*tr(W^2)W] with W ~ W_m^2(n,Sigma), represented as a LaTeX string:
writeLines(wishmom_sym(c(1, 1), 1, 1, latex=TRUE)) # iw = 1, for complex Wishart distribution
#> [(2n^2)p_{(3)}(\Sigma)+(n^3+2n)p_{(2, 1)}(\Sigma)+(n^2)p_{(1, 1, 1)}(\Sigma)]\Sigma
#> +[(n^2+2)p_{(2)}(\Sigma)+(3n)p_{(1, 1)}(\Sigma)]\Sigma^2
#> +(2n^2+4)p_{(1)}(\Sigma)\Sigma^3
#> +(6n)\Sigma^4

Moments of Inverse β-Wishart:

iwishmom()

The function iwishmom() computes $\mathbb{E}\left[\prod_{j=1}^r \mbox{tr}(W^{-j})^{f_j}W^{-iw}\right]$ numerically, where W ∼ Wmβ(n, Σ). When iw = 0, it computes $\mathbb{E}[\prod_{j=1}^r \mbox{tr}(W^{-j})^{f_j}]$.

Arguments

  • n: degrees of freedom of the β-Wishart distribution
  • S: covariance matrix of the β-Wishart distribution
  • f: a vector of nonnegative integers fj that represents the power for tr(W−j), j = 1, …, r
  • iw: Power of W−1
  • alpha: The type of Wishart distribution (α = 2/β):
    • 1/2: Quaternion Wishart
    • 1: Complex Wishart
    • 2: Real Wishart (default)

Output

When iw = 0, it returns $\mathbb{E}[\prod_{j=1}^r \mbox{tr}(W^{-j})^{f_j}]$. When iw ≠ 0, it returns $\mathbb{E}[\prod_{j=1}^r \mbox{tr}(W^{-j})^{f_j}W^{-iw}]$.

Examples

# Example 1: For E[tr(W^{-1})^2] with W ~ W_m^1(n,S), where n and S are defined below:
n <- 20
S <- matrix(c(25, 49,
              49, 109), nrow=2, ncol=2)
iwishmom(n, S, 2) # iw = 0, for real Wishart distribution
#> [1] 0.0006680892

# Example 2: For E[tr(W^{-1})^2*tr(W^{-3})W^{-2}] with W ~ W_m^1(n,S), where n and S are defined below:
n <- 20
S <- matrix(c(25, 49,
              49, 109), nrow=2, ncol=2)
iwishmom(n, S, c(2, 0, 1), 2, 2) # iw = 2, for real Wishart distribution
#>               [,1]          [,2]
#> [1,]  1.328434e-10 -6.101692e-11
#> [2,] -6.101692e-11  2.824292e-11

# Example 3: For E[tr(W^{-1})^2*tr(W^{-3})] with W ~ W_m^2(n,S), where n and S are defined below:
n <- 20
S <- matrix(c(25, 49 + 2i,
              49 - 2i, 109), nrow=2, ncol=2)
iwishmom(n, S, c(2, 0, 1), 0, 1) # iw = 0, for complex Wishart distribution
#> [1] 1.17985e-08

# Example 4: For E[tr(W^{-1})*tr(W^{-2})^2*tr(W^{-3})^2*W^{-1}] with W ~ W_m^2(n,S), where n and S are defined below:
n <- 30
S <- matrix(c(25, 49 + 2i,
              49 -2i, 109), nrow=2, ncol=2)
iwishmom(n, S, c(1, 2, 2), 1, 1) # iw = 1, for complex Wishart distribution
#>                             [,1]                        [,2]
#> [1,]  1.348928e-21+0.000000e+00i -6.116211e-22+2.496413e-23i
#> [2,] -6.116211e-22-2.496413e-23i  3.004350e-22+0.000000e+00i

iwishmom_sym()

The function iwishmom_sym() generates an analytical expression of $\mathbb{E}\left[\prod_{j=1}^r \mbox{tr}(W^{-j})^{f_j}W^{-iw}\right]$, where W ∼ Wmβ(n, Σ). When iw = 0, it generates an analytical expression of $\mathbb{E}[\prod_{j=1}^r \mbox{tr}(W^{-j})^{f_j}]$.

Arguments

  • f: a vector of nonnegative integers fj that represents the power for tr(W−j), j = 1, …, r
  • iw: Power of W−1
  • alpha: The type of Wishart distribution (α = 2/β):
    • 1/2: Quaternion Wishart
    • 1: Complex Wishart
    • 2: Real Wishart (default)
  • latex: The type of output
    • TRUE: LaTeX expression
    • FALSE: Dataframe (default)

Output

When iw = 0, the output is an analytical expression of $\mathbb{E}[\prod_{j=1}^r \mbox{tr}(W^{-j})^{f_j}]$. When iw ≠ 0, the output is an analytical expression of $\mathbb{E}[\prod_{j=1}^r \mbox{tr}(W^{-j})^{f_j}W^{-iw}]$. If latex = FALSE (default), the output is a data frame that stores the coefficients for the analytical expression. If latex = TRUE, the output is a $\LaTeX$ formatted string of the result in terms of ñ and Σ, where ñ = n − m + 1 − α and m is the dimension of the β-Wishart distribution.

Examples

# Example 1: For E[tr(W^{-1})^4] with W ~ W_m^1(n,Sigma), represented as a dataframe:
iwishmom_sym(4) # iw = 0, for real Wishart distribution
#> $dataframe
#>     kappa      h_kappa_numerator
#> 1       4              240n1-288
#> 2     3,1           64n1^2-256n1
#> 3     2,2        12n1^2-60n1+216
#> 4   2,1,1  12n1^3-72n1^2+36n1+72
#> 5 1,1,1,1 n1^4-7n1^3+n1^2+35n1-6
#> 
#> $denominator
#> [1] "n1^8-7n1^7-11n1^6+107n1^5+34n1^4-388n1^3-24n1^2+288n1"

# Example 2: For E[tr(W^{-1})*tr(W^{-2})W^{-1}] with W ~ W_m^1(n,Sigma), represented as a dataframe:
iwishmom_sym(c(1,1), 1) # iw = 1, for real Wishart distribution
#> $dataframe
#>   i   rho          c_numerator
#> 1 1     3 4n1^3-16n1^2-20n1+24
#> 2 1   2,1 n1^4-6n1^3+3n1^2+6n1
#> 3 1 1,1,1     n1^3-6n1^2+3n1+6
#> 4 2     2    2n1^3-10n1^2+36n1
#> 5 2   1,1       10n1^2-42n1+36
#> 6 3     1 4n1^3-16n1^2+60n1-72
#> 7 4  <NA>          40n1^2-48n1
#> 
#> $denominator
#> [1] "n1^8-7n1^7-11n1^6+107n1^5+34n1^4-388n1^3-24n1^2+288n1"

# Example 3: For E[tr(W^{-1})^4] with W ~ W_m^2(n,Sigma), represented as a LaTeX string:
writeLines(iwishmom_sym(4, 0, 1, latex=TRUE)) # iw = 0, for complex Wishart distribution
#> [(30\tilde{n})p_{(4)}(\Sigma^{-1})
#> +(16\tilde{n}^2-24)p_{(3,1)}(\Sigma^{-1})
#> +(3\tilde{n}^2+18)p_{(2,2)}(\Sigma^{-1})
#> +(6\tilde{n}^3-24\tilde{n})p_{(2,1,1)}(\Sigma^{-1})
#> +(\tilde{n}^4-8\tilde{n}^2+6)p_{(1,1,1,1)}(\Sigma^{-1})]
#> /(\tilde{n}^8-14\tilde{n}^6+49\tilde{n}^4-36\tilde{n}^2)

# Example 4: For E[tr(W^{-1})*tr(W^{-2})W^{-1}] with W ~ W_m^2(n,Sigma), represented as a LaTeX string:
writeLines(iwishmom_sym(c(1, 1), 1, 1, latex=TRUE)) # iw = 1, for complex Wishart distribution
#> [[(2\tilde{n}^3-8\tilde{n})p_{(3)}(\Sigma^{-1})+(\tilde{n}^4-4\tilde{n}^2)p_{(2,1)}(\Sigma^{-1})+(\tilde{n}^3-4\tilde{n})p_{(1,1,1)}(\Sigma^{-1})]\Sigma^{-1}
#> +[(\tilde{n}^3+6\tilde{n})p_{(2)}(\Sigma^{-1})+(5\tilde{n}^2)p_{(1,1)}(\Sigma^{-1})]\Sigma^{-2}
#> +(2\tilde{n}^3+12\tilde{n})p_{(1)}(\Sigma^{-1})\Sigma^{-3}
#> +(10\tilde{n}^2)\Sigma^{-4}]
#> /(\tilde{n}^8-14\tilde{n}^6+49\tilde{n}^4-36\tilde{n}^2)


Auxiliary Functions

Below is a list of auxiliary functions that are called by wishmom, iwishmom, wishmom_sym, and iwishmom_sym.

ip_desc()

The function ip_desc() generates all integer partitions of a given integer k in a reverse lexicographical order.

Arguments

  • k: A positive integer to be partitioned

Output

A matrix where each row represents an integer partition of k, listed in a reverse lexicographical order.

Examples

# Example 1: List of integer partitions of 3
ip_desc(3)
#>      [,1] [,2] [,3]
#> [1,]    3    0    0
#> [2,]    2    1    0
#> [3,]    1    1    1

# Example 2: List of integer partitions of 5
ip_desc(5)
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    5    0    0    0    0
#> [2,]    4    1    0    0    0
#> [3,]    3    2    0    0    0
#> [4,]    3    1    1    0    0
#> [5,]    2    2    1    0    0
#> [6,]    2    1    1    1    0
#> [7,]    1    1    1    1    1

dkmap()

The function dkmap() computes the mapping matrix Dk discussed in Appendix B of Hillier and Kan (2024), modified for the general β-Wishart case. The returned matrix is Dk but with n in the diagonal elements removed.

Arguments

  • k: The order of the mapping matrix Dk (a positive integer)
  • alpha: The type of β-Wishart distribution (α = 2/β):
    • 1/2: Quaternion Wishart
    • 1: Complex Wishart
    • 2: Real Wishart (default)

Output

The mapping matrix Dk but with n removed from its diagonal.

Examples

# Example 1: Compute the mapping matrix for k = 2, real Wishart
dkmap(2)
#>      [,1] [,2] [,3] [,4]
#> [1,]    2    1    1    0
#> [2,]    2    1    0    1
#> [3,]    4    0    0    0
#> [4,]    0    4    0    0

# Example 2: Compute the mapping matrix for k = 1, complex Wishart
dkmap(1, 1)
#>      [,1] [,2]
#> [1,]    0    1
#> [2,]    1    0

# Example 3: Compute the mapping matrix for k = 2, quaternion Wishart
dkmap(2, 1/2)
#>      [,1] [,2] [,3] [,4]
#> [1,] -1.0  1.0    1    0
#> [2,]  0.5 -0.5    0    1
#> [3,]  1.0  0.0    0    0
#> [4,]  0.0  1.0    0    0

denpoly()

The function denpoly() computes the coefficients of the denominator polynomial for the elements $\tilde{\mathcal{H}}_k$ and $\tilde{\mathcal{C}}_k$. The function returns a vector containing the coefficients in descending powers of ñ, with the last element being the coefficient of ñ.

Arguments

  • k: The order of the polynomial (a positive integer)
  • alpha: The type of β-Wishart distribution (α = 2/β):
    • 1/2: Quaternion Wishart
    • 1: Complex Wishart
    • 2: Real Wishart (default)

Output

A vector containing the coefficients of the denominator polynomial in descending powers of ñ for the elements of ℋ̃k and $\tilde{\mathcal{C}}_k$.

Examples

# Example 1: Compute the denominator polynomial for k = 3 and alpha = 2
# Output corresponds to the polynomial n1^5-3n1^4-8n1^3+12n1^2+16n1,
# where n1 is \eqn{\tilde{n}}
denpoly(3)
#> [1]  1 -3 -8 12 16

# Example 2: Compute the denominator polynomial for k = 2 and alpha = 1
# Output corresponds to the polynomial n1^3-n1, where n1 is \eqn{\tilde{n}}
denpoly(2, alpha = 1)
#> [1]  1  0 -1

qk_coeff()

The function qk_coeff() computes the coefficient matrix 𝒞k, which is obtained based on Corollary 1 of Hillier and Kan (2024), after a modification for the general β-Wishart case. 𝒞k is represented as a 3-dimensional array where each slice along the third dimension represents a coefficient matrix of the polynomial in descending powers of n.

Arguments

  • k: The order of the 𝒞k matrix
  • alpha: The type of Wishart distribution (α = 2/β):
    • 1/2: Quaternion Wishart
    • 1: Complex Wishart
    • 2: Real Wishart (default)

Output

A 3-dimensional array representing 𝒞k, a matrix of constants that allow us to obtain 𝔼[pλ(W)Wr], where r + |λ| = k and W ∼ Wmβ(n, Σ).

Examples

# Example 1:
qk_coeff(2) # For real Wishart distribution with k = 2
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2]
#> [1,]    1    1
#> [2,]    2    0

# Example 2:
qk_coeff(3, 1) # For complex Wishart distribution with k = 3
#> , , 1
#> 
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    0    0    0
#> [2,]    0    1    0    0
#> [3,]    0    0    1    0
#> [4,]    0    0    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2] [,3] [,4]
#> [1,]    0    2    1    0
#> [2,]    2    0    0    1
#> [3,]    2    0    0    1
#> [4,]    0    2    1    0
#> 
#> , , 3
#> 
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    0    0    1
#> [2,]    0    1    1    0
#> [3,]    0    2    0    0
#> [4,]    2    0    0    0

# Example 3:
qk_coeff(2, 1/2) # For quaternion Wishart distribution with k = 2
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2]
#> [1,] -0.5    1
#> [2,]  0.5    0

wish_ps()

The function wish_ps() computes the coefficient matrix ℋk that allows us to compute 𝔼[pκ(W)], which is obtained based on Proposition 5 of Hillier and Kan (2024), after a modification for the general β-Wishart case. ℋk is represented as a 3-dimensional array where each slice along the third dimension represents a coefficient matrix of the polynomial in descending powers of n.

Arguments

  • k: The order of the ℋk matrix
  • alpha: The type of Wishart distribution (α = 2/β):
    • 1/2: Quaternion Wishart
    • 1: Complex Wishart
    • 2: Real Wishart (default)

Output

A 3-dimensional array representing ℋk, a matrix of constants that allows us to obtain 𝔼[pκ(W)], where |κ| = k and W ∼ Wmβ(n, Σ).

Examples

# Example 1:
wish_ps(3) # For real Wishart distribution with k = 3
#> , , 1
#> 
#>      [,1] [,2] [,3]
#> [1,]    1    0    0
#> [2,]    0    1    0
#> [3,]    0    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2] [,3]
#> [1,]    3    3    0
#> [2,]    4    1    1
#> [3,]    0    6    0
#> 
#> , , 3
#> 
#>      [,1] [,2] [,3]
#> [1,]    4    3    1
#> [2,]    4    4    0
#> [3,]    8    0    0

# Example 2:
wish_ps(4, 1) # For complex Wishart distribution with k = 4
#> , , 1
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    1    0    0    0    0
#> [2,]    0    1    0    0    0
#> [3,]    0    0    1    0    0
#> [4,]    0    0    0    1    0
#> [5,]    0    0    0    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0    4    2    0    0
#> [2,]    3    0    0    3    0
#> [3,]    4    0    0    2    0
#> [4,]    0    4    1    0    1
#> [5,]    0    0    0    6    0
#> 
#> , , 3
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    5    0    0    6    0
#> [2,]    0    7    3    0    1
#> [3,]    0    8    2    0    1
#> [4,]    6    0    0    5    0
#> [5,]    0    8    3    0    0
#> 
#> , , 4
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0    4    1    0    1
#> [2,]    3    0    0    3    0
#> [3,]    2    0    0    4    0
#> [4,]    0    4    2    0    0
#> [5,]    6    0    0    0    0

# Example 3:
wish_ps(2, 1/2) # For quaternion Wishart distribution with k = 2
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2]
#> [1,] -0.5    1
#> [2,]  0.5    0

qkn_coeff()

The function qkn_coeff() computes the inverse of the coefficient matrix $\tilde{\mathcal{C}}_k$, which is obtained based on Corollary 2 of Hillier and Kan (2024), after a modification for the general β-Wishart case. $\tilde{\mathcal{C}}_k^{-1}$ is represented as a 3-dimensional array where each slice along the third dimension represents a coefficient matrix of the polynomial in descending powers of ñ.

Arguments

  • k: The order of the $\tilde{\mathcal{C}}_k$ matrix
  • alpha: The type of Wishart distribution (α = 2/β):
    • 1/2: Quaternion Wishart
    • 1: Complex Wishart
    • 2: Real Wishart (default)

Output

A 3-dimensional array representing $\tilde{\mathcal{C}}_k^{-1}$, a matrix of constants that allow us to obtain 𝔼[pλ(W−1)W−r], where r + |λ| = k and W ∼ Wmβ(n, Σ).

Examples

# Example 1:
qkn_coeff(2) # For real Wishart distribution with k = 2
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2]
#> [1,]   -1   -1
#> [2,]   -2    0

# Example 2:
qkn_coeff(3, 1) # For complex Wishart distribution with k = 3
#> , , 1
#> 
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    0    0    0
#> [2,]    0    1    0    0
#> [3,]    0    0    1    0
#> [4,]    0    0    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2] [,3] [,4]
#> [1,]    0   -2   -1    0
#> [2,]   -2    0    0   -1
#> [3,]   -2    0    0   -1
#> [4,]    0   -2   -1    0
#> 
#> , , 3
#> 
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    0    0    1
#> [2,]    0    1    1    0
#> [3,]    0    2    0    0
#> [4,]    2    0    0    0

# Example 3:
qkn_coeff(2, 1/2) # For quaternion Wishart distribution with k = 2
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2]
#> [1,]  0.5   -1
#> [2,] -0.5    0

iwish_ps()

The function iwish_ps() computes the inverse of the coefficient matrix $\tilde{\mathcal{H}}_k$ that allows us to compute 𝔼[pκ(W−1)], which is obtained based on Eq.(82) of Hillier and Kan (2024), after a modification for the general β-Wishart case. $\tilde{\mathcal{H}}_k^{-1}$ is represented as a 3-dimensional array where each slice along the third dimension represents a coefficient matrix of the polynomial in descending powers of ñ.

Arguments

  • k: The order of the $\tilde{\mathcal{H}}_k$ matrix
  • alpha: The type of Wishart distribution (α = 2/β):
    • 1/2: Quaternion Wishart
    • 1: Complex Wishart
    • 2: Real Wishart (default)

Output

A 3-dimensional array representing $\tilde{\mathcal{H}}_k^{-1}$, a matrix of constants that allows us to obtain 𝔼[pκ(W−1)], where |κ| = k and W ∼ Wmβ(n, Σ).

Examples

# Example 1:
iwish_ps(3) # For real Wishart distribution with k = 3
#> , , 1
#> 
#>      [,1] [,2] [,3]
#> [1,]    1    0    0
#> [2,]    0    1    0
#> [3,]    0    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2] [,3]
#> [1,]   -3   -3    0
#> [2,]   -4   -1   -1
#> [3,]    0   -6    0
#> 
#> , , 3
#> 
#>      [,1] [,2] [,3]
#> [1,]    4    3    1
#> [2,]    4    4    0
#> [3,]    8    0    0

# Example 2:
iwish_ps(4, 1) # For complex Wishart distribution with k = 4
#> , , 1
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    1    0    0    0    0
#> [2,]    0    1    0    0    0
#> [3,]    0    0    1    0    0
#> [4,]    0    0    0    1    0
#> [5,]    0    0    0    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0   -4   -2    0    0
#> [2,]   -3    0    0   -3    0
#> [3,]   -4    0    0   -2    0
#> [4,]    0   -4   -1    0   -1
#> [5,]    0    0    0   -6    0
#> 
#> , , 3
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    5    0    0    6    0
#> [2,]    0    7    3    0    1
#> [3,]    0    8    2    0    1
#> [4,]    6    0    0    5    0
#> [5,]    0    8    3    0    0
#> 
#> , , 4
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0   -4   -1    0   -1
#> [2,]   -3    0    0   -3    0
#> [3,]   -2    0    0   -4    0
#> [4,]    0   -4   -2    0    0
#> [5,]   -6    0    0    0    0

# Example 3:
iwish_ps(2, 1/2) # For quaternion Wishart distribution with k = 2
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2]
#> [1,]  0.5   -1
#> [2,] -0.5    0

qkn_coeffr()

The function qkn_coeffr() computes the coefficient matrix $\tilde{\mathcal{C}}_k$ for the general β-Wishart case. Elements of $\tilde{\mathcal{C}}_k$ are rational polynomials of ñ. The output contains two components: c and den. c is a 3-dimensional array where each slice along the third dimension represents a coefficient matrix of the numerator polynomial in descending powers of ñ, and den is a vector that represents the coefficients of the denominator polynomial in descending power of ñ.

Arguments

  • k: The order of the $\tilde{\mathcal{C}}_k$ matrix
  • alpha: The type of Wishart distribution (α = 2/β):
    • 1/2: Quaternion Wishart
    • 1: Complex Wishart
    • 2: Real Wishart (default)

Output

The output has two components: c and den. c is a 3-dimensional array representing the numerator polynomial of $\tilde{\mathcal{C}}_k$, and den is a vector representing the denominator polynomial of $\tilde{\mathcal{C}}_k$, where $\tilde{\mathcal{C}}_k$ is a matrix of constants that allow us to obtain 𝔼[pλ(W−1)W−r], where r + |λ| = k and W ∼ Wmβ(n, Σ).

Examples

# Example 1:
qkn_coeffr(2) # For real Wishart distribution with k = 2
#> $c
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2]
#> [1,]    0    1
#> [2,]    2   -1
#> 
#> 
#> $den
#> [1]  1 -1 -2
 
# Example 2:
qkn_coeffr(3, 1) # For complex Wishart distribution with k = 3
#> $c
#> , , 1
#> 
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    0    0    0
#> [2,]    0    1    0    0
#> [3,]    0    0    1    0
#> [4,]    0    0    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2] [,3] [,4]
#> [1,]    0    2    1    0
#> [2,]    2    0    0    1
#> [3,]    2    0    0    1
#> [4,]    0    2    1    0
#> 
#> , , 3
#> 
#>      [,1] [,2] [,3] [,4]
#> [1,]    0    0    0    2
#> [2,]    0    0    2    0
#> [3,]    0    4   -2    0
#> [4,]    4    0    0   -2
#> 
#> 
#> $den
#> [1]  1  0 -5  0  4

# Example 3:
qkn_coeffr(2, 1/2) # For quaternion Wishart distribution with k = 2
#> $c
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2]
#> [1,]  0.0  1.0
#> [2,]  0.5  0.5
#> 
#> 
#> $den
#> [1]  1.0  0.5 -0.5

iwish_psr()

The function iwish_psr() computes the coefficient matrix $\tilde{\mathcal{H}}_k$ for the general β-Wishart case. Elements of $\tilde{\mathcal{H}}_k$ are rational polynomials of ñ. The output contains two components: c and den. c is a 3-dimensional array where each slice along the third dimension represents a coefficient matrix of the numerator polynomial in descending powers of ñ, and den is a vector that represents the coefficients of the denominator polynomial in descending power of ñ.

Arguments

  • k: The order of the $\tilde{\mathcal{C}}_k$ matrix
  • alpha: The type of Wishart distribution (α = 2/β):
    • 1/2: Quaternion Wishart
    • 1: Complex Wishart
    • 2: Real Wishart (default)

Output

The output has two components: c and den. c is a 3-dimensional array representing the numerator polynomial of $\tilde{\mathcal{H}}_k$, and den is a vector representing the denominator polynomial of $\tilde{\mathcal{H}}_k$, where $\tilde{\mathcal{H}}_k$ is a matrix of constants that allow us to obtain 𝔼[pλ(W−1)], where |λ| = k and W ∼ Wmβ(n, Σ).

Examples

# Example 1:
iwsih_psr(2) # For real Wishart distribution with k = 2
#> $c
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2]
#> [1,]    0    1
#> [2,]    2   -1
#> 
#> 
#> $den
#> [1]  1 -1 -2
 
# Example 2:
iwish_psr(4, 1) # For complex Wishart distribution with k = 4
#> $c
#> , , 1
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    1    0    0    0    0
#> [2,]    0    1    0    0    0
#> [3,]    0    0    1    0    0
#> [4,]    0    0    0    1    0
#> [5,]    0    0    0    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0    4    2    0    0
#> [2,]    3    0    0    3    0
#> [3,]    4    0    0    2    0
#> [4,]    0    4    1    0    1
#> [5,]    0    0    0    6    0
#> 
#> , , 3
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    1    0    0   10    0
#> [2,]    0    3    6    0    2
#> [3,]    0   16   -6    0    1
#> [4,]   10    0    0    1    0
#> [5,]    0   16    3    0   -8
#> 
#> , , 4
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0    4   -3    0    5
#> [2,]    3    0    0    3    0
#> [3,]   -6    0    0   12    0
#> [4,]    0    4    6    0   -4
#> [5,]   30    0    0  -24    0
#> 
#> , , 5
#> 
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0    0    0    0    0
#> [2,]    0   12   -9    0   -3
#> [3,]    0  -24   18    0    6
#> [4,]    0    0    0    0    0
#> [5,]    0  -24   18    0    6
#> 
#> 
#> $den
#> [1]   1   0 -14   0  49   0 -36   0

# Example 3:
iwish_psr(2, 1/2) # For quaternion Wishart distribution with k = 2
#> $c
#> , , 1
#> 
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> 
#> , , 2
#> 
#>      [,1] [,2]
#> [1,]  0.0  1.0
#> [2,]  0.5  0.5
#> 
#> 
#> $den
#> [1]  1.0  0.5 -0.5


References

Díaz-García, José and Gutiérrez-Jáimez, Ramón (2011). On Wishart distribution: som extension. Linear Algebra and its Applications, 435, 1296-1310.

Drensky, Vesselin, Edelman, Alan, Genoar, Tierney, Kan, Raymond, and Koev, Plamen (2021). The Densities and Distributions of the Largest Eigenvalue and the Trace of a Beta-Wishart Matrix. Random Matrices: Theory and Applications, 10(1).

Letac, Gérard, and Massam, Héelène (2004). All invariant moments of the Wishart distribution. Scandinavian Journal of Statistics, 31, 295-318.

Hillier, Grant, and Kan, Raymond (2024). On the expectations of equivariant matrix-valued functions of Wishart and inverse Wishart Matrices. Scandinavian Journal of Statistics, 51, 697-723.