Power fit to some experimental data

I have to fit some data to a power law

$$F=\alpha q^{\beta}$$

being $q$ and $F$ the experimental data points. What I usually do is taking logs so that

$$\ln(F) = \beta \ln(q)+\ln(\alpha)$$

and apply least squares with uncertainties in both variables. How may I approach this in the case that some of the values of $q$ (and eventually $F$) are negative?. If it helps $\beta$ should be compatible with a value of $1$, but I want a general approach.

Edit:

The law to be fitted is Coulomb's law, $F$ is the force between to charged spheres. One of them has fixed charge, the other is a variable we call $q$. The proportionality constant $\alpha$ is related to $\varepsilon_0$ with a known expression in powers of $\beta=\frac{R}{d}$ where $d$ is the distance between the centers of the spheres and $R$ is their radius. So, I want to determine, from experiment, both $\alpha$ and $\beta$ and compare them with $\varepsilon_0$ and $1$.

The actual values are

$$\begin{array}{|c|c|} \hline q / \mathrm{nC} & F / \mathrm{mN} \\ \hline -26,8 \pm 0.8 & -1.5 \\ \hline -18.8 \pm 0.5 & -1.1 \\ \hline -9.9 \pm 0.2 & -0.6 \\ \hline 9.9 \pm 0.2 & 0.5 \\ \hline 18.8 \pm 0.5 & 1.0 \\ \hline 26.8 \pm 0.8 & 1.5 \\ \hline 35 \pm 2 & 2.1 \\ \hline \end{array}$$

with the uncertainty in the last decimal place for the forces.

-
I encourage you to seriously consider the answers besides the accepted one. Yes, taking logs is sometimes inappropriate, but sometimes it is the right thing to do. Moreover, your problems have nothing to do with normal vs. log-normal uncertainties; your problem is that you posit $F$ is a power law in $q$, and your data emphatically tells you only a complex number can be a power of $q$, and really you should have $\lvert F\rvert \propto \lvert q\rvert^\beta$ as your model. – Chris White Jun 3 '14 at 22:05
@ChrisWhite thanks. I did considered them. The absolute value model however assumes the force to be charge independent, which might obscure the case in which the experimental setup it is not that way. – Jorge Jun 4 '14 at 13:07

Being a mathematician, I feel that I should point out a common misconception here. Don't feel too bad, plenty of mathematicians (including myself) have fallen into this trap. Basically, if you want to fit data to power law using least squares methods, then you should NOT fit a straight line in log log space. You should fit a power law using non-linear least squares to the original data in linear space without any kind of transformation.

The basic intuition behind this is this. Least squares method assume that the error distribution is normal meaning that there is one true (unknown) value which you are trying to measure and the measurements you are taking are above/below the true value staying close to the true value so the distribution is gaussian. The problem is that log is not a linear function so when you take log log of the data, numbers bigger than one are pushed together and numbers less than one are spread apart and namely, the error distribution is not normal any more in log log space....meaning that least squares is not guaranteed to converge to the true values any more. In fact it has been proven that fitting in log-log space consistently gives you biased results.

Fitting a straight line in log-log was fashionable back in the day before GHz processors. Now whatever you are using to do the computation, most likely has the ability to do non-linear least squares power law fit to the original data so that is the one you should do. Since power-law is so prevalent in science, there are many packages and techniques for doing them efficiently, correctly, and fast.

I refer you to these published results and a rant here. This is the my question which sparked this discussion.

So by OP's request here is an actual example. The actual function is $y=f(x)=2x^{-4}=ax^b$. I take a bunch of points from $x\in[1,2]$ and $x=10$. Note that $f(10)=0.0002$ and I change this $y$-value to $y=0.00002$ while keeping all others the same and then fitting them using least squares to estimate $a$ and $b$.

Using non-linear least squares fit to the functional form $y=ax^b$ gives us

\begin{eqnarray} a&=&2\\ b&=&-4\\ r^2&=&0.9999. \end{eqnarray}

Using linear least squares fit to the functional form $y=mx+n$ in log-log space gives us

\begin{eqnarray} a &=& exp(n) = 2.46\\ b &=& m = -4.57\\ r^2 &=& 0.9831. \end{eqnarray}

