wolffd@0: function M = sample_gaussian(mu, Sigma, N)
wolffd@0: % SAMPLE_GAUSSIAN Draw N random row vectors from a Gaussian distribution
wolffd@0: % samples = sample_gaussian(mean, cov, N)
wolffd@0: 
wolffd@0: if nargin==2
wolffd@0:   N = 1;
wolffd@0: end
wolffd@0: 
wolffd@0: % If Y = CX, Var(Y) = C Var(X) C'.
wolffd@0: % So if Var(X)=I, and we want Var(Y)=Sigma, we need to find C. s.t. Sigma = C C'.
wolffd@0: % Since Sigma is psd, we have Sigma = U D U' = (U D^0.5) (D'^0.5 U').
wolffd@0: 
wolffd@0: mu = mu(:);
wolffd@0: n=length(mu);
wolffd@0: [U,D,V] = svd(Sigma);
wolffd@0: M = randn(n,N);
wolffd@0: M = (U*sqrt(D))*M + mu*ones(1,N); % transform each column
wolffd@0: M = M';
wolffd@0: