Clear Filters
Clear Filters

How could I fit a mixture of gaussians to 1D data?

31 views (last 30 days)
I have a data and i want to fit it by a mixture of gaussian, but I didn't know the existing number of gaussians.
I use this code:
[ndata text alldata] = xlsread('battery2.xlsx','li ion');
[R2,C2]=size(alldata);
Life = ndata(:,72);
A=round(Life);
[valu,idxC, idxV] = unique(A);
n = accumarray(idxV,1);
m= [n idxC] ;
y = [valu n];
binWidth = 1; %2 1
binCtrs = 0:binWidth:172;
figure
count= hist(Life,binCtrs);
hist(Life,binCtrs,'k');
%%rescale
counts = hist(Life,binCtrs);
nu = length(Life);
prob = counts / (nu * binWidth);
figure
plot(binCtrs,prob,'ko');hold on
k=2
paramEsts= gmdistribution.fit(Life,k)
xgrids = linspace(0,173,100);
y1 = gaussian1D(xgrids, paramEsts.mu(1), paramEsts.Sigma(1));
y2 = gaussian1D(xgrids, paramEsts.mu(2), paramEsts.Sigma(2));
% line(xgrids,y1,'color','r', 'LineWidth',2)
axis([-10 175 0 0.02]);
C=5;
w= [paramEsts.PComponents(1) paramEsts.PComponents(2)];
x=xgrids;
values=binCtrs;
[mu_est, sigma_est, w_est, counter, difference] = gaussian_mixture_model(values, C, 1.0e-3);
mu_est'
sigma_est'
w_est'
% compare empirical data to estimated distribution
p1_est = w_est(1) * norm_density(x, mu_est(1), sigma_est(1));
p2_est = w_est(2) * norm_density(x, mu_est(2), sigma_est(2));
p3_est = w_est(3) * norm_density(x, mu_est(3), sigma_est(3));
p4_est = w_est(4) * norm_density(x, mu_est(4), sigma_est(4));
p5_est = w_est(5) * norm_density(x, mu_est(5), sigma_est(5));
plot2 = plot(x, p1_est+p2_est+p3_est+p4_est+p5_est, 'r--', 'linewidth', 2);
but it didn't fit, how i could change the code to fit the data?
  2 Comments
Amr Hashem
Amr Hashem on 6 Mar 2017
I try this also:
k=2;
paramEsts= gmdistribution.fit(Life,k)
xgrids = linspace(0,173,100);
y1 = gaussian1D(xgrids, paramEsts.mu(1), paramEsts.Sigma(1));
y2 = gaussian1D(xgrids, paramEsts.mu(2), paramEsts.Sigma(2));
line(xgrids,y1,'color','r', 'LineWidth',2)
plot(xgrids, y1, 'r-');
plot(xgrids, y2, 'r-');
But, it also didn't work.

Sign in to comment.

Accepted Answer

Amr Hashem
Amr Hashem on 26 Mar 2017
Edited: Amr Hashem on 26 Mar 2017
Finally, I got it :)
clear all
clc
close all
tic
%%Read the raw data
[ndata text alldata] = xlsread('battery2.xlsx','li ion'); % open excel sheet 101 93
[R2,C2]=size(alldata);
% Col 72 stores for each event after how many months the event occurred (starting from the production date)
Life = ndata(:,72); % 72 coulumn of device age in months
binWidth = 1; %width of the bar equal 1
binCtrs = 0:binWidth:172; % the bars start from 0 to 172. d=729 w=104 93 m=24
figure
count= hist(Life,binCtrs);
hist(Life,binCtrs,'k'); % b k m %hold on; plot(n,count,'ro')
% grid;
set(gca,'XTick',[0:24:144]);
xlabel('Time (Months)');
ylabel('Failure count');
title('Failure count of Lithium-ion (Li-ion) batteries from 2010 through 2014');
AAY=legend({'Total=5618'});
axis([-5 144 0 100])
%%normalize the failure values to be from 0 to 1.
counts = hist(Life,binCtrs);
nu = length(Life);
prob = counts / (nu * binWidth);
figure
plot(binCtrs,prob,'ko');hold on
%%finding the number of components
%using minimum value of AIC
AIC = zeros(1,5);
obj = cell(1,5);
for Kk = 1:5
obj{Kk} = gmdistribution.fit(Life,Kk);
AIC(Kk)= obj{Kk}.AIC;
end
[minAIC,numComponents] = min(AIC);
numComponents
numComponents=4;
paramEsts= gmdistribution.fit(Life,numComponents)
%mu of mixture of gaussians of 4 components
MU=[paramEsts.mu(1);paramEsts.mu(2);paramEsts.mu(3);paramEsts.mu(4)];
SIGMA=cat(3,[paramEsts.Sigma(1)],[paramEsts.Sigma(2)],[paramEsts.Sigma(3)],[paramEsts.Sigma(4)]);
PPp=[paramEsts.PComponents(1),paramEsts.PComponents(2),paramEsts.PComponents(3),paramEsts.PComponents(4)];
objA = gmdistribution(MU,SIGMA,PPp);
xgridss=transpose(linspace(0,173,100));
plot(xgridss,pdf(objA,xgridss),'r-','linewidth',2)
toc

More Answers (2)

Amr Hashem
Amr Hashem on 25 Mar 2017
Edited: Amr Hashem on 25 Mar 2017
I find a code which help to determine the number of components:
%using minimum value of AIC
AIC = zeros(1,7);
obj = cell(1,7);
for Kk = 1:7
obj{Kk} = gmdistribution.fit(Data,Kk);
AIC(Kk)= obj{Kk}.AIC;
end
[minAIC,numComponents] = min(AIC);
numComponents
https://stackoverflow.com/questions/10359866/gaussian-mixture-model-components/18530779#18530779?newreg=93f752c3d2e94032bf9d30433218f245

Amr Hashem
Amr Hashem on 9 Apr 2017
Edited: Amr Hashem on 9 Apr 2017
Another authorized answer
https://www.mathworks.com/help/stats/examples/fitting-custom-univariate-distributions.html

Community Treasure Hunt

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

Start Hunting!