OccupancyDetection/GMM_EM.m at master · shounakgithub/OccupancyDetection · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
function [ im ] = GMM_EM( A, M )
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here

clc
    clear all
    close all
M=2
A = [-1.83939462300000;-2.68764102500000;-0.777492961000000;-0.0197987390000000;-3.18466871700000;-1.04872330300000;0.108827286100000;4.05436342690000;2.42814115150000;3.97864172220000;0.500954996100000;0.716370675900000;0.655569174100000;0.0579992681000000;1.70449330260000;1.93268780800000;-0.0133452090000000;0.821949999500000;2.09859394240000;1.27716471270000;-1.44246710600000;-0.846637941000000;-1.19439465400000;0.702298623100000;-5.05093694200000;-3.75297322400000;-4.27430058300000;-3.30623626100000;-3.01960620500000;-2.56127845600000;3.68724341220000;3.49602522030000;5.28212707270000;3.95694410450000;2.44718638520000;-1.47612602200000;-1.55373247000000;-0.847133235000000;-1.40562192200000;-2.38363250300000;-0.174278648000000;-1.38829996100000;1.83097486020000;3.94853902420000;-2.01786597900000;0.570787939700000;-0.711407157000000;-2.86440678100000;0.560318583600000;-1.34504596200000;-1.66149633700000;-2.56052512300000;-1.19076544700000;-1.49338616100000;-2.04630789900000;-2.61374415800000;-2.04813183300000;-1.07456719900000;-0.517888851000000;1.69896802710000;-2.03365621700000;-0.789543109000000;0.707898775500000;-0.0469421340000000;2.70751081170000;-2.59667787300000;-0.855045858000000;-0.524784310000000;-1.53658678900000;-1.73409792600000;-0.377194127000000;1.11727121300000;-1.83536559600000;-1.77077017000000;-0.137102969000000;0.862157619300000;-0.297303975000000;0.0925337075000000;1.30374914350000;2.49182502410000;0.125706723600000;2.46630033060000;2.83497399130000;3.62266768200000;1.45046233330000;-0.666693667000000;1.08729966490000;2.19428489390000;2.61389727210000;1.30882213710000;-0.224822874000000;1.93275437530000;0.376196612600000;-0.0766632450000000;2.19642364500000;1.88760376610000];


[m n dim]=size(A);

PC1 = A(:,1);

minPC1=min(PC1);
shiftPC1=PC1+abs(minPC1);

mPC1=max(PC1);
cPC1=std(PC1)^2;

% PC2 = A(:,2);
%
% minPC2=min(PC2);
% shiftPC2=PC2+abs(minPC2);
%
% mPC2=max(PC2);
% cPC2=std(PC2)^2;
%
% xi = [PC1,PC2];
xi =PC1;

[N dim]=size(A);

Pi(1,1:M)=1/M;

for i=1:M
    %mu(i,:) = [mPC1*rand(1) mPC2*rand(1)];
    mu(i)=mPC1*rand(1);
    %covar(:,:,i) = [cPC1 0; 0 cPC2];
    covar(:,:,i)=cPC1;
    p(i)=rand(1);
end
P=sum(p);
p=p./P;


for i=1:1
    tot(:,1)=zeros(1,N);
    for j=1:M
        %tot(:,1)=tot(:,1)+(mvnpdf(xi,mu(j,:),covar(:,:,j))*p(j));
        tot(:,1)=tot(:,1)+(mvnpdf(xi,mu(j),covar(:,:,j))*p(j));
    end

   for j=1:M
       %w(j,:) = (mvnpdf(xi,mu(j,:),covar(:,:,j))*p(j))./tot;
        w(j,:) = (mvnpdf(xi,mu(j),covar(:,:,j))*p(j))./tot; %bi --> Posterior
        p(j)=sum(w(j,:))/length(w(j,:));% Prior
        mu(j,:) = (w(j,:)*xi/sum(w(j,:)));
        %covar(:,:,j) = diag([cov(PC1.*w(j,:)'),cov(PC2.*w(j,:)')]);
        covar(:,:,j) = diag([cov(PC1.*w(j,:)')]);
        Prob(i,j)=p(j);
        if i > 1
            if Prob(i,j) == Prob(i-1)
                a(j)=1;
            end
        end

   end
   if i > 1
       p(M)=1-sum(p(1):p(M-1));
       if sum(a) == M
           ans=i-1
           break
       end
   end

end


variable=sum(w,1);
for j=1:length(w)
    for i=1:M
       var1=w(i,j);
       var2=variable(i)-var1;
       if var1*(M-1) > var2
          im(j)= (i-1)*(1/(M-1));
       end
    end
end

end