Code covered by the BSD License  

Highlights from
Gaussian Mixture Model

Gaussian Mixture Model

by

 

The code finds out the parameters of a gaussian mixture model by Expectation-Maximization Algorithm.

GMmodel(x,no_gaus_distr)
%Developed by: Ravi Shankar, 3rd year B.tech, IIT Guwahati.
%------------------------------------------------------------------//-------------------------------------------------------------------------------------%
% Notation used :-
% P(no_gaus_distr) - number of gaussian distribution
% N(num_rows) - number of training data
% D(dimension) - dimension of each data

% function to find the mixing ratio, covariance matrices and means of the
% gaussian distribution using reiterative procedure or Expectation 
% Maximization algorithm.


function [Wt,Meanmat,Cov] = GMmodel(x,no_gaus_distr)

% the above expression is the definition of the function GMmodel which 
% takes in the training data and the number of gaussian distribution 
% to model the data as input and the output consist of the three matrices, 
% first one is mixing ratios whch is P * N matrix and contains the 
% fraction of datapoints belonging to each gaussian distribution, the
% second one is mean matrix that contains the means of each gaussian
% distribution in P * D matix and the third one is covariance matrix which
% is D * D * P matrix and has covariance of each distribution in D * D
% matrix.

%-----------------------------------------------------//---------------------------------------------------------------------------------------------------%
         % initialization of the Meanmat, Wt and Cov is done in this block
         % of code. Initially K-mean clustering technique is being applied
         % to find the weights, means and covariance matrices.

         sizemat = size(x);
         tot_rows = sizemat(1);
         dimension = sizemat(2);
         Wt = zeros(no_gaus_distr,1);
         Cov = zeros(dimension,dimension,no_gaus_distr);
         Meanmat = kmclust(x,no_gaus_distr);
         
         % kmclust is a function that does the K-mean clustering out of
         % many data sets and returns the mean of each cluster.
         
         Eff_no_pnts = zeros(no_gaus_distr,1);
         indexcol = zeros(tot_rows,1);
         
         % the for loop used next is for indexing each row to the mean from
         % which its distance is minimum and the output is stored in
         % indexcol matrix.
    
         
         for currow = 1:tot_rows
                threshold = 10^20;
                for curmean = 1:no_gaus_distr
                    sqdist = 0;
                    for curdim = 1:dimension
                        sqdist = sqdist + ((x(currow,curdim) - Meanmat(curmean,curdim))^2);
                    end
                    dist = sqdist^(0.5);
                    if (dist <= threshold)
                       indexcol(currow,1) = curmean;
                       threshold = dist;
                    end
                end
         end
         
         % the covariance matrix is initialised first with values that are
         % found using the datapoints which are belonging to the same
         % cluster.
         
         for loopvar1 = 1:no_gaus_distr
             count = 0;
             tempcov = zeros(dimension,dimension);
             for loopvar2 = 1:tot_rows
                    if (indexcol(loopvar2,1) == loopvar1)
                        count = count + 1;
                        tempcov = tempcov + ((x(loopvar2,:)-Meanmat(loopvar1,:))')*(x(loopvar2,:)-Meanmat(loopvar1,:));
                    end
             end
             Wt(loopvar1,1) = count/tot_rows;
             tempcov = tempcov./count;
             %det(tempcov)
             Cov(:,:,loopvar1) = tempcov;
             %if (loopvar1 == 1)
             %    Cov(:,:) = tempcov(:,:);
             %else
             %    Cov = [Cov;tempcov];
             %end
         end
         %save('covariance first.mat','Cov');
         %exit (0);
         %tempvar2 = 0;
         %for currow = 1:tot_rows
         %   tempvar = 0;
             %cov_var = 1;
         %   for curdistr = 1:no_gaus_distr
         %        %tempcov = Cov(cov_var:(dimension*curdistr),:);
         %        tempcov = Cov(:,:,curdistr);
                 %cov_var = cov_var + dimension; 
         %        exp_part = exp(-0.5*(x(currow,:)-Meanmat(curdistr,:))*(1/tempcov)*((x(currow,:)-Meanmat(curdistr,:))'));
         %        tempvar = tempvar + ((Wt(curdistr,1))*(1/((det(tempcov))^(0.5)))*exp_part);
         %    end
         %    tempvar2 = tempvar2 + log(tempvar);
         %end
         %likelihoodnew = -((tot_rows*dimension)/2)*log(2*pi) + tempvar2;
         %flag_new = exp(likelihoodnew*(10^-5));
         likelihoodold = -10^30;
         likelihoodnew = -10^25;
         flag = 10^(-10);
         %flag_new = 10000;
         %flag_old = 1000;
         iteration = 0;
         
         % responsibilities matrix is nothing but the responsibiliets of each
         % gaussian to take each point and it is a N * P matrix one row for
         % each data set and cloumn for each gaussian distribution.
         
         responsibilities = zeros(tot_rows,no_gaus_distr);
         
         
%------------------------------------------------------//-------------------------------------------------------------------------------------------%
         % loop starts here 

         % the loop runs until the difference between successive iteration
         % becomes less than 0.001 times the value of log of the
         % probability of observing the training set of data
         % while (abs(likelihoodold-likelihoodnew)>0.001 && iteration <=500)
         while (((likelihoodnew-likelihoodold)>flag || iteration<25) && iteration<200)
               iteration = iteration + 1 %abs(likelihoodold-likelihoodnew)
               likelihoodold = likelihoodnew;
               %flag_old = flag_new;
 
               
               
%----------------------------------------------------//----------------------------------------------------------------------------------------------%               
               % this chunk involves finding the responsibilities of each
               % gaussian distribution for each point in N*P matrix.
               % inv is a predefined function of matlab that finds the
               % inverse of the matrix given in its argument. This is
               % essentially the most time consuming part if the dimension
               % of feature is very large.
               
               for currow = 1:tot_rows 
                   %cov_var = 1;
                   evidence = 0;
                   for loopvar1 = 1:no_gaus_distr
                       %tempcov = Cov(cov_var:(dimension*loopvar1),:);
                       tempcov = Cov(:,:,loopvar1);
                       %cov_var = cov_var + dimension;
                       exp_part = exp(-0.5*(x(currow,:)-Meanmat(loopvar1,:))*(inv(tempcov))*((x(currow,:)-Meanmat(loopvar1,:))'));
                       denominator = (Wt(loopvar1,1))*(1/((det(tempcov))^(0.5)))*exp_part;
                       responsibilities(currow,loopvar1) = denominator;
                       evidence = evidence + denominator;
                   end
                   responsibilities(currow,:) = responsibilities(currow,:)./evidence;
               end
               
               %str = 'finished first job'
               
%-------------------------------------------------//--------------------------------------------------------------------------------------------------%



               % this chunk calculates new means, weights and effective 
               % number of points.   

               for loopvar1 = 1:no_gaus_distr
                      tempvar = zeros(1,dimension);
                      tempvar2 = 0;
                      for currow = 1:tot_rows
                          tempvar = tempvar + (responsibilities(currow,loopvar1).*x(currow,:));
                          tempvar2 = tempvar2 + responsibilities(currow,loopvar1);
                      end
                      Meanmat(loopvar1,:) = tempvar./tempvar2;
                      Eff_no_pnts(loopvar1,1) = tempvar2;
                      Wt(loopvar1,1) = Eff_no_pnts(loopvar1,1)/tot_rows;
               end
               
              % str = 'finished second job'
              
%-----------------------------------------------//------------------------------------------------------------------------------------------------------%



              % this chunk of code is for recalculation of covariance
              % matrix. By recalculation it means the covariance matrix is
              % filled with new data sets that makes use of the
              % responsibilities calculated earlier.
              
              %cov_var = 1;
              for loopvar1 = 1:no_gaus_distr
                  tempcov = zeros(dimension,dimension);
                  tempvar2 = 0;
                  for currow = 1:tot_rows
                         tempcov = tempcov + responsibilities(currow,loopvar1)*(((x(currow,:))')*(x(currow,:)));
                         tempvar2 = tempvar2 + responsibilities(currow,loopvar1);
                  end
                  Cov(:,:,loopvar1) = (tempcov/tempvar2) - (Meanmat(loopvar1,:)')*(Meanmat(loopvar1,:));
                  %Cov(cov_var:(dimension*loopvar1),:) = tempcov./Eff_no_pnts(loopvar1,1);
                  %cov_var = cov_var + dimension;
              end
              
              % str = 'finished third job'
              
%---------------------------------------------//---------------------------------------------------------------------------------------------------------%



              % this chunk calculates the log likelihood in order to 
              % compare the previous value with the current value. 
              
              tempvar2 = 0;
              for currow = 1:tot_rows
                     tempvar = 0;
                     %cov_var = 1;
                     for curdistr = 1:no_gaus_distr
                         tempcov = Cov(:,:,curdistr);
                         %tempcov = Cov(cov_var:(dimension*curdistr),:);
                         %cov_var = cov_var + dimension;
                         exp_part = exp(-0.5*(x(currow,:)-Meanmat(curdistr,:))*(inv(tempcov))*((x(currow,:)-Meanmat(curdistr,:))'));
                         tempvar = tempvar + ((Wt(curdistr,1))*(1/((det(tempcov))^(0.5)))*exp_part);
                     end
                     tempvar2 = tempvar2 + log(tempvar);
              end
       
              
              % str = 'finished fourth job'
                
              likelihoodnew = -((tot_rows*dimension)/2)*log(2*pi) + tempvar2
              if (iteration==3)
                  flag = -1*(likelihoodnew/(10^(5)))
              end
              likelihoodnew-likelihoodold
              %flag_new = exp(likelihoodnew*(10^-5));
              %exp(abs(flag_old-flag_new))
         end   
%-------------------------------------------//----------------------------------------------------------------------------------------------------------%         

Contact us