How do I use fitgmdist with a set of data that is not in a histogram format

6 visualizaciones (últimos 30 días)
This set of data below gives me a trimodal distributions and I would like to run a statistical analysis on it. When I run fitgmdist it takes the data and turns it into a histogram and then tells me mixing prop and mean. Is there a way to run my data as it looks without Matlab turning it into a histogram and causing a large shift between 0-1?
X = [0.91,3.55,1.20,2.35,1.58,0.75,6.65,10.22,15.30,19.38,6.98,7.79,3.43,3.68,2.04,1.93,1.03,1.14,2.16,5.65,0.48,0.75,0.28,0.20,0.15,0.20,0.11,0.07,0.05];
bar(X, 'stacked')
figure(1)
  3 comentarios
Max Mortensen
Max Mortensen el 19 de En. de 2021
For now I am looking for the mean, and the std deviation of each mode. I found that fitgmdist could give the mean of each mode but when I use fitgmdist matlab changes the data to a histagram, so the modes are no longer centered around 2, 10, and 20 but now they are centered around 0.5, and 5 after that.
Yes I am expecting there to be multiple modes within each distribution I plot using a bar graph.
Adam Danz
Adam Danz el 19 de En. de 2021
With so few data and only 1 "repetition" of data for each x-value, curve fitting will be difficult.
An alternative would be to split the bars into 3 groups, one for each mode. That could be done with differentiation and smoothing or perhaps with findpeaks. Then fit each mode; fitdist could come in handy for that.

Iniciar sesión para comentar.

Respuesta aceptada

Adam Danz
Adam Danz el 19 de En. de 2021
Editada: Adam Danz el 20 de En. de 2021
You could fit the distribution to a tri-modal mixture of gaussians but very little data and a very narrow peaks (1-2 bars), you'll need to set bounds to the parameters. Even then you won't get great fits.
The demo below walks through the use of fmincon and fits the data to a tri-modal mixture of gaussians. Its loosely based on a fitting proceedure I recently used to fit a bi-modal distribution (with much more data) and you'll likely need to tweek it. The results are highly sensitive to the bounds. See inline comments for details.
% Original data to fit
X = [0.91,3.55,1.20,2.35,1.58,0.75,6.65,10.22,15.30,19.38,6.98,7.79,3.43,3.68,2.04,1.93,1.03,1.14,2.16,5.65,0.48,0.75,0.28,0.20,0.15,0.20,0.11,0.07,0.05];
xvals = 1:numel(X);
% Define a tri-modal gaussian with 3 independent peaks with 11 parameters:
% b(1) vertical off
% b(2) slope
% b(3) amp1
% b(4) mu1
% b(5) std1
% b(6) amp2
% b(7) mu2
% b(8) std2
% b(9) amp3
% b(10) mu3
% b(11) std3
mixGausMdl = @(b,x) b(1) + b(2) * x(:, 1) + b(3) * exp( -(x(:, 1) - b(4)).^2/(2*b(5)^2)) + b(6) * exp(-(x(:, 1) - b(7)).^2/(2*b(8)^2)) + b(9) * exp(-(x(:, 1) - b(10)).^2/(2*b(11)^2));
% Set initial guesses, lower and upper bounds to each parameter
initGuess = [0, 0.1, 3, 2, 1.06, 19, 10, 5.88 4 20 1];
lowbounds = [0 -inf 0 0 .1 0 0 0.1 0 20 .1];
highbounds = [0 inf 10 10 10 inf 25 10 6 30 10];
% Poisson log-likelihood; fmincon finds the minimum of this function. Thus,
% will find the max of the poisson loglikelihood.
poissonLogLikelihood = @(r,lambda) r(:)'*log(lambda(:)) - sum(lambda);
% loss functions for fitting
nllFcn = @(params) -poissonLogLikelihood(X, mixGausMdl(params, xvals(:)));
% Do fitting
opts = optimset('LargeScale', 'off', 'MaxIter', 50000, 'Display', 'off');
[paramEst, nll, exitflag] = fmincon(nllFcn, initGuess, [], [], [], [], lowbounds, highbounds,[],opts);
% show results
T = table((1:numel(paramEst))',paramEst(:), initGuess(:), lowbounds(:), highbounds(:), ...
'VariableNames', {'idx','Estimates','initGuess','lowBound','HiBound'}, ...
'RowNames',{'VertOff','slope','amp1','mu1','std1','amp2','mu2','std2','amp3','mu3','std3'});
disp(T)
idx Estimates initGuess lowBound HiBound ___ _________ _________ ________ _______ VertOff 1 0 0 0 0 slope 2 0.0073667 0.1 -Inf Inf amp1 3 2.5688 3 0 10 mu1 4 2.3787 2 0 10 std1 5 1.1271 1.06 0.1 10 amp2 6 13.232 19 0 Inf mu2 7 9.9742 10 0 25 std2 8 2.4188 5.88 0.1 10 amp3 9 3.6307 4 0 6 mu3 10 20 20 20 30 std3 11 1.0583 1 0.1 10
% Plot results
clf()
bar(xvals, X, 'FaceColor',[.8 .8 .8])
hold on
% Plot fit
xx = linspace(xvals(1),xvals(end),100);
yy = mixGausMdl(paramEst,xx(:));
plot(xx,yy,'r-','LineWidth',5)
% Show each guassian separately
gaus = @(x,mu,sig,amp,vo)amp*exp(-(((x-mu).^2)/(2*sig.^2)))+vo;
y1 = gaus(xx, paramEst(4), paramEst(5), paramEst(3), paramEst(1));
y2 = gaus(xx, paramEst(7), paramEst(8), paramEst(6), paramEst(1));
y3 = gaus(xx, paramEst(10), paramEst(11), paramEst(9), paramEst(1));
plot(xx,y1, 'c:', 'LineWidth',3)
plot(xx,y2, 'k:', 'LineWidth',3)
plot(xx,y3, 'b:', 'LineWidth',3)

Más respuestas (0)

Community Treasure Hunt

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

Start Hunting!

Translated by