Estimating parameter uncertainty with fminsearch

I have four sets of data that I would like to fit to four lines, but with certain constraints (e.g. the four lines should all have the same slope, etc.). Some of the data sets overlap in x, but not y, i.e. they can't be fit to a single function. Currently I'm fitting all four simultaneously by calculating the sum of squares for each set/line, adding them up, and using fminsearch on the total sum of squares function. It works well, but fminsearch does not return any uncertainty on the fit parameters. Is there a way to either calculate the uncertainty on the parameters returned by fminsearch, or do this fit (with the constraints) without using fminsearch?

I found this post that seems to suggest I can use the Hessian of my sum of squares function: https://stackoverflow.com/questions/38175179/uncertainties-errors-for-non-continuous-minimizer-in-octave However, I'm not very well-versed in statistics and couldn't quite understand whether this is the right approach. Any help is appreciated! Thanks in advance!

[John D'Errico]

Is all you are trying to do is fit a model with a single common slope, but where the intercept term varies for each subset of data?

If so, then using fminsearch is wildly the wrong thing to do for many reasons. But I'm not positive of that, because I wonder if that is really the entire problem. Your description is too vague to be sure here.

If I am correct, then it is also true that computing what you need is also very easy, even arguably for the parameter uncertainty.

Other questions remain however. Even if all that is true, what toolboxes do you have? Two that would be pertinent in this respect would be the stats tool box and the curve fitting toolbox.

[A]

Thanks for your reply, John. I should have been a bit more descriptive in my original question. The four sets of data look roughly like the four branches of an X shape. I want all four to have the same slope (up to a sign). If the top two branches intersect at (xt,yt) and the bottom two branches intersect at (xb,yb), I want xt and xb to be the same (yt and yb should be different). Currently I'm fitting to four parameters: overall x and y offsets, the distance between yt and yb, and the slope of all four lines. What is a good way to fit all four data sets at once and get uncertainty estimates on the four parameters? I have both the Stats and Curve Fitting Toolboxes. I appreciate your help!

Answer 1

[John D'Errico]

I think I am starting to understand what you are askng. One very good option is to use abs. I've been told I could use better abs myself. oops. wrong kind of abs. :)

That is, as long as the two upper segments seem to intersect at a point, as well as the two lower segments, then for example, you could use abs to fit each half. Thus the model essentially becomes

y_up = yt + S*abs(x - xt)

and

y_down = yb - S*abs(x - xt)

I could show you how to combine both pieces into one model. As you can see there, xt is the same for both sub-models, as is S, the slope of those lines. The nice thing is, abs gets you part of the way there. Keeping yb and yt separate alows you to have a separation between the two halves vertically.

The next question then becomes, do you know for each point, which branch of the shape these points fall on? Most important would be a vector that would identify the upper and lower branches. For example, suppose you had a boolean vector V that is either 1 for any point that lives on the upper branch, and 0 if the point lives along the lower branch.

N = 100;
x = rand(N,1)*10;
V = +(rand(N,1) > 0.5); % + makes V a double

% the true values of these parameters so I can generate some data.
% Don't tell anyone what they really are!
yt0 = 4;
yb0 = 2;
xt0 = 5;
S0 = 2.5;

y = (V*2 - 1).*abs(x - xt0)*S0 + yt0*V + yb0*(1-V) + randn(N,1)*.5;

plot(x,y,'o')

Yes, I know ground truth for all those parameters, thus [xt, yb, yt, S]. But unless I have some fake data, I won't be able to fit a model. I'll use the curve fitting toolbox, since you have it, and it will very conveniently predict the parameter uncertainties.

myfitfun = @(xt,yb,yt,S,x,V) (V*2 - 1).*abs(x - xt)*S + yt*V + yb*(1-V);
ft = fittype(myfitfun,'indep',{'x','V'})
ft = 
     General model:
     ft(xt,yb,yt,S,x,V) = (V*2-1).*abs(x-xt)*S+yt*V+yb*(1-V)

All seems good there. I wonder if fit will manage to survive without even giving it starting values?

mdl = fit([x,V],y,ft)
     General model:
     mdl(x,V) = (V*2-1).*abs(x-xt)*S+yt*V+yb*(1-V)
     Coefficients (with 95% confidence bounds):
       xt =        4.99  (4.954, 5.026)
       yb =       2.082  (1.9, 2.263)
       yt =       3.938  (3.752, 4.124)
       S =       2.535  (2.474, 2.596)

Note that as far as fit is concerned, this is a two independent variable model, with x and V being independent variables, and with 4 unknown parameters to estimate. The CFTB actually thinks this is a surface I am fitting, not that I really care if I might hurt its feelings. :)

As it is, the parameter estimates are pretty decent, with the bounds nicely enclosing the ground truth values.

You should recognize that as I built this model, the slopes of all 4 lines are the same, except for a sign swap on two of the segments. As well, the two halves are split top and bottom, but they join at a common value of xt for both halves of the X.

ypred = myfitfun(mdl.xt,mdl.yb,mdl.yt,mdl.S,x,V);
plot(x,y,'ko',x,ypred,'xr')
legend('data','model')

The hardest thing about this was making up the fake data. All it really requires for the fit is to know which branch of the X any point lies on. That is, top or bottom? It might have been smarter if I had provided some decent guesses for the parameters to fit, but I seem to have gotten away with it. That need not always be the case though.

[A]

This worked perfectly! Very neat to use V as an additional independent variable – this should help with some other fitting problems I've been working on as well (and yes, I do know which points belong to which branch). Thank you so much for your help, I really appreciate it!