rmatrixbeta {matrixsampling}

Matrix Beta sampler

Description

Samples a matrix Beta (type I) distribution.

Usage

rmatrixbeta(n, p, a, b, Theta1 = NULL, Theta2 = NULL, def = 1,
  checkSymmetry = TRUE)

Arguments

n

sample size, a positive integer

p

dimension, a positive integer

a, b

parameters of the distribution, positive numbers with constraints given in Details

Theta1

numerator noncentrality parameter, a positive semidefinite real matrix of order p; setting it to NULL (default) is equivalent to setting it to the zero matrix

Theta2

denominator noncentrality parameter, a positive semidefinite real matrix of order p; setting it to NULL (default) is equivalent to setting it to the zero matrix

def

1 or 2, the definition used; see Details

checkSymmetry

logical, whether to check the symmetry of Theta1 and Theta2

Details

A Beta random matrix U is defined as follows. Take two independent Wishart random matrices S₁ ~ W_p(2a,I_p,Θ₁) and S₂ ~ W_p(2b,I_p,Θ₂).

• definition 1: U = (S₁+S₂)^-½S₁(S₁+S₂)^-½

• definition 2: U = S₁^½(S₁+S₂)^-1S₁^½

In the central case, the two definitions yield the same distribution. Under definition 2, the Beta distribution is related to the Beta type II distribution by U ~ V(I+V)^-1.

Parameters a and b are positive numbers that satisfy the following constraints:

• if both Theta1 and Theta2 are the null matrix, a+b > (p-1)/2; if a < (p-1)/2, it must be half an integer; if b < (p-1)/2, it must be half an integer

• if Theta1 is not the null matrix, a >= (p-1)/2; if b < (p-1)/2, it must be half an integer

• if Theta2 is not the null matrix, b >= (p-1)/2; if a < (p-1)/2, it must be half an integer

Value

A numeric three-dimensional array; simulations are stacked along the third dimension (see example).

Warning

Definition 2 requires the calculation of the square root of S₁ ~ W_p(2a,I_p,Θ₁) (see Details). While S₁ is always positive semidefinite in theory, it could happen that the simulation of S₁ is not positive semidefinite, especially when a is small. In this case the calculation of the square root will return NaN.

Note

The matrix variate Beta distribution is usually defined only for a > (p-1)/2 and b > (p-1)/2. In this case, a random matrix U generated from this distribution satisfies 0 < U < I. For an half integer a \le (p-1)/2, it satisfies 0 \le U < I and rank(U)=2a. For an half integer b \le (p-1)/2, it satisfies 0 < U \le I and rank(I-U)=2b.

Examples

Bsims <- rmatrixbeta(10000, 3, 1, 1)
dim(Bsims) # 3 3 10000