function log_y = log_mvnpdf(X, Mu, Sigma)
%LOG_MVNPDF Log of multivariate normal probability density function (pdf).
%   LOG_Y = LOG_MVNPDF(X) returns the n-by-1 vector Y, containing the log probability
%   density of the multivariate normal distribution with zero mean and
%   identity covariance matrix, evaluated at each row of the n-by-d matrix
%   X. Rows of X correspond to observations and columns correspond to
%   variables or coordinates.
%
%   LOG_Y = LOG_MVNPDF(X,MU) returns the log density of the multivariate normal
%   distribution with mean MU and identity covariance matrix, evaluated
%   at each row of X.  MU is a 1-by-d vector, or an n-by-d matrix, in which
%   case the density is evaluated for each row of X with the corresponding
%   row of MU.  MU can also be a scalar value, which MVNPDF replicates to
%   match the size of X.
%
%   LOG_Y = LOG_MVNPDF(X,MU,SIGMA) returns the log density of the multivariate normal
%   distribution with mean MU and covariance SIGMA, evaluated at each row
%   of X.  SIGMA is a d-by-d matrix, or an d-by-d-by-n array, in which case
%   the density is evaluated for each row of X with the corresponding page
%   of SIGMA, i.e., MVNPDF computes Y(I) using X(I,:) and SIGMA(:,:,I).
%   Pass in the empty matrix for MU to use its default value when you want
%   to only specify SIGMA.
%
%   If X is a 1-by-d vector, MVNPDF replicates it to match the leading
%   dimension of MU or the trailing dimension of SIGMA.
%
%   Example:
%
%      mu = [1 -1];
%      Sigma = [.9 .4; .4 .3];
%      X = mvnrnd(mu, Sigma, 10);
%      log_p = log_mvnpdf(X, mu, Sigma);
%
%   See also MVNPDF, MVNRND, NORMPDF.
  
  if nargin < 1 | isempty(X)
    error('Requires the input argument X.');
  elseif ndims(X) > 2
    error('X must be a matrix.');
  end

  % Get size of data.  Column vectors provisionally interpreted as multiple scalar data.
  [n,d] = size(X);

  % Assume zero mean, data are already centered
  if nargin < 2 | isempty(Mu)
    X0 = X;

    % Get scalar mean, and use it to center data
  elseif prod(size(Mu)) == 1
    X0 = X - Mu;

    % Get vector mean, and use it to center data
  elseif ndims(Mu) == 2
    [n2,d2] = size(Mu);
    if d2 ~= d % has to have same number of coords as X
      error('X and MU must have the same number of columns.');
    elseif n2 == n % lengths match
      X0 = X - Mu;
    elseif n2 == 1 % mean is a single row, rep it out to match data
      X0 = X - repmat(Mu,n,1);
    elseif n == 1 % data is a single row, rep it out to match mean
      n = n2;
      X0 = repmat(X,n2,1) - Mu;
    else % sizes don't match
      error('X or MU must be a row vector, or X and MU must have the same number of rows.');
    end

  else
    error('MU must be a matrix.');
  end

  % Assume identity covariance, data are already standardized
  if nargin < 3 | isempty(Sigma)
    % Special case: if Sigma isn't supplied, then interpret X
    % and Mu as row vectors if they were both column vectors
    if d == 1 & prod(size(X)) > 1
      X0 = X0';
      [n,d] = size(X0);
    end
    xRinv = X0;
    sqrtInvDetSigma = 1;

    % Single covariance matrix
  elseif ndims(Sigma) == 2
    % Special case: if Sigma is supplied, then use it to try to interpret
    % X and Mu as row vectors if they were both column vectors.
    if (d == 1 & prod(size(X)) > 1) & size(Sigma,1) == n
      X0 = X0';
      [n,d] = size(X0);
    end

    % Make sure Sigma is the right size
    if size(Sigma,1) ~= d | size(Sigma,2) ~= d
      error('SIGMA must be a square matrix with size equal to the number of columns in X.');
    else
      % Make sure Sigma is a valid covariance matrix
      [spd,R] = isspd(Sigma);
      if spd
	% Create array of standardized data, vector of inverse det
	xRinv = X0 / R;
	sqrtInvDetSigma = 1 / prod(diag(R));
      else
	error('SIGMA must be symmetric and positive definite.');
      end
    end

    % Multiple covariance matrices
  elseif ndims(Sigma) == 3
    % Special case: if Sigma is supplied, then use it to try to interpret
    % X and Mu as row vectors if they were both column vectors.
    if (d == 1 & prod(size(X)) > 1) & size(Sigma,1) == n
      X0 = X0';
      [n,d] = size(X0);
    end

    % Data and mean are a single row, rep them out to match covariance
    if n == 1 % already know size(Sigma,3) > 1
      n = size(Sigma,3);
      X0 = repmat(X0,n,1); % rep centered data out to match cov
    end

    % Make sure Sigma is the right size
    if size(Sigma,1) ~= d | size(Sigma,2) ~= d
      error('Each page of SIGMA must be a square matrix with size equal to the number of columns in X.');
    elseif size(Sigma,3) ~= n
      error('SIGMA must have one page for each row of X.');
    else

      % Create array of standardized data, vector of inverse det
      xRinv = zeros(n,d);
      sqrtInvDetSigma = zeros(n,1);
      for i = 1:n
	% Make sure Sigma is a valid covariance matrix
	[spd,R] = isspd(Sigma(:,:,i));
	if spd
	  xRinv(i,:) = X0(i,:) / R;
	  sqrtInvDetSigma(i) = 1 / prod(diag(R));
	else
	  error('SIGMA must be symmetric and positive definite.');
	end
      end
    end

  elseif ndims(Sigma) > 3
    error('SIGMA must be a matrix or a 3 dimensional array.');
  end

  % Exponents in pdf are the inner products of the standardized data
  quadform = sum(xRinv.^2, 2);
  
  log_y = -d/2*log(2*pi) + log(sqrtInvDetSigma) -0.5*quadform;


function [t,R] = isspd(Sigma)
%ISPDS Test if a matrix is positive definite symmetric
% T = ISPDS(SIGMA) returns a logical indicating whether the matrix SIGMA is
% square, symmetric, and positive definite, i.e., it is a valid full rank
% covariance matrix.
%
% [T,R] = ISPDS(SIGMA) returns the cholesky factor of SIGMA in R.  If SIGMA
% is not square symmetric, ISPDS returns [] in R.

% Test for square, symmetric
  [n,m] = size(Sigma);
  if (n == m) & all(all(abs(Sigma - Sigma') < 10*eps*max(abs(diag(Sigma)))))
    % Test for positive definiteness
    [R,p] = chol(Sigma);
    if p == 0
      t = true;
    else
      t = false;
    end
  else
    R = [];
    t = false;
  end