How do I generate a pdf from some known percentile values

Hi Matlab community,

I have percentile values that describe a distribution of possible sea level rise magnitudes. I would like to be able to generate a probability density function that closely approximates the actual distribution from which the percentile values were generated. Can anyone suggest how to achieve this please?

Example data:

prctiles = [5 10 30 33 50 67 70 90 95];

SLR = [3.2760 3.5265 4.1286 4.2013 4.5566 4.9151 4.9836 5.6045 5.9105];

Answer 1

[Star Strider]

A slightly different approach —

prctiles = [5    10    30    33    50    67    70    90    95];
SLR = [3.2760    3.5265    4.1286    4.2013    4.5566    4.9151    4.9836    5.6045    5.9105];

B = fminsearch(@(b) norm(prctiles/100 - cdf('Normal',SLR,b(1),b(2))), [SLR(prctiles==50);rand])

B = 2×1

4.5589 0.8095

SLRv = linspace(min(SLR), max(SLR));
yfit = cdf('Normal', SLRv, B(1), B(2));

figure
plot(SLR, prctiles/100, 'p')
hold on
plot(SLRv, yfit, '-r')
hold off
grid
title(sprintf('$p = N(%.2f, %.3f)$',B), 'Interpreter','latex')
legend('SLR','Fitted Noprmal Distribution', 'Location','NW')

Experiment to get different results.

.

[Image Analyst]

Yes, that could be. I assumed the SLR was the sea level rise and he had thousands of measurements for how much it rose or fell from one time point to another, and he gave just a few samples from that and we needed to find the distribution. But you assume that the SLR points are the points from the CDF and so you find the equation that would give that CDF. I think that maybe you are closer to what @Christopher Stokes meant.

[Star Strider]

It appears that ‘SLR’ are the values that correspond to the percentiles given by ‘prctiles’, making the parameter estimates straightforward with any nonlinear parameter estimation technique (although nlmfit would likely be preferable to fminsearch since it also produces statistics on the parameter estimates and the fit).

The PDF uses the same parameter arguments as the CDF, so that calculation is then straightforward.

I have no idea what the actual ‘SLR’ data represent. That was not stated, so I make no assumptions about them, other than that they apparently correspond to the ‘prctiles’ values.

.

[Christopher Stokes]

This is a great solution, thank you. You are right that the SLR values are estimates of future sea level rise rates, and each value corresponds to a percentile value that describes it's probability of occuring based on many climate simulations.

Your solution works really well for the example I gave previously, which is approximately a normal distribution. The further into the future the SLR projections go, the more skewed the distribution becomes, and so the solution starts to fall down. Is there any way the solution could be generalised to include other probability distribution families (a bit like allfitdist from file exchange)? I know this is a much bigger ask!

Here's some example data for SLR in the year 2100:

SLR = [7.4280 8.2520 10.1210 10.3470 11.6500 13.0120 13.2370 16.3701 18.5500];

prctiles = [5 10 30 33 50 67 70 90 95];

[Star Strider]

As always, my pleasure!

The only way I can think of is to test them against all appropriate distributions, and in my trials, the Generalised Extreme Value distirbution provided the best fit, and fitted the tails appropriately. (I added the fitnlm calls to get statistics on the parameters and the fit. Using only fitnlm would be appropriate. The fminsearch function can be more robust in estimating the parameters, since it doesn’t use graident-descent, so I usually start with it for well-posed optimisations such as this.)

SLR = [7.4280    8.2520   10.1210   10.3470   11.6500   13.0120   13.2370   16.3701   18.5500];
prctiles = [5    10    30    33    50    67    70    90    95];

B = fminsearch(@(b) norm(prctiles/100 - cdf('Generalized Extreme Value',SLR,b(1),b(2),b(3))), [0; rand; SLR(prctiles==50)])

B = 3×1

-0.0557 2.7225 10.6097

GEVmdl = fitnlm(SLR(:), prctiles(:)/100, @(b,x)cdf('Generalized Extreme Value',x,b(1),b(2),b(3)), B)

GEVmdl = 
Nonlinear regression model:
    y ~ F(b,x)

Estimated Coefficients:
          Estimate        SE        tStat       pValue  
          _________    ________    _______    __________

    b1    -0.055765    0.023704    -2.3525      0.056861
    b2       2.7225    0.043684     62.323    1.1473e-09
    b3        10.61    0.027091     391.64    1.8705e-14


