3
$\begingroup$

I'm trying to fit a Cole-Davidson formula:

CD = a(Sin[b1*ArcTan[2 \[Pi]*f*t1]]/(2 \[Pi]*f (1 + (2 \[Pi]*f*t1)^2)^{(b1/2))}+4 (Sin[b1*ArcTan[4 \[Pi]*f*t1]]/(4 \[Pi]*f (1 + (4 \[Pi]*f*t1)^2)^{(b1/2)}))), 

(a, b1 and t1 are parameters) to experimental data set (in the form: {Amplitude, Frequency}:

data={{40003000, 2.0216}, {34219000, 2.7139}, {29275000, 
  3.5391}, {25042000, 4.6468}, {21420000, 6.2962}, {18326000, 
  8.5129}, {15675000, 12.868}, {13411000, 18.275}, {11470000, 
  21.079}, {9809700, 28.439}, {8394500, 41.372}, {7179000, 
  59.07}, {6141300, 69.911}, {5253900, 105.59}, {4496500, 
  134.92}, {3844800, 180.5}, {3290300, 201.62}, {2813700, 
  402.62}, {2408900, 379.57}, {2060600, 750.51}, {1762800, 
  628.69}, {1507800, 928.53}, {1290700, 999.34}, {1103900, 
  1432.5}, {1000700, 1207.9}, {944060, 1796.1}, {826930, 
  1529.6}, {682810, 1599.4}, {591110, 2049.9}, {564410, 
  2235.5}, {506030, 2156.5}, {467050, 2008}, {385810, 1855}, {319510, 
  2476.7}, {316810, 2650.9}, {270960, 1634.3}, {263540, 
  1825.8}, {218140, 2599.8}, {180400, 2007.3}}

To fit model to data set I used:

FindFit[data, CD, {a, t1, b1}, f, Method -> "X"]

where X denotes methods a had tried (Newton, Lavenberg-Marquardt etc). And here is my problem: I cannot get a good fit to the experimental data. It's not even close fit - the curve and the data points dont match at any point. I tried to change max iterations, put constraints on parameters - no luck. Any thoughts about whats wrong with my simple code are appreciated.

$\endgroup$
  • 1
    $\begingroup$ Your CD is not syntactically valid. $\endgroup$ – David G. Stork Feb 4 '16 at 23:42
  • $\begingroup$ Check that your parentheses are balanced in CD. Note that in Mathematica curly braces {} are used to identify lists, and not to modify operation precedence. You should only use parentheses () for that. In all honesty, I am not sure how you even got a fit at all. $\endgroup$ – MarcoB Feb 5 '16 at 0:42
5
$\begingroup$

Re-writing you function:

fun[a_, b_, w_, t_] := 
 a (Sin[b*ArcTan[w t]]/(w (1 + (w t)^2)^(b/2)) + 
    4 (Sin[b*ArcTan[2 w t]]/(2 w (1 + (2 w t)^2)^(b/2))))

Modifying data:

datam = {2 Pi #2, #1} & @@@ data;

After playing with manipulate tentative starting guesses:

nlm = NonlinearModelFit[datam, 
  fun[a, b, w, t], {{a, 10^13}, {b, 0.1}, {t, 0.1}}, w]

Parameter estimates:

nlm["BestFitParameters"]

yields: {a -> -5.74097*10^7, b -> -0.462728, t -> 1.24964}

Rescaling function (from rad/s to Hz) and visualizing:

rf[f_] := Normal[nlm] /. w -> 2 Pi f
Show[LogLogPlot[rf[x], {x, 1, 2100}], ListLogLogPlot[Reverse /@ data]]

enter image description here

I apologize if I have misinterpreted but perhaps this might be useful. (noting: (i) parameter estimates may violate desired constraints)

$\endgroup$
  • $\begingroup$ Thanks a lot! Your comments were crucial for solving the problem. It seems that there is at least one more physical process "behind" the data, and the simple CD formula is not enough to fitting the curve properly to the data set. $\endgroup$ – user19388 Feb 5 '16 at 8:32
4
$\begingroup$

This is an extended comment on @ubpdqn 's answer. First, I am completely ignorant about the Cole-Davidson formula. But if one plots all of the data using

Show[ListLogLogPlot[Reverse /@ data, PlotRange -> Full], 
 LogLogPlot[rf[x], {x, Min[data[[All, 2]]], Max[data[[All, 2]]]}]]

one sees a lack of fit:

Lack of fit figure

I ask: Is that lack of fit in that one section of the curve important or not?

On the original scale the lack of fit is less noticeable:

Show[ListPlot[Reverse /@ data, PlotRange -> Full], 
 Plot[rf[x], {x, Min[data[[All, 2]]], Max[data[[All, 2]]]}]]

Figure on original scale

I'd argue that the answer about the apparent lack of fit depends on the objective associated with the predictions rather than anything intrinsic in the data.

$\endgroup$
  • $\begingroup$ I agree entirely with your comment and think it is very important. A power law springs to mind and you are right about extreme value loss of fit. +1. My principal aim was to act as a catalyst...Thanks :) $\endgroup$ – ubpdqn Feb 5 '16 at 6:29

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.