Random variate of a singular Wishart distribution with non-integral degrees of freedom

Question

Let's say that I have a random variable $X$ that follows a singular Wishart distribution with $\nu$ degrees of freedom and a shape matrix $\Sigma$ that is $p\times p$. We can write "symbolically": $$X\sim W_p(\nu,\Sigma)$$ where $\nu\le p-1$ for the singular case.

If $\nu$ is integral then I can get a random variate by simply taking $\nu$ independent random normal variables $Y_i=\mathcal N(0,\Sigma)$, where $\mathcal N(0,\Sigma)$ is the normal distribution with no mean and a covariance of $\Sigma$. Then my random variate will be: $X=\sum_{i=1}^\nu Y_iY_i'$. If $\nu$ is less or equal to $p-1$, this random variate will be non-invertible with probability of one.

Let's say that I would like to get a random variate for a non-integral value of $\nu$, for example $\nu=\frac12$. In this situation, I obviously can't use the summation of normal variables. What would be the appropriate way to get a random variate for such value of $\nu$?

Also, I am not sure it makes perfect sense to have a random variate for $\nu=\frac12$ as if I sum two independent variates of this kind then I can't ensure that the rank of my random variate will be $1$ as what it should be normally ($1=\nu+\nu$).

Answer 1

The Wishart Distribution is defined on the manifold $\mathcal{M}(p)$ of all positive-definite symmetric (psd) $p\times p$ matrices. In the coordinate system $(x_{ij}, 1\le j \le i \le p)$ (which identifies this space of matrices with a subset of $\mathbb{R}^{p(p+1)/2}$) the density of the "standard" Wishart distribution with $n\gt p-1$ "degrees of freedom" at a matrix $x$ is given by

$$f_n(x) = C_n \det(x)^{(n-p-1)/2}\, e^{-\operatorname{tr}(x)/2}$$

