Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Optimizing Gaussian mix fit

Subject: Optimizing Gaussian mix fit

From: RAJATH

Date: 11 Feb, 2010 14:13:02

Message: 1 of 5

Hello all,

  The code here is for Gaussian mixture fit using EM algorithm. The code was originally written by Mr Ohad Gal.
  Since this function has been used many times in the program, computation time has increased a lot.
  Please do help me in optimizing it.
  
----------------------------------------------------------------------------------------------------------------
Minor changes like prior multiplication of matrix multiplication are of very little help.
The major culprit here are the two eqns:

  P = C ./ (Ic*sqrt(sig2)) .* exp( -((X*Ir - Ic*u).^2)./(2*Ic*sig2) );
  for m = 1:M
        Z(:,m) = (P(:,m)*t(m))./(P*t(:));
  end
-------------------------------------------------------------------------------------------------------------------
Thanks in advance.

----------------------------------------- Code ------------------------------------------------------------------
% format: [u,sig,t,iter] = fit_mix_gaussian( X,M )
% input: X - input samples, Nx1 vector
% M - number of gaussians which are assumed to compose the distribution
% output: u - fitted mean for each gaussian
% sig - fitted standard deviation for each gaussian
% t - probability of each gaussian in the complete distribution
% iter- number of iterations done by the function

% initialize and initial guesses
N = length( X );
Z = ones(N,M) * 1/M; % indicators vector
P = zeros(N,M); % probabilities vector for each sample and each model
t = ones(1,M) * 1/M; % distribution of the gaussian models in the samples
 u = linspace(min(X),max(X),M); % mean vector
 sig2 = ones(1,M) * var(X) / sqrt(M); % variance vector

 m = 0:1:M-1;
 u = (m+0.5)*max(X)/M;
 sig2 = [1 1 1 1] * max(X).^2/M^2;


C = 1/sqrt(2*pi); % just a constant
Ic = ones(N,1); % - enable a row replication by the * operator
Ir = ones(1,M); % - enable a column replication by the * operator
Q = zeros(N,M); % user variable to determine when we have converged to a steady solution
thresh = 1e-2;
step = N;
last_step = inf;
iter = 0;
min_iter = 10;
MAX_ITER = 20;

% main convergence loop, assume gaussians are 1D

while ((( abs((step/last_step)-1) > thresh) & (step>(N*eps)) ) | (iter<min_iter) )
    
    % E step
    % ========
    Q = Z;
    P = C ./ (Ic*sqrt(sig2)) .* exp( -((X*Ir - Ic*u).^2)./(2*Ic*sig2) );
    for m = 1:M
        Z(:,m) = (P(:,m)*t(m))./(P*t(:));
    end
        
    % estimate convergence step size and update iteration number
    prog_text = sprintf(repmat( '\b',1,(iter>0)*12+ceil(log10(iter+1)) ));
    iter = iter + 1;
    last_step = step * (1 + eps) + eps;
    step = sum(sum(abs(Q-Z)));
