2
$\begingroup$

I am fitting a power law $Y = aX^b$ to a data set. There are two ways to do this:

  • One is on log-scale, which means fitting $\log(Y)$ to $a_1 + b_1\log(X)$.
  • In other case, using linear scale to fit $Y = a_2X^{b_2}$.

Mathematically, both should be same, with $a_1 = \log (a_2)$ and $b_1 = b_2$. However, I am getting different fitting parameters by two methods (first looking visually poor!).

Here is a working example:

Y = {{1., 4.69246}, {2., 61.9172}, {3., 428.652}, {4., 2687.97}, {5., 10869.6}};

nlm1 = NonlinearModelFit[Log[Y], a1 + b1 x, {a1, b1}, x];
nlm1["ParameterTable"]
Show[ListPlot[Y], Plot[Exp[nlm1[Log[x]]], {x, 0.5, 6}]]

nlm2 = NonlinearModelFit[Y, a2*x^b2, {a2, b2}, x];
nlm2["ParameterTable"]
Show[ListPlot[Y], Plot[nlm2[x], {x, 0.5, 6}]]
Improve this question
$\endgroup$
1

2 Answers 2

Reset to default
11
$\begingroup$

The difference is in the way you treat the residuals. The real model is usually:

Y == a X^b + error

where error follows some distribution (usually assumed to be a normal distribution).

When you log both sides, you get:

Log[Y] == Log[a X^b + error]

Note how you can't actually simplify the log on the right any more because of the + under the Log. (edit note also how there is no guarantee that a X^b + error is actually positive, so it's not even clear you're allowed to take the log at all!)

However, NonLinearModelFit implicitly assumes a model of the form y == f[x] + error, so the two formulations aren't equivalent.

Of course, you could start from the other side and assume that the real model is:

Log[Y] == Log[a X^b] + error == Log[a] + b Log[X] + error

This formulation is often very useful or even preferable over the previous one (many physical phenomena play out on log scales rather than linear scales), but then the corresponding model in normal space would be:

Y == Exp[error] * (a X^b)

Again: when you assume an additive error model one way, you will end up with a non-additive error model on the other side. That's why these two models don't fit in the same way.

My rule of thumb is: if you can find a simple transformation (log-log, log-linear, linear-log etc.) that makes your data linear, you should probably go with that. There's usually a good justification for it and it's a lot easier to judge the quality of your fit when it's just a straight line.

Improve this answer
$\endgroup$
6
  • 2
    $\begingroup$ +1 for explaining different ways that error might be incorporated. It is my biased option that some folks at this site forget (and even seem unaware) that the error structure is an essential component of a "model". The next step (which I know practice) is the examination of residuals for any deviations from a proposed error structure. But with so few observations here, only very gross departures can be observed. $\endgroup$
    – JimB
    Aug 11 at 14:54
  • $\begingroup$ Log[a X^b] expands to Log[a] + b Log[X] rather than Log[a] * b Log[X], unless I missed something. Your overall point still stands of course. $\endgroup$
    – Mad Physicist
    Aug 12 at 3:02
  • $\begingroup$ thanks for an elaborative explanation. I get your point. However, the data I provided, gives a better fit for ax^b, rather than a + bLog(x) (which makes data linear !!) $\endgroup$
    – user49535
    Aug 12 at 12:23
  • $\begingroup$ @MadPhysicist Yes, thanks for catching the typo. $\endgroup$
    – Sjoerd Smit
    Aug 14 at 9:43
  • $\begingroup$ @user49535 I don't know if I agree with you. The fit LinearModelFit[Log[Y], {1, x}, x] is not perfect, but neither is NonlinearModelFit[Y, a x^b, {a, b}, x], which has significantly greater residuals for the 2nd and 3rd data points than the others. $\endgroup$
    – Sjoerd Smit
    Aug 14 at 9:48
2
$\begingroup$

Log[Y! will not only logarithmize the function value, but also the argument.

If you fit the logarithms, the errors are equally distributed. However, if you then take the exponent, the first points will have small errors and the last points large errors. This is clearly visible in:

Y = {{1., 4.69246}, {2., 61.9172}, {3., 428.652}, {4., 2687.97}, {5., 
    10869.6}};
Ylog = {#[[1]], Log[#[[2]]]} & /@ Y;
nlm1 = NonlinearModelFit[Ylog, a1 + b1 x, {a1, b1}, x];
nlm1["ParameterTable"]
Show[ListPlot[Y], Plot[Exp[nlm1[x]], {x, 0.5, 6}]]

enter image description here

Improve this answer
$\endgroup$

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.