gmmeval_new.m

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

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: GMM parameter
%			gmmParam.mix(i).mu: dim x gaussianNum matrix where each column is a mean vector
%			gmmParam.mix(i).sigma: 1 x gaussianNum vector where each element represents the covariance of a gaussian
%			gmmParam.mix(i).w: 1 x gaussianNum vector where each element is a weighting coefficient
%		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:
%		mixNum=3;
%		gmmParam.mix(1).mu=-5; gmmParam.mix(1).sigma=1; gmmParam.mix(1).w=0.1;
%		gmmParam.mix(2).mu= 0; gmmParam.mix(2).sigma=4; gmmParam.mix(2).w=0.5;
%		gmmParam.mix(3).mu= 5; gmmParam.mix(3).sigma=3; gmmParam.mix(3).w=0.4;
%		x=linspace(-10, 10, 101);
%		logProb = gmmEval(x, gmmParam);
%		prob=exp(logProb);
%		plot(x, prob, '.-');
%		line(x, gmmParam.mix(1).w*gaussian(x, gmmParam.mix(1).mu, gmmParam.mix(1).sigma), 'color', 'r');
%		line(x, gmmParam.mix(2).w*gaussian(x, gmmParam.mix(2).mu, gmmParam.mix(2).sigma), 'color', 'm');
%		line(x, gmmParam.mix(3).w*gaussian(x, gmmParam.mix(3).mu, gmmParam.mix(3).sigma), 'color', 'g');
%
%	Another example, to plot 3D GMM:
%		gmmParam.mix(1).mu=[-3; 3]; gmmParam.mix(1).sigma=1; gmmParam.mix(1).w=0.2;
%		gmmParam.mix(2).mu=[ 3;-3]; gmmParam.mix(2).sigma=3; gmmParam.mix(2).w=0.3;
%		gmmParam.mix(3).mu=[ 3; 3]; gmmParam.mix(3).sigma=2; gmmParam.mix(3).w=0.5;
%		bound = 8;
%		pointNum = 31;
%		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);
%		mesh(xx, yy, zz); axis tight; box on

%	Roger Jang, 20000602

if nargin==0, selfdemo; return; end

mu   =[gmmParam.mix.mu];
sigma=[gmmParam.mix.sigma];
w    =[gmmParam.mix.w];

if length(w)==1		% Reduce to a single gaussian!
	if w~=1
		fprintf('Error in the case of single gaussian where w is not equal to 1!\n');
	end
	logProb=sum(gaussianLog(data, mu, sigma));	% No need to have the second output
	return
end

[dim, dataNum]=size(data);
gaussianNum=length(sigma);
log2pi=log(2*pi);
logw=log(w);
logProb=zeros(1, dataNum);
logGaussianProb=zeros(gaussianNum, dataNum);
for i=1:gaussianNum
	dataMinusMu = data-mu(:,i)*ones(1, dataNum);
	logGaussianProb(i,:) = (-sum(dataMinusMu.*dataMinusMu, 1)/sigma(i)-dim*(log2pi+log(sigma(i))))/2;
end

for i=1:dataNum
%	logProb(i)=mixLogSum(logw(:)+logGaussianProb(:,i));
	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
% == 1D example
mixNum=3;
gmmParam.mix(1).mu=-5; gmmParam.mix(1).sigma=1; gmmParam.mix(1).w=0.1;
gmmParam.mix(2).mu= 0; gmmParam.mix(2).sigma=4; gmmParam.mix(2).w=0.5;
gmmParam.mix(3).mu= 5; gmmParam.mix(3).sigma=3; gmmParam.mix(3).w=0.4;
x=linspace(-10, 10, 101);
logProb = feval(mfilename, x, gmmParam);
prob=exp(logProb);
figure; plot(x, prob, '.-');
line(x, gmmParam.mix(1).w*gaussian(x, gmmParam.mix(1).mu, gmmParam.mix(1).sigma), 'color', 'r');
line(x, gmmParam.mix(2).w*gaussian(x, gmmParam.mix(2).mu, gmmParam.mix(2).sigma), 'color', 'm');
line(x, gmmParam.mix(3).w*gaussian(x, gmmParam.mix(3).mu, gmmParam.mix(3).sigma), 'color', 'g');
% === 2D example
gmmParam.mix(1).mu=[-3; 3]; gmmParam.mix(1).sigma=1; gmmParam.mix(1).w=0.2;
gmmParam.mix(2).mu=[ 3;-3]; gmmParam.mix(2).sigma=3; gmmParam.mix(2).w=0.3;
gmmParam.mix(3).mu=[ 3; 3]; gmmParam.mix(3).sigma=2; gmmParam.mix(3).w=0.5;
bound = 8;
pointNum = 31;
x = linspace(-bound, bound, pointNum);
y = linspace(-bound, bound, pointNum);
[xx, yy] = meshgrid(x, y);
data = [xx(:), yy(:)]';
logProb = feval(mfilename, data, gmmParam);
zz = reshape(exp(logProb), pointNum, pointNum);
figure; mesh(xx, yy, zz); axis tight; box on