The $95\%$ bounds for $b$ are $(-4.689, -4.451)$ which doesn't even include the true value for $b$. Basically log pushes numbers bigger than one together and spreads apart numbers less than one. So in linear space changing $y=0.0001$ to $y=0.00001$ is inconsequential. The algorithm doesn't care too much about it. The change in the residual is very tiny. In log-log space the difference is rather huge (an entire order of magnitude) and the point has become an "outlier" from the nice linear trend. The change in residual is now large. Since the algorithm is trying to minimize the sum of the squares, the algorithm can't ignore the deviation and must try to accommodate the outlier point so the line is skewed messing up the slope and the intercept.

Here is the MATLAB code and the graphs. You can easily reproduce this and play with it yourself. Blue are the data points. Black is the power fit in linear space. Red is the fitted line from log-log least squares fitting. The top panel shows all three in linear space and the bottom panel shows all three in log-log space to emphasize the differences. The exponent is off by more than a half between the two methods. But how see how good the second fit it ($r^2=0.98$) and the confidence interval is nowhere near the true value of $b$. The algorithm is very confident that $b\neq-4$.

The same effect will work in reverse with large changes in large numbers. For example, changing a data point from 10000 to 20000 in linear space will cause a huge change but in log-log space the change is not a big deal at all so the algorithm will again give misleading results.

close all
clc
clear all

x = [1:0.01:2 10];
y = 2*x.^(-4);
y(end)
y(end) = 0.00002;

