How to find a inflexion point with only discrete datas?

Hi everyone,

I'm currently working on a fonction that need to find the inflexion point of a curve that I drew with only discrete datas (I have a abs vector with the abscissa and a vector y with ordinate). The typical shape of the curve look like this :

I don't use matlab a lot so maybe there's a very easy way to get it.

Here is an example for this plot :

x = [ 0 0.05 0.1 0.17 0.25 0.37 0.5 0.72 1 1.42 2 2.83 4 5.67 8 11.32 15 21.22 30 42.44 60 84.87 120 169.72 339.42 480 720 960 1200 1440];

y = [ 17.836 17.742 17.716 17.693 17.669 17.642 17.615 17.578 17.536 17.484 17.422 17.347 17.26 17.159 17.048 16.934 16.849 16.765 16.706 16.669 16.644 16.626 16.612 16.6 16.58 16.572 16.563 16.557 16.552 16.547];

The abcissa are first computed into their log10 form before being plotted?

I hope someone can help me but anyway, have a good day !

PS : I don't know if it's useful but I use matlab R2012a.

Answer 1

[Azzi Abdelmalek]

x = [ 0  0.05  0.1  0.17  0.25  0.37  0.5  0.72  1  1.42  2  2.83  4  5.67  8  11.32  15  21.22  30  42.44  60  84.87  120  169.72 339.42  480  720  960  1200  1440];
y = [ 17.836  17.742  17.716  17.693  17.669  17.642  17.615  17.578  17.536  17.484  17.422  17.347  17.26  17.159 17.048  16.934  16.849  16.765  16.706  16.669  16.644  16.626  16.612  16.6  16.58  16.572  16.563  16.557  16.552 16.547];
xx=log10(x)
df1=gradient(y,xx)
df2=gradient(df1,xx)
id=sign(df2)
idx=strfind(id,[-1 1])
inflexionP=x(idx+1)

[Justin Miroir]

