Plotting a two peak model as a fitted line on a probplot.

7 views (last 30 days)
Hi,
I am trying to plot a two peak distribution model on a probplot.
and then using this example from MathWorks: Probability plots - MATLAB probplot - MathWorks United Kingdom
However, as my distribution uses a custom mle function is doesn't seem to work as desired. This is the example code I am working with.
%example from fist link above
clear; clc;
rng(10) % For reproducibility
x3 = [trnd(20,1,50) trnd(4,1,100)+3];
histogram(x3)
pdf_normmixture = @(x3,p,mu1,mu2,sigma1,sigma2) ...
p*normpdf(x3,mu1,sigma1) + (1-p)*normpdf(x3,mu2,sigma2);
pStart = .5;
muStart = quantile(x3,[.25 .75]);
sigmaStart = sqrt(var(x3) - .25*diff(muStart).^2);
start = [pStart muStart sigmaStart sigmaStart];
lb = [0 -Inf -Inf 0 0];
ub = [1 Inf Inf Inf Inf];
options = statset('MaxIter',300,'MaxFunEvals',600);
paramEsts = mle(x3,'pdf',pdf_normmixture,'Start',start, ...
'LowerBound',lb,'UpperBound',ub,'Options',options);
xgrid = linspace(1.1*min(x3),1.1*max(x3),200);
pdfgrid = pdf_normmixture(xgrid, ...
paramEsts(1),paramEsts(2),paramEsts(3),paramEsts(4),paramEsts(5));
figure(1)
histogram(x3,'Normalization','pdf')
hold on
plot(xgrid,pdfgrid,'-')
hold off
xlabel('x3')
ylabel('Probability Density')
legend('Sample Data','Fitted pdf','Location','best')
%end of first example
%example from second link above
figure(2)
probplot(x3,'noref')
h = probplot(gca,pdf_normmixture,paramEsts);
h.Color = 'r';
h.LineStyle = '-';
%end of second example
Is there a way to display something like this on a probplot?
  3 Comments
dpb
dpb on 7 Oct 2022
Edited: dpb on 7 Oct 2022
Agreed. This penchant of having multiple disparate objects that are similar but not the same is a pain. It seems the Statistics Toolbox has been particularly bad about introducing stuff before it's fully ripe (the ill-fated datatset comes to mind, too.)
Looks like you'll have to build the probablllity plot for the custom distribution by evaluating it directly without the benefit of the prepared probplot.
I've not looked inside it so see just how much one might be able to cop as a starting point or if it is exceedingly complex making use of the zillion undocumented behind-the-sceenes pieces.

Sign in to comment.

Accepted Answer

dpb
dpb on 7 Oct 2022
Edited: dpb on 7 Oct 2022
Actually, when I had a few minutes to think about it, it's not too hard; you just have to define the associated CDF function to go w with your PDF and use it instead...
%example from fist link above
rng(10) % For reproducibility
x3 = [trnd(20,1,50) trnd(4,1,100)+3];
%histogram(x3)
pdf_normmixture = @(x3,p,mu1,mu2,sigma1,sigma2) ...
p*normpdf(x3,mu1,sigma1) + (1-p)*normpdf(x3,mu2,sigma2);
pStart = .5;
muStart = quantile(x3,[.25 .75]);
sigmaStart = sqrt(var(x3) - .25*diff(muStart).^2);
start = [pStart muStart sigmaStart sigmaStart];
lb = [0 -Inf -Inf 0 0];
ub = [1 Inf Inf Inf Inf];
options = statset('MaxIter',300,'MaxFunEvals',600);
paramEsts = mle(x3,'pdf',pdf_normmixture,'Start',start, ...
'LowerBound',lb,'UpperBound',ub,'Options',options);
xgrid = linspace(1.1*min(x3),1.1*max(x3),200);
pdfgrid = pdf_normmixture(xgrid, ...
paramEsts(1),paramEsts(2),paramEsts(3),paramEsts(4),paramEsts(5));
probplot(x3,'noref')
% define a cdf function to match...
cdf_normmixture = @(x3,p,mu1,mu2,sigma1,sigma2) ...
p*normcdf(x3,mu1,sigma1) + (1-p)*normcdf(x3,mu2,sigma2);
h = probplot(gca,cdf_normmixture,paramEsts);
h.Color = 'r';
h.LineStyle = '-';
title('Probability plot for Mixed Normal Distribution')
  2 Comments
dpb
dpb on 7 Oct 2022
Glad to help -- I knew that was what needed to be able to add; just wasn't thinking clearly how simple the CDF function was to write until came to me in a realization was over-complicating things...

Sign in to comment.

More Answers (0)

Products


Release

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by