How to generate random symmetric positive definite matrices using MATLAB?

The algorithm I described in the comments is elaborated below. I will use $\tt{MATLAB}$ notation.

function A = generateSPDmatrix(n)
% Generate a dense n x n symmetric, positive definite matrix

A = rand(n,n); % generate a random n x n matrix

% construct a symmetric matrix using either
A = 0.5*(A+A'); OR
A = A*A';
% The first is significantly faster: O(n^2) compared to O(n^3)

% since A(i,j) < 1 by construction and a symmetric diagonally dominant matrix
%   is symmetric positive definite, which can be ensured by adding nI
A = A + n*eye(n);

end

Several changes are able to be used in the case of a sparse matrix.

function A = generatesparseSPDmatrix(n,density)
% Generate a sparse n x n symmetric, positive definite matrix with
%   approximately density*n*n non zeros

A = sprandsym(n,density); % generate a random n x n matrix

% since A(i,j) < 1 by construction and a symmetric diagonally dominant matrix
%   is symmetric positive definite, which can be ensured by adding nI
A = A + n*speye(n);

end

In fact, if the desired eigenvalues of the random matrix are known and stored in the vector rc, then the command

A = sprandsym(n,density,rc);

will construct the desired matrix. (Source: MATLAB sprandsym website)

A usual way in Bayesian statistics is to sample from a probability measure on real symmetric positive-definite matrices such as Wishart (or Inverse-Wishart).

I don't use Matlab but a quick check on Google gives this command (available in the Statistics toolbox):

W = wishrnd(Sigma,df)

where Sigma is some user-fixed positive definite matrix such as the identity and df are degrees of freedom.

While Daryl's answer is great, it gives symmetric positive definite matrices with very high probability , but that probability is not 1. This method gives a random matrix being symmetric positive definite matrix with probability 1.

My answer relies on the fact that a positive definite matrix has positive eigenvalues.

The matrix symmetric positive definite matrix A can be written as , A = Q'DQ , where Q is a random matrix and D is a diagonal matrix with positive diagonal elements. The elements of Q and D can be randomly chosen to make a random A.

The matlab code below does exactly that

function A = random_cov(n)

Q = randn(n,n);

eigen_mean = 2; 
% can be made anything, even zero 
% used to shift the mode of the distribution

A = Q' * diag(abs(eigen_mean+randn(n,1))) * Q;

return