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

7 views (last 30 days)
Andrew Feenan on 5 Oct 2022
Commented: dpb on 7 Oct 2022
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.
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
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?
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.

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...
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 CommentsShowHide 1 older comment
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...

### Categories

Find more on Half-Normal Distribution in Help Center and File Exchange

R2022b

### Community Treasure Hunt

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

Start Hunting!

Translated by