how to find the inflection point from data contained in an array?

Dear all,

I have an array of size (360*1 cell) where in each cell I have a set of variables an data as the one that I attach here (data.xlsx). If I plot the variables (Volume, Price) I obtain the following scatter plot:

As you can see, the set of points shown an 'inverse S-Shape'.

My question is, how can I find the two 'inflection points' that I marked in red? I want to find the two points using just numerical data instead of using some analytical equation, because so far all my analytical fits have gone wrong (so far, I tried poly 3 and smoothing spline).

Again @Star strider I tried to use your code from question https://fr.mathworks.com/matlabcentral/answers/295156-how-to-find-the-inflection-point-of-a-curve-in-matlab

but is not working with my data. Any help would be really welcome :)

Also, how can I automatize this process so it works for all my cells contained in my array?

If this is really complicated and its better to use an analytical equation. Then, how do I fit 360 curves without doing one by one? and how do I find the inflection points?

Thanks in advance!

[darova]

I found this on the internet

The picture shows something else instead of yours

[Angelavtc]

Yes, I know @darova, actually you helped me to find those point on my graphs https://fr.mathworks.com/matlabcentral/answers/515523-second-derivative-from-a-smoothing-spline-fit , the problem is that when I use the curve fit I find several points where my graph change from Concave to convex, something like this graph:

and as it was suggested by @JohnD'Errico this could be due to the bad fit that I used. Therefore I want to avoid a fit curve and just calculate those points using my numerical data, as I am sure, my excel file follows an 'inverse S-Shape' as it represents the supply of electricity during one hour (also called the Merit Order Curve). I am just interested in the 2 red points that I marked before. Thank you!

[darova]

Looks you are looking for points of maximum curvature (small radius), not inflection points

[Angelavtc]

Right, well if I come back to an spline approximation I tried to use the SLM file:

x = Volume; 
y = Price;
pp = spline(x,y);% 
n=300;
xx =linspace(x(1), x(end),n);
yy= ppval(pp,xx); 

x_max = slmpar(pp,'maxfpp')

I got the following error because there are duplicated x values, how can I solve that?

Error using chckxy (line 50)
The first input must contain unique values.

[darova]

What about this?

A = xlsread('data.xlsx');
x = A(:,1);
y = A(:,2);
%%

ii = [1:4 5:10:length(x)];
x1 = x(ii);
y1 = y(ii);
t1 = [0; cumsum(hypot(diff(x1),diff(y1)))];
t2 = linspace(0,t1(end),100);
x2 = spline(t1,x1,t2);
y2 = spline(t1,y1,t2);

dy = diff(y2)./diff(x2);                       % first derivative
d2y = diff(y2,2)./diff(x2(2:end)).^2;          % second derivative
kk = d2y./(1+dy(2:end).^2).^(3/2);             % curvature (1/r)

plot(x,y,'.b')
hold on
plot(x2,y2,'.-r')
plot(x2(2:end-1),kk*1e6,'.-g')
hold off
legend('original data','spline','curvature','location','northwest')

[Angelavtc]

@darova thanks! this is really good in graph terms! in addition how can I just create a table that contains the pairs (x2,y2) that correspond to the max curvature? I mean the points that correspond to the red circles?

And how can I automatize this procedure given the fact that I have to do the same for 360 curves that arise from 360 timetables contained into a cell array?

Thanks again!

[darova]

What about max and min? Are those peaks always so high?

[Angelavtc]

yep :), @darova I dont know why I cant accept your answer, as you reply really soon to my bad question

Answer 1

[darova]

Please accept my answer

A = xlsread('data.xlsx');
x = A(:,1);
y = A(:,2);
%%

ii = [1:4 5:10:length(x)];
x1 = x(ii);
y1 = y(ii);
t1 = [0; cumsum(hypot(diff(x1),diff(y1)))];    % arc-length (parameter for interpolation)
t2 = linspace(0,t1(end),100);                  % new arc-length
x2 = spline(t1,x1,t2);
y2 = spline(t1,y1,t2);

dy = diff(y2)./diff(x2);                       % first derivative
d2y = diff(y2,2)./diff(x2(2:end)).^2;          % second derivative
kk = d2y./(1+dy(2:end).^2).^(3/2);             % curvature (1/r)

plot(x,y,'.b')
hold on
plot(x2,y2,'.-r')
plot(x2(2:end-1),kk*1e6,'.-g')
hold off
legend('original data','spline','curvature','location','northwest')

Answer 2

[John D'Errico]

YOU ESSENTIALLY JUST ASKED THIS QUESTION. How many times will you ask it?

However, there are two important points I need to explain in what you are doing.

First: replicates. If you have replicated x values, then no, you cannot use spline (OR pchip) directly. Those are interpolating tools. You have what, 16 data points at a volume of 15670? What prediction do you expect an INTERPOLATING spline to produce at that Volume?

plot(Vol,Price,'o')

The simplest way to deal with replicates is to average the values of y, for each replicated x. You can download my consolidator tool from the File exchange. It is designed to do exactly that.

[Vave,Pave] = consolidator(Vol,Price);
plot(Vave,Pave,'o')

xlabel 'Average reduced Volume'
ylabel 'Average reduced Price'
grid on

Still, are you serious? I can count, what roughly at least 5 different inflection points. I see various points where the relationship is possibly discontinuous. Exactly what do you realistically expect to find here? Do you seriously expect to be able to identify second derivative changes, to see where that second derivative changes sign, and to have even the slightest amount of confidence in the result?

What happens when you try to fit this data with a spline?

spl = spline(Vave,Pave);
fnplt(spl)
hold on
plot(Vave,Pave,'ro')
xlabel 'Average reduced Volume'
ylabel 'Average reduced Price'
grid on

There is an old saying about computers, and the use thereof. Garbage in, garbage out. As you can see, the spline has produced just what I would have expected - garbage as a result.

You can do better in terms of a plot, with the use of pchip.

spl = pchip(Vave,Pave);
fnplt(spl)
hold on
plot(Vave,Pave,'ro')
xlabel 'Average reduced Volume'
ylabel 'Average reduced Price'
grid on

While this seems to be more realistic, using a pchip interpolant to find second derivative changes on crappy data is literally insane, because pchip is less able to estimate those second derivatives. It is not as highly differentiable as would be a spline.

Here is the second derivative of the pchip interpolant.

spl = pchip(Vave,Pave);
spl2 = fnder(spl,2);
fnplt(spl2)

Sorry, but the result is exactly what I would expect on this data. Useless.

Hoping to use any method to accurately find an inflection point on that data is almost a laughable idea. I'm sorry, but it is. At the very least, there would be multiple inflection points. I'm sorry, but you are kidding yourself in this task.

[Angelavtc]

yes, @JohnD´errico you are right, I have been kidding myself with all this problem. In any case, thank you!