0
$\begingroup$

I'm trying to fit a Voigt distribution to some measurements, but no matter what I try Mathematica simply never finishes, even if I do MaxIterations -> 1 and very low precision goals. For example you can take this data:

DD = 100; theta2B = 13.742; Lamda = 0.999047433;
ϕ[DD_, theta2_] := 2 Pi DD (Sin[theta2/360*Pi] - Sin[theta2B/360*Pi])/Lamda;
profile = 1/ϕ[DD, theta2]^2 - Sin[2 ϕ[DD, theta2]]/ϕ[DD, theta2]^3 + (1 - 
  Cos[2 ϕ[DD, theta2]])/(2 ϕ[DD, theta2]^4);
data = Transpose[{Table[i, {i, 7, 20, 0.01}], Table[profile, {theta2, 7, 20, 0.01}]}];
ListPlot[data, PlotRange -> All]

And function definition:

voigt[I_, δ_, σ_, x_] := I*PDF[VoigtDistribution[δ, σ], {x}];

And then fitting using:

Clear[a, x, δ, σ];
model3 = voigt[a, δ, σ, x];
fit3 = NonlinearModelFit[data, 
  model3, {{a, 1}, {δ, 1}, {σ, 1}}, x, 
  Method -> NMinimize];
a = fit3[[1, 2, 1, 2]];
δ = fit3[[1, 2, 2, 2]];
σ = fit3[[1, 2, 3, 2]];
fitplot2 = LogPlot[model3, {x, 0, 8.5}]

The only error message that I sometimes get is that the output is unreal. The output from voigt[] is on the form {y + .0i} (I don't know why it's in {}), so I've tried to only take the real part.

$\endgroup$
  • 2
    $\begingroup$ The VoigtDistribution is peaked around zero (Median[VoigtDistribution[a, b]] = 0; play also with Manipulate[ Plot[voigt[a, b, c, x], {x, -5, 5}], {a, 1, 20}, {b, 1, 20}, {c, 1, 20}]) and your data is peaked (rather narrowly) around 14. $\endgroup$ – corey979 Nov 29 '16 at 19:23
  • $\begingroup$ Oops, I accidently gave you the non-zero-centered data... $\endgroup$ – bjorn Nov 29 '16 at 20:12
  • 2
    $\begingroup$ In your definition of voigt you use capital I, which is the complex number I. I suggest you use different variable (E, D, I are symbols). Also, PDF[VoigtDistribution[δ, σ], {x}] threads PDF over the list {x} - if you want value of PDF at a point, use PDF[VoigtDistribution[δ, σ], x]. $\endgroup$ – Gosia Nov 29 '16 at 20:17
3
$\begingroup$

First, the VoigtDistribution is centered around zero, and the OP's data has a narrow peak at around 14. Let's shift it to zero (crudely):

pos = Max@Position[(Transpose@data)[[2]], (Transpose@data)[[2]] // Max]

675

peak = (Transpose@data)[[1, pos]]

13.74

data = Transpose[Transpose@data - {peak, 0}];
plot = ListPlot[data, PlotRange -> All]

enter image description here

Second - and important - it is a matter of good starting values for the parameters for NonlinearModelFit. Playing with

Manipulate[
 Plot[voigt[a, b, c, x], {x, -5, 5}, PlotRange -> {{-6, 6}, {0, 1}}],
  {a, 0.1, 1, 0.1}, {b, 0.1, 1, 0.1}, {c, 0.1, 1, 0.1}]

I choose {1, 0.2, 0.2}.

model3 = voigt[a, b, c, x];
fit3 = NonlinearModelFit[data, model3, {{a, 0.8}, {b, 0.2}, {c, 0.2}}, x]

which with no errors gives

Normal@fit3

enter image description here

a = fit3[[1, 2, 1, 2]]
b = fit3[[1, 2, 2, 2]]
c = fit3[[1, 2, 3, 2]]
fitplot2 = Plot[model3, {x, -10, 10}, PlotRange -> All, PlotStyle -> Red];

0.742524

0.0809207

0.223623

and

Show[fitplot2, plot]

enter image description here

$\endgroup$
  • $\begingroup$ Thank you! I accidentally gave you the wrong data (it should of course be zero-centered). With your method it's working for me now, but when I used NMinimize method before submitting here on this relatively simple data it still hadn't finished after one hour. $\endgroup$ – bjorn Nov 29 '16 at 20:25
  • 1
    $\begingroup$ First, it's Method -> "NMinimize", not Method -> NMinimize as in your code (see the docs). Second, "NMinimize" indeed runs long enough to not wait for the result. On the other hand, "Newton" works (although takes twice as long as the automatic methods). You can also play with other methods (listed in the docs). "LevenbergMarquardt" also works in a time even slightly shorter than the automatic method. $\endgroup$ – corey979 Nov 29 '16 at 22:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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