Thank you for this answer but I think there is a problem (correct me if I'm wrong) : this code is OK if the inflexion point is exactly at one of the data points but that's not necessarily the case. In fact this give you the interval where you can find it, right?

But I think that's the way to take. I don't know how to do it but I should get a function with the data points and find the exact point where the slope of the tangent is the steepest. Is it possible?

Again thank you for your answer.

[Azzi Abdelmalek]

What exact point means? You will never find it numerically, because all data are stored in memory with a certain error. You can insert other points by interpolation then find your point with a certain approximation.

[Justin Miroir]

Is it not possible to get a curve with (let's say) 10000 who fits the original data points and then find the inflexion point?

[Azzi Abdelmalek]

Like I said, you can interpolate to find your 10000 points, you can use interp1 function

x = [ 0  0.05  0.1  0.17  0.25  0.37  0.5  0.72  1  1.42  2  2.83  4  5.67  8  11.32  15  21.22  30  42.44  60  84.87  120  169.72 339.42  480  720  960  1200  1440];
y = [ 17.836  17.742  17.716  17.693  17.669  17.642  17.615  17.578  17.536  17.484  17.422  17.347  17.26  17.159 17.048  16.934  16.849  16.765  16.706  16.669  16.644  16.626  16.612  16.6  16.58  16.572  16.563  16.557  16.552 16.547];
xi=linspace(min(x),max(x),10000);
yi=interp1(x,y,xi)

[Justin Miroir]

Sorry I just saw your comment. Thank you for your answer but I used the method from John, it seems to work :) Thank you anyway!

Answer 2

[John D'Errico]

You have a nice smooth curve. So you can use an interpolating spline.

x = [ 0  0.05  0.1  0.17  0.25  0.37  0.5  0.72  1  1.42  2  2.83  4  5.67  8  11.32  15  21.22  30  42.44  60 84.87  120  169.72  339.42  480  720  960  1200  1440];

y = [ 17.836  17.742  17.716  17.693  17.669  17.642  17.615  17.578  17.536  17.484  17.422 17.347  17.26  17.159  17.048  16.934  16.849  16.765  16.706  16.669  16.644  16.626  16.612 16.6  16.58  16.572  16.563  16.557  16.552  16.547];

xlog = log10(x);
pp = spline(xlog(2:end),y(2:end));
fnplt(pp)
hold on
plot(xlog,y,'o')

Note that I fit the curve from the second point on, because the log function is singular at 0, approaching -inf as a limit.

The problem is, finding an inflection point can be difficult to pin down exactly, because it is equivalent to estimating the second derivative of a curve, only given data. This is one variety of something called an ill-posed problem. Differentiation of data tends to amplify any noise in the data, and do so greatly. Thus, a good reason for why is called ill-posed.

So, we can now differentiate that cubic spline twice, then look for where the second derivative crosses zero.

ppd2 = fnder(pp,2);
fnplt(ppd2)

Ok, as I said, it gets a bit nasty looking at the ends, because as I said, differentiation amplifies noise. Two derivatives make it worse. We can find all locations where the curve crosses zero easily enough though.

Here, I'll just employ a tool from my SLM toolbox.

d2zerolocs = slmsolve(ppd2,0)
d2zerolocs =
   -1.0711  -0.63129  -0.58148  0.93178  2.6653  2.7129  2.929

slmsolve uses interpolation on the actual function to find all locations on the curve that satisfy the given y value. In this case, I'm using it to find all locations that satisfy the given second derivative value.

We know from the plot that the location of interest is near 1, on the log scale, so log(x) is approximately 0.93. (In fact, I'd not even trust that location to more than about 1 significant digit.) If you want that in the original x domain, just exponentiate using 10^().

10.^d2zerolocs(4)
ans =
       8.5462

I suppose, if you are going to use my slmtoolbox anyway, to do the solve, then you could afford to use it to do the fit too.

slm = slmengine(xlog(2:end),y(2:end),'plot','on','decreasing','on','knots',8,'result','pp');

So maybe a little smoother.

ppd2 = fnder(pp,2);
fnplt(ppd2)
grid on

ppd2 = fnder(slm,2);
fnplt(ppd2)

d2zerolocs = slmsolve(ppd2,0)
d2zerolocs =
    -1.0979      0.87456       2.4792       2.6467

This time, it came out a wee bit below 0.9. As I said before, I'd not trust ANY estimator of the location of an inflection point to more than about ONE significant digit. So be careful. Computers return 16 digits. But here, almost all of those digits are meaningless, virtually random garbage.

10.^d2zerolocs(2)
ans =
       7.4914

The SLM toolbox is found here :

http://www.mathworks.com/matlabcentral/fileexchange/24443-slm-shape-language-modeling

I just noticed that the posted version did not have slmsolve in it yet. So I've posted a new release, that does include slmsolve. You may need to wait a few minutes before it becomes visible to the world.

[Justin Miroir]

Thank you for your answer (and your time), it helped me a lot to explain it step by step. Just a little thing : given that the scale use log10, I find the answer using 10^(ans) and not exp(ans).

I'll try it on my computer and I'll tell you if I have any problem.

[John D'Errico]

DOH! Sorry. You said log10. In fact, I even used log10 there in the example. But then I exponentiated using exp. Of course, you realize, I did that as a test, to see if you were awake and still reading. Yeah. Right. I'll modify the answer to be correct.

[Justin Miroir]

Haha right :) Unfortunatelly, I encountered a problem : I cannot use fnplt or fnder on matlab... Do you have an idea why?

it gives me : Undefined function 'fnplt' for input arguments of type 'struct'.

For the beginning I just copied/pasted your code but anyway it does not work :(

[Justin Miroir]

And you told me the problem was the fact that you get a lot of noise problems when you derivate the function twice. Is it not possible to get with the spline the tangent of each point and then keep the one that has the steepest slope? The inflection point should be at this point.

[John D'Errico]

fnplt and fnder are part of Curve fitting toolbox. Some years ago, they were found in the old splines TB, which then got merged into the curve fitting TB.

You can use plotslm to plot a spline in pp form instead of using fnplt.

fnder is easily done by a simple hack. I've attached to this comment a file that does differentiation of a spline in pp form, called ppder.

Yes, you have a good iea, if you want to find the location automatically. You could probably automatically choose the solution that happens at the point of steepest slope.

pp = slmengine(xlog(2:end),y(2:end),'plot','on','decreasing','on','knots',8,'result','pp');
>> ppd1 = ppder(pp,1);
>> ppd2 = ppder(pp,2);
>> d2zerolocs = slmsolve(ppd2,0)
d2zerolocs =
      -1.0979      0.87456       2.4792       2.6467
>> slopes = ppval(d2zerolocs,ppd1)
slopes =
    -0.085348     -0.76269    -0.056551    -0.057712

This should work with an interpolating spline too, but you need to be careful. If the curve is at all noisy, things will go to hell.

pp = spline(xlog(2:end),y(2:end));
ppd1 = ppder(pp,1);
ppd2 = ppder(pp,2);
d2zerolocs = slmsolve(ppd2,0)
d2zerolocs =
      -1.0711     -0.63129     -0.58148      0.93178       2.6653       2.7129        2.929

slopes = ppval(d2zerolocs,ppd1)
slopes =
    -0.079017     -0.14944     -0.14775     -0.76369    -0.051691    -0.051895    -0.047643

So the simple algorithm of looking for the root that occurs at the point of maximum slope seems to work. BE CAREFUL THOUGH. If your curve has noise in it, then you will be back here, asking why it got the wrong root. :) My answer will be to use SLM for the fit. :)

[Justin Miroir]

Ok thank you for everything. I downloaded the last version of matlab so I can use the other toolboxes :). I'll create a function that find it automatically and that plot the corresponding vertical line to verify if the point is well evaluated.

Still thank you for everything, you helped me a lot. Now I have to implement 3 or 4 more functions :) (I'm working on my thesis).

[John D'Errico]

Good luck on your thesis.

[RAYDEL RAMOS]

The "slmsolve" function is available in matlab 2019a?

I would like to get all the inflexion points as explained in several comments but "slmsolve" is not defined in matlab and then I have not been able to get the desired results. Whenever I run the script I get the following alert message

"Indefinite or variable function 'slmsolve'"

[RAYDEL RAMOS]

Help, I'm a little desperate !!!!