bilinear curvefit

I'm trying to perform a bilinear curve fit. I tried to create a custom eqn in the curve fitting toolbox (x<a)*(p1*x+p2)+(x>=a)*(p3) but it won't iterate on 'a' as I was hoping. I've also been playing with >> f1 = fittype('(x<a)*(p1*x+p2)+(x>=a)*(p3)', 'independent', 'x', 'coefficients', {'p1', 'p2','p3','a'}) >> [fitobj,gof1] = fit(WavgMay, kskcMay, f1)

p3 could be replaced with p1*a+p2 but I wasn't sure how matlab was handling it.

Any suggestions for the best way to do this?

one more thing, the fit is giving really low values of 'a' which i assume is because it's start with a random number between 0 and 1. How do I incorporate a better starting point?

Ultimately I was hoping to adapt code from http://www.nbb.cornell.edu/neurobio/land/PROJECTS/MKG23curvefit/index.html to test the sensitivity of the initial guess to help narrow down the threshold value of 'a'.

Thanks for your help. ds

Answer 1

[John D'Errico]

This is essentially just a linear spline, so a pair of piecewise linear segments, continuous at the join point.

The problem is, locating the break point is somewhat more difficult. First, this is an old post, and no data was even provided. So I need to make up some data to show how it works.

x = rand(100,1)*10;
y = 2*x.*(x < 4) + (12 - x).*(x >= 4) + randn(size(x))/5;
plot(x,y,'.')

Now, if the data lived on a nice uniform spacing, we could use the findchangepts utility. (It was introduced long after this post was made though. But even so, the scattered spacing in x would preclude the use of findchangepts.)

Ok, next, how should we implement the piecewise linear spline? Since we have only two segments, so a single break, the answer is easy. Use the absolute value function! Since the question was written in terms of the curve fitting toolbox, I'll do the same.

mdl = fittype('a0 + (slope1+slope2)/2*x + (slope2-slope1)/2*abs(x - b)','indep','x');

In general, it is hugely important to have good starting values for the parameters. But we can be pretty lazy here. The break point is probably most important. And even there, I'll just pick a point in the middle of the data. I could probably have used polyfit to generate better estimates of the slopes, but this will probably suffice.

coeff0 = [mean(y),1,1,mean(x)]

coeff0 = 1×4

4.7663 1.0000 1.0000 5.0621

fittedmdl = fit(x,y,mdl,'start',coeff0)

fittedmdl = 
     General model:
     fittedmdl(x) = a0 + (slope1+slope2)/2*x + (slope2-slope1)/2*abs(x - b)
     Coefficients (with 95% confidence bounds):
       a0 =       6.088  (5.983, 6.193)
       b =       4.024  (3.977, 4.071)
       slope1 =       1.972  (1.917, 2.026)
       slope2 =      -1.014  (-1.037, -0.9904)

plot(fittedmdl,x,y,'o')

I've parameterized it so that the two segments have corresponding slopes of slope1 and slope2. The breakpoint is b. As you can see, slope1 was estimated to have value 1.972, when the original curve would have had a slope of 2. Likewise, the second segment was also pretty well estimated at near -1.

[Grace]

Hi, I have a similar task with the original question and tried your code excactly in R2022a. It generated this result:

General model:

fittedmdl(x) = a0 + (slope1+slope2)/2*x + (slope2-slope1)/2*abs(x - b)

Coefficients (with 95% confidence bounds):

a0 = 5.403 (-2.471e+07, 2.471e+07)

b = -5.083 (-7.585e+07, 7.585e+07)

slope1 = 0.427 (-9.731e+06, 9.731e+06)

slope2 = 0.01379 (-0.1184, 0.146)

Similarly, I tried the same approach to my data and had this outcome:

What I need to find are the slopes of the two straight lines and the intersection location. Thanks so much!