function [Da,Ds]=repsnmd3(w,sigma,mu,k,ss,n)

% REPSNMD3    Repeated sampling from a tri-variate normal mixture density.
%
%             Example: load nmd3
%                      [Da,Ds]=repsnmd3(wlav,sigmalav,mulav,'knorm',50,10);
% 
%             REPSNMD3(W,SIGMA,MU,K,SS,N) takes N samples of size SS from the 
%             tri-variate normal mixture density in W, SIGMA, MU.
%
%             For each sample a measure of the difference between the KDE of the
%             sample and the true density of the whole dataset is taken.  KDEs are 
%             based upon the kernel K.  
%
%             ** The type of algorithm used for generating the smoothing **
%             ** parameter H is defined within the routine.              **
%
%             Christian C. Beardah 1996

inc=32;

for i=1:length(w),
  xmin(i)=mu(1,i)-3.5*sqrt(sigma(1,3*i-2));
  xmax(i)=mu(1,i)+3.5*sqrt(sigma(1,3*i-2));

  ymin(i)=mu(2,i)-3.5*sqrt(sigma(2,3*i-1));
  ymax(i)=mu(2,i)+3.5*sqrt(sigma(2,3*i-1));

  zmin(i)=mu(3,i)-3.5*sqrt(sigma(3,3*i));
  zmax(i)=mu(3,i)+3.5*sqrt(sigma(3,3*i));
end;

xpop=linspace(min(xmin),max(xmax),inc);
ypop=linspace(min(ymin),max(ymax),inc);
zpop=linspace(min(zmin),max(zmax),inc);

[X,Y,Z]=meshgrid(xpop,ypop,zpop);

fpop=nmd3plot(X,Y,Z,w,sigma,mu,0);

deltax=xpop(2)-xpop(1);
deltay=ypop(2)-ypop(1);
deltaz=zpop(2)-zpop(1);

for i=1:n,

  s=snmd3(ss,w,sigma,mu);

% H selection algorithm:

  h=[hns(s(:,1),k) hns(s(:,2),k) hns(s(:,3),k)];

  fsam=kdest3d(X,Y,Z,s,k,h);

  figure;
  subplot(121);
  [Tx,Ty,Tz,Xd,Yd,Zd,f,d]=pccont3(40,s,'knorm',inc,h);
  colormap(white);
  subplot(122);
  ht=pclevel(60,d); 
  const=ht(1); 
  surfcv4(Xd,Yd,Zd,f,inc,const,'w');
  colormap(white);

% Absolute difference

  Da(i)=deltax*deltay*deltaz*sum(sum(abs(fsam-fpop)));

% Squared difference

  Ds(i)=deltax*deltay*deltaz*sum(sum((fsam-fpop).^2));

end;