1
$\begingroup$

I am working with a fairly large dataset, and at one point I need it to be fit to a complicated model (along the lines of 

y00 + y01 y + y01p5 y^(3/2) + ay02 y^2 + ay02p5 y^(5/2) + ay03 y^3 + ay03p5 
y^(7/2) + ay04 y^4 + ay04p5 y^(9/2) + ay05 y^5 + ay05p5 y^(11/2) + ay06 y^6 +
ay06p5 y^(13/2) + ay07 y^7 + ay07p5 y^(15/2) + ay03L y^3 Log[y] + ay04L y^4 
Log[y] + ay04p5L y^(9/2) Log[y] + ay05L y^5 Log[y] + ay05p5L y^(11/2) Log[y] + 
ay06L y^6 Log[y] + ay06p5L y^(13/2) Log[y] + ay07L y^7 Log[y] + ay07p5L  
y^(15/2) Log[y] + ay08L y^8 Log[y] + ay06L2 y^6 Log[y]^2 + ay07L2 y^7 Log[y]^2  
+ ay07p5L2 y^(15/2) Log[y]^2 + ay08L2 y^8 Log[y]^2 + ay08p5L2 y^(17/2) 
Log[y]^2 + ay09L2 y^9 Log[y]^2 + ay09p5L2 y^(19/2) Log[y]^2

etc. for more terms). Anyway, for a while, we had been using 

NonlinearModelFit[{yvalues, modelvalues}, form, {parameterlist}, y]

with great success. It gave accurate results without more than .001 "standard error" for any term. However, I recently tried to go back and do the same thing with 

LinearModelFit[{yvalues, modelvalues}, {yterms (functions of y) in the model}, y]

thinking it might even improve the error and/or speed (since the model is linear). Instead, this route gave way more error -- up to 10^10 in some cases.

The lower order terms match up pretty well between the two (though again, NLM is more accurate, based on known results), but at the highest orders they diverge wildly.

Does anyone know what's going on?

Edit: To give a better sense of the specifics, here is the (precision-truncated) LMFit output model:

1.0000 - 3.7113 y + 12.566 y^(3/2) - 4.9285 y^2 - 38.293 y^(5/2) + 115.73 y^3 
- 101.51 y^(7/2) - 117.50 y^4 + 719.13 y^(9/2) - 1216.9 y^5 + 958.93 y^(11/2) 
+ 2034.8 y^6 - 7782.0 y^(13/2) + 15384. y^7 - 12116. y^(15/2) - 1.9534*10^7 
y^8 + 2.7183*10^11 y^(17/2) - 8.1524 y^3 Log[y] + 26.372 y^4 Log[y] - 102.45 
y^(9/2) Log[y] + 58.320 y^5 Log[y] + 236.81 y^(11/2) Log[y] - 950.36 y^6 
Log[y] + 1080.4 y^(13/2) Log[y] + 133.71 y^7 Log[y] - 5213.2 y^(15/2) Log[y] - 
1.5259*10^6 y^8 Log[y] + 3.1658*10^10 y^(17/2) Log[y] + 33.231 y^6 Log[y]^2 - 
85.076 y^7 Log[y]^2 + 417.63 y^(15/2) Log[y]^2 - 31961. y^8 Log[y]^2 + 
9.7743*10^8 y^(17/2) Log[y]^2

Here is the NLMFit model:

1.0000 - 3.7113 y + 12.566 y^(3/2) - 4.9285 y^2 - 38.293 y^(5/2) + 115.73 y^3 
- 101.51 y^(7/2) - 117.50 y^4 + 719.13 y^(9/2) - 1216.9 y^5 + 958.93 y^(11/2) 
+ 2034.8 y^6 - 7782.0 y^(13/2) + 15384. y^7 - 12154. y^(15/2) - 14424. y^8 + 
93631. y^(17/2) - 8.1524 y^3 Log[y] + 26.372 y^4 Log[y] - 102.45 y^(9/2) 
Log[y] + 58.320 y^5 Log[y] + 236.81 y^(11/2) Log[y] - 950.36 y^6 Log[y] + 
1080.4 y^(13/2) Log[y] + 133.71 y^7 Log[y] - 5215.7 y^(15/2) Log[y] + 10899.  
y^8 Log[y] - 14225. y^(17/2) Log[y] + 33.231 y^6 Log[y]^2 - 85.076 y^7 
Log[y]^2 + 417.59 y^(15/2) Log[y]^2 - 314.80 y^8 Log[y]^2 - 530.32 y^(17/2) 
Log[y]^2

