200个论坛币紧急求助利用MATLAB实现GMM模型的修改

1. function varargout = gmm(X, K_or_centroids)
   
2. % ============================================================
   
3. % Expectation-Maximization iteration implementation of
   
4. % Gaussian Mixture Model.
   
5. %
   
6. % PX = GMM(X, K_OR_CENTROIDS)
   
7. % [PX MODEL] = GMM(X, K_OR_CENTROIDS)
   
8. %
   
9. %  - X: N-by-D data matrix.
   
10. %  - K_OR_CENTROIDS: either K indicating the number of
    
11. %       components or a K-by-D matrix indicating the
    
12. %       choosing of the initial K centroids.
    
13. %
    
14. %  - PX: N-by-K matrix indicating the probability of each
    
15. %       component generating each point.
    
16. %  - MODEL: a structure containing the parameters for a GMM:
    
17. %       MODEL.Miu: a K-by-D matrix.
    
18. %       MODEL.Sigma: a D-by-D-by-K matrix.
    
19. %       MODEL.Pi: a 1-by-K vector.
    
20. % ============================================================
    
21. 
    
22.     threshold = 1e-15;
    
23.     [N, D] = size(X);
    
24. 
    
25.     if isscalar(K_or_centroids)
    
26.         K = K_or_centroids;
    
27.         % randomly pick centroids
    
28.         rndp = randperm(N);
    
29.         centroids = X(rndp(1:K), :);
    
30.     else
    
31.         K = size(K_or_centroids, 1);
    
32.         centroids = K_or_centroids;
    
33.     end
    
34. 
    
35.     % initial values
    
36.     [pMiu pPi pSigma] = init_params();
    
37. 
    
38.     Lprev = -inf;
    
39.     while true
    
40.         Px = calc_prob();
    
41. 
    
42.         % new value for pGamma
    
43.         pGamma = Px .* repmat(pPi, N, 1);
    
44.         pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);
    
45. 
    
46.         % new value for parameters of each Component
    
47.         Nk = sum(pGamma, 1);
    
48.         pMiu = diag(1./Nk) * pGamma' * X;
    
49.         pPi = Nk/N;
    
50.         for kk = 1:K
    
51.             Xshift = X-repmat(pMiu(kk, :), N, 1);
    
52.             pSigma(:, :, kk) = (Xshift' * ...
    
53.                 (diag(pGamma(:, kk)) * Xshift)) / Nk(kk);
    
54.         end
    
55. 
    
56.         % check for convergence
    
57.         L = sum(log(Px*pPi'));
    
58.         if L-Lprev < threshold
    
59.             break;
    
60.         end
    
61.         Lprev = L;
    
62.     end
    
63. 
    
64.     if nargout == 1
    
65.         varargout = {Px};
    
66.     else
    
67.         model = [];
    
68.         model.Miu = pMiu;
    
69.         model.Sigma = pSigma;
    
70.         model.Pi = pPi;
    
71.         varargout = {Px, model};
    
72.     end
    
73. 
    
74.     function [pMiu pPi pSigma] = init_params()
    
75.         pMiu = centroids;
    
76.         pPi = zeros(1, K);
    
77.         pSigma = zeros(D, D, K);
    
78. 
    
79.         % hard assign x to each centroids
    
80.         distmat = repmat(sum(X.*X, 2), 1, K) + ...
    
81.             repmat(sum(pMiu.*pMiu, 2)', N, 1) - ...
    
82.             2*X*pMiu';
    
83.         [dummy labels] = min(distmat, [], 2);
    
84. 
    
85.         for k=1:K
    
86.             Xk = X(labels == k, :);
    
87.             pPi(k) = size(Xk, 1)/N;
    
88.             pSigma(:, :, k) = cov(Xk);
    
89.         end
    
90.     end
    
91. 
    
92.     function Px = calc_prob()
    
93.         Px = zeros(N, K);
    
94.         for k = 1:K
    
95.             Xshift = X-repmat(pMiu(k, :), N, 1);
    
96.             inv_pSigma = inv(pSigma(:, :, k)+eye(D));
    
97.             tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);
    
98.             coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma));
    
99.             Px(:, k) = coef * exp(-0.5*tmp);
    
100.         end
     
101.     end
     
102. end

复制代码

MATLAB代码如上，第96行在执行时会出现
警告: 矩阵接近奇异值，或者缩放错误。结果可能不准确。RCOND =  1.718855e-21。
> In gmm>calc_prob at 96
  In gmm at 40
  In gmm_accuracy at 4


代码是从这个网页获取的：http://blog.pluskid.org/?p=39&cpage=1#comments
关于这个问题也有相关的分析和解释：http://freemind.pluskid.org/mach ... ariance-estimation/

我不要求搞懂原理，只需要代码能够运行就好！
200个论坛币悬赏

楼主你好~我目前也遇到了这个问题。经过改错后已经可以运行：
问题原因：代码中没有处理singular问题
改正： 加几句代码，处理singular问题
把你第96行的function换成这个就好
function Px = calc_prob()
        Px = zeros(N, K);
        for k = 1:K
            Xshift = X-repmat(pMiu(k, : ), N, 1);%x-u
            lemda=1e-5;
            conv=pSigma(:, :, k)+lemda*diag(diag(ones(D)));%这里处理singular问题，为协方差矩阵加上一个很小lemda*I
            inv_pSigma = inv(conv);%协方差的逆
            tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);%(X-U_k)sigma.*(X-U_k),tmp是个N*1的向量
            coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma));%前面的参数
            Px(:, k) = coef * exp(-0.5*tmp);%把数据点 x 带入到 Gaussian model 里得到的值
        end
    end
end

200。。请拿来啊~~

改成pinv
就行了
给我两百吧！

magicsun 发表于 2015-2-8 13:52
改成pinv
就行了
给我两百吧！

试过，
错误使用 svd
SVD 的输入不能包含 NaN 或 Inf。

出错 pinv (line 18)
[U,S,V] = svd(A,'econ');

出错 gmm/calc_prob (line 96)
            inv_pSigma = pinv(pSigma(:, :, k)+eye(D));

出错 gmm (line 40)
        Px = calc_prob();

悬赏还可以再增加，只要能解决问题！

我跑DSGE的过程中遇到了同样的问题，不知道是MATLAB版本的问题还是咋滴啊

我在跑这个代码的时候也遇到了这个文题，不知道解决了没有啊