Limit divergence of polynomial fitting

I am approximating a set of data with a polynomial to later find the cumulative density function. I am using a polynomial since I do not know the distribution of the set of data, that may even vary, and because it should be fully automatized. The problem is thatat the end of my set of data, once the polynomial reaches the last data, divergest to infinity. Since it should approximate a probability density function it shouldn't happen, because like so it influence all the results. How can I fix that?

      p=polyfit(X,N,21);
      x=linspace(0,4000,1000);
      y=polyval(p,x);

Answer 1

[John D'Errico]

The simple answer for this is to ask a question. WHY IN THE NAME OF GOD AND LITTLE GREEN APPLES ARE YOU USING A POLYNOMIAL????? Ok, maybe I'm being too hard on you. Polyfit is there for you to use. Nothing in the documentation tells you it does silly things, or that this is the fundamental nature of polynomials. Somewhere teachers should teach you, right after they tell you that you CAN use a polynomial, that you probably, usually, should not do so. And then explain why not.

Seriously, DON'T USE A POLYNOMIAL. And NEVER use a high order polynomial. They go to hell faster than a low order one. A polynomial order of 21? WAY TOO HIGH. But then you should not use a polynomial at all here.

Polynomials do exactly what you saw. They go to infinity. Think about it. what other behavior can you possibly expect? Can a polynomial approximate a function that does non-trivial (non-constant) stuff within some region, then goes constant as x goes to plus or minus infinity? Think about it. A polynomial can never be a model for a process that behaves like a CDF. You have this really high order polynomial, with powers as high as 21. X is as large as 4000?

What is 4000^21? How about other numbers near there?

4000^21
ans =
    4.398e+75

3950^21
ans =
   3.3771e+75

3500^21
ans =
   2.6634e+74

The highest power of X with a non-zero coefficient will dominate the behavior of a polynomial as X goes to +/-inf. So your high order polynomial produces complete and utter garbage out there. You know this is the behavior that you see. And simple mathematics tells you that you cannot expect any other behavior from a polyomial model. That is exactly the behavior you expect from a polynomial, the innate, fundamental behavior that exists in any polynomial, and something you cannot change.

The obvious answer to my original question is you have no idea what else you might use. Sorry, but that is clearly true.

And, the answer to that is to use a spline. There are several choices you might use there. The simplest is pchip, which will interpolate your data. (Of course, you could have just tried interp1, which can use pchip, or various options as you would direct it to use.)

pp = pchip(X,N);
x = linspace(0,4000,1000);
y = ppval(pp,x);

DON'T use the function spline, as it will show oscillatory behavior in the tails. If you use interp1 instead, then the 'linear', 'pchip', or 'cubic' methods should all be ok. While 'linear' will just do piecewise linear interpolation, that is probably entirely adequate for this purpose.

IF you feel that you needed to smooth the curve instead of interpolate it, then download my SLM toolbox from the file exchange. It will require the optimization toolbox. But pchip should be entirely adequate here.

[NF]

Yes, I don't have any idea to what else I could use. However the numbers were just put as test, nothing specific, I only took the last run I did. I'm going now to try the pchip function, hoping this one works as a probability density function that can be than integradet for cdf. Thanks!

[NF]

I have used this spline, however as I thought, being this an interpolation and starting from a bar chart to find a probability density function, the curveis not what I was looking for. In this picture there is what matlab gives me.

I should obtain a curve that has a "classical" distribution, like a Weibull or Lognormal distribution.

Is there some tool that may help to do so? As I said I don't know which distribution the data have and even if I discover it it may change in other situations. Thanks!!

[John D'Errico]

First of all, don't use pchip on an empirically posed hstogram for the PDF.

Better is to use it on the CDF. Then differentiate it to get an approximation for the PDF.

Let me give you an example. I don't have your data, and I'm too lazy to make some Weibull data.

x = exp(randn(1,100));
hist(x,100)

Yes. You and I know this is lognormally distributed data. SHH. Don't tell anyone else. It is a secret.

x = sort(x);
p = ones(size(x));p = p/numel(p);

So an empirical cdf would simply be:

n = numel(x);
cp = linspace(0,1,n);
plot(x,cp,'-')
grid on

A little bumpy. But if all you are willing to get is an empirical CDF, then how much can you ask for?

We can fit it with a spline. But again, you cannot possibly think a polynomial would have a chance in any universe that it would fit that.

Again, polynomials simply do NOT have the property that they are flat as x goes to infinity.

We can build a pchip spline fit. It can be evaluated, as you see.

pp = pchip(x,cp);
ppval(pp,.5)
ans =
      0.27979

xinterp = linspace(0,max(x),1000);
plot(x,cp,'bo',xinterp,ppval(pp,xinterp),'r-')
title 'Interpolated CDF'

A little bumpy, as you might expect an interpolated CDF to be. Yes, I could differentiate that if I really, desperately wanted an approximation to the PDF, but it would be not terribly useful.

We could use a regression spline, coupled with a set of intelligent constraints on the shape. My SLM toolbox does exactly that. It can be found on the file exchange.

slm = slmengine(x,cp,'knots',[0,exp(linspace(log(.1),log(max(x)),6))],'leftvalue',0,'rightvalue',1,'rightslope',0,'increasing','on','concavedown',[1,max(x)],'plot','on');

Differentiate that function to get an empirical PDF. Because the CDF as generated has the property that it goes to 1.0 exactly at the top end, we know the integral of the PDF will be exactly 1.

The derivative curve an an empirical PDF. In blue, I overlaid the known lognormal PDF. Given a relatively small sample of only 100 data points, this seems eminently reasonable.

That works reasonably well, if you have no idea what the distribution would be.

Of course, I could have used tools from the stats toolbox too, to fit a distribution to the same data. For example:

PD = fitdist(x','lognorm')
PD = 
  LognormalDistribution

    Lognormal distribution
         mu = -0.157192   [-0.36937, 0.054986]
      sigma =   1.06933   [0.938878, 1.24221]

As I generated the data, the correct distribution would have had mu=0 & sigma=1, which is entirely consistent with the fitted distribution parameters.

So the red curve is the empirical spline PDF. Blue is the estimate from fitdist, and green is the actual distribution that I sampled the data from. 100 data points is not a lot to estimate a distribution though.

[John D'Errico]

Finally, what can you do what you have absolutely no clue what your distribution is?

FITDIST is pretty good here.

PD = fitdist(x','kernel','support','positive')
PD = 
KernelDistribution

    Kernel = normal
    Bandwidth = 0.514084
    Support = positive

Here, I told it nothing more than that the distribution has support on the positive half of the real line, thus [0,inf].

plot(xinterp,pdf(PD,xinterp))
hold on
plot(xinterp,lognpdf(xinterp,0,1),'g-')

As you can see, the result (the blue curve) is very consistent with all the other fits I did before.

[NF]

Thanks, the approximation works perfectly for my needs.

I just have one problem with that. The data that I get consists of thousands of values and it happens that some are the equal. In such a situation Matlab returns me the error "Error using chckxy (line 51) The data sites should be distinct."

I found that this error is due to the fact that the same x-value has two or more different y-values and then Matlab does not know how to interpolate. How should I treat these doubles? Thank you

[John D'Errico]

So the pchip spline is the one that gets upset at the reps? It would do so. And chckxy would be part of that code.

There are two simple solutions. One is to use a tool like my consolidator to remove the replicates.

x = rand(10,1);
x = [x;x];
y = exp(x);

spl = pchip(x,y)
Error using chckxy (line 51)
The data sites should be distinct.
Error in pchip (line 58)
[x,y,sizey] = chckxy(x,y);

Which of course fails.

[xu,yu] = consolidator(x,y,@mean);

spl = pchip(xu,yu);
fnplt(spl)

You could have done something similar using uniquetol in combination with accumarray.

Consolidator is available for download from the file exchange.

https://www.mathworks.com/matlabcentral/fileexchange/8354-consolidator

[NF]

Now it works perfectly. Thank you so much for all the help!!