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

Hi,

I am trying to plot a two peak distribution model on a probplot.

For example, Fit Custom Distributions - MATLAB & Simulink Example - MathWorks United Kingdom

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?

[dpb]

I've not poked at the new objects in the Statistics TB much; I quit the active consulting gig long enough ago they weren't yet available.

But, to use probplot looks like you have to build your custom distrubtion via the DistributionFitter App with a custom function. It seems to have a bunch of needs including building a m-file that describes it an putting it in a custom-created +prob directory on your MATLABPATH. That's pretty deep...I thought there should be some way with makedist to take custom mle object and convert it, but doesn't seem to be; you can supply specific values of coefficients to the know distributions, but there's not the ability to reference a custom one that I can find there.

I gotta' run at the moment, but look at the following to see the way they say to go about it --

<Custom-distributions-using-the-distribution-fitting-tool>

[Andrew Feenan]

Thanks for your response.

The distribution fitter app isn't particularly suitable for my use case but thanks for looking into it.

It' frustrating this isn’t possible on a probplot like it is on a histogram etc.

[dpb]

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.

Answer 1

[dpb]

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')

[Andrew Feenan]

Thank you very much. You have saved me so much time. I was in the middle of building the probability plot for my distribution function.

[dpb]

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...