Bilinear fit
theigormagician
Sat, 09/03/2011 - 11:38 am
Sat, 09/03/2011 - 11:38 am
This is described as an example in the Igor manual. In Igor 6.22A, you can read about this by executing
DisplayHelpTopic "Conditionals"So how about fitting it with a quadratic function?
September 3, 2011 at 06:49 pm - Permalink
Wave w
Variable x
//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
//CurveFitDialog/ Equation:
//CurveFitDialog/ f(x) = (b0+m0*x)*(x<xlim) + ((b0+m0*xlim) + m1*(x-xlim)) * (x>=xlim)
//CurveFitDialog/ End of Equation
//CurveFitDialog/ Independent Variables 1
//CurveFitDialog/ x
//CurveFitDialog/ Coefficients 4
//CurveFitDialog/ w[0] = b0
//CurveFitDialog/ w[1] = m0
//CurveFitDialog/ w[2] = xlim
//CurveFitDialog/ w[3] = m1
return (w[0]+w[1]*x)*(x<w[2]) + ((w[0]+w[1]*w[2]) + w[3]*(x-w[2])) * (x>=w[2])
End
September 6, 2011 at 09:01 am - Permalink
That example was concocted to illustrate how to enter conditionals into the New Fit Function dialog. As such, it is a pretty minimal example. Among other things, there is nothing that requires the two line segments to have the same Y value at the point where they join.
One thing to watch out for with fitting functions like this is that the coefficient that sets the break point between the two pieces of the line is problematic when it comes to computing derivatives of the function with respect to the coefficients. It will require an epsilon wave.
Here is my standard explanation of what the epsilon wave does:
-------------------
Curve fitting uses partial derivatives of your fit function with respect to the fit coefficients in order to find the gradient of the chi-square surface. It then solves a linearized estimate of the chi-square surface to find the next estimate of the solution (the minimum in the chi-square surface).
Since Igor doesn't know anything about your fit function, it (he?) must use a numerical approximation of the derivatives, which are calculated by finite differences. That is, a model value is calculated at the present estimate of the fit coefficients, then each coefficient is perturbed by a small amount and the derivative is calculated from the difference.
The epsilon wave sets the size of the perturbation used to calculate the derivative. If you don't use an epsilon wave, Igor sets epsilon to:
if your coefficient wave is single precision,
eps[i] = 1e-4*coef[i], unless coef[i] is zero, in which case eps[i] = 1e-4
if your coefficient wave is double precision,
eps[i] = 1e-10*coef[i], unless coef[i] is zero, in which case eps[i] = 1e-10
There are a couple of reasons you might need to explicitly set the epsilon. One is if your fit function is insensitive to a coefficient. That is, perturbing the coefficient makes a very small change in the model value. You can get a situation where the dependence of the model is so small that floating-point truncation results in no change in the model. You have to set epsilon to a sufficiently large value that the model actually changes when the perturbation is applied.
Another case where you need to set epsilon is where the model is discrete or noisy. This can happen if the model involves a table look-up or a series solution of some sort. In the case of table look-up, epsilon needs to be large enough to make sure you get two different values out of the table.
In the case of a series solution, you have to stop summing terms in the series at some point. If the truncation of the series results in less than full floating-point resolution of the series, you need to make sure epsilon is large enough that the change in the model is larger than the resolution of the series. A series might include something like a numerical solution of an ODE (using IntegrateODE). It could also involve FindRoots or Optimize, each of which gives you an approximate result, and run faster if you don't demand high precision.
-------------------
In the case of the break point in the bilinear fit function, you can see that the derivatives won't change much as the break point varies within the space between two X values. It will change a lot when it crosses over an X value. So if the delta for estimating the derivative is small compared to the X increment, it is very likely to fall within the space between two X values. The solution is to use an epsilon value larger than the X increment. Here is an example:
I first pasted the example referenced in the first posting into the Procedure window. Then I made some "data" to fit:
make junk= 1e-3*x^2+gnoise(1)
display junk
Then I tried a fit and it failed with a Singular Matrix error:
W_coef[0] = {0,.1,-20,.2,60}
FuncFit/NTHR=0 PieceWiseLineFit W_coef junk /D
**** Singular matrix error during curve fitting ****
There may be no dependence on these parameters:
W_coef[4]
So I went back to the Curve Fitting dialog and added Epsilon values. You do that on the Coefficient tab: first select New Wave from the Epsilon menu, then you can edit the boxes in the Epsilon column in the Coefficients list. This is the result:
Make/D/O/N=5 eps
W_coef[0] = {0,.1,-20,.2,60}
eps[0] = {1e-6,1e-6,1e-6,1e-6,1}
FuncFit/NTHR=0 PieceWiseLineFit W_coef junk /D /E=eps
Fit converged properly
fit_junk= PieceWiseLineFit(W_coef,x)
W_coef={-0.78659,0.066702,-9.5764,0.19855,68.164}
V_chisq= 154.609;V_npnts= 128;V_numNaNs= 0;V_numINFs= 0;
V_startRow= 0;V_endRow= 127;
W_sigma={0.267,0.00678,0.879,0.00879,3.77}
Coefficient values ± one standard deviation
w_0 =-0.78659 ± 0.267
w_1 =0.066702 ± 0.00678
w_2 =-9.5764 ± 0.879
w_3 =0.19855 ± 0.00879
w_4 =68.164 ± 3.77
The actual values of epsilon aren't particularly critical. I used 1e-6 for the non-critical ones; they should simply be "small" compared to the actual coefficients. The last value I set to 1 because my X values are spaced by 1.0.
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
September 6, 2011 at 10:30 am - Permalink
This reminds me; for instance the nonlinear fitting routines in the GNU Scientific Library (https://www.gnu.org/s/gsl/) require the user to provide a function that calculates the partial derivatives for each of the parameters, at arbitrary locations. Initially I found that annoying, but it turns out that fairly often I can derive analytical expressions for the derivatives, and if not then the reasoning is that the user can simply call the GSL derivative estimators.
I haven't run into issue with this, but I wonder if there is value in allowing the specification of a user-defined function that calculates the partial derivatives explicitly.
September 6, 2011 at 12:41 pm - Permalink
Having said that, if I can find a situation where it provided real benefit I would consider it. I haven't seen a situation where it really was compelling. Some say that analytic derivatives will be calculated faster, but in fact, most derivatives involve many more operations than the basic fit function itself. I have also been told that analytic derivatives are more accurate. That is true, but again, in my experience, the derivatives don't need to be extremely accurate. The most important aspect is that they reliably point downhill in the function gradient.
I may be worth pointing out that in this particular case, it would be very hard to write an analytic derivative for the slope-break coefficient!
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
September 6, 2011 at 05:19 pm - Permalink