function z=mvkdest(x,A,k,h)

% MVKDEST     Multivariate Kernel Density Estimator.
%
%             MVKDEST(X,A,K,H) evaluates a multivariate Kernel Density 
%             Estimate with kernel given by the function in K.
%             'KMVN' is the multivariate normal distribution.
%
%             X contains the point(s) for which the density 
%             estimate is required.  
%
%             The multivariate data points are stored row by row in A.
%
%             The optional parameter H represents the window width.
%             If omitted a simplistic 'optimal' value is used.
%
%             Christian C. Beardah 1993

[n d]=size(A); % Establish dimensions of A, 
               % n=number of data points,
               % d=dimension of the data. 

s=sqrt((sum(sqrt(diag(cov(A)))))/d);

[nx d]=size(x); % Establish number of data points in matrix x.

if nargin==3,  % If no value of h has been specified by the user then 
               % calculate one.
  if k=='kmvn',  % Normal kernel
    a=(4/(2*d+1))^(1/(d+4));
  end; 

  h=a*s*n^(-1/(d+4)); % Calculate 'optimum' value of h.
end;

for j=1:nx,
  B(j,:)=A(1,:);
end;

sum=feval(k,(x-B)/h);

for i=2:n,     % Sum over all data points
  for j=1:nx,
    B(j,:)=A(i,:);
  end; 
  sum=sum+feval(k,(x-B)/h);
end;

z=sum/(n*h^d); % z holds the kernel density estimates for the points in x
 