3
$\begingroup$

I'm trying to fit formula:

    modelx[d1_?NumberQ, dyf1_?NumberQ, f_] := 6.62*10^-21 (NIntegrate[72/d1^3 (u^2/(81 + 9 u^2 - 2 u^4 + u^6))*((u^2 (d1^2/dyf1))/(u^4 + (2*Pi*f*(d1^2/dyf1))^2)) + 288/d1^3 (u^2/(81 + 9 u^2 - 2 u^4 + u^6))*((u^2 (d1^2/dyf1))/(u^4 + (4*Pi*f*(d1^2/dyf1))^2)), {u, 0, \[Infinity]}])

to experimentral data: 

data={{19998000, 13.068}, {17123000, 13.62}, {14666000, 14.265}, 
{12556000,14.887}, {10756000, 15.245}, {9206200, 16.122}, {7883900, 
17.043}, {7883900, 16.644}, {6750700, 16.843}, {5781800, 
17.551}, {4951600, 17.722}, {4241300, 18.256}, {3631400, 
18.826}, {3108000, 19.284}, {2661600, 19.299}, {2278600, 
19.923}, {1952000, 20.209}, {1671100, 20.375}, {1431100, 
20.604}, {1225500, 21.004}, {1048900, 21.313}, {899510, 
21.358}, {769990, 21.644}, {659570, 22.018}, {564850, 
22.265}, {483000, 22.34}, {413680, 22.488}, {354450, 
22.65}, {303300, 22.701}, {260180, 22.772}, {222610, 
22.985}, {190750, 23.041}, {163280, 23.228}, {139720, 
23.202}, {119640, 23.453}, {102650, 23.539}, {87900, 
23.537}, {75165, 23.584}, {64234, 23.753}, {54904, 23.607}, {47340, 
23.842}, {40338, 24.074}, {34484, 23.853}, {29634, 23.833}, {25406, 
23.931}, {21739, 23.896}, {18573, 23.842}, {15857, 24.066}, {13725, 
24.008}, {11746, 23.97}, {9924.4, 24.282}}

to estimate the values of the parameter d1 and dyf1. For this purpose, I've been trying to use NonlinearModelFit:

nlm = NonlinearModelFit[data, modelx[d1, dyf1, f] , {{d1, 10^-10}, {dyf1, 10^-10}}, f, Method -> "LevenbergMarquardt", MaxIterations -> 300]

The problem: I cannot get a good fit to the experimental data. It's not even close fit - from the parameter table it may be seen, that standard errors of estimated values are grater that the estimates. Also, visually the fit seems to be poor in the log-log scale:

Show[ListLogLogPlot[data], LogLogPlot[nlm[f], {f, 5000, 50000000}]]

I know that the function may be well fitted to experimental data; the "good" values for parameters are:

d1=2.33*10^-10; dyf1=1.08*10^-11;

As you may see, in the nlm formula I used rather nice initial conditions for the fitting procedure, but no luck. Any help appreciated!

Edit

After suggestions (thanks Jim and Michael!) I tried new code to find an impact of factor "6.62*10^-21" on values of d1 and dyf1, that is

model = Integrate[#, {u, 0, \[Infinity]}, 
Assumptions -> a > 0 && t > 0 && f > 0] & /@
Apart[a*(u^2/(81 + 9 u^2 - 2 u^4 + 
u^6))*((u^2 )/(u^4 + (2*Pi*f*(t))^2)) + 
4 a (u^2/(81 + 9 u^2 - 2 u^4 + 
u^6))*((u^2 )/(u^4 + (4*Pi*f*(t))^2))];

where a==72(6.62*10^-21*b/d1^3) and b==d1^2/dyf1. Now, when I'm trying to run NonlinearModelFit, I get "The function value is not a list of real numbers [...]". Any help appreciated!

Improve this question
$\endgroup$
7
  • 3
    $\begingroup$ Have you tried rescaling the units of your data and your model? $\endgroup$
    – J. M.'s eventual burnout ♦
    Feb 17, 2016 at 17:41
  • 1
    $\begingroup$ When I run this with Mathematica 10.2 (Windows 7) I get "The integrand...has evaluated to non-numerical values for all sampling points in the region with boundaries {{[Infinity],0.}}." That would suggest that there might be some issues with the function modelx. Did you get such errors/warnings? $\endgroup$
    – JimB
    Feb 17, 2016 at 19:52
  • $\begingroup$ J.M. - I have not tried rescaling... does Mathematica work better in (0,1) regime? Jim Baldwin - yes, I get such warnings, but also the values of the parameter are produced. Maybe something is wrong with modelx, but for now I dont see what it might be. $\endgroup$
    – user19388
    Feb 17, 2016 at 20:16
  • $\begingroup$ Is this a statistical model of some kind? It's a monster. Have you considered using a simpler model? Linearizing the model to some order around the values of d1 and dyf1 you expect? $\endgroup$
    – Searke
    Feb 17, 2016 at 21:17
  • 1
    $\begingroup$ @Searke. For whatever it's worth the same situation holds for the great increase in the use of linear mixed models (not to mention nonlinear mixed models) where there are multiple levels and types of random effects. $\endgroup$
    – JimB
    Feb 18, 2016 at 20:03