where $C_n$ is the normalizing constant. Given any specified positive-definite $p\times p$ matrix $V,$ another Wishart distribution can be generated upon multiplying a standard Wishart $x$ to produce $Vx.$ (This is the multidimensional analog of a change of scale: see "Higher Dimensional Geometries" at https://stats.stackexchange.com/a/564628/919.)

One important property of this distribution is that the diagonal elements of $x$ will all have scaled chi-squared distributions with $n$ degrees of freedom.

The Bartlett Decomposition, as described in Wikipedia, generates $x$ in the form $L^\prime AA^\prime L$ where $V = L^\prime L$ is the Cholesky decomposition of $V$ ($L$ is lower triangular), the diagonal of $A$ is a sequence of $\chi(n), \chi(n-1), \chi(n-p+1)$ variables (that is, square roots of chi-squared variables), and below the diagonal the components are standard Normal; all $p(p+1)/2$ of these variables are independent. This gives a practical and fairly simple way to generate Wishart variables for any $V$ and $n$ that are mathematically possible. This R function gives one implementation directly modeled on the Wikipedia article (using its notation).

#
# Random Wishart matrices.
# Returns a rank-3 array indexed by the third variable.
#
rWishart <- function(n, V, df) {
  d <- dim(V)[1]
  L <- chol(V)
  X <- replicate(n, {
    A <- diag(sqrt(rchisq(d, df + 1 - seq_len(d))), d)
    A[lower.tri(A)] <- rnorm(d*(d-1)/2)
    Y <- crossprod(L, A)
    tcrossprod(Y)
  })
  array(X, dim = c(d, d, n))
}

As an example of its use, I set $$V = \pmatrix{2 & -3 \\ -3 & 6}$$ (whence $p=2$ and the implied correlation coefficient is extremely negative), specified the degrees of freedom in a variable df, and generated $20,000$ realizations from this Wishart distribution as

X <- rWishart(2e4, V, df)

The result is a rank-three array. That is, X has three subscripts. Its last subscript determines one $p\times p$ matrix. For instance, the first realization is X[,,1] and the last realization is X[,,2e4].

Here are histograms of the components of these realizations for two different values of $n,$ $1 + 1/\pi \approx 1.32$ (close to the minimum of $1$) and $20$ (fairly large). Because symmetry implies $x_{12}=x_{21}$ for any such realization, only three histograms are needed. Notice the differing value scales among the plots.

On top of the histograms I have plotted the chi-squared distributions (left and middle histograms) and the "variance-Gamma" distribution for the off-diagonal $x_{12}$ component. The agreement is good, indicating the formulas in the Wikipedia article are correct (and that the code has correctly implemented them).

Of course the components of a Wishart random matrix are not independent. Here is a scatterplot matrix (of the first two thousand) of the previous results. I show all four components to demonstrate the realizations are indeed all symmetric. The color coding is the same as before.

The R code to produce the figures might be helpful, so it is reproduced below.

#
# Example.
# This generates one set of realizations for each parameter value in `df,
# storing the results in a list `lX`.
#
set.seed(17)
V <- matrix(c(2,-3,-3,6), 2) # Must be psd
df <- c(1 + 1/pi, 20)
lX <- lapply(df, function(df) rWishart(2e4, V, df = df))
#
# Density of the off-diagonal elements.
# (Works only for 2D matrices!)
#
dvgamma <- function(x, V, df) {
  s <- prod(sqrt(diag(V)))
  rho <- V[1,2] / s
  y <- x / (s * (1 - rho^2))
  besselK(abs(y), (df-1)/2) * exp(rho * y) * abs(x)^((df-1)/2) / 
    (gamma(df/2) * sqrt(2^(df-1) * pi * (1 - rho^2) * s^(df + 1)))
}
#
# Plot histograms of the components.
#
sigma2 <- diag(V)
colors <- hsv(seq(0, 2/3, length.out=3), .25, .8)
par(mfrow=c(2,3))
for (i in seq_along(lX)) {
  X <- lX[[i]]
  n <- df[i]
  hist(X[1,1,], freq=FALSE, breaks=102, col=colors[1], 
       main=bquote(df==.(1 + signif(n-1,2))), xlab=expression(X["1,1"]))
  curve(dchisq(x/sigma2[1], n)/sigma2[1], lwd=2, add=TRUE)
hist(X[2,2,], freq=FALSE, breaks=102, col=colors[2], 
       main=bquote(df==.(1 + signif(n-1,2))), xlab=expression(X["2,2"]))
  curve(dchisq(x/sigma2[2], n)/sigma2[2], lwd=2, add=TRUE)
hist(X[1,2,], freq=FALSE, breaks=102, col=colors[3], 
       main=bquote(df==.(1 + signif(n-1,2))), xlab=expression(X["1,2"]))
  curve(dvgamma(x, V, n), add=TRUE, lwd=2)
}
par(mfrow=c(1,1))

Plot the scatterplot matrix for one of the sets of realizations.

X <- lX[[1]]
panel.hist <- function(x, col) { # Modified from the man page for pairs
  COLOR <<- COLOR + 1
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(usr[1:2], 0, 1.5) )
  h <- hist(x, plot = FALSE, breaks=20)
  breaks <- h$breaks; nB <- length(breaks)
  y <- h$counts; y <- y/max(y)
  rect(breaks[-nB], 0, breaks[-1], y, col=colors[COLOR])
}
colors <- colors[c(1,2,2,3)]
COLOR <- 0
e <- c(expression(X["1,1"]), expression(X["2,1"]), expression(X["1,2"]), expression(X["2,2"]))
pairs(t(matrix(X, 4))[1:2000, ], col="#00000020", diag.panel = panel.hist, labels=e)

Answer 2

The answer to this question is negative, we must have an integral number of degrees of freedom $\nu$ if $\nu\le p-1$. To prove it, we can look at the singular Wishart, defined as: $$W_p(\nu,\Sigma)=\sum_{i=1}^nY_iY_i'$$ with $Y_i\sim\mathcal N(0,\Sigma)$, and $\Sigma$ a positive definite matrix.

The rank of the random variate will be $\nu$ almost surely if $\nu\le p-1$ and $\nu$ is integral. In order to keep the property that $W_p(\nu_1,\Sigma)+W_p(\nu_2,\Sigma)=W_p(\nu_1+\nu_2,\Sigma)$ for all $\nu_1$ and $\nu_2$, we must ensure that $\nu_1$ and $\nu_2$ are integral numbers. This can be seen from the definition of the singular Wishart above.

We can try to prove that there exists a way to define $W_p(\frac12,\Sigma)$ so adding two of these independent random variates will create a matrix of rank $1$ but that requires the two variates to share the same (degenerate) eigenspace. Unless we constraint it explicitly then we can't have this property. So we need to have $\nu$ integral until we reach full rank.