cf_LogRV_WilksLambdaNC: Characteristic function of a linear combination of...

In gajdosandrej/CharFunToolR: Numerical Computation Cumulative Distribution Function and Probability Density Function from Characteristic Function

Characteristic function of a linear combination of independent LOG-TRANSFORMED WILK's LAMBDA distributed random variables

Description

cf_LogRV_WilksLambdaNC(t, p, n, q, delta, coef, MAX) evaluates characteristic function of a linear combination (resp. convolution) of independent LOG-TRANSFORMED WILK's LAMBDA distributed random variables, with non-central distributions specified by the parameters p, n, q, the non-centrality parameters delta, and the coefficients coef.

That is, cf_LogRV_WilksLambdaNC evaluates the characteristic function of a random variable Y = coef_1*W_1 +...+ coef_N*W_N, such that cf_Y(t) = cf_W_1(coef_1*t) *...* cf_W_N(coef_N*t), where cf_W_i(t) is CF of W_i = log(Lambda_i), and each Lambda_i has non-central WILK's LAMBDA distribution, Lambda_i ~ Lambda(p_i,m_i,n_i,delta_i), for i = 1,...,N.

In particular, \Lambda_i = det(E_i)/det(E_i + H_i), with E_i and H _i being independent random matrices with Wishart distributions, E_i ~ Wp_i(m_i,\Sigma_i) with central Wishart distribution, and H_i ~ Wp_i(n_i,\Sigma_i,delta) with non-central Wishart distribution, with n_i >= p_i and \Sigma_i > 0 (an unknown positive definite symmetric covariance matrix), for all i = 1,...,N.

The non-central distribution of MINUS LOG-TRANSFORMED WILK's LAMBDA STATISTIC, say L ~ Lambda(p,n,q,delta), with L in (0,1), is specified by its characteristic function

cf_{-log(L)}(t) = cf_LogRV_Beta(-t,n/2,q/2) * ... * cf_LogRV_Beta(-t,(n+1-p)/2,q/2) * Hypergeom1F1Mat(-i*t,-i*t+(n+q)/2,-delta/2)

, i.e. the characteristic function of the non-central distribution differs from the central distribution by a factor specified by the confluent hypergeometric function 1F1(a;b;X) of a matrix argument X, with the complex (vector) parameters specified as a = -i*t, b = -i*t+(n+q)/2, and the matrix argument X = -delta/2, where delta is the non-centrality (matrix) parameter of the distribution.

Usage

cf_LogRV_WilksLambdaNC(t, p, n, q, delta, coef, MAX)

Arguments

t	vector or array of real values, where the CF is evaluated.
p	vector of the dimension parameters p = (p_1,...,p_N). If empty, default value is p = (1,...,1).
n	vector of degrees of freedom of the Wishart matrices E_i, n = (n_1,...,n_N). If empty, default value is n = (1,...,1).
q	vector of degrees of freedom of the Wishart matrices H_i, q = (q_1,...,q_N). If empty, default value is q = (1,...,1).
delta	p.s.d. matrix or vector of nonnegative eigenvalues or array of matrices or vectors of eigenvalues of the non-centrality parameters, delta = (delta_1,...,delta_N). Default value is an empty array.
coef	vector of the coefficients of the linear combination of the log-transformed random variables. If coef is scalar, it is assumed that all coefficients are equal. If empty, default value is coef = 1.
MAX	the maximum number of partitions used for computing the hypergeometric 1F1 function with matrix argument, for more details see HypergeompFqMat. If MAX is empty, default value is set to MAX = 0. If MAX = 0, then the algorithm uses the fast (approximate) method for computing the confluent hypergeometric function 1F1(a;b;X), see Hypergeom1F1MatApprox, otherwise the algorithm uses the algorithm Hypergeom1F1Mat based on computing the truncated expansion series of 1F1(a;b;X) with MAX number of partitions.

Details

Computing CF of the LOG-TRANSFORMED NON-CENTARL WILK's LAMBDA random variable depends on computing the generalized hypergeometric function of matrix argument. cf_LogRV_WilksLambdaNC uses modified implementation for computing the truncated hypergeometric function, see the algorithm HypergeompFqMat, as originaly suggested in Koev and Edelman (2006). The truncation of the hypergeometric series is controled by the parameter MAX (start with MAX = 20). If necessary, the parameter MAX should be set to sufficiently large value but this can be numerically intractable. If MAX = 0 (now it is the default value), the algorithm uses the fast (approximate) method for computing the confluent hypergeometric function 1F1(a;b;X), see Hypergeom1F1MatApprox, otherwise the algorithm uses the algorithm Hypergeom1F1Mat based on computing the truncated expansion series of 1F1(a;b;X) with MAX number of partitions. The approximate alternative method (if MAX = 0) is based on heuristic arguments (not formally proven), and the value of the confluent hypergeometric function 1F1(a,b,X) of a matrix argument is calculated (approximately) as 1F1(a;b;X) ~ 1F1(a;b;x(1)) * ... * 1F1(a;b;x(p)), where 1F1(a;b;x(1)) is the scalar value confluent hypergeometric function 1F1(a,b,x(i)) with (x[1],...,x[p]) = eigen(X).

By default (or if MAX is set to value MAX = 0), cf_LogRV_WilksLambdaNC uses this alternative method for computing 1F1(a;b;X).

Value

Characteristic function cf(t) of a linear combination of independent LOG-TRANSFORMED WILK's LAMBDA distributed random variables.

Note

