Rachel-Zhang 提供的源码。
高斯混合模型
没有输入参数判断，没有协方差是否可逆验证。
我要用语音处理的，电脑卡死机，逆矩阵不是所有的都有的。
或者用文库里面的代码
function varargout = gmm(X, K_or_centroids)
% ============================================================
% Expectation-Maximization iteration implementation of
% Gaussian Mixture Model.
%
% PX = GMM(X, K_OR_CENTROIDS)
% [PX MODEL] = GMM(X, K_OR_CENTROIDS)
%
% - X: N-by-D data matrix.----------NxD的矩阵
% - K_OR_CENTROIDS: either K indicating the number of--------单个数字K/[K] 或者 KxD矩阵的聚类
% components or a K-by-D matrix indicating the
% choosing of the initial K centroids.
%
% - PX: N-by-K matrix indicating the probability of each--------NxK矩阵，第N个数据点占第K个单一高斯概率密度
% component generating each point.
% - MODEL: a structure containing the parameters for a GMM:
% MODEL.Miu: a K-by-D matrix.-------------KxD矩阵，初始化聚类样本，后面循环为每个数据点概率归一化后再聚类概率归一化后的均值矩阵
% MODEL.Sigma: a D-by-D-by-K matrix.------DxDxK矩阵，初始化数据点对于聚类的方差[聚类等概率]，后面循环为均值矩阵改变以后的方差
% MODEL.Pi: a 1-by-K vector.-----------1xK矩阵，初始化数据点使用聚类的概率分布,后面循环高斯混合概率系数归一化的分母Nk/N,高斯混合加权系数向量
% ============================================================
threshold = 1e-15;%阈值
[N, D] = size(X);%矩阵X的行N，列D
if isscalar(K_or_centroids)%判断值是否为1x1矩阵即单个数字
K = K_or_centroids;%取[k]的k
% randomly pick centroids
rndp = randperm(N);%返回一行由1到N个的整数无序排列的向量
centroids = X(rndp(1:K), :);%取出X矩阵行打乱后的矩阵的前K行
else
K = size(K_or_centroids, 1);%取矩阵K_or_centroids的行数
centroids = K_or_centroids;%取矩阵K_or_centroids矩阵
end
% initial values
[pMiu pPi pSigma] = init_params();
%初始化 嵌套函数后面可见。KxD矩阵pMiu聚类采样点，1*K向量pPi使用同一个聚类采样点出现概率，D*D*K的矩阵pSigma矩阵X的列项j对于聚类采样点的协方差
Lprev = -inf; %inf表示正无究大，-inf表示为负无究大
while true
Px = calc_prob();%NxK矩阵Px存放聚类点k（共有聚类点K个）对于全部数据点的正态分布生成样本的概率密度
% new value for pGamma
pGamma = Px .* repmat(pPi, N, 1);%NxK矩阵pGamma在使用聚类采样点k的条件下某个数据点n生成样本概率密度（条件概率密度）
pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K); %NxK矩阵pGamma对于使用数据点条件概率密度行向归一化。
%求每个样本由第K个聚类，也叫“component“生成的概率）
% new value for parameters of each Component
Nk = sum(pGamma, 1);%1xK矩阵Nk第k个聚类点被数据点使用的概率总和
pMiu = diag(1./Nk) * pGamma' * X;
%KxD矩阵 重新计算每个component的均值 维数变化KxK*KxN*NxD=KxD 数据点进行 聚类概率归一化.条件概率密度.数据点=得到均值（期望）
%均值=Σ概率Pi*数据点某一属性Xi；而这里还多了个聚类概率Nk
pPi = Nk/N; %更新混合高斯的加权系数
for kk = 1:K %重新计算每个component的协方差矩阵
Xshift = X-repmat(pMiu(kk, :), N, 1);%NxD矩阵Xshift在某一个聚类点的情况下，每个属性在这个聚类下的对均值（期望）差数(X-μi)
pSigma(:, :, kk) = (Xshift' * ...
(diag(pGamma(:, kk)) * Xshift)) / Nk(kk);%DxD矩阵pSigma(::,i) 概率Pi 聚类概率=1/Nk(i)
%第i个方差矩阵= Σ(X-μi)转置*概率Pi*(X-μi)*第i个聚类概率
end
% check for convergence
L = sum(log(Px*pPi')); %求混合高斯分布的似然函数
if L-Lprev < threshold %随着迭代次数的增加，似然函数越来越大，直至不变
break; %似然函数收敛则退出
end
Lprev = L;
end
if nargout == 1 %如果返回是一个参数的话，那么varargout=Px;
varargout = {Px};
else %否则，返回[Px model],其中model是结构体
model = [];
model.Miu = pMiu;
model.Sigma = pSigma;
model.Pi = pPi;
varargout = {Px, model};
end
function [pMiu pPi pSigma] = init_params()
pMiu = centroids;%得X矩阵中的任意K行，KxD矩阵 聚类点
pPi = zeros(1, K);%获取K维零向量[0 0 ...0] 加权系数（每行数据与聚类中点最小距离的概率）
pSigma = zeros(D, D, K);%获取K个DxD的零矩阵
%distmat为D维距离差平方和
% hard assign x to each centroids %X有NxD；sum(X.*X, 2)为Nx1； repmat(sum(X.*X, 2), 1, K)行向整体1倍，列向整体K倍；结果NxK
distmat = repmat(sum(X.*X, 2), 1, K) + ... %distmat第j行的第i个元素表示第j个数据与第i个聚类点的距离，如果数据有4个，聚类2个，那么distmat就是4*2矩阵
repmat(sum(pMiu.*pMiu, 2)', N, 1) - 2*X*pMiu'; %sum(A，2)结果为每行求和列向量，第i个元素是第i行的求和；
[dummy labels] = min(distmat, [], 2); %返回列向量dummy和labels，dummy向量记录distmat的每行的最小值，labels向量记录每行最小值的列号（多个取第一个），即是第几个聚类，labels是N×1列向量，N为样本数
for k=1:K
Xk = X(labels == k, :); %把标志为同一个聚类的样本组合起来
pPi(k) = size(Xk, 1)/N; %求混合高斯模型的加权系数，pPi为1*K的向量
pSigma(:, :, k) = cov(Xk);
%如果X为Nx1数组，那么cov(Xk)求单个高斯模型的协方差矩阵
%如果X为NxD(D>1)的矩阵，那么cov(Xk)求聚类样本的协方差矩阵
%cov()求出的为方阵--《概率论与数理统计》-多维随机变量的数字特征，且是对称矩阵（上三角和下三角对称）--《线性代数》
%pSigma为D*D*K的矩阵
end
end
function Px = calc_prob()
Px = zeros(N, K);%NxK零矩阵
for k = 1:K
Xshift = X-repmat(pMiu(k, :), N, 1); %NxD矩阵Xshift表示为对于一个1xD聚类点向量行向增N倍的样本矩阵-Uk，第i行表示xi-uk
inv_pSigma = inv(pSigma(:, :, k)); %求协方差矩阵的逆，pSigmaD*D*K的矩阵， inv_pSigmaD*D的矩阵,Σ^-1
tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);
%tmp为N*1矩阵，第i行表示（xi-uk)^T*Sigma^-1*(xi-uk) 即-(x-μ)转置x(Σ^-1).(x-μ) ----矩阵有叉乘（矩阵乘）和点乘
coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma));
%求多维正态分布中指数前面的系数，（2π)^(-d/2)*|(Σ^-(1/2))|
Px(:, k) = coef * exp(-0.5*tmp); %NxK矩阵求单独一个D维正态分布生成样本的概率密度或贡献
end
end
end
%聚类 数据挖掘
%矩阵算法 线性代数
%概率密度，近似然 数理统计
被我稍微修改备注的
function model = gmmEM(data,Kk,option)
%
% K 为model数分为几个聚类
% Reference: PRML p438-439
tic
if nargin < 3
option.eps = 1e-12;%收敛设置
option.maxiter = 1000;%
end
global num_data K;
K=Kk;
X = data'; %X为D*N型数据，跟PRML对样本数据描述相反
[dim,num_data] = size(X);
%Initialize
%-------------------------------
%K = numel(unique(data.y));
[inx, C,~] = kmeans(data,K);%[inx, C,~] = kmeans(X',K);把对象分为K个聚类
mu = C';%聚类点均值转置
%inx D个数据项的聚类编号
pai = zeros(1,K);%加权系数清零
E = zeros(dim,dim,K);
for k=1:K
pai(k) = sum(inx==k); %inx==index 统一聚类的个数统计
E(:,:,k) = eye(dim);%E单位矩阵 初始化协方差矩阵
end
pai = pai/num_data;%各个聚类的加权系数
iter = 0;%游标、统计迭代次数
log_val = logGMM(X,mu,E,pai);
try
iter
while iter<option.maxiter
iter = iter+1;
% E step
Yz = compu_Yz(X,mu,E,pai);%聚类k归一化后的多维条件正态概率
%前面使用了数据点使用k聚类归一化（纵向归一化）而后造成没有横向归一化
% M step 高斯迭代（横向归一化需要进行均值，协方差矩阵更新）
NK = sum(Yz);%每个聚类被使用的比例总和存放1Xk, 用于后面的横向归一化
for k = 1:K
mu(:,k) = 1/NK(k)*X*Yz(:,k); %均值矩阵更新 DXN NXK
%均值=Σ概率Pi*数据点某一属性Xi；而这里还多了个聚类概率Nk
Ek = zeros(dim,dim);%协方差矩阵清零
%n=1:num_data;
for n=1:num_data
Ek = Ek+Yz(n,k)*(X(:,n)-mu(:,k))*(X(:,n)-mu(:,k))';
end
E(:,:,k) = 1/NK(k)*Ek;
%第i个方差矩阵= Σ(X-μi)转置*概率Pi*(X-μi)*第i个聚类概率
end
pai = NK/num_data;
%检查是否收敛
log_val_new = logGMM(X,mu,E,pai);
if abs(log_val_new-log_val) < option.eps
model.Yz = Yz;
model.mu = mu;
model.E = E;
model.iter = iter;
return;
end
log_val = log_val_new;
if mod(iter,10)==0
disp(['-----进行第' num2str(iter) '迭代...'])
end
end
catch
fprintf('出错!!\n已经迭代次数:%d\n最大迭代次数\n',iter,option.maxiter);
model.Yz = Yz;
model.mu = mu;
model.E = E;
model.iter = iter;
%return;
end
if iter==option.maxiter
fprintf('达到最大迭代次数%d',maxiter)
model.Yz = Yz;
model.mu = mu;
model.E = E;
model.iter = iter;
end
toc
%model.usedTime = toc-tic;
end
%------------------------------------
function val = logGMM(X,mu,E,pai)
%样本数据X（行数据项，列属性）。均值矩阵。协方差矩阵（方阵，对角对称）。加权系数
global num_data K
val = 0;
%N = size(X,2);
%K = size(mu,2);
for n=1:num_data
tmp = 0;
for k=1:K
p = mvnpdf(X(:,n),mu(:,k),E(:,:,k));
tmp = tmp+pai(k)*p;
end
val = val + log(tmp);
end
end
%---------------------------------------------------------------
function Yz = compu_Yz(X,mu,E,pai)
%样本数据X（行数据项，列属性）。均值矩阵。协方差矩阵。加权系数
global num_data K
Yz = zeros(num_data,K);%清零
for n = 1:num_data
Y_nK = zeros(1,K);%清零 第k个聚类的
for k=1:K
Y_nK(k) = pai(k)*mvnpdf(X(:,n),mu(:,k),E(:,:,k));
%库内mvnpdf 返回DX1数组Y_nk (DX1;DX1;DXD)
%mvnpdf 单独一个N维正态分布生成样本的概率密度（离散此处值当作密度）
%pai(k)*正态概率密度=条件概率密度（在使用某个聚类的前提下）
end
Yz(n,:) = Y_nK/sum(Y_nK);%概率归一化 （一个数据使用k个聚类的比例）
%求第n个数据分别在聚类中的单独一个N维正态分布生成样本的概率（密度）
end
end
matlab自带mvnpdf的源代码
function y = mvnpdf(X, Mu, Sigma)
%MVNPDF Multivariate normal probability density function (pdf).
% Y = MVNPDF(X) returns the probability density of the multivariate normal
% distribution with zero mean and identity covariance matrix, evaluated at
% each row of X. Rows of the N-by-D matrix X correspond to observations or
% points, and columns correspond to variables or coordinates. Y is an
% N-by-1 vector.
%
% Y = MVNPDF(X,MU) returns the density of the multivariate normal
% distribution with mean MU and identity covariance matrix, evaluated
% at each row of X. MU is a 1-by-D vector, or an N-by-D matrix, in which
% case the density is evaluated for each row of X with the corresponding
% row of MU. MU can also be a scalar value, which MVNPDF replicates to
% match the size of X.
%
% Y = MVNPDF(X,MU,SIGMA) returns the density of the multivariate normal
% distribution with mean MU and covariance SIGMA, evaluated at each row
% of X. SIGMA is a D-by-D matrix, or an D-by-D-by-N array, in which case
% the density is evaluated for each row of X with the corresponding page
% of SIGMA, i.e., MVNPDF computes Y(I) using X(I,:) and SIGMA(:,:,I).
% If the covariance matrix is diagonal, containing variances along the
% diagonal and zero covariances off the diagonal, SIGMA may also be
% specified as a 1-by-D matrix or a 1-by-D-by-N array, containing
% just the diagonal. Pass in the empty matrix for MU to use its default
% value when you want to only specify SIGMA.
%
% If X is a 1-by-D vector, MVNPDF replicates it to match the leading
% dimension of MU or the trailing dimension of SIGMA.
%
% Example:
%
% mu = [1 -1]; Sigma = [.9 .4; .4 .3];
% [X1,X2] = meshgrid(linspace(-1,3,25)', linspace(-3,1,25)');
% X = [X1(:) X2(:)];
% p = mvnpdf(X, mu, Sigma);
% surf(X1,X2,reshape(p,25,25));
%
% See also MVTPDF, MVNCDF, MVNRND, NORMPDF.
% Copyright 1993-2011 The MathWorks, Inc.
if nargin<1
error(message('stats:mvnpdf:TooFewInputs'));
elseif ndims(X)~=2
error(message('stats:mvnpdf:InvalidData'));
end
% Get size of data. Column vectors provisionally interpreted as multiple scalar data.
[n,d] = size(X);
if d<1
error(message('stats:mvnpdf:TooFewDimensions'));
end
% Assume zero mean, data are already centered
if nargin < 2 || isempty(Mu)
X0 = X;
% Get scalar mean, and use it to center data
elseif numel(Mu) == 1
X0 = X - Mu;
% Get vector mean, and use it to center data
elseif ndims(Mu) == 2
[n2,d2] = size(Mu);
if d2 ~= d % has to have same number of coords as X
error(message('stats:mvnpdf:ColSizeMismatch'));
elseif n2 == n % lengths match
X0 = X - Mu;
elseif n2 == 1 % mean is a single row, rep it out to match data
X0 = bsxfun(@minus,X,Mu);
elseif n == 1 % data is a single row, rep it out to match mean
n = n2;
X0 = bsxfun(@minus,X,Mu);
else % sizes don't match
error(message('stats:mvnpdf:RowSizeMismatch'));
end
else
error(message('stats:mvnpdf:BadMu'));
end
% Assume identity covariance, data are already standardized
if nargin < 3 || isempty(Sigma)
% Special case: if Sigma isn't supplied, then interpret X
% and Mu as row vectors if they were both column vectors
if (d == 1) && (numel(X) > 1)
X0 = X0';
d = size(X0,2);
end
xRinv = X0;
logSqrtDetSigma = 0;
% Single covariance matrix
elseif ndims(Sigma) == 2
sz = size(Sigma);
if sz(1)==1 && sz(2)>1
% Just the diagonal of Sigma has been passed in.
sz(1) = sz(2);
sigmaIsDiag = true;
else
sigmaIsDiag = false;
end
% Special case: if Sigma is supplied, then use it to try to interpret
% X and Mu as row vectors if they were both column vectors.
if (d == 1) && (numel(X) > 1) && (sz(1) == n)
X0 = X0';
d = size(X0,2);
end
%Check that sigma is the right size
if sz(1) ~= sz(2)
error(message('stats:mvnpdf:BadCovariance'));
elseif ~isequal(sz, [d d])
error(message('stats:mvnpdf:CovSizeMismatch'));
else
if sigmaIsDiag
if any(Sigma<=0)
error(message('stats:mvnpdf:BadDiagSigma'));
end
R = sqrt(Sigma);
xRinv = bsxfun(@rdivide,X0,R);
logSqrtDetSigma = sum(log(R));
else
% Make sure Sigma is a valid covariance matrix
[R,err] = cholcov(Sigma,0);
if err ~= 0
error(message('stats:mvnpdf:BadMatrixSigma'));
end
% Create array of standardized data, and compute log(sqrt(det(Sigma)))
xRinv = X0 / R;
logSqrtDetSigma = sum(log(diag(R)));
end
end
% Multiple covariance matrices
elseif ndims(Sigma) == 3
sz = size(Sigma);
if sz(1)==1 && sz(2)>1
% Just the diagonal of Sigma has been passed in.
sz(1) = sz(2);
Sigma = reshape(Sigma,sz(2),sz(3))';
sigmaIsDiag = true;
else
sigmaIsDiag = false;
end
% Special case: if Sigma is supplied, then use it to try to interpret
% X and Mu as row vectors if they were both column vectors.
if (d == 1) && (numel(X) > 1) && (sz(1) == n)
X0 = X0';
[n,d] = size(X0);
end
% Data and mean are a single row, rep them out to match covariance
if n == 1 % already know size(Sigma,3) > 1
n = sz(3);
X0 = repmat(X0,n,1); % rep centered data out to match cov
end
% Make sure Sigma is the right size
if sz(1) ~= sz(2)
error(message('stats:mvnpdf:BadCovarianceMultiple'));
elseif (sz(1) ~= d) || (sz(2) ~= d) % Sigma is a stack of dxd matrices
error(message('stats:mvnpdf:CovSizeMismatchMultiple'));
elseif sz(3) ~= n
error(message('stats:mvnpdf:CovSizeMismatchPages'));
else
if sigmaIsDiag
if any(any(Sigma<=0))
error(message('stats:mvnpdf:BadDiagSigma'));
end
R = sqrt(Sigma);
xRinv = X0./R;
logSqrtDetSigma = sum(log(R),2);
else
% Create array of standardized data, and vector of log(sqrt(det(Sigma)))
xRinv = zeros(n,d,superiorfloat(X0,Sigma));
logSqrtDetSigma = zeros(n,1,class(Sigma));
for i = 1:n
% Make sure Sigma is a valid covariance matrix
[R,err] = cholcov(Sigma(:,:,i),0);
if err ~= 0
error(message('stats:mvnpdf:BadMatrixSigmaMultiple'));
end
xRinv(i,:) = X0(i,:) / R;
logSqrtDetSigma(i) = sum(log(diag(R)));
end
end
end
elseif ndims(Sigma) > 3
error(message('stats:mvnpdf:BadCovariance'));
end
% The quadratic form is the inner products of the standardized data
quadform = sum(xRinv.^2, 2);
y = exp(-0.5*quadform - logSqrtDetSigma - d*log(2*pi)/2);