Note the wild disparity in the terms with y^8 or higher. Most importantly, the NLM is overwhelmingly accurate for these higher terms, while the LM is not.

Improve this question
$\endgroup$
9
  • 3
    $\begingroup$ Fitting high order polynomials (even as low an order as 4 or 5) can be problematic to fit because of the high correlation among the coefficient estimators and the loss of significant digits when raising a number to a large power. General advice would be (1) Consider a different basis than a polynomial, (2) standardize both the dependent and independent variable values, and (3) increase the working precision. Specific advice would require specifics from you such as some data and a specific model. $\endgroup$
    – JimB
    Commented Aug 2, 2016 at 18:13
  • $\begingroup$ Giving the estimated coefficients helps (unless what you've given is just a truncated model) but you really need to give more specifics. What is the range for the values of y and what's the approximate sample size? What's the mean square error from the NonlinearModelFit? I'd prefer more information but at least with those items one could simulate some data to see what might be going on. $\endgroup$
    – JimB
    Commented Aug 2, 2016 at 20:11
  • $\begingroup$ The y values are all of order 10^-10 through 10^-35, and the sample size is 33. I'll come back and give some more info when I get a chance. $\endgroup$
    – cmm0052
    Commented Aug 2, 2016 at 20:38
  • 2
    $\begingroup$ 33 ? With 31 parameters to fit (plus an error variance) ? But your first sentence says you have a fairly large dataset. What am I missing? $\endgroup$
    – JimB
    Commented Aug 2, 2016 at 20:55
  • 2
    $\begingroup$ I think a lot of folks here would like to solve this issue for you. But I'm not seeing enough information given (with responses bordering on excuses) to help us help you. @Feyre's suggestion about scaling fits perfectly (no pun intended) with your equation. Log[y] = Log[y*10^15] - Log[10^15], for example. Mathematica can take care of the extra (but known) terms that might be created. And "33 values of y are coupled to 51 values of b": When did b come to the party? $\endgroup$
    – JimB
    Commented Aug 3, 2016 at 2:21

1 Answer 1

Reset to default
2
$\begingroup$

This is an extended comment. When plotting the "NLMFit" model using the range of y values you give I get the following:

f[y_] := 1.0000 - 3.7113 y + 12.566 y^(3/2) - 4.9285 y^2 - 38.293 y^(5/2) + 115.73 y^3 -
  101.51 y^(7/2) - 117.50 y^4 + 719.13 y^(9/2) - 1216.9 y^5 + 958.93 y^(11/2) +
  2034.8 y^6 - 7782.0 y^(13/2) + 15384. y^7 - 12154. y^(15/2) - 14424. y^8 +
  93631. y^(17/2) - 8.1524 y^3 Log[y] + 26.372 y^4 Log[y] - 102.45 y^(9/2) Log[y] +
  58.320 y^5 Log[y] + 236.81 y^(11/2) Log[y] - 950.36 y^6 Log[y] +
  1080.4 y^(13/2) Log[y] + 133.71 y^7 Log[y] - 5215.7 y^(15/2) Log[y] +
  10899. y^8 Log[y] - 14225. y^(17/2) Log[y] + 33.231 y^6 Log[y]^2 -
  85.076 y^7 Log[y]^2 + 417.59 y^(15/2) Log[y]^2 - 314.80 y^8 Log[y]^2 -
  530.32 y^(17/2) Log[y]^2

Plot[f[y], {y, 10^-35, 10^-10}]

Fit with default scaling

If I include zero in the vertical axis, I get

Plot[f[y], {y, 10^-35, 10^-10}, PlotRange -> {All, {0, 1}}]

Fit with zero in vertical axis

It appears that in this case, none of the terms higher order than 1 have any meaning. If so, then it doesn't really matter that LinearModelFit and NonlinearModelFit give different coefficients (at least with respect to predictions). But if you are then using the estimated coefficients as if they were "data" in another fitting procedure, I'd would strongly discourage that and do all of the fitting in a single model - or just use the interpolation functions available in Mathematica.

Again, the general advice given above by both me and @Feyre is just that: general advice. For your particular issue I think you need specific advice based on your dataset. That means supplying the model to be fit along with the data. If it's confidential data, then make it impossible to obtain the exact data by adding some constant only known to you to the dependent variable (the same constant to all values). But we need a complete example.

Improve this answer
$\endgroup$

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.