2 Answers 2

Reset to default
6
$\begingroup$

Following @J.M. 's suggestion to scale the data I also reparameterized the model:

(* Reparameterize model *)
(* b == d1^2/dyf1 *)
(* a == 6.62*10^-21*12*b/d1^3 *)

modelx[a_?NumericQ, b_?NumericQ, f_?NumericQ] := 
(NIntegrate[ ((6 a u^4)/( (4 b^2 f^2 \[Pi]^2 + u^4) (81 + 9 u^2 - 2 u^4 + u^6)) + 
    (19 a u^4)/( (16 b^2 f^2 \[Pi]^2 + u^4) (81 + 9 u^2 - 2 u^4 + u^6))), {u, 0, ∞}])

(* Scale the predictor values *)
data2 = data;
data2[[All, 1]] = data2[[All, 1]]/1000000;

(* Fit model *)
nlm = NonlinearModelFit[data2, modelx[a, b, f], {{a, 36}, {b, 0.0045}}, f, 
  Method -> "LevenbergMarquardt", MaxIterations -> 300]

(* Plot results *)
Show[ListLogLogPlot[data], LogLogPlot[nlm[f/1000000], {f, 5000, 50000000}]]

(* Convert fitted parameters to original parameters *)
{d1 -> (4.2987918004513846`*^-7 b^(1/3))/a^(1/3), 
  dyf1 -> 1.847961094362806`*^-13/(a^(2/3) b^(1/3))} /.nlm["BestFitParameters"]
(* {d1 -> 2.1174835472049714`*^-8,dyf1 -> 9.980071510204195`*^-14} *)

Curve fit

Update

Below is a modified version that temporarily reparameterizes the model but also follows @MichaelE2 's use of Integrate to obtain estimates of the original parameters along with their estimated standard errors.

data = {{19998000, 13.068}, {17123000, 13.62}, {14666000, 
    14.265}, {12556000, 14.887}, {10756000, 15.245}, {9206200, 
    16.122}, {7883900, 17.043}, {7883900, 16.644}, {6750700, 
    16.843}, {5781800, 17.551}, {4951600, 17.722}, {4241300, 
    18.256}, {3631400, 18.826}, {3108000, 19.284}, {2661600, 
    19.299}, {2278600, 19.923}, {1952000, 20.209}, {1671100, 
    20.375}, {1431100, 20.604}, {1225500, 21.004}, {1048900, 
    21.313}, {899510, 21.358}, {769990, 21.644}, {659570, 
    22.018}, {564850, 22.265}, {483000, 22.34}, {413680, 
    22.488}, {354450, 22.65}, {303300, 22.701}, {260180, 
    22.772}, {222610, 22.985}, {190750, 23.041}, {163280, 
    23.228}, {139720, 23.202}, {119640, 23.453}, {102650, 
    23.539}, {87900, 23.537}, {75165, 23.584}, {64234, 
    23.753}, {54904, 23.607}, {47340, 23.842}, {40338, 
    24.074}, {34484, 23.853}, {29634, 23.833}, {25406, 
    23.931}, {21739, 23.896}, {18573, 23.842}, {15857, 
    24.066}, {13725, 24.008}, {11746, 23.97}, {9924.4, 24.282}};
(* Scale down the predictor variable *)
data2 = data;
data2[[All, 1]] = data2[[All, 1]]/1000000;

(* Convert d1 and dyf1 such that a = d1^2/dyf1 and b = 6.62*10^(-21)a/d1^3 = 6.62*10^(-21)/(d1 dyf1) *)
model = Integrate[#, {u, 0, ∞},
     Assumptions -> b > 0 && a > 0 && f > 0] & /@
   Apart[(72 b) (u^2/(81 + 9 u^2 - 2 u^4 + u^6))*((u^2 )/(u^4 + (2*Pi*f*a)^2)) +
     (288 b) (u^2/(81 + 9 u^2 - 2 u^4 + u^6))*((u^2 )/(u^4 + (4*Pi*f*a)^2))];

(* Estimate transformed parameters *)
nlm = NonlinearModelFit[data2, model, {{b, 1}, {a, 4.4 10^(-9)}}, f,
   Method -> "LevenbergMarquardt", MaxIterations -> 300, 
   WorkingPrecision -> $MachinePrecision];

(* Plot results *)
Show[ListLogLogPlot[data], 
 LogLogPlot[nlm[f/1000000], {f, 5000, 50000000}]]

(* Convert back to original parameters *)
convert = Solve[{b == 6.62*10^(-21)/(d1 dyf1), a == d1^2/dyf1}, {d1, dyf1}][[3]]
(* {d1->a^(1/3)/b^(1/3),dyf1->1/(a^(1/3) b^(2/3))} *)
sol = nlm["BestFitParameters"]
(* {b->2.61041953398039418124946423422022403882,a->0.00437789081254762809545490882186499512} *)
{d1, dyf1} /. convert /. sol
(* {2.2308525602769456*10^(-8),1.1367810114017232*10^(-13)} *)

(* Construct covariance matrix for the estimators of d1 and dyf1 
   using the Delta method *)
cov = nlm["CovarianceMatrix"]
(* {{0.00001818912969975069104554318298941637,1.
923673480190859470043920077963200094427339451309835*10^-7},
{1.923673480190859470043920077963200095293828391530791*10^(-7),
5.1499953674006162567008452940009580161083650398944*10^(-9)}} *)
g = Grad[{d1, dyf1} /. convert, {b, a}];
var = ((g.cov.Transpose[g]) /. sol)

(* Standard error for estimator of d1 *)
var[[1, 1]]^0.5
(* 0.0006105953679914561 *)

(* Standard error for estimator of dyf1 *)
var[[2, 2]]^0.5
(* 0.020012556571506964 *)

(* Correlation of estimators *)
var[[1, 2]]/(var[[1, 1]] var[[2, 2]])^0.5
(* -0.9760212710024363 *)

Note that the correlation of the estimators of the original coefficients is extremely high. Such high correlations can also make it difficult obtain convergence.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Jim, can you check estimates and standard errors nlm["ParameterTable"] in your nb? It seems that SE are still greater than estimates (at least i got such results). Plot looks nice, but the values of the parameters may be akward. $\endgroup$
    – user19388
    Feb 18, 2016 at 8:57
  • 1
    $\begingroup$ Yes, I didn't show how to convert the standard errors for the re-parameterized model to that of the original model. I'll do that shortly. $\endgroup$
    – JimB
    Feb 18, 2016 at 15:50
  • 1
    $\begingroup$ I'm sad that I cannot give a second upvote. Nicely done! $\endgroup$
    – J. M.'s eventual burnout ♦
    Feb 19, 2016 at 7:04
5
$\begingroup$

You could calculate with arbitrary precision, using the exact integral instead of a numeric one. The integration takes a few seconds, but it's not too long.

model = Integrate[#, {u, 0, ∞}, 
     Assumptions -> d1 > 0 && dyf1 > 0 && f > 0] & /@ 
   Apart[72/
       d1^3 (u^2/(81 + 9 u^2 - 2 u^4 + 
          u^6))*((u^2 (d1^2/dyf1))/(u^4 + (2*Pi*f*(d1^2/dyf1))^2)) + 
     288/d1^3 (u^2/(81 + 9 u^2 - 2 u^4 + 
          u^6))*((u^2 (d1^2/dyf1))/(u^4 + (4*Pi*f*(d1^2/dyf1))^2))];

nlm = NonlinearModelFit[data, model, {{d1, 10^-10}, {dyf1, 10^-10}}, 
  f, Method -> "LevenbergMarquardt", MaxIterations -> 300, 
  WorkingPrecision -> $MachinePrecision]

Show[ListLogLogPlot[data], LogLogPlot[nlm[f], {f, 5000, 50000000}]]

Mathematica graphics

There's a warning about the precision of data, but it's unimportant in this case.

Improve this answer
$\endgroup$
4
  • 2
    $\begingroup$ It seems the lesson for me is that all of the standardizing and reparameterization of models that I've done for years in SAS and R to get convergence can be replaced by learning how to take advantage of the precision options in Mathematica. $\endgroup$
    – JimB
    Feb 18, 2016 at 5:57
  • 2
    $\begingroup$ @JimBaldwin From my experience, it is not always possible to cope vastly different scales in $x$ and $y$ by increasing precision. I strongly suspect that at least built-in "LevenbergMarquardt" has hard-coded minimal step size (in the sense of options "StartingScaledStepSize" and "MaxScaledStepSize") equals approximately to 10^-10 which becomes insufficient when the scales differ too much (for example, magnitude 10^12 or more, depends on the number of variables). $\endgroup$
    – Alexey Popkov
    Feb 18, 2016 at 8:00
  • 2
    $\begingroup$ @JimBaldwin Yes, what Alexey says seems right. OTOH, the precision noise that sometimes occurs in using NIntegrate, or the round-off in plain machine precision, can be a significant factor, too. Scaling is a standard strategy. Increasing working precision is another, in Mathematica at least; it could be looked at as a lazy way, in that it's easy to try without having to think too hard. $\endgroup$
    – Michael E2
    Feb 18, 2016 at 11:42
  • $\begingroup$ @Michael E2 - thanks for your help! With your suggestion I got a good value for ratio d1^2 / dyf (~10^-9, physicists are smilling). Now I'm trying to get the exact values of d1 and dyf. I have noticed, that adding a constant before the integral as small as 10^-21 makes some new problems in notebook, so I decided to get the best of two worlds and mix yours and Jim Baldwin's (reparametrizing only) ideas. For now it ended with "The function value is not a list of real numbers with dimensions [...]"... Looks like another sleepless night $\endgroup$
    – user19388
    Feb 18, 2016 at 18:42

Your Answer

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

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