gmmeval.m

来自「一个关于数据聚类和模式识别的程序,在生物化学,化学中因该都可以用到.希望对大家有」· M 代码 · 共 98 行

function [logProb, gaussianProb] = gmmEval(data, gmmParam);
% gmmEval: Evaluation of a GMM (Gaussian mixture model)
%	Usage: [logProb, gaussianProb] = gmmEval(data, gmmParam);
%		data: dim x dataNum matrix where each column is a data point
%		gmmParam(i): Parameters for Gaussian component i
%			gmmParam(i).mu: a mean vector of dim x 1
%			gmmParam(i).sigma: a covariance matrix of 3 possible dimensions:
%				1 x 1: identity covariance matrix times a constant for each Gaussian
%				dim x 1: diagonal covariance matrix for each Gaussian
%				dim x dim: full covariance matrix for each Gaussian
%			gmmParam(i).w: a weighting factor
%		logProb: 1 x dataNum vector of output probabilities
%		gaussianProb(i,j) is the probability of data(:,j) to the i-th Gaussian (This is for gmmTrain.m only)
%
%	For example, to plot 2D GMM:
%
%		data=linspace(-10, 10, 101);
%		gmmParam(1).mu = -5; gmmParam(1).sigma = 1; gmmParam(1).w = 0.1;
%		gmmParam(2).mu =  0; gmmParam(2).sigma = 4; gmmParam(2).w = 0.5;
%		gmmParam(3).mu =  5; gmmParam(3).sigma = 3; gmmParam(3).w = 0.4;
%		logProb = gmmEval(data, gmmParam);
%		prob=exp(logProb);
%		figure; plot(data, prob, '.-');
%		line(data, gmmParam(1).w*gaussian(data, gmmParam(1)), 'color', 'r');
%		line(data, gmmParam(2).w*gaussian(data, gmmParam(2)), 'color', 'm');
%		line(data, gmmParam(3).w*gaussian(data, gmmParam(3)), 'color', 'g');
%
%	Another example to plot 3D GMM:
%
%		gmmParam(1).mu = [-3, 3]'; gmmParam(1).sigma = [5, 2]'; gmmParam(1).w = 0.3;
%		gmmParam(2).mu = [3, -3]'; gmmParam(2).sigma = [4, 1]'; gmmParam(2).w = 0.3;
%		gmmParam(3).mu =  [3, 3]'; gmmParam(3).sigma = [1, 4]'; gmmParam(3).w = 0.4;
%		bound = 8;
%		pointNum = 51;
%		x = linspace(-bound, bound, pointNum);
%		y = linspace(-bound, bound, pointNum);
%		[xx, yy] = meshgrid(x, y);
%		data = [xx(:), yy(:)]';
%		logProb = gmmEval(data, gmmParam);
%		zz = reshape(exp(logProb), pointNum, pointNum);
%		subplot(2,2,1);
%		mesh(xx, yy, zz); axis tight; box on
%		subplot(2,2,2);
%		contour(xx, yy, zz, 30); axis image; box on

%	Roger Jang, 20000602, 20080726

if nargin<1, selfdemo; return; end

[dim, dataNum]=size(data);
gaussianNum=length(gmmParam);
log2pi=log(2*pi);
logProb=zeros(1, dataNum);
logGaussianProb=zeros(gaussianNum, dataNum);

if prod(size(gmmParam(1).sigma))==1	% identity covariance matrix times a constant for each Gaussian
	for i=1:gaussianNum
		dataMinusMu = data-repmat(gmmParam(i).mu, 1, dataNum);
		logGaussianProb(i,:) = (-sum(dataMinusMu.*dataMinusMu, 1)/gmmParam(i).sigma-dim*(log2pi+log(gmmParam(i).sigma)))/2;
	end
elseif prod(size(gmmParam(1).sigma))==dim	% diagonal covariance matrix for each Gaussian
	for i=1:gaussianNum
		dataMinusMu = data-repmat(gmmParam(i).mu, 1, dataNum);
		logGaussianProb(i,:) = (-sum(dataMinusMu.*dataMinusMu./repmat(gmmParam(i).sigma, 1, dataNum), 1)-dim*log2pi-log(prod(gmmParam(i).sigma)))/2;
	end
else	% full covariance matrix for each Gaussian
	for i=1:gaussianNum
		dataMinusMu = data-repmat(gmmParam(i).mu, 1, dataNum);
		logGaussianProb(i,:) = (-sum((inv(gmmParam(i).sigma)'*dataMinusMu).*dataMinusMu, 1)-dim*log2pi-log(det(gmmParam(i).sigma)))/2;
	end
end

logw=log([gmmParam.w]');
for i=1:dataNum
	logProb(i)=mixLogSumMex(logw(:)+logGaussianProb(:,i));
end

if nargout>1
	gaussianProb=exp(logGaussianProb);	% This output is necessary for gmmTrain.m!
end

% ====== Self demo
function selfdemo
gmmParam(1).mu = [-3, 3]'; gmmParam(1).sigma = [5, 2]'; gmmParam(1).w = 0.3;
gmmParam(2).mu = [3, -3]'; gmmParam(2).sigma = [4, 1]'; gmmParam(2).w = 0.3;
gmmParam(3).mu =  [3, 3]'; gmmParam(3).sigma = [1, 4]'; gmmParam(3).w = 0.4;
bound = 8;
pointNum = 51;
x = linspace(-bound, bound, pointNum);
y = linspace(-bound, bound, pointNum);
[xx, yy] = meshgrid(x, y);
data = [xx(:), yy(:)]';
logProb = gmmEval(data, gmmParam);
zz = reshape(exp(logProb), pointNum, pointNum);
subplot(2,2,1);
mesh(xx, yy, zz); axis tight; box on
subplot(2,2,2);
contour(xx, yy, zz, 30); axis image; box on