Nonlinear Fitting Constraints

Share icon

Tue, 11/01/2011 - 05:04 pm

I know nonlinear fitting constraints are not allowed, but does anyone have a clever way around this problem?

Here is what I would like to do:
Curevefit with dblexp_XOffset with nonlinear constraints of "k1*k2 - k3*k4 >0", "k1*k2 - k3*k4 < 0" (really I want k1*k2=k3*k4 )
I have tried holding one variable and performing the curvefit within a while loop until my condition (within a modest tolerance) is met.
I change the variable by a small amount each time the loop runs.

Display/K=1 life
Make/O/T/N=2 const = {"k2 >0", "k4 > k2", "k4 > 0"}
Variable V_Fiterror=0
Curvefit/Q/NTHR=0/TBOX=768/W=2 dblexp_XOffset life(9.7E-9, rightx(life)*0.95)/C=const/NWOK
Variable dum
if (V_FitError==0)
do
    if (k3*k4>k1*k2)
        dum = k4*0.99
        const[2] =  "k4 <" + num2str(dum)
    else
        dum = 1.01*k4
        const[2]= "k4 >" +  num2str(dum)
    endif
   
    Curvefit/Q/NTHR=0/TBOX=768/W=2/G dblexp_XOffset life(9.7E-9, rightx(life)*0.95)/C=const
    if (V_FitError!=0)
        break
    endif
while (abs(k1*k2 - k3*k4) > 0.05*(k3*k4))
endif


This works with spotty success - I figure there must be a better way.
FYI: "Life" is a 1D wave containing lifetime data, which is scaled in nanoseconds (using Igor 6.2.2.0).
There isn't a good way to do that with curve fitting.

You might be able to come up with an implementation of a penalty function that would do it. You would write a function that evaluates the model, plus the constraint equation, and computes a penalty score. The penalty score might be the basic sum of squared errors plus something that measures the amount by which your constraint is violated. Then you would use the Optimize operation to minimize the penalty function.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

November 1, 2011 at 05:16 pm - Permalink

johnweeks wrote:

You might be able to come up with an implementation of a penalty function that would do it. You would write a function that evaluates the model, plus the constraint equation, and computes a penalty score. The penalty score might be the basic sum of squared errors plus something that measures the amount by which your constraint is violated. Then you would use the Optimize operation to minimize the penalty function.


Gencurvefit,cannot do constraints, but it does apply upper and lower bounds for a given parameter. In addition, you can specify a cost function for it to use whilst it does the curvefit. By default this is an internal Chi2 implementation but you can do what John suggests and apply a penalty.

Function Chi2(coefs, y_obs, y_calc, s_obs)
    Wave coefs, y_obs, y_calc, s_obs
    make/n=(numpnts(y_obs))/free/d diff
    diff = ((y_obs-y_calc)/s_obs)^2
    return sum(diff)
end

Function mypenalty(coefs, y_obs, y_calc, s_obs)
    Wave coefs, y_obs, y_calc, s_obs
    variable penaltyconst = 10
    variable penaltypower = -1
    make/n=(numpnts(y_obs))/free/d diff
    diff = ((y_obs-y_calc)/s_obs)^2
    return sum(diff) + penaltyconst * (k1*k2 - k3*k4)^penaltypower
End

November 1, 2011 at 05:47 pm - Permalink

To make the fitting faster, you can write your own fitting function. You only need three parameters rather than four.

Use fitting parameters K1, K2 and K3; and replace K4 by K1*k2/K3.

Also set the constraint to make sure K3 != 0.




BMangum wrote:
I know nonlinear fitting constraints are not allowed, but does anyone have a clever way around this problem?

Here is what I would like to do:
Curevefit with dblexp_XOffset with nonlinear constraints of "k1*k2 - k3*k4 >0", "k1*k2 - k3*k4 < 0" (really I want k1*k2=k3*k4 )





Yaohua Liu
Neutron & X-ray scattering/ Materials Science/ Argonne National Lab

November 2, 2011 at 08:35 am - Permalink

Thanks for all the advice.
I think Yaohua made a great point, in my case I can simply re-write the fitting function such that my nonlinear constraints are met by are in the fit function itself.
This gaurantees that my condition is satisfied.

In the parlance of the double exponential fit I now use A2 = T1*A1/T2 (k3 = k1*k2/k4 for dblexp_XOffset)

I also decided to add an extra variable to allow for a little tolerance and slightly better fits.
Here is my new fit function:
Function dbl_exp_mod(w,x):fitfunc
    Wave w
    variable x
    return w[0]+w[1]*exp(-((-w[5]+x)/w[2])) + w[4]*(w[1]*w[2]/w[3])*exp(-((-w[5]+x)/w[3]))
End

In this new fit function I now place the following constraints:
"k4 < 1+tol"
"k4 > 1-tol"
It is working fine now.
Thanks again!

November 2, 2011 at 09:05 am - Permalink