function[q] = mvgaussian(p,mn,cv)

%dimension of input data
d = size(mn,2); 
incv = inv(cv);

%remove mean
for i = 1:d
    subm(:,i) = p(:,i) - mn(i);
end

submi = subm * incv;
%given x (column of row vectors ) ( x * inv(cv) * x' )
mcvmt = dot( submi , subm, 2); 

dnorm = 1.0 / ((2*pi)^(d(1)/2) * norm(cv)^0.5);

q = dnorm * exp ( -0.5 * mcvmt );