Curve fitting to NaN for certain x-values
Hello Igor users, I am trying to fit a user-built three-exponential function to a data set to extract kinetic information. The curve fit dialog is below:
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) = A3*((ERFEXP(KB,D,(1/(K1+K2+K3+(K1*4))),X)*(K1+K2+K3+(K1*4)))*K1*K2*(ERFEXP(KB,D,(1/(-K2-K3)),X)*(K2-K3)+ERFEXP(KB,D,(1/(-K1-K2-(K1*4))),X)*(K1-K2+(K1*4))-ERFEXP(KB,D,(1/(-K1-K3-(K1*4))),X)*(K1-K3+(K1*4)))/((K2-K3)*(K2-K1*(1+4))*(K3-K1*(1+4))))+\
//CurveFitDialog/ OS
//CurveFitDialog/ End of Equation
//CurveFitDialog/ Independent Variables 1
//CurveFitDialog/ X
//CurveFitDialog/ Coefficients 7
//CurveFitDialog/ w[0] = KB
//CurveFitDialog/ w[1] = D
//CurveFitDialog/ w[2] = OS
//CurveFitDialog/ w[3] = A3
//CurveFitDialog/ w[4] = K1
//CurveFitDialog/ w[5] = K2
//CurveFitDialog/ w[6] = K3
return w[3]*((ERFEXP(w[0],w[1],(1/(w[4]+w[5]+w[6]+(w[4]*4))),X)*(w[4]+w[5]+w[6]+(w[4]*4)))*w[4]*w[5]*(ERFEXP(w[0],w[1],(1/(-w[5]-w[6])),X)*(w[5]-w[6])+ERFEXP(w[0],w[1],(1/(-w[4]-w[5]-(w[4]*4))),X)*(w[4]-w[5]+(w[4]*4))-ERFEXP(w[0],w[1],(1/(-w[4]-w[6]-(w[4]*4))),X)*(w[4]-w[6]+(w[4]*4)))/((w[5]-w[6])*(w[5]-w[4]*(1+4))*(w[6]-w[4]*(1+4))))+w[2]
End
The x-data runs from -50 to 4000 (time) and the function yields "NaN" at about 250 (note the x-data is not evenly spaced, the steps are close together close to zero and get farther apart farther from zero). I do not think I am generating 0/0 in the fit and I can't quite figure out where my bug is.
Note that "erfexp" is a user-designed exponential function that is the convolution of an error function with a Gaussian, KB is the gaussian FWHM, D is "time zero", A3 is an amplitude, and the k's are rate constants. Do you have any debugging suggestions? Many thanks!
I would add a debugging section to catch the NaN in ERFEXP, like this:
if( numtype(w0) != 0 || numtype(w1) != 0 || numtype(w2) != 0 || numtype(xx) != 0 )
Debugger
endif
Variable result = ... your code here
if( numtype(result) != 0 )
Debugger
endif
return result
End
Then look at the coefficients in the calling routine (THIRDSTATE) to see who's to blame.
December 9, 2021 at 12:48 pm - Permalink
It's hard for my eyes and brain to understand such a long complicated expression. I would be tempted to break it up into multiple terms assigned to variables, and then combine those terms in a shorter, simpler expression. It looks to me like you have a situation where you are multiplying one ERFEXP() by another and another. Your description sounds like they should be added? But as I said, I can't really parse that expression without a lot of pain :)
I would also add a similar catcher in THIRDSTATE, like this, just in case something in that expression is at fault.
if (numtype(result))
Debugger
endif
return result
end
The Curve Fitting dialog won't care about the extra code (though you might see it when it shows you the function).
December 9, 2021 at 04:15 pm - Permalink
Both erf() and Gauss() are well-behaved for large values of X (in both directions), but you don't show us what's in ERFEXP(), so we can"t tell if erf() and gauss() might be combined in a bad way.
The exp() function itself will overflow at fairly modest values of X.
December 9, 2021 at 04:17 pm - Permalink