Ver.: 05-Oct-2018 11:46:29 (consistent with Matlab CharFunTool v1.3.0, 10-Aug-2018 15:46:49).

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] Witkovsky, V., 2018. Exact distribution of selected multivariate test criteria by numerical inversion of their characteristic functions. arXiv preprint arXiv:1801.02248.

See Also

For more details see WIKIPEDIA: https://en.wikipedia.org/wiki/Wilks

Other Continuous Probability Distribution: cfS_Arcsine(), cfS_Beta(), cfS_Gaussian(), cfS_Laplace(), cfS_Rectangular(), cfS_Student(), cfS_TSP(), cfS_Trapezoidal(), cfS_Triangular(), cfS_Wigner(), cfX_ChiSquare(), cfX_Exponential(), cfX_FisherSnedecor(), cfX_Gamma(), cfX_InverseGamma(), cfX_LogNormal(), cf_ArcsineSymmetric(), cf_BetaNC(), cf_BetaSymmetric(), cf_Beta(), cf_ChiSquare(), cf_Exponential(), cf_FisherSnedecorNC(), cf_FisherSnedecor(), cf_Gamma(), cf_InverseGamma(), cf_Laplace(), cf_LogRV_BetaNC(), cf_LogRV_Beta(), cf_LogRV_ChiSquareNC(), cf_LogRV_ChiSquare(), cf_LogRV_FisherSnedecorNC(), cf_LogRV_FisherSnedecor(), cf_LogRV_MeansRatioW(), cf_LogRV_MeansRatio(), cf_LogRV_WilksLambda(), cf_Normal(), cf_RectangularSymmetric(), cf_Student(), cf_TSPSymmetric(), cf_TrapezoidalSymmetric(), cf_TriangularSymmetric(), cf_vonMises()

Examples

## EXAMPLE 1
# CF of log of noncentral Wilks Lambda RV distribution Lambda(p,n,q,delta)
p <- 5
n <- 10 # d.f. of within SS&P
q <- 3  # d.f. of between SS&P
delta <- sort(runif(p))
coef <- 1
t <- seq(-10, 10, length.out = 201)
# MAX <- 0
cf1 <- function(t) cf_LogRV_WilksLambdaNC(t, p, n, q, delta, coef, 0)
# MAX <- 20
cf2 <- function(t) cf_LogRV_WilksLambdaNC(t, p, n, q, delta, coef, 20)
plotReIm2(list(cf1, cf2), list(t, t), title = 'CF of log of non-central Wilks Lambda RV')

## EXAMPLE 2
# CF of a weighted linear combination of minus log Wilks Lambda RVs
p <- c(5, 5, 5)
n <- c(10, 15, 20)
q <- c(3, 2, 1)
# delta <- c(sort(runif(p[1])), sort(runif(p[2])), sort(runif(p[3])))
delta <- matrix(list(runif(p[1]), runif(p[2]), runif(p[3])), 1, 3)
coef <- -c(10, 15, 20) / 45
t <- seq(-20, 20, length.out = 201)
plotReIm(function(t) {cf_LogRV_WilksLambdaNC(t, p, n, q, delta, coef)}, t,
         title = 'CF of a weighted linear combination of -log Wilks Lambda RVs')

## EXAMPLE 3
# PDF/CDF of minus log Wilks Lambda RV (p=10, n=30, q=5) from its CF
p <- 10
n <- 30
q <- 5
delta <- c(1, 2, 3, 10, 30)
coef <- -1
cf0 <- function(t) cf_LogRV_WilksLambdaNC(t, p, n, q, coef = coef)
cf <- function(t) cf_LogRV_WilksLambdaNC(t, p, n, q, delta, coef)
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)
matplot(cbind(result0$x, result$x), cbind(result0$cdf, result$cdf),
        xlab = 'x', ylab = 'CDF',
        main = expression(paste('CDFs of -log(',Lambda,') under null and alternative hypothesis')))
matplot(cbind(result0$x, result$x), cbind(result0$pdf, result$pdf),
        xlab = 'x', ylab = 'CDF',
        main = expression(paste('PDFs of -log(',Lambda,') under null and alternative hypothesis')))

## EXAMPLE 4 (Compare exact vs. simulated non-central Wilk's distribution)
p <- 10 # p - length of the sample vectors (dimensionality)
n <- 30 # n - sample size / degrees of freedon
q <- 5  # q - degrees of freedom due to hypothesis model
N <- 10000  # N - number of simulation samples
M <- replicate(q, runif(p)) # M - the (true) mean (p x q)-matrix
delta <- eigen(M %*% t(M))$values # delta - the eigenvalues of non-centrality matrix
L <- rep(0, N)
for(i in 1:N) {
        X <- replicate(n, rnorm(p))
        E <- X %*% t(X)
        Y <- replicate(q, rnorm(p)) + M
        H <- Y %*% t(Y)
        L[i] <-  -log(det(E)/det(E + H))
}
# Exact and the empirical CDF of -log(L)
cf <- function(t) cf_LogRV_WilksLambdaNC(t, p, n, q, delta, -1)
options <- list()
options$xMin <- 0
options$SixSigmaRule <- 6
prob <- c(0.9, 0.95, 0.99)
result <- cf2DistGP(cf, prob = prob, options = options)
Fn <- ecdf(L)
matplot(cbind(result$x, result$x), cbind(result$cdf, Fn(result$x)),
        xlab = "x", ylab = "CDF / ECDF",
main = "Exact vs. empirical CDF of the non-central distribution")

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