Curve fitting with functional variables

Share icon

Wed, 09/02/2009 - 01:23 pm

I have some curve fitting that I have been doing in Mathematica, but would like to do in Igor. I have attempted this in the past, but ran into problems because of the complexity and length of the function I am trying to fit to (it would easily take over a page to print out). Basically, the expression is simplified by grouping of terms into functional variables. In Mathematica, I can define these variables as needed, and they allow a concise expression of the final function to which I am fitting. In Igor, I don't see how I can accomplish the same thing. It looks like I have to write it all out without any simplification. I can do this, but it is a painstakingly slow process. In addition, when I have done it in the past, Igor seems to choke on the length of the function.

To give a short example of what I am talking about:

a = (1/b+1/c)^-1

d= (Exp[-w^2 a]

f= Log[sqrt[cosh^2[d]]]

This is anything close to the actual function I am dealing with, but those details are (I believe) unimportant. What is important is a final function ("f") that is dependent on other "functional" variables. How does one deal with them within Igor? What suggestions does the community have for this type of situation? This is fairly straightforward to do in Mathematica, but inconvenient as one has to export the data just right, and set up the directory for Mathematica to import, etc. In addition, graphs are VERY difficult to manipulate in Mathematica, so I'd rather not have to resort to using it if possible. I'm really not sure where to go from here.
patarroyo wrote:
I have some curve fitting that I have been doing in Mathematica, but would like to do in Igor. .... when I have done it in the past, Igor seems to choke on the length of the function.

To give a short example of what I am talking about:

a = (1/b+1/c)^-1

d= (Exp[-w^2 a]

f= Log[sqrt[cosh^2[d]]]


You do not say what your independent parameter is. All you give are fitting coefficients. I am going to assume that w in your expression is the independent parameter. Also, why you have a sqrt of a squared value is a mystery to me. I assume this is an illustration and will make it just sqrt of cosh. Finally, the only two coefficients that I see in this are b and c.

// ww[0] = b
// ww[1] = c

Function MyFitFunction(ww,xx)
     wave ww
     variable xx

     variable a, d, f

     a = (1/ww[0]) + (1/ww[1])
     d = exp(-xx^2/a)
     f = log(sqrt(cosh(d)))

    return f

end


Basically, you can define the parameters within your fit function to collapse all complexities in step-by-step fashion, until you are returning from the fit function the value of just one "internal function".

--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAH

September 2, 2009 at 03:52 pm - Permalink

jjweimer wrote:

Also, why you have a sqrt of a squared value is a mystery to me.


Don't forget that the square root of a squared real number is always positive. Or, he could've just used abs()

September 2, 2009 at 10:07 pm - Permalink

Thanks for the help. I just made up a function off the top of my head, hence the cosh^2 sqrt mistake. What I have are basically 5 dependent variables, one independent, and a messy function that relates them all together. I am going to try to do what JJ suggested, though I still have my concerns about the length of this function... only one way to find out!

Update: Here is the code I am trying to implement:

#pragma rtGlobals=1     // Use modern global access method.
//ww[0]=deltaW
//ww[1]=t_a
//ww[2]=t_b
//ww[3]=R_a
//ww[4]=R_b


Function ChemExFit(ww,tcp)
     wave ww
     variable tcp
 
     variable aminus, aplus, psi, zeta,xi,eta,Dplus,Dminus,lambda1,R21
     
     aminus=ww[3]-ww[4]+1/ww[1]-1/ww[2]
     aplus=ww[3]+ww[4]+1/ww[1]+1/ww[2]
     psi=aminus^2-ww[0]^2+4/(ww[1]*ww[2])
     zeta=2*ww[0]*aminus
     xi=tcp*sqrt((psi+sqrt(psi^2+zeta^2)))/(2*sqrt(2))
     eta=tcp*sqrt((-psi+sqrt(psi^2+zeta^2)))/(2*sqrt(2))
     Dplus=0.5*(1+(psi+2*ww[0]^2)/(sqrt(psi^2+zeta^2)))
     Dminus=0.5*(-1+(psi+2*ww[0]^2)/(sqrt(psi^2+zeta^2)))  
     lambda1=Exp(-tcp*aplus/2+ln(sqrt(Dplus*(Cosh(xi))^2-Dminus*(Cos(eta))^2)+sqrt(Dplus*(Sinh(xi))^2+Dminus*(Sin(eta))^2)))
     R21=-ln(lambda1)/tcp
 
   return R21
end


Basically, R21 is the function that I am trying to fit my data to. tcp is my independent variable. The wave are the fit parameters. The last time I tried this I actually tried typing in the expanded "R21" as a custom fit with no luck. I get the feeling that a lot of the concepts here are pretty basic for Igor, but I have just had no reason to deal with such a complicated function in Igor before, hence my ignorance. I guess i am kind of stuck of where to go from here though as far as using Igor's curve fitting functionality.

September 3, 2009 at 09:26 am - Permalink

I think I got it! I saw another JJ reply to a question regarding curve fitting. I just added ": FuncFit" after my function name, and that caused it to appear in the curve fit dialog! Most excellent!

-Mike

September 3, 2009 at 11:10 am - Permalink

Problem: curve fit works great, but my fit abruptly drops to negative infinity past a certain value of the independent variable. I believe this clipping is caused by a precision error of some kind. Is there a way for me to get around this type of issue?

September 14, 2009 at 03:02 pm - Permalink

patarroyo wrote:
Problem: curve fit works great, but my fit abruptly drops to negative infinity past a certain value of the independent variable. I believe this clipping is caused by a precision error of some kind. Is there a way for me to get around this type of issue?


Well, cosh(x) and sinh(x) both become very large at rather modest values of x. Is that perhaps part of the problem?

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

September 15, 2009 at 03:58 pm - Permalink

Yes I believe it is. The question then becomes is there a way around that? It would seem that I am being hit by rounding errors here. Suggestions?

September 16, 2009 at 11:09 am - Permalink

patarroyo wrote:
Yes I believe it is. The question then becomes is there a way around that? It would seem that I am being hit by rounding errors here. Suggestions?


At some point you simply have to accept the limitations of floating-point arithmetic. Even though floating-point numbers often give the impression of being continuous, they actually have finite resolution.

Having said that, there are sometiimes ways to get around such problems by careful ordering of operations (like do divisions before mutliplications if that will tend to keep the results nearer to 1). You may be able to scale the problem to keep the magnitudes nearer to 1d.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

September 17, 2009 at 03:21 pm - Permalink