simulate.m

function [x]          = simulate(k,p,t0,s0,S0,numberofdatapoints)
% p-dimensional, k components, number of points
pi=[];

for i=1:k
    pi=gamrnd(1+repmat(0,1,1),ones(k,1),k,1);  % gammas to make a Dirichlet
end

%    under uniform Dir prior
pi=pi/sum(pi);                              % creates a Dirichlet
[pi,i]=sort(pi,'descend');
pi(k)=0.005;%make sure one component is very low probability
pi=pi./sum(pi);


for j=1:k
    %mus and Sigmas for each component
    Sigma(:,:,j) = iwishrnd(S0,s0+p-1);
    mu(:,j) = sqrt(t0)*chol(Sigma(:,:,j))'*randn(p,1);
end

z=zeros(numberofdatapoints,1);
for j=1:numberofdatapoints
    z(j)=1+sum(cumsum(pi)<=rand(1,1));
    i=z(j);
    x(:,j)=mu(:,i)+chol(Sigma(:,:,i))'*randn(p,1);
end

nz=sum(outer(z,1:k,'=='))';

Project Homepage: