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$
3
  • 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, 2016 at 19:23
  • $\begingroup$ Oops, I accidently gave you the non-zero-centered data... $\endgroup$
    – a20
    Nov 29, 2016 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, 2016 at 20:17

1 Answer 1

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$
2
  • $\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$
    – a20
    Nov 29, 2016 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, 2016 at 22:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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