polynomial curve fitting without regression (power law)

298 Views Asked by At

I have a set of points which seem to follow the curve of a polynomial function of the form $x\mapsto x^a$ where $a$ is well chosen.

Is there an easier one to compute $a$ than using a polynomial regression ?

2

There are 2 best solutions below

4
On BEST ANSWER

@caverac gave the answer but I would like to add a point.

The model being $y=x^a$ and using $n$ data points $(x_i,y_i)$ making the transform $$\log(y)=a\log(x)\to v=a w$$ and using linear regression means that you minimize $$SSQ_1=\sum_{i=1}^n(a \log(x_i)-\log(y_i))^2$$ The problem is that what is measured is $y$ and not any of its possible transforms. This means that, after this first step which will provide a good estimate of $a$, you should use nonlinear regression in order to minimize $$SSQ_2=\sum_{i=1}^n(x_i^a-y_i)^2$$ writing the derivative with respect to $a$ and setting it equal to $0$ then leads to $$\sum_{i=1}^n x_i^{2a}\log(x_i)-\sum_{i=1}^n x_i^{a}y_i\log(x_i)=0$$ which can be solved using Newton method for example. If you use it, the iterates would be given by $$a_{n+1}=a_n-\frac{\sum _{i=1}^n x_i^{a_n} \log (x_i) \left(x_i^{a_n}-y_i\right)}{\sum _{i=1}^n x_i^{a_n}\log^2(x_i) \left(2 x_i^{a_n}-y_i\right)}$$ which can easily be done using Excel; for sure the iterations must start with $a_0$ given by the first step.

Depending on the level of noise, the parameters could be quite different and the extrapolation could be quite different too.

If the errors are small, you could easily show that stopping at the first step is almost equivalent to the minimization of the sum of squares of relative errors.

For illustration purposes, let us use the following data $$\left( \begin{array}{cc} x_i & y_i \\ 10 & 20000 \\ 11 & 40000 \\ 12 & 80000 \\ 13 & 120000 \\ 14 & 160000 \\ 15 & 220000 \end{array} \right)$$ The first step will lead to $a=4.49304$. Starting Newton iterations $$\left( \begin{array}{cc} n & a_n \\ 0 & 4.49304 \\ 1 & 4.55058 \\ 2 & 4.54037\\ 3 & 4.53994 \\ 4 & 4.53993 \end{array} \right)$$ Now, let us suppose that you want to extrapolate the model to $x=20$; the first one would lead to $y=700779$ while the second one would give $y=806476$ which is quite different ($15$% relative difference) even if the parameter just changed by $1$%.

If you have to use the model for interpolation instead, it would be a good idea to minimize the sum of squares of relative errors and a very similar process would lead to $a=4.45794$.

Please laso notice that, using this trick, you would have something intermediate between the two models with $$a=\frac{\sum _{i=1}^n y_i \log (x_i) \log (y_i)}{\sum _{i=1}^n y_i \log^2 (x_i)}$$ Applied to the data used for illustration, this would lead to $a=4.53282$ which is quite close to the value obtaind by the rigourous procedure.

0
On

Linear regression, if $y = x^a$ then $\ln y = a \ln x$, call $v=\ln y$ and $u = \ln x$, then you have a set of points of the form $(u, v)$ that will follow a linear relation with slope $a$ and intercept $0$.

Here is a way of finding $a$: The residual sum of squares for the fit is

$$ E(a) = \frac{1}{2}\sum_i(v_i - au_i)^2 $$

Taking the derivative with respect to $a$ we get

$$ \frac{dE}{da} =-\sum_i u_i(v_i - a u_i) $$

Because we want to find the minimum of $E$ we set this to zero and find

$$ a = \frac{\sum_i u_iv_i}{\sum_i u_i^2} $$