Two-phase and Multiphase linear Regression

Here's a user-defined fitting function that will fit an arbitrary number of segments, along with the X values at the intersections of the segments:

Function MultiLineFit(w, x) : fitfunc
    Wave w
    Variable x
    
    Variable nSegs = NumPnts(w)/2
    // format of coefficient wave is
    // w[0] = a1
    // w[1] = b1
    // w[2] = x1
    // w[3] = b2
    // w[4] = x2    (if necessary)
    // w[5] = b3    (only used if x2 is present
    
    Variable result = w[0]
    if (x < w[2])
        return result + w[1]*x
    endif
    
    Variable i
    Variable a2 = w[0]
    for (i = 1; i < nSegs; i += 1)
        Variable xIndex = 2*i
        a2 = (w[xIndex-1] - w[xIndex + 1])*w[xIndex] + a2
        if ( (i == nSegs-1) || (x < w[2*(i+1)]) )
            return a2 + w[xIndex+1]*x
        endif
    endfor
    
    // if we get here, something went wrong
    return NaN
end

The function figures out how many segments you are fitting by the number of points in the coefficient wave (w input). The first coefficient is the y intercept of the first segment. The second is the slope and the third is the X value at the joint between the first and second segments.

After that, the y intercept of each segment is calculated from the slopes and joint values of the previous segments. The fourth is the slope of the second segment and the fifth is the next segment starting X. The sixth would be the next slope, etc. Note that you should not put in the last X value of the data set as the end of the final segment; only the X values between the segments are required.

I tried it on a fake data set created by adding a bit of noise to the output of the fitting function. It worked well for my contrived example. I could imagine it might not work well if the slopes between segments are too similar or if the segment join points don't leave enough points in a segment. It would also fail if one of the X values get somehow set to a higher value than the next. If that happens, you could try a constraint like k2 < .9*k4. The .9 means, don't get too close. You will probably have to adjust that on a case-by-case basis.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

June 3, 2008 at 01:58 pm - Permalink