Fitting implicit function or implicit wave?

Share icon

Sat, 11/29/2014 - 04:02 am

Hello Forum Viewers,

I seem to have a different situation from "Fitting Implicit Functions" and "All-At-Once Fitting Functions", or maybe I didn't understand these examples...

My data are:
data_x[0], data_x[1], ...., data_x[n_data-1]
data_y[0], data_y[1], ...., data_y[n_data-1]
cw[0], cw[1], ...cw[n_data-1] // cw is a constant array with each cw[i] corresponding to data_x[i]

My goal is to fit data_y vs data_x using function my(w, x, cw) as follows (of course this does not have the required form, just to show you guys the idea):
//
Function my(w, x, cw): FitFunc
     wave w; variable x; wave cw // not good for FitFunc of course
     variable v1=w[0]
     variable v2=w[1]
//
     variable angle1=0   // angle1 has a constant value
     variable angle2=cw[i]  // angle2 changes with data_x where i=0, 1, ...., n_data-1
//
     variable temp1= my_func1(angle1, v1)    
     variable temp2= my_func1(angle2, v1)    
     variable temp3= temp1*temp2+1 
     variable temp4= my_func2(x, v2)
     return temp3*temp4
End
//
Function my_func1(angle, param1)
variable angle, param1
return (1+cos(angle))/(1+2*cos(angle)*param1)
End
//
Function my_func2(x, param2)
variable x, param2
return (1-param2)/(1+cos(x))
End

In short, my fitting_function can not be expressed as f(w, x) and although data_x and cw do have some functional relationship, an explicit expression would be difficult.

Any suggestions? Thank you very much!

zh  
I don't think you have an implicit function, you just have a function with auxiliary data.

First- I think you need an all-at-once function because you need to enforce one-to-one correspondence between cw[i] and your data[i].

Second- you could simply look up your auxiliary wave inside the fit function:
Function fitfunc(pw, yw, xw)
     Wave pw, yw, xw

     Wave cw          // look up auxiliary wave (must be named exactly "cw")
     ... do your computations ...
     yw = ... whatever ...
end

If you're going to have a variety of cw's depending on circumstances, you may want to encapsulate all that in a structure fit function. But using a structure fit function is hard as you need a driver function to set up the structure and call FuncFit.

A note: if you have an auxiliary wave that requires point-for-point correspondence between your data and the elements in the auxiliary wave, there are some restrictions on what you can do:

Don't use subranges of your input data. You must use all the data in your input waves. If you need to fit a subrange, extract the required data into another wave.

You can't use the /D auto-destination feature of FuncFit. You will have to use /D=. That way you get point-for-point destination values.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

December 1, 2014 at 10:33 am - Permalink

Thanks John for all the details. Now my fitting works by using fitfunc(pw, yx, xw) and by looking up the auxiliary wave inside the fit function. But now it encounters "9 iterations with no decrease in chi square" much more often than before (when using fitfunc(w, x)). In the all-at-once fit function I need to call other user functions which are defined in the form user_func(w, x). Do you think this could cause problems, as when user_func(w, x) is called by all-at-once fit function the variable x is treated as a wave? Thanks.

zh

December 2, 2014 at 08:28 am - Permalink

That message means that the initial guesses are as good as they can be. If you know that's not true, there are other things that can cause it.

Make sure that all your waves are double-precision. Curve fitting is very sensitive to floating-point roundoff errors. Your coefficient wave is used directly in the fit, so it is important that it be double-precision.

If the look-up in your cw wave involves values of limited precision, it's possible that the derivatives aren't being computed accurately. They don't need to extremely accurate, they certainly need to be of the correct sign and non-zero. You may need to use an epsilon wave to get good derivatives. Here's a description of epsilon waves I've sent to other customers:
-------------------------------------
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.
-------------------------------------

If this doesn't help, you might try posting an experiment file containing the data you want to fit, your initial guesses and the code for your fitting function. I'll see if there's anything I can do to help.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

December 2, 2014 at 09:22 am - Permalink

I think I found the problem. It is caused by the complex nature of my fitting function. When the fitting function is strongly dependent on initial guess, it is very likely to get into no convergence situation. Thanks a lot!

December 3, 2014 at 02:59 am - Permalink