meancov.m

The pattern recognition matlab toolbox

%MEANCOV Estimation of the means and covariances from multiclass data
% 
%   [U,G] = MEANCOV(A,N)
% 
%	INPUT
% 	A		Dataset
%	N		Normalization to use for calculating covariances: by M, the number
%					of samples in A (N = 1) or by M-1 (default, unbiased, N = 0).
%
% OUTPUT
% 	U		Mean vectors
% 	G		Covariance matrices
%
% DESCRIPTION  
%	Computation of a set of mean vectors U and a set of covariance matrices G
%	of the C classes in the dataset A. The covariance matrices are stored as a
%	3-dimensional matrix G of the size K x K x C, the class mean vectors as a
%	labeled dataset U of the size C x K.
%
% The use of soft labels or target labels is supported.
% 
%	SEE ALSO 
%	DATASETS, GAUSS, NBAYESC, DISTMAHA

% Copyright: R.P.W. Duin, [e-mail]
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: meancov.m,v 1.4 2007/04/13 09:30:54 duin Exp $

function [U,G] = meancov(a,n)

	prtrace(mfilename);

	% N determines whether the covariances are normalized by M (N = 1) or by 
	% M-1 (unbiased, N = 0), where M is the number of objects.

	if (nargin < 2 | isempty(n))
		prwarning(4,'normalisation not specified, assuming by M-1');
		n = 0;
	end

	if (n ~= 1) & (n ~= 0)
		error('Second parameter should be either 0 or 1.')
	end
	
	if (isdouble(a))			% A is a matrix: compute mean and covariances
		U = mean(a);							% 	in the usual way.
		if nargout > 1
			G = cov(a,n); 
		end
    
  elseif (isdatafile(a))
    
		if nargout  < 2
    	U = meancov_datafile(a,n);
		else
      [U,G] = meancov_datafile(a,n);
    end
    
  elseif (isdataset(a))
    
		[m,k,c] = getsize(a);
		if nargout == 2
			G = zeros(k,k,c);
		end
  	if (islabtype(a,'crisp'))
			
			if (c==0)  % special solution if all data is unlabeled
				U = mean(+a);
				if nargout > 1
					G = cov(+a,n); 
				end
			else
				for i = 1:c     
					J = findnlab(a,i);
					if isempty(J)
						U(i,:) = repmat(NaN,1,k);
						if nargout > 1
							G(:,:,i) = repmat(NaN,k,k);
						end
					else
						U(i,:) = mean(a(J,:),1);
						if (nargout > 1)
							G(:,:,i) = covm(a(J,:),n);	
						end
					end
				end
			end
			labu = getlablist(a);
  	elseif (islabtype(a,'soft'))
  		problab = gettargets(a);
			% Here we also have to be careful for unlabeled data
			if (c==0)
				prwarning(2,'The dataset has soft labels but no targets defined: using targets 1');
				U = mean(+a);
				if nargout > 1
					G = cov(+a,n);
				end
			else
				U = zeros(c,k);
				for i = 1:c

					% Calculate relative weights for the means.
					g = problab(:,i); nn = sum(g); g = g/mean(g); 

					U(i,:) = mean(a.*repmat(g,1,k));	% Weighted mean vectors	

					if (nargout > 1)

						u  = mean(a.*repmat(sqrt(g),1,k));

						% this appears to be needed to weight cov terms properly
						G(:,:,i) = covm(a.*repmat(sqrt(g),1,k),1) - U(i,:)'*U(i,:) + u'*u;

						% Re-normalise by M-1 if requested.
						if (n == 0)
							G(:,:,i) = m*G(:,:,i)/(m-1);
						end
					end
				end
			end
			labu = getlablist(a);
  	else
			% Default action.
  		U = mean(a);
			if nargout > 1
  			G = covm(a,n);
			end
			labu = [];
  	end

		% Add attributes of A to U.
		U = dataset(U,labu,'featlab',getfeatlab(a), ...
                                'featsize',getfeatsize(a));
		if (~islabtype(a,'targets'))
		%	p = getprior(a);
			U = setprior(U,a.prior); 
		end
    
  else
    error('Illegal datatype')
    
	end

return



%MEANCOV Datafile overload

function [u,g] = meancov_datafile(a,n)

	if nargin < 2, n = []; end
	
  [m,k,c] = getsize(a);
  k = 0; % just for detecting first loop
	next = 1;
	while next > 0
		[b,next] = readdatafile(a,next);
    if k == 0
      k = size(b,2);
      if c == 0
        u = zeros(1,k);
      else
        u = zeros(c,k);
      end
      if nargout > 1
        if c == 0
          g = zeros(k,k);
        else
          g = zeros(k,k,c);
        end
      end
    end
		%compute for each call to readdatafile contributions to mean and cov
		bb = + b;
		if c == 0
			u = u + sum(bb,1);
		else
			nlab = getnlab(b);
			for j=1:size(b,1)
				if nlab(j) > 0
					u(nlab(j),:) = u(nlab(j),:) + bb(j,:);
					if nargout > 1
						g(:,:,nlab(j)) = g(:,:,nlab(j)) + bb(j,:)'*bb(j,:);
					end
				end
			end
		end
	end

	f = classsizes(a);
	u = u ./ repmat(f',1,k);

	if nargout == 2
		for j=1:c
			g(:,:,j) = g(:,:,j)/f(j) - u(j,:)'*u(j,:);
			if isempty(n) | n == 0
				g(:,:,j) = (f(j)/(f(j)-1))*g(:,:,j);
			end
		end
	end
	u = dataset(u,getlablist(b));
	u.prior = b.prior;
return