gauss.m

The pattern recognition matlab toolbox

%GAUSS Generation of a multivariate Gaussian dataset
% 
% 	A = GAUSS(N,U,G,LABTYPE)
%
% INPUT
%	  N		     Array of number of objects to generate for each class
%	  U		     Dataset with means, labels and priors of classes (default: 
%              zeros(C,K), where C is the length of array N and K = 1)
%	  G		     Covariance matrix for each class (default: eye(K))
%	  LABTYPE	 Label type (default: 'crisp')
%
% OUTPUT
%	  A		     Dataset containing multivariate Gaussian data
%
% DESCRIPTION
% Generation of N K-dimensional Gaussian distributed samples for C classes.
% The covariance matrices should be specified in G (size K*K*C) and the
% means, labels and prior probabilities can be defined by the dataset U with
% size (C*K). If U is not a dataset, it should be a C*K matrix and A will
% be a dataset with C classes.
%
% If N is a vector, exactly N(I) objects are generated for class I, I = 1..C.
% 
% EXAMPLES
% 1. Generation of 50 points for each of two 1-dimensional distributions with
%    mean -1 and 1 and with variances 1 and 2:
%
%	     GAUSS([50 50],[-1;1],CAT(3,1,2))
%
%   Note that the two 1-dimensional class means should be given as a column
%   vector [1;-1], as [1 -1] defines a single 2-dimensional mean. Note that
%   the 1-dimensional covariance matrices degenerate to scalar variances,
%   but have still to be combined into a collection of square matrices using
%   the CAT(3,....) function.
%
% 2. Generation of 300 points for 3 classes with means [0 0], [0 1] and 
%    [1 1] and covariance matrices [2 1; 1 4], EYE(2) and EYE(2):
%
%      GAUSS(300,[0 0; 0 1; 1 1]*3,CAT(3,[2 1; 1 4],EYE(2),EYE(2)))
%
% SEE ALSO 
% DATASETS, PRDATASETS

% 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: gauss.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function a = gauss(n,u,g,labtype)

	prtrace(mfilename);

	if (nargin < 1)
		prwarning (2,'number of samples not specified, assuming N = 50'); 	
		n = 50;
		end
		cn = length(n);
	if (nargin < 2)
		prwarning (2,'means not specified; assuming one dimension, mean zero');
		u = zeros(cn,1); 
	end;
	if (nargin < 3)
		prwarning (2,'covariances not specified, assuming unity');
	 	g = eye(size(u,2)); 
	end
	if (nargin < 4)
		prwarning (3,'label type not specified, assuming crisp');
		labtype = 'crisp'; 
	end

	% Return an empty dataset if the number of samples requested is 0.

	if (length(n) == 1) & (n == 0)
		a = dataset([]); 
		return
	end

	% Find C, desired number of classes based on U and K, the number of 
	% dimensions. Make sure U is a dataset containing the means.

	if (isa(u,'dataset'))
		[m,k,c] = getsize(u);
		lablist = getlablist(u);
		p = getprior(u);
	else
		[m,k] = size(u); 		
		c = m;
		lablist = genlab(ones(c,1));
		u = dataset(u,lablist);
		p = ones(1,c)/c;
	end

	if (cn ~= c) & (cn ~= 1)
		error('The number of classes specified by N and U does not match');
	end

 	% Generate a class frequency distribution according to the desired priors.
	n = genclass(n,p);

	% Find CG, the number of classes according to G. 
	% Make sure G is not a dataset.

	if (isempty(g))
		g = eye(k); 
		cg = 1;
	else
		g = +g; [k1,k2,cg] = size(g);
		if (k1 ~= k) | (k2 ~= k)
			error('The number of dimensions of the means U and covariance matrices G do not match');
		end
		if (cg ~= m & cg ~= 1)
			error('The number of classes specified by the means U and covariance matrices G do not match');
		end
	end

	% Create the data A by rotating and scaling standard normal distributions 
	% using the eigenvectors of the specified covariance matrices, and adding
	% the means.

	a = [];
	for i = 1:m
		j = min(i,cg);						% Just in case CG = 1 (if G was not specified).

		% Sanity check: user can pass non-positive definite G.		
		[V,D] = eig(g(:,:,j)); V = real(V); D = real(D);	
		a = [a; randn(n(i),k)*sqrt(D)*V' + repmat(+u(i,:),n(i),1)];
	end

	% Convert A to dataset by adding labels and priors.

	labels = genlab(n,lablist);
	a = dataset(a,labels,'lablist',lablist,'prior',p);

	% If non-crisp labels are requested, use output of Bayes classifier.
	switch (labtype)
		case 'crisp'
			;
		case 'soft'
			w = nbayesc(u,g); 		
			targets = a*w*classc;
			a = setlabtype(a,'soft',targets);
		otherwise
			error(['Label type ' labtype ' not supported'])
	end

	a = setname(a,'Gaussian Data');

return