# How to sample from a Wishart distribution?

From Wikipedia, we know that $n$, the degrees of freedom, should be larger than $p-1$ where $p$ is the dimension of the scale matrix.

Also, from the bottom part of the same article, we see "Bartlett decomposition" to sample Wishart distribution where the diagonal element is generated from $\chi^{2}(n - i + 1)$.

However, if you type the following commands in Matlab 2011b:

Sigma = eye(2)
nu = 1
wishrnd(Sigma,nu)


It will return something! Here, nu (degrees of freedom) is no longer greater than the dimension of the scale matrix, right?

Also, if you look close in to the implementation of wishrnd by issuing type wishrnd, we see the following line:

a = diag(sqrt(chi2rnd(df-(0:n-1))));


Clearly, the implementation generates $\chi^{2}(n-i)$ instead of $\chi^{2}(n-i+1)$.

Something incompatible here?

• See the first answer. This is a special case. But, I do doubt that the parameterization for this function differs from the normal definition. – Liangjie Hong Mar 8 '12 at 0:13

## 1 Answer

I don't see any incompatibility here; in fact, I think the Wikipedia article you cite provides your answer:

In fact the above definition can be extended to any real $n > p − 1$. If $n \leq p − 2$, then the Wishart no longer has a density—instead it represents a singular distribution.

Having $n<p$ doesn't mean the distribution is undefined; clearly from the definition of the Wishart distribution we can generate a sample by taking the (scaled) sample covariance matrix of $n$ i.i.d. multivariate Gaussian random variables. The result will be singular in the case of $n<p$, but that doesn't mean the distribution is undefined or that sampling from it should cause an error.

EDIT: The above assumes you're talking about integer values of $n$; otherwise you do get an error as expected:

>> wishrnd(eye(3), 1.5)
Error using wishrnd (line 50)
Degrees of freedom must be a positive integer, or a non-integer value larger than
the dimension of SIGMA minus 1.

• Thanks for the insights about the case. However, this still does not related to how it is generated, right? – Liangjie Hong Mar 8 '12 at 0:12
• First: notice that the code you're referring to is only used for large df (df>81+p). In the test case you give, it's just generating random normals using the code in the preceding block. Second: the code matches the Wikipedia page; df-(0:n-1)) is the same as df-1+(1:n)). If the line read df-(1:n) then that would correspond to $\chi^2(n-i)$ like you claim. – bnaul Mar 8 '12 at 1:18