CODEWITHSUNDEEP
matlabdataset
Rosalie
Rosalie

Reputation: 15

Find the inflection point of a discrete dataset in Matlab

I have a dataset (an array with x and y values) to which I eventually want to fit a nonlinear model. In order to increase the chances that the fit will be good and that the model will converge, I want to supply initial parameters as close to the solution as possible. One of the parameter of the model I want to fit corresponds to the inflection point of the curve. Is there a way I can estimate the inflection point of a discrete dataset (not a continuous function) in Matlab? I don't want to have to curvefit beforehand since this is intended to accelerate the curvefit and increase its probability of being successful.

Upvotes: 0

Views: 3404

Answers (1)

user2271770
user2271770

Reputation:

The idea is to work in 2 steps:

  1. First define what is the noise level, and "simplify" the (x, y) polygonal line to some "smooth" subset (xs, ys). After the simplification, every change in curve features will be considered significant.

  2. Look for the changes in convexity of the polygonal line (xs, ys).

For first point you may use the Douglas-Peucker algorithm, implemented here.

The second point is based on the signed area constructed from 3 consecutive points of the polyline:

      | x0   y0   1 |
    1 |             |
A = - | x1   y1   1 |
    2 |             |
      | x2   y2   1 |

The polyline is "convex" - i.e. positive curvature - where these areas are positive, and "concave" - i.e. negative curvature - when the area is negative. The changes in curvature are the places where the inflection points are located.

After downloading the dpsimplify function from File Exchange, you could run the following code (suppose that x and y vectors are already existing):

%// Part 1.
[ps,~] = dpsimplify([x(:),y(:)], 1e-3);  %// adjust "tol" above noise level
xs = ps(:,1);
ys = ps(:,2);

%// Part 2
I = 1:numel(xs)-2;
sgnA = sign( ...
    xs(I+0).*ys(I+1) ...
  + xs(I+1).*ys(I+2) ...
  + xs(I+2).*ys(I+0) ...
  - xs(I+0).*ys(I+2) ...
  - xs(I+1).*ys(I+0) ...
  - xs(I+2).*ys(I+1) ...
);
k_inflex = find(2 == abs(diff(sgnA)));
x_inflex = xs(k_inflex);
y_inflex = ys(k_inflex);

Upvotes: 2

Related Questions