Experimental Physics – How to Fit Power to Experimental Data

Question

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.

Answer 1

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.

---

Addendum:

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')

---

Second Addendum:

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.