2
$\begingroup$

I am struggling finding pseudo Code for a non-linear fit of the following function:

$y = a\, x^b$

Package NLS in R does perform well, but utilizing external software is not practicable in my program anymore. I would like to avoid this by implementing a short class by myself. Literature links including the solution are also welcome.

Update: I already have implemented the log transform with linear square fit. That is not the solution I am looking for. I need to implement NLS fit.

Update 2, giving reason why the log transformed method is not sufficient:

My original data:

Allometry function related to forestry

The log transformed data with a linear model line in red, one can see that the log transformed residuals are a heavy influence:

And its log transformed counterpart, with fitted line in red

The original data again, the red line is the power function with parameters derived from the log transformed method, the blue line is the power function predicted with package NLS2 in R:

The original data with power function derived from the log transform method in red and the better function derived with R-package NLS2 in blue

Clearly the log-transformed residuals have too much influence, the NLS method forces the prediction on the stable part of the model.

$\endgroup$
  • 3
    $\begingroup$ Does this help? mathworld.wolfram.com/LeastSquaresFittingPowerLaw.html $\endgroup$ – Christian Clason Apr 26 '16 at 16:08
  • 2
    $\begingroup$ One can simply do this by first performing the log of your function, resulting in: $log(y) = log(a) + b*log(x)$. Using this representation, the problem becomes a Linear Least Square to find $log(a)$ and $b$. Obviously then you can obtain $a$ alone. $\endgroup$ – spektr Apr 26 '16 at 18:50
  • $\begingroup$ @choward yes I am absolutely aware of this technique and have it already implemented. But I need a NLS fit as I do not want to minimize log transformed residuals... $\endgroup$ – Jan Hackenberg Apr 26 '16 at 22:24
  • $\begingroup$ Question: what in your mind is the problem with minimizing the log transformed residuals? Do you think some finite precision issue by computing a bunch of logs will make the result less accurate? $\endgroup$ – spektr Apr 29 '16 at 15:38
  • 1
    $\begingroup$ @choward That is exactly the point, added three R plots to my question which should clearify. Gauss Newton prediction seems to solve efficiently the problem in my C++ code. $\endgroup$ – Jan Hackenberg Apr 29 '16 at 17:44
5
$\begingroup$

Chapter 10 of Numerical Analysis (Richard L. Burden, J. Douglas Faires) gives good readable pseudo code for Newtons method. The start parameters are taken from the solution of the linear problem of the log transformed data. (see comment of Christian Clason giving a link (Linear Solution, Mathworld Wolfram) )

while not reached convergence criteria
      calucate vector of residuals F(x) and Jacobian matrix J(x)
      solve J(x)y = -F(x)
      set x = x + y

Here $x = a,b$

$-F(x)$ is for $n$ points a vector with length $n$:

the $ith$ element is $- a\, x_i^b + y_i$

and the Jacobian $J(x)$ is a $n,2$ Matrix with the $ith$ row:

$J_{i,1} = x_i^b$

and

$J_{i,2} = a\,log(x_i)\,x_i^b $

$\endgroup$

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.