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

fit gaussian mixturegmdistribution.fitmixture of gaussian

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?

Best Answer

  • 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