Gaussian pdf(probability density function)

function [z] = gausspdf(x,mu,sigma);

% GAUSSPDF  Values of the Gaussian probability density
%           function
%
%    GAUSSPDF(X,MU,SIGMA) returns the likelihood of the point or set of points X with respect to a
%   Gaussian process with mean MU and covariance SIGMA. MU is a 1*D vector (where D is the
%   dimension of the process), SIGMA is a D*D matrix  and X is a N*D matrix where each row is a
%    D-dimensional point.
%
%    See also MEAN, COV, GLOGLIKE

[N,D] = size(x);
if (min(N,D) == 1), x=x(:)'; end;

invSig = inv(sigma); %  inverse of ()
mu = mu(:)';           %???

x = x-repmat(mu,N,1); % (x-mu)

z = sum( ((x*invSig).*x), 2 ); %(
z = exp(-0.5*z) / sqrt((2*pi)^D * det(sigma));

Note:

B = sum(A,dim) sums along the dimension of A specified by scalar dim. The dim input is an integer value from 1 to N, where N is the number of dimensions in A. Set dim to 1 to compute the sum of each column, 2 to sum rows, etc.

Y = inv(X) returns the inverse of the square matrix X. A warning message is printed if X is badly scaled or nearly singular.
In practice, it is seldom necessary to form the explicit inverse of a matrix. A frequent misuse of inv arises when solving the system of linear equations . One way to solve this is with x = inv(A)*b. A better way, from both an execution time and numerical accuracy standpoint, is to use the matrix division operator x = A\b. This produces the solution using Gaussian elimination, without forming the inverse.