Clear Filters
Clear Filters

How to create a Gaussian mixed model (GMM) with many components?

3 views (last 30 days)
Hi, I have a series of 1D data which is in the histogram. I want to use a GMM equivalent model for this data. I have written the following code, can you guide me in the continuation of this code?
clear all ; close all ; clc
warning off ;
load Data
% type 1: section [1,37,78,187,249,288,364,396,451,532,563,607]
X = 1:607 ;
D = D(X) ;
ft = fittype( 'smoothingspline' );
opts = fitoptions( 'Method', 'SmoothingSpline' );
opts.SmoothingParam = 0.015;
[SF, ~] = fit( X', D, ft, opts );
D = SF(X) ;
D = normalize(D-min(D),'norm') ;
%% Get Data from bar
figure ;
B = bar(X,D,'EdgeColor','none') ; hold on
Y = B.YData ;
%% peak clusters
[P data] = findpeaks(-Y+max(Y),'MinPeakWidth',4.5,'MinPeakDistance',10) ; % 'WidthReference','halfheight'
plot(data,P,'vg') ; data = data' ;
indexSep = [1 data' length(X)] ;
K = length(indexSep)-1 ;
index = zeros(size(X)) ;
for i=1:K
index(indexSep(i):indexSep(i+1)) = i ;
end
index = index' ;
%% Show Clusters
G = {'r*','b*','g*','k*','c*','y*','m*','ro','bo','go','ko','co','yo'} ;
for i = 1:K
plot(X(index==i),Y(index==i),[G{i}])
end
%% Does it need weight?
%% fit components
Model = 'a*exp(-((x-b)/c)^2)' ;
for i = 1:K
opts = fitoptions( 'Method', 'NonlinearLeastSquares','MaxIter', 500,'Robust','Bisquare');
needindex = find(maxk(Y(index==i),1)==Y(index==i)) ;
opts.StartPoint = [max(Y(index==i)) indexSep(i)+X(needindex(1)) (max(X(index==i))-min(X(index==i)))/2];
opts.Lower = [max(Y(index==i))-5 opts.StartPoint(2)-10 0] ;
opts.Upper = [max(Y(index==i))+5 opts.StartPoint(2)+10 Inf] ;
f = fit(X(index==i)',Y(index==i)',Model,opts) ;
AICoeff = coeffnames(f) ;
AIValue = cellstr(num2str(coeffvalues(f)')) ;
fs = replace(Model,AICoeff,AIValue) ;
fs = replace(fs,{'--','-+','+-','++','*','^'},{'+','-','-','+','.*','.^'}) ;
eval(['F' num2str(i) ' = fs ;'])
eval(['f' num2str(i) ' = inline(F' num2str(i) ') ;'])
% need fit with full sum(FK)
end
%% Show
newx = -400:1200 ;
for i = 1:K
eval(['newy(' num2str(i) ',:) = f' num2str(i) '(newx) ;'])
plot(newx,newy(i,:),'LineWidth',1,'Color','k')
end
axis([X(1)-10 X(end)+10 0 max(Y)])
%% In this section I want to see an equivalent GMM that matches this data
options = statset('MaxIter', 500);
gm = gmdistribution.fit(Y',K-1, 'Options', options);
gm = fitgmdist(Y',K) ;
y = pdf(gm,newx');
plot(newx,y,'LineWidth',2,'Color','r')
Indeed, I need a way to calculate GMM coefficients(correct gm.mu,gm.Sigma,\pi_{k}) and match each component.
  1 Comment
Mohamed Hassan
Mohamed Hassan on 22 Jun 2022
add this part to your code:
ComponentMeans = gm.mu;
ComponentCovariance = gm.Sigma;
MixtureProportions = gm.ComponentProportion
I have a question please. Could I apply your code to an image? I wanna show the gaussian distributions on the histogram.

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!