HypergeompFqMat: Hypergeometric function of a matrix argument
In gajdosandrej/CharFunToolR: Numerical Computation Cumulative Distribution Function and Probability Density Function from Characteristic Function

HypergeompFqMat

R Documentation

Hypergeometric function of a matrix argument

Description

HypergeompFqMat(a, b, x, y, alpha, MAX, lam) computes the truncated hypergeometric function pFq^alpha(a;b;x;y) with parameters a and b and of two matrix arguments, x and y.

Usage

HypergeompFqMat(a, b, x, y, alpha = 2, MAX = 20, lam)

Arguments

`a`	parameter of the hypergeometric function `pFq^alpha(a;b;x;y)`.
`b`	parameter of the hypergeometric function `pFq^alpha(a;b;x;y)`.
`x`	matrix argument (given as vector of eigenvalues).
`y`	optional second matrix argument (vector of eigenvalues).
`alpha`	parameter of the hypergeometric function `pFq^alpha(a;b;x;y)`, default value is `alpha = 2`.
`MAX`	maximum number of partitions, `\|kappa\|<=MAX`, default value is `MAX = 20`.
`lam`	optional parameter, `kappa<=lam`.

Details

Here, hypergeometric function pFq^alpha(a;b;x;y) is defined as

pFq^alpha(a;b;x;y) = pFq^alpha(a1,...,ap;b1,...,bp;x;y) = sum_k sum_kappa [((a1)_kappa^(alpha) ... (a1)_kappa^(alpha)) / (k!(b1)_kappa^(alpha) ... (b1)_kappa^(alpha)] C_kappa^alpha(X),

where kappa = (kappa1,kappa2,...) is a partition of k and (ak)_kappa^(alpha), (bk)_kappa^(alpha) denote the generalized Pochhammer symbols, and C_kappa^alpha(X) is the Jack function.

For statistical applications considered here we need to use the confluent hypergeometric function 1F1(a;b;X) and the Gauss hypergeometric function 2F1(a,b;c;X) in matrix argument. These can be expressed as special cases of the generalized hypergeometric function pFq^alpha(a;b;X), with the parameter alpha = 2 (the case with zonal polynomials).

For more details and definition of the hypergeometric function pFq^alpha(a;b;x;y) with matrix argument see, e.g., Koev and Edelman (2006) or Muirhead (2009).

Value

List including the following items:

`s`	hypergeometric sum, `pFq^alpha(a;b;x;y)`,
`ss`	partial sums.

Note

Ver.: 10-Dec-2018 17:35:21 (consistent with Matlab CharFunTool v1.3.0, Ver.: 23-Oct-2017 11:47:05).

Author(s)

MHG LICENCE:
This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

If you use MHG in any program or publication, please acknowledge its author by adding a reference.

This is R version of modified (by Viktor Witkovsky) version of the original MATLAB code hg.m by Plamen Koev

References

[1] Koev, P. and Edelman, A., 2006. The efficient evaluation of the hypergeometric function of a matrix argument. Mathematics of Computation, 75(254), 833-846.

[2] Muirhead RJ. Aspects of multivariate statistical theory. John Wiley & Sons; 2009 Sep 25.

Examples

## EXAMPLE 1
# Evaluate 2F3^9([3 4];[5 6 7];[1 2];[8,9]) & kappa<=(4,3), |kappa|<=6
a <- t(c(3, 4))
b <- t(c(5, 6, 7))
x <- t(c(1, 2))
y <- t(c(8, 9))
alpha <- 9
MAX <- 6
lam <- c(4, 3)
result <- HypergeompFqMat(a, b, x, y, alpha, MAX, lam)

## EXAMPLE 2
# CF of minus log of noncentral Wilks Lambda RV distribution
## cf_LogRV_WilksLambdaNC is using HypergeompFqMat
p <- 10
m <- 30 # elsewhere it is denoted as n (d.f. of within SS&P)
n <- 5  # elsewhere it is denoted as q (d.f. of between SS&P)
X <- c(0, 0, 0.1, 1, 20)
coef <- -1
MAX <- 25
t <- seq(-10, 10, length.out = 201)
plotReIm(function(t)
        cf_LogRV_WilksLambdaNC(t, p, m, n, X, coef, MAX),
        t,
        title = 'CF of log of noncentral Wilks Lambda RV')

## EXAMPLE 3
# PDF/CDF of minus log Wilks Lambda RV (p=10, m=30, n=5) from its CF
p <- 10
m <- 30
n <- 5
X <- c(0, 0, 0.1, 1, 20)
MAX <- 30
coef <- -1
cf0 <- function(t) cf_LogRV_WilksLambdaNC(t, p, m, n, coef = coef, MAX = MAX)
cf <- function(t) cf_LogRV_WilksLambdaNC(t, p, m, n, X, coef, MAX)
prob <- c(0.9, 0.95, 0.99)
options <- list()
options$xMin <- 0
result0 <- cf2DistGP(cf0, prob = prob, options = options)
result <- cf2DistGP(cf, prob = prob, options = options)
print(result)
plot(result0$x, result0$cdf, type="l", col="orange",
     xlab = 'x', ylab = 'CDF',
     main = expression(paste('CDFs of -log(',Lambda,') under null and alternative hypothesis')))
lines(result$x, result$cdf, col="green")
plot(result0$x, result0$pdf, type="l", col="yellow",
     xlab = 'x', ylab = 'PDF',
     main = expression(paste('PDFs of -log(',Lambda,') under null and alternative hypothesis')))
lines(result$x, result$pdf, col="pink")

## EXAMPLE 4 (!!!LONG CALCULATION of Hypergeom2F1Mat!!!)
# Non-null distribution of log-transformed T statistic, W = -log(T)
# Test statistic for testing equality of 2 covariance matrices
p <- 5          # p - length of the sample vectors (dimensionality)
n1 <- 15        # n1 - sample size from the first normal poulation
n2 <- 20        # n1 - sample size from the second normal poulation
nu1 <- n1 - 1   # nu1 - degrees of freedom
nu2 <- n2 - 1   # nu1 - degrees of freedom
nu <- nu1 + nu2 # nu1 - degrees of freedom
# delta - eigenvalues of the non-centrality matrix Delta
delta <- c(1.1, 1.325, 1.55, 1.775, 2.0)
# Create the characteristic function of null distribution
c <- GammaMultiLog(nu / 2, p) - GammaMultiLog(nu1 / 2, p) - GammaMultiLog(nu2 / 2, p)
cfH0 <- function(t) {exp(c + GammaMultiLog((1 - 1i * t) * nu1 / 2, p) +
GammaMultiLog((1 - 1i * t) * nu2 / 2, p) - GammaMultiLog((1 - 1i * t) * nu / 2, p))}
# Create the characteristic function of non-null distribution
MAX <- 50
cfHA <- function(t) {cfH0(t) * prod(delta)^(-nu1 / 2) *
Hypergeom2F1Mat(nu / 2, (1 - 1i * t) * nu1 / 2, (1 - 1i * t) * nu / 2, 1 - 1 / delta, MAX)[[1]]}
# Evaluate PDF/CDF and selected quantiles of W = -log(T)
x <- seq(50, 90, length.out = 100)
prob <- c(0.9, 0.95, 0.975, 0.99, 0.999)
result <- cf2DistGP(cfHA, x, prob)

gajdosandrej/CharFunToolR documentation built on June 3, 2024, 7:46 p.m.