I have to fit two data sets to two functions(eps1,eps2) with same parameters (w0,f,g)

Initially, I know only one function

eps2[w]=f*w0^2*g*w/((w0^2 - w^2)^2 + γ^2*w^2)

and with Kramers-Kronig relation:

eps1[w]=epsinf + 1/Pi*NIntegrate[eps2[W]/(W-w),
  {W,-Infinity,Infinity},Method -> "PrincipalValue", Exclusions -> Automatic];

So i have two functions eps1 and eps2.


model[index_, w_?NumericQ, epsinf_?NumericQ, g_?NumericQ, 
  w0_?NumericQ, f_?NumericQ] :=KroneckerDelta[index - 1]*eps1 +  
  KroneckerDelta[index - 2] *eps2

allData = Join[{1, Sequence @@ #} & /@ re, {2, Sequence @@ #} & /@ im]; 

% re and im are the data sets

nlm = NonlinearModelFit[allData, model[index, w, epsinf, g, w0, 
    f], {{epsinf, 1}, {w0, 3}, {f, 5}, {g, 1}}, {index, w}, MaxIterations -> 1000];

...It works well but the problem is that it's too slow . I 'am wondering if I could make something to run faster. It takes 20min and if I add more parameters it takes much more.

  • $\begingroup$ Your code doesn't run. It has errors in how the functions are defined. Without some example data, it's not possible to execute the fit. $\endgroup$ – Jens Mar 18 '15 at 23:07
  • $\begingroup$ I have spell-checked the code to the best of my knowledge in the field and can provide a semi-working example (because even reducing iterations still has the NonlinearModelFit work excruciatingly slowly and throw a bunch of errors), as well as generate a set of fake data. I'm still adamant, that OP should take all possible steps to eliminate the KK-transform. Most terms used for fitting optical data do have a closed-form expression for their KK-transform. $\endgroup$ – LLlAMnYP Mar 20 '15 at 11:30

Are we colleagues, by any chance?

Anyway, there's a well known closed-form expression for eps1. Rather than doing the very costly NIntegrate to get eps1 from eps2 (which is reevaluated many, many times during the execution of NonLinearModelFit, I suggest you write down the analytical expression for eps1.

For the form of eps2 that you show here, you can define eps == eps1 + I eps2, where eps[w] == epsinf + w0^2 * f/(w0^2-w^2-I g w) (I am assuming, that you call what conventionally is known as "delta-epsilon" the variable f, and also that the greek gamma in your definition of eps2 is a typo and should be g or vice-versa).

Having eliminated NIntegrate, the code should run much faster, but a complex function will require proper adjustment of the NormFunction in the fitting procedure (I'm not quite sure, you might have to use FindFit instead of NonLinearModelFit).


Even with restricted datasets (I imagine, this would be some reflectivity data) the general idea holds. In the unfortunate case that you have only eps2 available, there is no use for a Kramers-Kronig transformation as you only have the imaginary part to fit to. In that case just fit as normal, model function with some parameters to known dataset of just eps2. Without NIntegrate things will go much faster.

Usually, if you have eps2 in some range, you'll also be able to do a Kramers-Kronig transformation to get eps1 in a somewhat narrower range. Extra data never hurts, although it will be somewhat ambiguous (high-energy oscillators can contribute an absolute offset, for example). Let eps_infinity be an extra fitting parameter, but bear in mind, that it will be unrelated to the "true" value of eps_infinity.

Even if you are using functions other than harmonic oscillators, most typical terms for fitting optical data will have a closed-form Kramers-Kronig transformation. If you permit yourself to use extremely arbitrary terms, you might as well just abandon fitting altogether.

Finally, in case you have a term for which you have a closed-form expression for the imaginary part, but not for the real part, there is a trick you should try to employ before resorting to calculating the Kramers-Kronig transformation over and over:

  • Let your term (imaginary part) be defined as term2[x,a,b,c...] where x is frequency, a is amplitude, b,c... are the remaining parameters.
  • Get the Kramers-Kronig transformation at specific values of parameters. That is, define the real part as term1[x,a0,b0,c0...] == KKTransform[term2[x,a0,b0,c0...].
  • Kramers-Kronig transformation gives you the real part at a fixed frequency, so you'll have to calculate it for several frequencies, get a list of datapoints of the form of {x1,term1[x1,a0,b0,c0...]}. Construct an InterpolationFunction from these. Now you have a Function object for term1[x[a0,b0,c0...].
  • Now you need to study, how it behaves with varying parameters. The first one, amplitude, is easy. You can easily state that term1[x,a,b0,c0...] == a*term1[x,a0,b0,c0...]/a0.
  • For the remaining parameters grab a pencil and paper and look closely at how the KK-transform integral changes with changing parameters. That will depend on the exact form of your term.
  • With some luck you'll arrive at a closed-form or InterpolationFunction-form of the KK-transform of your term. That's already much better than NIntegrate.
| improve this answer | |
  • $\begingroup$ thank you for answering me ! I am undergraduate physics student ! $\endgroup$ – Christos Mar 18 '15 at 20:37
  • $\begingroup$ Yes, i understand this ,but sometimes only eps2 is given . I dont have always harmonic oscillators. $\endgroup$ – Christos Mar 18 '15 at 20:43

Not the answer you're looking for? Browse other questions tagged or ask your own question.