Least squares fit to multiple differently-sized data sets simultaneously

I can use lsqcurvefit to simultaneously fit a system of non-linear equations using matrix arguments. However, I think the matrix requirement to have every element defined means that all data sets must be the same size? I can't find a way around that.

To make it clearer, if I make up some simple data sets which are the same size:

x=[1 3 5 7;     5 7 9 11;    9 11 13 15];
y=[10 11 12 13; 10 11 12 13; 10 11 12 13];

and want to fit each with a straight line, where they all share the same gradient (and so the solver must simulatenously fit that across all three data sets, but fit the individual offset to each set separately):

%variables v(1:4) are m, c1, c2, c3
fn = @(v,xdata)[v(1).*xdata(1,:) + v(2); v(1).*xdata(2,:) + v(3); v(1).*xdata(3,:) + v(4)];

I can set up a guess and solve using:

x0 = [1; 10; 9; 8];
fitvars = lsqcurvefit(fn,x0,x,y);

The problem comes when the data sets don't have the same number of points - if the second row of xdata and ydata only had three points, how do I pass the data to lsqcurvefit?

I've tried padding the matrices with NaN but then the objective function returns undefined numbers and the solver complains. I tried converting the data and the function to cell arrays to allow different lengths, but lsqcurvefit won't take it.

I'm aware this particular example can be solved trivially, but the real functions and data sets are non-linear (which explains why I'm using lsqcurvefit) and the system has a minimum of five data sets to solve with one shared variable between each function, as well as several unique variables and constants to each.

Thanks for any help.

Answer 1

[John D'Errico]

Yes. There is a way. Simply store the data in a cell array. Make the vectors COLUMN vectors. I was lazy here when I did that, using transposes.

x={[1 3 5 7]' [5 7 9 11]' [9 11 13 15]'};
y={[10 11 12 13]' [10 11 12 13]' [10 11 12 13]'};
yy = vertcat(y{;});

%variables v(1:4) are m, c1, c2, c3
fn = @(v,xdata) [v(1).*xdata{1} + v(2); v(1).*xdata{2} + v(3); v(1).*xdata(3,:) + v(4)];

x0 = [1; 10; 9; 8];
fitvars = lsqcurvefit(fn,x0,x,yy);

[Richard Cobley]

Thanks John - this approach worked.

I tidied up a couple of typos, dropped out one point from the data to check it works, and wrote a (terrible) loop approach to re-shape the output from the function to match the original cell array format. The working solution is below for anyone else who has this problem.

Thanks a lot,

Richard.

x={[1 3 5 7]' [5 7 9 11]' [9 11 15]'};
y={[10 11 12 13]' [10 11 12 13]' [10 11 13]'};
yy = vertcat(y{:});

%variables v(1:4) are m, c1, c2, c3
fn = @(v,xdata) [v(1).*xdata{1} + v(2); v(1).*xdata{2} + v(3); v(1).*xdata{3} + v(4)];

x0 = [1; 10; 9; 8];
fitvars = lsqcurvefit(fn,x0,x,yy);

yfit = testfn(fitvars,x);

%re-shape fitted y points back to original format
yfitcell=cell(size(y));
yyc = yy;
for counter1 = 1 : size(y,2)
   yfitcell{counter1} = yyc(1:length(y{counter1}));
   yyc(1:length(y{counter1}))=[];
end

figure; hold on;
cellfun(@plot,x,y,{'.' '.' '.'})
cellfun(@plot,x,yfitcell,{'--' '--' '--'})

[Vipultomar]

Hi. Thanks for sharing the solution. So overall all the ydata (of different dimensions) have been concatenated in a column and are being compared to theoretical ydata values from the function (which in turn should also be in the same order). Right?

[Riccardo Marinelli]

Hi, thank you for your discussion. What could I do if I wanded to perform a fit weighted on some vectors associated to the ydata?

To give an example, I would like to fit

x1 = [2 4 9]; y1 = [5 13 26]; x2 = [0 2 4 9]; y2 = [1 6 10 28]

but to y1 and y2 I associate uncertainty vectors

err1 = [0.1 0.2 0.5]; err2 =[0.2 0.4 0.1 0.6];

and as a weight I'd like to take

weight1 = 1./(err1.^1); weight2 = 1./(err2.^2);

[John D'Errico]

@Riccardo Marinelli

There are so many unspoken things in your question that I cannot easily answer you. Is this a LINEAR regression? Or a nonlinear regression, as the original question asked?

Do you really mean to compute the weights as you did?

weight1 = 1./(err1.^1); weight2 = 1./(err2.^2); 

Note that weight1 only uses err1.^1, yet you are squaring the weights for weight2. Was that intentional, or a typo?

The original question was a complicated one, where some of the parameters are shared between models, but you did not seem to indicate that.

In fact, it looks like all you are asking to do is a weighted linear regression, perhaps of the form y=a*x+b, and to do it several times. As such, you could just use a loop, calling the function lscov each time, for each subproblem. (lscov is perhaps the simplest tool to perform a weighted regression, that I know you will have because it is not part of any toolbox.) Could you do that for all of the unknowns at once, in one big call? Well, yes, but that is surely more effort than you need, and unless you understand how to use sparse matrices, it would be inefficient.

[Riccardo Marinelli]

Hi @John D'Errico, yes, err1.^1 was a typo.

My situation was this: I needed to perform a simultaneous fit of seven curves with shared parameters since I didn't possess enough data to do many fits. I also desired to weight the fit as it is usually done on the squared inverse of the uncertainties

I found the aweseome custom function nlinmultifit() that allows me also to obtain easily the confidence interval of the fits, so yesterday I stuck with this function.

If you are interested, here is what I obtained with some made up data to make a test:

Answer 2

[Qu Cao]

The data sets should have the same number of points. The 'xdata' and 'ydata' should be well defined vetors or matrices.

[Richard Cobley]

The data sets I need to fit do not have the same number of points - that is the problem I am trying to solve.