% The power law fit on the original data
ft = fittype( 'power1' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
[fitresult1, gof1] = fit(x',y', ft, opts )

% The linear fit in log-log space
logx = log(x);
logy = log(y);
[fitresult2, gof2] = fit(logx',logy',fittype( 'poly1' ))

% The plots
x1 = min(x):0.01:max(x);
y1 = fitresult1.a*x1.^fitresult1.b;
y2 = exp(fitresult2.p1*log(x1)+fitresult2.p2);
figure
subplot(2,1,1)
plot(x,y,'o',x1,y1,'k',x1,y2,'r')
xlim([min(x),max(x)+1])
subplot(2,1,2)
loglog(x,y,'o',x1,y1,'k',x1,y2,'r')
xlim([min(x),max(x)+1])
legend('Data','Power Fit','Linear Fit')


Does non linear least squares manage a set of points with experimental/numerical errors?

I am assuming by manage you mean if NLLSF will give us an unbiased estimate of the parameters. The answer to that is, least squares fitting works if and only if the errors are gaussian. So if the errors are random errors (with the expected value and arithmetic mean being zero) then least squares gives you an unbiased estimate of the true parameter values. If you have any type of systematic error then least squares is not guaranteed to work. We often assume/want-to-think that we only have random error unless there is good evidence to the contrary.

And also does give the uncertainties in both a and b based on those parameters?

Again assuming errors being normally distrubuted, there are standard methods for estimating the error deviation. The mean is already assumed to be zero. We estimate the standard error for each parameter and then we use it to compute the confidence intervals. So yes, if I ONLY have experimental data and I have no idea what the real values of the parameters are, then I can still compute the confidence intervals.

Now an introductory bibliographic reference other than wikipedia would be very helpful.

I don't know of a good stats/regression book I can recommend. I would say just pick your favorite intro to stats book. And there are those papers by Clauset that I linked to both available on Arxiv for free.

TL;DR

1. Physicist think everything is a power law. It isn't. Even if the data looks like its power law, it probably isn't.

2. If you do legitimately have a power law, don't fit a straight line in log-log space. Just use NLLSF to fit a power law on the original data. If there is still any doubt, I can easily make up a numerical example and show you concretely how different the two answers can be.

-
It would be nice, for the sake of reference and having a concrete example. I'm kind of interested in misuses and misconceptions in the usage of mathematics in the contexts of physics, so thank you. – Jorge May 30 '14 at 19:50
NLLS is only (strictly) appropriate if your measurements have additive normally distributed errors. If the errors are mostly multiplicative and log-normal, you do want to use LLS in log space; if they're something more complicated, you either need more advanced methods, or you just pretend your errors are nice and normal / log-normal, even if they really aren't, and hope that your data is clean enough to still give you a decent fit. The real problem isn't so much finding the right fitting method, but that the right method depends on the error distribution, and you usually don't know that. – Ilmari Karonen May 30 '14 at 21:47
@IlmariKaronen Yes, that's why I said any kind of least squares (linear or non-linear) should only be done if the error distribution is gaussian. If you don't know or don't want to/can't assume anything about the error distribution then fitting in log-log space is just as invalid and we move into non-parametric methods. All I am saying is usually people want to do power-law-NLLSF and think that straight-line-LLSF in log-log space is equivalent while it isn't. – Fixed Point May 30 '14 at 22:50
@Jorge I have added an example. Please take a look. – Fixed Point Jun 1 '14 at 1:07
Thank you so much for this answer. However I think that there is a typo in the definition of $y=f(x)=2x^4$, shouldn't be $y=f(x)=2x^{-4}$. Does non linear least squares manage a set of points with experimental/numerical errors? And also does give the uncertainties in both $a$ and $b$ based on those parameters? Now an introductory bibliographic reference other than wikipedia would be very helpful. Thanks again. – Jorge Jun 1 '14 at 9:34

Are you sure that it is $q$ taking the negative value and not $\alpha$?

Non-integer powers of negative numbers are complex, so this surely is not what your experiment returned as values for $F$.

If you already know that $\beta$ can only take integer values, you could fit $$F_1 = \alpha q^1$$ $$F_2 = \alpha q^2$$ and $$F_3 = \alpha q^3$$ and show that $\beta = 1$ is the best fit.

If indeed it is $\alpha$ taking negative values along with $F$, then just fit $$-F = -\alpha q^\beta$$

If F only sometimes takes negative values, it might be shifted, i.e. $$F = \alpha q^\beta - \gamma$$ Then you can fit $$F + \gamma = \alpha q^\beta$$ where $\gamma$ is the minimum of $F$.

In your case, you can decide whether the charges are same sign or alternate sign. Write up the observation that alternate signs attract, while same signs repel and fit the absolute value of force over the absolute value of charge, finding a slope of $1/4\pi\epsilon_0$ and a power of 1. Of course, for the absolute values you can use the $\log$ thing, which (edit) apperently is not (/edit) very elegant!

-
Or the model is incorrect if $F$ can take negative values. – Bernhard May 30 '14 at 14:47
@Bernhard Very true! Still I inferred from the OP's statement that the model is given and therefore assumed to be true. – Neuneck May 30 '14 at 14:50
I was thinking something along this shifting line, I'm editing the question with more info. – Jorge May 30 '14 at 14:51
@Jorge I updated my post (the last paragraph is new). – Neuneck May 30 '14 at 15:23
@Neuneck I also thought about it, taking absolute values, but how may I justify that step? The shifting one is much more obvious but harder to implement. – Jorge May 30 '14 at 15:26

Taking the logarithms of negative numbers isn't a viable solution.

When you postulate $F=\alpha q^{\beta}$ you are implicitly supposing that is F is a real number, q is positive.

Just make a least square fit using the absolute values of $F$, $\alpha$ and $q$:

If $F=\alpha q^{\beta}$ then $|F|=|\alpha| |q|^{\beta}$, so $\ln(|F|) = \beta \ln(|q|)+\ln(|\alpha|)$.

That certainly will work.

Next study the sign: it is evident from data that with negative charge you get negative force. Then you can conclude that $\alpha$ is positive.

-
Could you please explain the implicit assumption? I don't see it. I don't see neither that the absolute values trick is allowed, but I thought about it. – Jorge May 30 '14 at 15:28
This is definitely the way to go, since Coulomb law is symmetrical in the sign of the charge (there is your implicit assumption). – Davidmh May 30 '14 at 15:32
Well, if you admit a power log as $q^{\beta}$ with beta real you need to admit $\beta$ real, otherwise you would need F to be a complex number. – Lelesquiz May 30 '14 at 15:33

This sort of problem is reasonably straightforward to solve numerically using Gauss-Newton or similar algorithms and your favourite programming language. Matlab even has an entire curve fitting toolbox which does everything very nicely.

These non-linear least squares techniques should work for any (differentiable) function, so have quite a lot of applications.

-