% fprintf( '%s%d iterations\n',prog_text,iter );

    % M step
    % ========
    Zm = sum(Z); % sum each column
    Zm(find(Zm==0)) = eps; % avoid dvision by zero
    u = (X')*Z ./ Zm;
    sig2 = sum(((X*Ir - Ic*u).^2).*Z) ./ Zm;
    t = Zm/N;

    if(iter > MAX_ITER) break; end;

end

% plot the fitted distribution
% =============================
sig = sqrt( sig2 );

[y, jy] = plot_mix_gaussian( u,sig,t,X);

Subject: Optimizing Gaussian mix fit

From: Bruno

Date: 12 Apr, 2010 17:54:04

Message: 2 of 5

The code is giving:


??? Error using ==> rdivide
Matrix dimensions must agree.

Error in ==> fit_mix_gaussian at 41
    P = C ./ (Ic*sqrt(sig2)) .* exp( -((X*Ir - Ic*u).^2)./(2*Ic*sig2) );



Any thoughts?
Regards,
Marchesi

Subject: Optimizing Gaussian mix fit

From: RAJATH

Date: 26 Apr, 2010 11:20:25

Message: 3 of 5

"Bruno " <bruno.marchesi@gmail.com> wrote in message <hpvmns$89d$1@fred.mathworks.com>...
> The code is giving:
>
>
> ??? Error using ==> rdivide
> Matrix dimensions must agree.
>
> Error in ==> fit_mix_gaussian at 41
> P = C ./ (Ic*sqrt(sig2)) .* exp( -((X*Ir - Ic*u).^2)./(2*Ic*sig2) );
>
>
>
> Any thoughts?
> Regards,
> Marchesi

Its fine. Code is working properly.

Rajath

Subject: Optimizing Gaussian mix fit

From: Greg Heath

Date: 27 Apr, 2010 05:31:07

Message: 4 of 5

On Apr 26, 7:20 am, "RAJATH " <rajathb...@gmail.com> wrote:
> "Bruno " <bruno.march...@gmail.com> wrote in message <hpvmns$89...@fred.mathworks.com>...
> > The code is giving:
>
> > ??? Error using ==> rdivide
> > Matrix dimensions must agree.
>
> > Error in ==> fit_mix_gaussian at 41
> >     P = C ./ (Ic*sqrt(sig2)) .* exp( -((X*Ir - Ic*u).^2)./(2*Ic*sig2) );
>
> > Any thoughts?
> > Regards,
> > Marchesi
>
> Its fine. Code is working properly.
>
> Rajath

No. I agree with Bruno

??? Error using ==> ./
Matrix dimensions must agree.

Error in ==> C:\MATLAB6p5p1\work\GaussianMixture\gaussmix.m
On line 93 ==> P = C ./ (Ic*sqrt(sig2)) .* exp( -((X*Ir -
Ic*u).^2)./(2*Ic*sig2) );

Please compare the code you are running with the code you posted.

Hope this helps.

Greg

Subject: Optimizing Gaussian mix fit

From: Greg Heath

Date: 28 Apr, 2010 12:03:16

Message: 5 of 5

On Apr 27, 1:31 am, Greg Heath <he...@alumni.brown.edu> wrote:
> On Apr 26, 7:20 am, "RAJATH " <rajathb...@gmail.com> wrote:
>
>
>
>
>
> > "Bruno " <bruno.march...@gmail.com> wrote in message <hpvmns$89...@fred.mathworks.com>...
> > > The code is giving:
>
> > > ??? Error using ==> rdivide
> > > Matrix dimensions must agree.
>
> > > Error in ==> fit_mix_gaussian at 41
> > >     P = C ./ (Ic*sqrt(sig2)) .* exp( -((X*Ir - Ic*u).^2)./(2*Ic*sig2) );
>
> > > Any thoughts?
> > > Regards,
> > > Marchesi
>
> > Its fine. Code is working properly.
>
> > Rajath
>
> No. I agree with Bruno
>
> ??? Error using ==> ./
> Matrix dimensions must agree.
>
> Error in ==> C:\MATLAB6p5p1\work\GaussianMixture\gaussmix.m
> On line 93  ==>     P   = C ./ (Ic*sqrt(sig2)) .* exp( -((X*Ir -
> Ic*u).^2)./(2*Ic*sig2) );
>
> Please compare the code you are running with the code you posted.

This is the culprit:

 sig2 = [1 1 1 1] * max(X).^2/M^2;

instead of

sig2 = ones(1,M) * max(X).^2/M^2;


However, I don't find this to be a reasonable initialization.

Try this (assumes a uniform distribution):

disp('INITIAL CENTERS VECTOR')

du = (max(X)-min(X))/(M+1);
u = linspace(min(X) +0.5*du,max(X)-0.5*du,M);

disp('INITIAL VARIANCE VECTOR')

sig2 = (du^2/12)*ones(1,M);

Hope this helps.

Greg

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us