Re: Fast vs. Slow NonlinearModelFit models
- To: mathgroup at smc.vnet.net
- Subject: [mg123990] Re: Fast vs. Slow NonlinearModelFit models
- From: "Dan O'Brien" <danobrie at gmail.com>
- Date: Wed, 4 Jan 2012 05:02:48 -0500 (EST)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <201112211152.GAA02211@smc.vnet.net> <201112220926.EAA13844@smc.vnet.net>
Ok this idea occurred to me, but I can't expand it so that it's a good
approximation over a large enough domain. Any other thoughts?
On 12/22/2011 3:26 AM, DrMajorBob wrote:
> You might fit using a short Series expansion in place of Erfc, then use
> the fitted parameters as initial guess in a second fit for the exact model.
>
> Bobby
>
> On Wed, 21 Dec 2011 05:52:34 -0600, Dan O'Brien<danobrie at gmail.com> wrote:
>
>> Hello,
>>
>> I'm wondering if anyone has a suggestion for speeding up the
>> NonlinearModelFit computation for model2 below. The model used to fit
>> my experimental data can be a simple complex function who's magnitude
>> squared is a lorentzian (model1), or, some may argue a more complete
>> model is the same function convolved with a gaussian which gives a
>> function that involves the complementary error function (model2), the
>> Voigt profile. As indicated below with AbsoluteTiming using these two
>> different models show VERY different computational times. I imagine it
>> has something to do with mathematica's internal implementation of Erfc
>>
>> In[1]:= $Version
>> (*Lorentzian Lineshape*)
>> \[Chi]R[\[Omega]_, {A_, \[Phi]A_, \[CapitalGamma]_, \[Omega]v_}] := (
>> A E^(I \[Phi]A))/(-(\[Omega] - \[Omega]v) - I \[CapitalGamma]);
>>
>> (*Voigt Lineshape*)
>> z[\[Sigma]_, \[CapitalGamma]_, \[Omega]v_, \[Omega]_] := (\[Omega] - \
>> \[Omega]v + I \[CapitalGamma])/(Sqrt[2] \[Sigma])
>> FadeevaF[z_] := Exp[-z^2] Erfc[-I z]
>> IH\[Chi]R[\[Omega]_, {\[Sigma]_,
>> A_, \[Phi]A_, \[CapitalGamma]_, \[Omega]v_}] :=
>> I Sqrt[\[Pi]/2] (
>> A E^(I \[Phi]A))/\[Sigma] FadeevaF[
>> z[\[Sigma], \[CapitalGamma], \[Omega]v, \[Omega]]]
>> (*Redefine functions so A is peak amplitude of the function magnitude \
>> squared*)
>> A\[Chi]R[\[Omega]_, {\[Sigma]_,
>> A_, \[Phi]A_, \[CapitalGamma]_, \[Omega]v_}] := \[Chi]R[ \[Omega], \
>> {Sqrt[Abs[A]] Abs[\[CapitalGamma]], \[Phi]A,
>> Abs[\[CapitalGamma]], \[Omega]v}]
>> AIH\[Chi]R[\[Omega]_, {\[Sigma]_,
>> A_, \[Phi]A_, \[CapitalGamma]_, \[Omega]v_}] :=
>> IH\[Chi]R[\[Omega], {\[Sigma], (
>> Sqrt[Abs[A]] E^(-(1/2) \[CapitalGamma]^2/\[Sigma]^2) Sqrt[2/\[Pi]]
>> Abs[\[Sigma]])/
>> Erfc[Abs[\[CapitalGamma]]/(
>> Sqrt[2] Abs[\[Sigma]])], \[Phi]A, \[CapitalGamma], \[Omega]v}]
>> (*Generate data with noise*)
>> dataplot =
>> ListPlot[data =
>> Table[{\[Omega],
>> Abs[AIH\[Chi]R[\[Omega], {4, 1.1, 0, 2.1, 0}] +
>> AIH\[Chi]R[\[Omega], {4, .6, \[Pi], 2.1, 18}]]^2 +
>> RandomReal[{-0.05, 0.05}]}, {\[Omega], -30, 50, .5}],
>> Joined -> True]
>> (*Models*)
>> model1 = Abs[
>> A\[Chi]R[#, {4, a1, 0, \[CapitalGamma]1, \[Omega]vo1}] +
>> A\[Chi]R[#, {4, a2, \[Pi], \[CapitalGamma]2, \[Omega]vo2}]]^2&;
>> model2 =
>> Abs[AIH\[Chi]R[#, {4, a1, 0, \[CapitalGamma]1, \[Omega]vo1}] +
>> AIH\[Chi]R[#, {4, a2, \[Pi], \[CapitalGamma]2, \[Omega]vo2}]]^2&;
>>
>> (*Fit*)
>> AbsoluteTiming[
>> mod1 = NonlinearModelFit[data,
>> model1[\[Omega]], {{a1, 1}, {\[CapitalGamma]1, 2}, {\[Omega]vo1,
>> 1}, {a2, .8}, {\[CapitalGamma]2, 2}, {\[Omega]vo2,
>> 15}}, \[Omega]]][[1]]
>> AbsoluteTiming[
>> mod2 = NonlinearModelFit[data,
>> model2[\[Omega]], {{a1, 1}, {\[CapitalGamma]1, 2}, {\[Omega]vo1,
>> 1}, {a2, .8}, {\[CapitalGamma]2, 2}, {\[Omega]vo2,
>> 15}}, \[Omega]]][[1]]
>>
>> Show[{dataplot,
>> Plot[{Normal[mod1[\[Omega]]],
>> Normal[mod2[\[Omega]]]}, {\[Omega], -30, 50}]}]
>>
>> Out[1]= "8.0 for Microsoft Windows (32-bit) (October 6, 2011)"
>>
>> Out[8]=****Plot Snipped
>>
>> Out[11]= 0.0390625 *Here is AbsoluteTiming for model1*
>>
>> Out[12]= 75.0078125 *Here is AbsoluteTiming for model2*
>>
>> Out[13]=****Plot Snipped
>