# Low quality of estimated parameters with NonlinearModelFit

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!

• Have you tried rescaling the units of your data and your model? – J. M.'s torpor Feb 17 '16 at 17:41
• 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? – JimB Feb 17 '16 at 19:52
• 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. – user19388 Feb 17 '16 at 20:16
• 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? – Searke Feb 17 '16 at 21:17
• @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. – JimB Feb 18 '16 at 20:03

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}]]  There's a warning about the precision of data, but it's unimportant in this case. • 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. – JimB Feb 18 '16 at 5:57 • @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). – Alexey Popkov Feb 18 '16 at 8:00 • @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. – Michael E2 Feb 18 '16 at 11:42 • @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 – user19388 Feb 18 '16 at 18:42 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} *)  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.

• 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. – user19388 Feb 18 '16 at 8:57
• 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. – JimB Feb 18 '16 at 15:50
• I'm sad that I cannot give a second upvote. Nicely done! – J. M.'s torpor Feb 19 '16 at 7:04