I have a problem of fitting individual gaussians after a Gaussian mixture model fitting. The means and sigmas for the individual gaussians are obtained from fitgmdist. Please advise.

My code:

if true
  % code
end
options=statset('MaxIter',1000);
GMModel = fitgmdist((log10(newstuff))',2,'Options',options, 'RegularizationValue',0.01);
[edgea centera]=hist(log10(newstuff),100);
gaussPdf= pdf(GMModel,centera');
A = sum(gaussPdf);
gaussPdf = gaussPdf/A;
storing_gausspdf(:,mm)=gaussPdf;
for n = 1:(GMModel.NumComponents),
  mu(n)          = GMModel.mu(n);
  sigma(n)       = sqrt(GMModel.Sigma(1,1,n));
  weight(n)      = GMModel.ComponentProportion(n);
  gaussPdfi(:,n) = weight(n)*normpdf(centera',mu(n),sigma(n))/A;
  ComponentFractions(1:2,mm)=GMModel.ComponentProportion;
  ComponentFractions(3:4,mm)=GMModel.mu;

end

    subplot(4,4,mm);
    plot(centera', gaussPdf, 'r', 'linewidth', 1.25);
    hold on;
    stairs(centera',(edgea./sum(edgea))');
    plot(centera', gaussPdfi(:,1),'-g');
    plot(centera', gaussPdfi(:,2),'-b');

    hold off;

end

The individual gaussians are in green and blue. The global gaussian is in red and the stair representation in blue is raw data. The rsquare for the global fit is 0.99. I also always observe that the second gaussian (in blue) is the one that is always underestimated in amplitude. Any pointers towards solving this will be of great help. Thanks!