# Matching X'X with Wishart Samples in R

$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.

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 :)