What's Wrong with Numerical Recipes' MEDFIT Routine?

The problem in MEDFIT is more in the theory than in the code. It's hard to point to one place in the code and say AhA! The author of MEDFIT assumed that the minimum of the residual will occur at a place where the derivative is zero. But for L1 approximation, that's not true: the approximating function is only C0 continuous, so there are places where the derivative doesn't exist, namely at every data point. The criterion to use is spelled out in section 4-4 of "The Approximation of Functions -- Vol 1: Linear Theory" by John R. Rice, Addison-Wesley, 1964 (102ff): a minimum occurs where an arbitrary perturbation increases the residual. In the case of a continuous and continuously differentiable approximant, this is the same as saying the derivative is zero, but it's not the case with the approximant for L1 fitting. So the first defect in MEDFIT is in assuming that the minimum occurs at a zero of the derivative of the approximating function.

The second defect arises in the transformation of the flawed theory into an algorithm. To understand this problem, consider first the problem of one parameter L1 fitting. In this case, what one wants is the median of the set. But if the size of the set is even, anywhere between the two media is an equally good solution. The algorithm in NR, however, considers only solutions that interpolate one or two data points. (In general, L1 fitting to N parameters will interpolate N data points). The result is that MEDFIT is always computing a derivative at a place where it doesn't exist. It turns out that what MEDFIT really needs is just a sign, and there's a 50-50 chance of getting that right, even if the derivative doesn't exist. But when there's just a little bit of data, and MEDFIT gets close to the solution, it will get alternating positive and negative derivatives, and just rattle back and forth between two points, one probably being a solution (but at which MEDFIT can't decide to quit).

The third kind of defect arises in the transformation of the flawed algorithm into code. Others have already pointed out these defects, which can be spotted even if one has no idea what the code is intended to do.

Incidentally, the Fortran and C versions could frequently have different behavior because the tests used in the C version in place of the Fortran SIGN function do not do the same as what the Fortran-77 standard says the SIGN function does.