Number of observations: 9, Error degrees of freedom: 6
Root Mean Squared Error: 0.00689
R-Squared: 1,  Adjusted R-Squared 1
F-statistic vs. zero model: 2.19e+04, p-value = 1.68e-12

SLRv = linspace(min(SLR), max(SLR));
yfit = cdf('Generalized Extreme Value', SLRv, B(1), B(2), B(3));

figure
plot(SLR, prctiles/100, 'p')
hold on
plot(SLRv, yfit, '-r')
hold off
grid
title(sprintf('$p = GEV(%.4f, %.3f, %.3f)$',B), 'Interpreter','latex')
legend('SLR','Fitted Generalized Extreme Value Distribution', 'Location','NW')

Testing this with the earlier data for comparison —

prctiles = [5    10    30    33    50    67    70    90    95];
SLR = [3.2760    3.5265    4.1286    4.2013    4.5566    4.9151    4.9836    5.6045    5.9105];

B = fminsearch(@(b) norm(prctiles/100 - cdf('Generalized Extreme Value',SLR,b(1),b(2),b(3))), [0; rand; SLR(prctiles==50)])

B = 3×1

-0.2704 0.7888 4.2785

GEVmdl = fitnlm(SLR, prctiles/100, @(b,x)cdf('Generalized Extreme Value',x,b(1),b(2),b(3)), B)

GEVmdl = 
Nonlinear regression model:
    y ~ F(b,x)

Estimated Coefficients:
          Estimate       SE         tStat       pValue  
          ________    _________    _______    __________

    b1    -0.27049    0.0075964    -35.607     3.271e-08
    b2     0.78886    0.0043945     179.51    2.0164e-12
    b3      4.2785    0.0026176     1634.5    3.5397e-18


Number of observations: 9, Error degrees of freedom: 6
Root Mean Squared Error: 0.00235
R-Squared: 1,  Adjusted R-Squared 1
F-statistic vs. zero model: 1.89e+05, p-value = 2.61e-15

SLRv = linspace(min(SLR), max(SLR));
yfit = cdf('Generalized Extreme Value', SLRv, B(1), B(2), B(3));

figure
plot(SLR, prctiles/100, 'p')
hold on
plot(SLRv, yfit, '-r')
hold off
grid
title(sprintf('$p = GEV(%.4f, %.3f, %.3f)$',B), 'Interpreter','latex')
legend('SLR','Fitted Generalized Extreme Value Distribution', 'Location','NW')

The Generalsed Extreme Value distribution would appear to work well for all the data, if the provided examples are representative of them.

.

[Christopher Stokes]

That's a great solution, thank you! Really appreciate your help and great to see how these optimisation functions can be applied.

[Star Strider]

As always, my pleasure!

It was an education for me as well!

.

Answer 2

[Image Analyst]

Have you seen fitdist() in the Stats toolbox?

Of course you'd be better off with much more data.

[Image Analyst]

Here's an example:

clc;    % Clear the command window.
fprintf('Beginning to run %s.m ...\n', mfilename);
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 17;
SLR = [3.2760    3.5265    4.1286    4.2013    4.5566    4.9151    4.9836    5.6045    5.9105];
% Plot data.
subplot(2, 1, 1);
bar(SLR)
grid on;
title('Original SLR Data', 'FontSize', fontSize);
xlabel('Index', 'FontSize', fontSize);
ylabel('SLR Value', 'FontSize', fontSize);
% Get distribution.
d = fitdist(SLR(:), 'Normal')
% Make curve, plot distribution.
% https://en.wikipedia.org/wiki/Normal_distribution
x = linspace(min(SLR), max(SLR), 1000);
amp = 1 / (d.sigma * sqrt(2*pi));
y = amp * exp(-(1/2) * ((x - d.mu) / d.sigma) .^ 2)
subplot(2, 1, 2);
plot(x, y, 'b-', 'LineWidth', 2);
grid on;
title('Estimated Distribution of SLR', 'FontSize', fontSize);
xlabel('SLR', 'FontSize', fontSize);
ylabel('PDF', 'FontSize', fontSize);

Pick the distribution that fits the theory of what the distribution should actually be. Hopefully you know this in advance. Actually you need to if you're going to model it. Otherwise just normalize your histogram and that is the actual PDF.

Answer 3

[Christopher Stokes]

Hi, I have been using that function quite a bit recently and the excellent wrapper function allfitdist() from the file exchange. My problem is that they rely on relatively large data samples from which to build a PDF when what I have is a small number of summary statistics that describe the distribution (i.e. not data samples). I imagine there is a way to use fitdist to do what I need to do, but I can't envisage how this would work.