1
$\begingroup$

$X'X \sim Wishart(\Sigma,n)$, however I'm having a tough time producing this in R.

Example:

data=cbind(rnorm(100,10,5),rnorm(100,5,2),rnorm(100,-4,3))

X=cbind(rnorm(1,10,5),rnorm(1,5,2),rnorm(1,-4,3))
t(X)%*%X
rWishart(10,99,cov(data))

The data generated from rWishart is not close to $X'X$. What am I doing wrong? The help documentation mentions $\Sigma$ should be a scaled matrix, however I'm unsure what this is.

$\endgroup$
4
$\begingroup$

a 'scale' matrix can be any matrix that is positive definite. Wishart distribution is often used in Bayesian hierarchical model to capture characteristics of a inverse covariance matrix.

Back to your problem, if you read the help page of rWishart carefully, it says:

If X1,...,Xm, Xi in R^p is a sample of m independent multivariate Gaussians with mean (vector) 0, and covariance matrix Σ, the distribution of M = X'X is W_p(Σ, m).

However in your toy example, you chose to sample $X_i$'s with different means and different variance, and the degree of freedom $p$ is predetermined by the size of your sample $X_i$, not randomly chosen.

A better example can be constructed as such:

require(MASS)
data = cbind(rnorm(100,0,5),rnorm(100,0,2),rnorm(100,0,3))
Sigma = cov(data) % this is a 3 by 3 matrix
eigen(Sigma) %check positiv definite

%construct X
X = mvtnorm(100,rep(0,3),Sigma)

%define d.f.
df = dim(X)[1]

%generate random wishart sample with df and Sigma
rWishart(10,df,Sigma)

% compute X'X
t(X)%*% X

You'll find that these random samples are roughly in the range of X'X, a more rigorous check can be done by looping said X'X a few times(~100,000) and take the empirical mean. In theory this should agree with the first moment of Wishart distribution m*Σ by law of large numbers.

You can certainly generate non-central Wishart distribution, a good reference on this topic (or in fact any matrix variate distribution) is to look at Matrix Variate Distributions by Gupta and Nagar.

Hope this helps :)

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.