EM アルゴリズムを使用して、特定のデータセットで 4 つのコンポーネントを含むガウス混合モデルをトレーニングしたいと考えています。このセットは 3 次元で、300 個のサンプルが含まれています。
rank(sigma) = 2問題は、EM アルゴリズムの約 6 ラウンド後、共分散行列 sigma が matlab ( 3 ではなく)に従って特異値に近くなることです。これは、ガウス分布を評価する複雑な値のような望ましくない結果につながりますgm(k,i)。
さらに、ガウスのログを使用してアンダーフローの問題を説明しました - E-step を参照してください。これが正しいかどうかはわかりませんが、責任 p(w_k | x^(i), theta) の exp を別の場所に持っていく必要があるかどうかはわかりません。
私の EM アルゴリズムの実装がこれまでのところ正しいかどうか教えてもらえますか? そして、特異共分散シグマに近い問題をどのように説明するのでしょうか?
これが EM アルゴリズムの私の実装です。
最初に、kmeans を使用してコンポーネントの平均と共分散を初期化しました。
load('data1.mat');
X = Data'; % 300x3 data set
D = size(X,2); % dimension
N = size(X,1); % number of samples
K = 4; % number of Gaussian Mixture components
% Initialization
p = [0.2, 0.3, 0.2, 0.3]; % arbitrary pi
[idx,mu] = kmeans(X,K); % initial means of the components
% compute the covariance of the components
sigma = zeros(D,D,K);
for k = 1:K
sigma(:,:,k) = cov(X(idx==k,:));
end
E-stepでは、次の式を使用して責任を計算しています。
w_k は k ガウス成分です。
x^(i) は単一のデータポイント (標本)
theta は、ガウス混合モデルのパラメータを表します: mu、Sigma、pi。
対応するコードは次のとおりです。
% variables for convergence
converged = 0;
prevLoglikelihood = Inf;
prevMu = mu;
prevSigma = sigma;
prevPi = p;
round = 0;
while (converged ~= 1)
round = round +1
gm = zeros(K,N); % gaussian component in the nominator
sumGM = zeros(N,1); % denominator of responsibilities
% E-step: Evaluate the responsibilities using the current parameters
% compute the nominator and denominator of the responsibilities
for k = 1:K
for i = 1:N
Xmu = X-mu;
% I am using log to prevent underflow of the gaussian distribution (exp("small value"))
logPdf = log(1/sqrt(det(sigma(:,:,k))*(2*pi)^D)) + (-0.5*Xmu*(sigma(:,:,k)\Xmu'));
gm(k,i) = log(p(k)) * logPdf;
sumGM(i) = sumGM(i) + gm(k,i);
end
end
% calculate responsibilities
res = zeros(K,N); % responsibilities
Nk = zeros(4,1);
for k = 1:K
for i = 1:N
% I tried to use the exp(gm(k,i)/sumGM(i)) to compute res but this leads to sum(pi) > 1.
res(k,i) = gm(k,i)/sumGM(i);
end
Nk(k) = sum(res(k,:));
end
Nk(k)は、M ステップで指定された式を使用して計算され、M ステップで新しい確率を計算するために使用されp(k)ます。
Mステップ
% M-step: Re-estimate the parameters using the current responsibilities
for k = 1:K
for i = 1:N
mu(k,:) = mu(k,:) + res(k,i).*X(k,:);
sigma(:,:,k) = sigma(:,:,k) + res(k,i).*(X(k,:)-mu(k,:))*(X(k,:)-mu(k,:))';
end
mu(k,:) = mu(k,:)./Nk(k);
sigma(:,:,k) = sigma(:,:,k)./Nk(k);
p(k) = Nk(k)/N;
end
収束を確認するために、次の式を使用して対数尤度が計算されます。
% Evaluate the log-likelihood and check for convergence of either
% the parameters or the log-likelihood. If not converged, go to E-step.
loglikelihood = 0;
for i = 1:N
loglikelihood = loglikelihood + log(sum(gm(:,i)));
end
% Check for convergence of parameters
errorLoglikelihood = abs(loglikelihood-prevLoglikelihood);
if (errorLoglikelihood <= eps)
converged = 1;
end
errorMu = abs(mu(:)-prevMu(:));
errorSigma = abs(sigma(:)-prevSigma(:));
errorPi = abs(p(:)-prevPi(:));
if (all(errorMu <= eps) && all(errorSigma <= eps) && all(errorPi <= eps))
converged = 1;
end
prevLoglikelihood = loglikelihood;
prevMu = mu;
prevSigma = sigma;
prevPi = p;
end % while
ガウス混合モデルの EM アルゴリズムの Matlab 実装に何か問題がありますか?
以前のトラブル:
問題は、 であるため、対数尤度を使用して収束を確認できないことです-Inf。これは、責任の式でガウスを評価する際に丸められたゼロ値に起因します (E ステップを参照)。
私の EM アルゴリズムの実装がこれまでのところ正しいかどうか教えてもらえますか? そして、丸められたゼロ値の問題をどのように説明するのでしょうか?
これが EM アルゴリズムの私の実装です。
最初に、kmeans を使用してコンポーネントの平均と共分散を初期化しました。
load('data1.mat');
X = Data'; % 300x3 data set
D = size(X,2); % dimension
N = size(X,1); % number of samples
K = 4; % number of Gaussian Mixture components
% Initialization
p = [0.2, 0.3, 0.2, 0.3]; % arbitrary pi
[idx,mu] = kmeans(X,K); % initial means of the components
% compute the covariance of the components
sigma = zeros(D,D,K);
for k = 1:K
sigma(:,:,k) = cov(X(idx==k,:));
end
E-ステップでは、次の式を使用して責任を計算しています
対応するコードは次のとおりです。
% variables for convergence
converged = 0;
prevLoglikelihood = Inf;
prevMu = mu;
prevSigma = sigma;
prevPi = p;
round = 0;
while (converged ~= 1)
round = round +1
gm = zeros(K,N); % gaussian component in the nominator -
% some values evaluate to zero
sumGM = zeros(N,1); % denominator of responsibilities
% E-step: Evaluate the responsibilities using the current parameters
% compute the nominator and denominator of the responsibilities
for k = 1:K
for i = 1:N
% HERE values evalute to zero e.g. exp(-746.6228) = -Inf
gm(k,i) = p(k)/sqrt(det(sigma(:,:,k))*(2*pi)^D)*exp(-0.5*(X(i,:)-mu(k,:))*inv(sigma(:,:,k))*(X(i,:)-mu(k,:))');
sumGM(i) = sumGM(i) + gm(k,i);
end
end
% calculate responsibilities
res = zeros(K,N); % responsibilities
Nk = zeros(4,1);
for k = 1:K
for i = 1:N
res(k,i) = gm(k,i)/sumGM(i);
end
Nk(k) = sum(res(k,:));
end
Nk(k)は、M ステップで指定された式を使用して計算されます。
Mステップ
% M-step: Re-estimate the parameters using the current responsibilities
mu = zeros(K,3);
for k = 1:K
for i = 1:N
mu(k,:) = mu(k,:) + res(k,i).*X(k,:);
sigma(:,:,k) = sigma(:,:,k) + res(k,i).*(X(k,:)-mu(k,:))*(X(k,:)-mu(k,:))';
end
mu(k,:) = mu(k,:)./Nk(k);
sigma(:,:,k) = sigma(:,:,k)./Nk(k);
p(k) = Nk(k)/N;
end
収束を確認するために、次の式を使用して対数尤度が計算されます。
% Evaluate the log-likelihood and check for convergence of either
% the parameters or the log-likelihood. If not converged, go to E-step.
loglikelihood = 0;
for i = 1:N
loglikelihood = loglikelihood + log(sum(gm(:,i)));
end
% Check for convergence of parameters
errorLoglikelihood = abs(loglikelihood-prevLoglikelihood);
if (errorLoglikelihood <= eps)
converged = 1;
end
errorMu = abs(mu(:)-prevMu(:));
errorSigma = abs(sigma(:)-prevSigma(:));
errorPi = abs(p(:)-prevPi(:));
if (all(errorMu <= eps) && all(errorSigma <= eps) && all(errorPi <= eps))
converged = 1;
end
prevLoglikelihood = loglikelihood;
prevMu = mu;
prevSigma = sigma;
prevPi = p;
end % while
最初のラウンドの後、loglikelihood約 700 です。2 番目のラウンドでは、E ステップの-Inf一部のgm(k,i)値がゼロであるためです。したがって、対数は明らかに負の無限大です。
ゼロの値もゼロにsumGM等しいため、行列 muand内のすべての NaN エントリにつながります。sigma
どうすればこの問題を解決できますか? 私の実装に何か問題があるかどうか教えてもらえますか? 何とかMatlabの精度を上げることで解決できるでしょうか？
編集：
gm(k,i) の exp() 項のスケーリングを追加しました。残念ながら、これはあまり役に立ちません。さらに数ラウンド後、まだアンダーフローの問題が発生します。
scale = zeros(N,D);
for i = 1:N
max = 0;
for k = 1:K
Xmu = X(i,:)-mu(k,:);
if (norm(scale(i,:) - Xmu) > max)
max = norm(scale(i,:) - Xmu);
scale(i,:) = Xmu;
end
end
end
for k = 1:K
for i = 1:N
Xmu = X(i,:)-mu(k,:);
% scale gm to prevent underflow
Xmu = Xmu - scale(i,:);
gm(k,i) = p(k)/sqrt(det(sigma(:,:,k))*(2*pi)^D)*exp(-0.5*Xmu*inv(sigma(:,:,k))*Xmu');
sumGM(i) = sumGM(i) + gm(k,i);
end
end
さらに、kmeans は、平均が M ステップで計算される次のラウンドとはまったく異なる平均を初期化することに気付きました。
kmeans:
mu = 13.500000000000000 0.026602138870044 0.062415945993735
88.500000000000000 -0.009869960132085 -0.075177888210981
39.000000000000000 -0.042569305020309 0.043402772876513
64.000000000000000 -0.024519281362918 -0.012586980924762
Mステップ後：
round = 2
mu = 1.000000000000000 0.077230046948357 0.024498886414254
2.000000000000000 0.074260118474053 0.026484346404660
3.000000000000002 0.070944016105476 0.029043085983168
4.000000000000000 0.067613431480832 0.031641849205021
次のラウンドでmuはまったく変わりません。ラウンド2と同じままです。
これは gm(k,i) のアンダーフローが原因ではないでしょうか? 私のスケーリングの実装が間違っているか、アルゴリズムの実装全体がどこかで間違っています:(
編集2
4 ラウンド後NaN、値を取得し、gm をさらに詳しく調べました。1 つのサンプルのみ (および 0.5 係数なし) を見るgmと、すべてのコンポーネントでゼロになります。matlab に入れgm(:,1) = [0 0 0 0]ます。これにより、sumGM はゼロに等しくなります -> NaN は、I がゼロで除算されるためです。で詳細を説明しました
round = 1
mu = 62.0000 -0.0298 -0.0078
37.0000 -0.0396 0.0481
87.5000 -0.0083 -0.0728
12.5000 0.0303 0.0614
gm(:,1) = [11.7488, 0.0000, 0.0000, 0.0000]
round = 2
mu = 1.0000 0.0772 0.0245
2.0000 0.0743 0.0265
3.0000 0.0709 0.0290
4.0000 0.0676 0.0316
gm(:,1) = [0.0000, 0.0000, 0.0000, 0.3128]
round = 3
mu = 1.0000 0.0772 0.0245
2.0000 0.0743 0.0265
3.0000 0.0709 0.0290
4.0000 0.0676 0.0316
gm(:,1) = [0, 0, 0.0000, 0.2867]
round = 4
mu = 1.0000 0.0772 0.0245
NaN NaN NaN
3.0000 0.0709 0.0290
4.0000 0.0676 0.0316
gm(:,1) = 1.0e-105 * [0, NaN, 0, 0.5375]
まず第一に、手段は変更されていないようで、kmeans の初期化と比較して完全に異なります。
そして、すべてのサンプル (ここのような最初のサンプルだけでなく) の出力によると、1 つのガウス成分のみに対応しますgm(:,1)。サンプルはすべてのガウス成分に「部分的に分散」されるべきではありませんか?
EDIT3：
したがって、 mu が変更されない問題は、M-step: の最初の行にあると思いますmu = zeros(K,3);。
アンダーフローの問題を説明するために、現在、ガウスのログを使用しようとしています:
function logPdf = logmvnpdf(X, mu, sigma, D)
Xmu = X-mu;
logPdf = log(1/sqrt(det(sigma)*(2*pi)^D)) + (-0.5*Xmu*inv(sigma)*Xmu');
end
新しい問題は、共分散行列シグマです。Matlab の主張: 警告: 行列が特異値に近いか、スケーリングが不適切です。結果が不正確になる場合があります。
6 ラウンド後、gm (ガウス分布) の虚数を取得します。
更新された E-Step は次のようになります。
gm = zeros(K,N); % gaussian component in the nominator
sumGM = zeros(N,1); % denominator of responsibilities
for k = 1:K
for i = 1:N
%gm(k,i) = p(k)/sqrt(det(sigma(:,:,k))*(2*pi)^D)*exp(-0.5*Xmu*inv(sigma(:,:,k))*Xmu');
%gm(k,i) = p(k)*mvnpdf(X(i,:),mu(k,:),sigma(:,:,k));
gm(k,i) = log(p(k)) + logmvnpdf(X(i,:), mu(k,:), sigma(:,:,k), D);
sumGM(i) = sumGM(i) + gm(k,i);
end
end