function [X, sample_mu, sample_covariance] = generate_mvn(mu, S, n)
  % this method does not relies on the statistics package
  % as it is not necessarily available on every machine
  % instead it makes use of the fact that the Cholesky decomposition
  % of the covariance matrix can be used to correlatce otherwise decorrelated points

  if (nargin == 0)
    mu=[0 0];
    S=[1 0; 0 1];
    n=500;
  elseif (nargin == 1)
    S=[1 0; 0 1];
    n=500;
  elseif (nargin == 2)
    n=500;
  end

  R = chol(S); % note that the result of the chol() method is an upper-triangular matrix
  X = randn(n, 2);
  X = (R'*X')'; % or equivalenty X=X*R
  X = bsxfun('plus', X, mu);
  sample_mu = mean(X);
  sample_covariance = cov(X);

  minimum = min(X);
  maximum = max(X);
  x1_interval = minimum(1)-.1:.05:maximum(1)+.1;
  x2_interval = minimum(2)-.1:.05:maximum(2)+.1;
  [A, B] = meshgrid(x1_interval, x2_interval);

% this is a rather ugly (yet working) way to compute multivariate normal pdfs for the bivariate case
% density = @(m, c, sample_matrix) 1/(2*pi)*det(c)^-.5*exp(-.5*sum((inv(c)*(sample_matrix-m)')'.*(sample_matrix-m), 2));
% Z = density(mu, S, [A(:) B(:)]);
  Z = 1/(2*pi) * exp (-sumsq ((bsxfun('minus', [A(:) B(:)], mu))/R, 2)/2) / prod (diag(R)); % again, we assume the bivariate case, in the more general case the first term should be (2*pi)^(-d/2), where d is the number of dimensions a data point consists of
  Z = reshape(Z, size(A));

  scatter(X(:,1), X(:,2));
  hold on;
  contour(A, B, Z);
  hold off;
  % make the data matrix X decorrelated somehow (HINT: think of the Cholesky decomposition of the covariance matrix)
end
