How to fit a curve by a series of segmented lines in Matlab?

Question

I have a simple loglog curve as above. Is there some function in Matlab which can fit this curve by segmented lines and show the starting and end points of these line segments ? I have checked the curve fitting toolbox in matlab. They seems to do curve fitting by either one line or some functions. I do not want to curve fitting by one line only.

If there is no direct function, any alternative to achieve the same goal is fine with me. My goal is to fit the curve by segmented lines and get locations of the end points of these segments .

Answer 1

First of all, your problem is not called curve fitting. Curve fitting is when you have data, and you find the best function that describes it, in some sense. You, on the other hand, want to create a piecewise linear approximation of your function.

I suggest the following strategy:

1. Split manually into sections. The section size should depend on the derivative, large derivative -> small section
2. Sample the function at the nodes between the sections
3. Find a linear interpolation that passes through the points mentioned above.

Here is an example of a code that does that. You can see that the red line (interpolation) is very close to the original function, despite the small amount of sections. This happens due to the adaptive section size.

function fitLogLog()
   x = 2:1000;
   y = log(log(x));

   %# Find section sizes, by using an inverse of the approximation of the derivative
   numOfSections = 20;
   indexes = round(linspace(1,numel(y),numOfSections));
   derivativeApprox = diff(y(indexes));
   inverseDerivative = 1./derivativeApprox;
   weightOfSection =  inverseDerivative/sum(inverseDerivative);   
   totalRange = max(x(:))-min(x(:));   
   sectionSize = weightOfSection.* totalRange;

   %# The relevant nodes
   xNodes = x(1) + [ 0 cumsum(sectionSize)];
   yNodes = log(log(xNodes));

   figure;plot(x,y);   
   hold on;
   plot (xNodes,yNodes,'r');
   scatter (xNodes,yNodes,'r');
   legend('log(log(x))','adaptive linear interpolation');
end

Answer 2

Andrey's adaptive solution provides a more accurate overall fit. If what you want is segments of a fixed length, however, then here is something that should work, using a method that also returns a complete set of all the fitted values. Could be vectorized if speed is needed.

Nsamp = 1000;     %number of data samples on x-axis
x = [1:Nsamp];    %this is your x-axis
Nlines = 5;       %number of lines to fit

fx = exp(-10*x/Nsamp);  %generate something like your current data, f(x)
gx = NaN(size(fx));     %this will hold your fitted lines, g(x)

joins = round(linspace(1, Nsamp, Nlines+1));  %define equally spaced breaks along the x-axis

dx = diff(x(joins));   %x-change
df = diff(fx(joins));  %f(x)-change

m = df./dx;   %gradient for each section

for i = 1:Nlines
   x1 = joins(i);   %start point
   x2 = joins(i+1); %end point
   gx(x1:x2) = fx(x1) + m(i)*(0:dx(i));   %compute line segment
end

subplot(2,1,1)
h(1,:) = plot(x, fx, 'b', x, gx, 'k', joins, gx(joins), 'ro');
title('Normal Plot')

subplot(2,1,2)
h(2,:) = loglog(x, fx, 'b', x, gx, 'k', joins, gx(joins), 'ro');
title('Log Log Plot')

for ip = 1:2
    subplot(2,1,ip)
    set(h(ip,:), 'LineWidth', 2)
    legend('Data', 'Piecewise Linear', 'Location', 'NorthEastOutside')
    legend boxoff
end

Answer 3

This is not an exact answer to this question, but since I arrived here based on a search, I'd like to answer the related question of how to create (not fit) a piecewise linear function that is intended to represent the mean (or median, or some other other function) of interval data in a scatter plot.

First, a related but more sophisticated alternative using regression, which apparently has some MATLAB code listed on the wikipedia page, is Multivariate adaptive regression splines.

The solution here is to just calculate the mean on overlapping intervals to get points

function [x, y] = intervalAggregate(Xdata, Ydata, aggFun, intStep, intOverlap)
% intOverlap in [0, 1); 0 for no overlap of intervals, etc.
% intStep    this is the size of the interval being aggregated.

minX = min(Xdata);
maxX = max(Xdata);

minY = min(Ydata);
maxY = max(Ydata);

intInc = intOverlap*intStep; %How far we advance each iteraction.
if intOverlap <= 0
   intInc = intStep;
end
nInt = ceil((maxX-minX)/intInc); %Number of aggregations

parfor i = 1:nInt
    xStart = minX + (i-1)*intInc;
    xEnd   = xStart + intStep;
    intervalIndices = find((Xdata >= xStart) & (Xdata <= xEnd));
    x(i) = aggFun(Xdata(intervalIndices));
    y(i) = aggFun(Ydata(intervalIndices));
end

For instance, to calculate the mean over some paired X and Y data I had handy with intervals of length 0.1 having roughly 1/3 overlap with each other (see scatter image):

[x,y] = intervalAggregate(Xdat, Ydat, @mean, 0.1, 0.333)

x =

Columns 1 through 8

0.0552    0.0868    0.1170    0.1475    0.1844    0.2173    0.2498    0.2834

Columns 9 through 15

0.3182    0.3561    0.3875    0.4178    0.4494    0.4671    0.4822

y =

Columns 1 through 8

0.9992    0.9983    0.9971    0.9955    0.9927    0.9905    0.9876    0.9846

Columns 9 through 15

0.9803    0.9750    0.9707    0.9653    0.9598    0.9560    0.9537

We see that as x increases, y tends to decrease slightly. From there, it is easy enough to draw line segments and/or perform some other kind of smoothing.

(Note that I did not attempt to vectorize this solution; a much faster version could be assumed if Xdata is sorted.)