Least squares fit versus averaging

Question

Suppose I want to find the spring constant of a spring from measured $F$ and $x$ data points. There are two basic ways to do this.

1. I could calculate the spring constant for each data point and average the results. $k = \frac{1}{N} \sum_i \frac{F_i}{x_i}$
2. I could find the least-squares fit of a linear function to the data, which would give me $k = \frac{\sum F_i x_i}{\sum x_i^2}$.

Is there a reason to prefer one of these methods over the other? I have been teaching my high school students to do the second as a simple introduction to data-fitting techniques, but I don't have a good answer for why they shouldn't just do the first, which feels more obvious to them.

Answer 1

That's a good question. The answer depends on how the data is collected. We will posit two simple models and furnish the best unbiased estimator for each model via the Gauss Markov theorem.

We will consider that the displacement $x$ is fixed in advance and that the force $F$ is then measured with error. Each model assumes the error is additive, mean zero, and independent. However, each model will differently describe the error variability.

---

First, suppose that the force $F$ is measured with error variance independent of the displacement $x$. Then, GM proves that $\frac{\sum_{i=1}^n x_iF_i}{\sum_{i=1}^n x_i^2}$ is the lowest variance unbiased estimator. Thus in this case you are correct that your students should use your estimator.

Second, suppose that the force $F$ is measured with error standard deviation proportional to $x$. Thus, larger forces are measured with more variability. Then, GM proves that $\frac{1}{n} \sum_{i=1}^n \frac{F_i}{x_i}$ is the lowest variance unbiased estimator. Thus in this case your students' intuition is dead-on.

---

In conclusion, which estimator is best depends on the extent the error variance depends on the displacement $x$. Do you think a given force and a force one hundred times larger are measured with equal precision? I suspect not. In practice, I would think the standard deviation is proportional to $x$ (and hence to $F$). Therefore I would use the estimator your students prefer.

Answer 2

Here is the quantitative (simulation) approach @Dr. Momo was mentioning: we can simulate some points for $x_i$ and $F_i$, give them some Gaussian errors, estimate $k$ with the two approaches and see how close the estimates are to the true parameter.

The assumptions underlying this simulation are likely a decent approximation for simple tabletop experiments, but one must be careful when choosing the method to apply to a specific problem to check how the uncertainties are distributed - not everything is Gaussian in reality, and figuring out these distributions is tricky business!

With that out of the way: these two estimators, while both unbiased, differ in their variance (@user551504 has discussed the reason why).

The result looks like this: the normalization is irrelevant, but one can clearly see the errors are on average smaller with the linear fit approach.

Here is the code (in python - it also requires the numpy and matplotlib packages):

import numpy as np
import matplotlib.pyplot as plt

# x and F errors' standard deviations
sigma_x = .1
sigma_F = .1

# true elastic parameter
k_true = 1.

# number of mock measurement points
n_measurement_points = 20

# number of trials to run
n_trials = 1000

def trial(seed: int):

    # generate N equally-spaced x points between 1 and 20, inclusive
    xi = np.linspace(1, 20, num=n_measurement_points)
    # generate the "true" forces, according to the linear law
    # (without the negative sign for simplicity, it doesn't change anything)
    Fi = k_true * xi

    # add Gaussian random noise to x and F
    rng = np.random.default_rng(seed)    
    xi += rng.normal(scale=sigma_x, size=xi.shape)
    Fi += rng.normal(scale=sigma_F, size=Fi.shape)

    # estimate the elastic parameter with the two approaches
    k_estimate_avg = 1 / n_measurement_points * np.sum(Fi / xi)
    k_estimate_fit =  np.sum(Fi * xi) / np.sum(xi**2)
    
    return k_estimate_avg - k_true, k_estimate_fit - k_true

# run the trial several times to get statistics for the variance
# of the estimator
avg_errs = []
fit_errs = []
for n in range(n_trials):
    # we seed the RNG with the index of the for loop - doesn't really matter,
    # it's useful for reproducibility
    err_avg, err_fit = trial(n)
    avg_errs.append(err_avg)
    fit_errs.append(err_fit)

# make a histogram of the results 
plt.hist(avg_errs, alpha=.5, density=True, label='average')
plt.hist(fit_errs, alpha=.5, density=True, label='fit')
plt.xlabel('Estimated minus true $k$')
plt.ylabel('Probability density')
plt.legend()

# to save the figure, one can do
# plt.savefig('linear_fit_vs_avg.png', dpi=150)

# to show it interactively, instead, do
# plt.show()

I'll also note that the result is qualitatively the same even if errors are only included for $x$ or for $F$.

Answer 3

As in all things statistical, it completely depends on your assumptions.

Why averages are useful: If you have a variable x which you can sample, and it comes from a normal distribution with some mean $\mu$ and standard deviation $\sigma$, it turns out that the maximum likelihood estimator for the true mean $\mu$ is the sample mean $\frac{1}{N} \sum x_{i}$: the sample mean is your best guess of the true mean. An example of this process is a very large gacha board (or Galton's board). If you don't know where the balls keep being dropped from, you can estimate it by measuring locations of where the balls end up. To briefly introduce the relevant math, when you try to maximize the product of Gaussians representing the probability of $\mu$ being the true mean: $\Pi_{i} \exp{(x_{i}-\mu)^{2}}$ and take a derivative = 0, you find $\sum (x_{i} - \mu) = 0$ so $\mu = \frac{1}{N}\sum x_{i}$ is the highest likelihood guess for the mean.

Why least-squares are useful: This one is in the name. It minimizes the least-squared error. This one becomes a maximum likelihood estimator when your residuals $\epsilon = F - kx$ are normally distributed: occurring with relative probability $P(\epsilon) \propto \exp{(-\epsilon^{2}/\sigma^{2})}$. This time, when trying to estimate $k$ to maximize the likelihood, we take derivatives with respect to $k$ which is inside of $\epsilon = F - kx$. Each term of the product will contribute a $2\epsilon \frac{d\epsilon}{dk} = 2\epsilon x = 2 Fx - 2 kxx$. Again, the sum of these will give us the best guess for k: $\sum F_{i}x_{i} - \sum k x_{i}^{2} = 0$ so $k = \frac{\sum F_{i}x_{i}}{\sum x_{i}^{2}}$

So the question is, do you think each time you sample a $F_{i}, x_{i}$ you are sampling a $k_{i} = F_{i}/x_{i}$ which comes from a normal distribution centered around the true constant k, or do you think the spring has a constant $k$ which will give you normally distributed residuals in $F_{i}$, $x_{i}$?

I think a great lesson for your students is not to just take averages or do least squares out of gut feeling, but to understand that these are commonly used mathematical machinery for a reason, derived by reason, and they come with assumptions.

Answer 4

The linearity approach has a few advantages:

1. It gives less weight to points near the origin. This is often a good thing, since those points generally have larger relative error bars.

2. It smooths out noise more effectively and is less impacted by outliers than the first approach.

You can play around with this in Excel and see how the two approaches differ. I'm sure that somebody can offer more quantitative insights, but that's my experience with the two approaches.

As a word of warning, though, if you're using the linearity approach, you should collect data at equal intervals along your domain. Even professionally, I often see linearity plots with 4-5 points clustered at x = 1 and a single point at x = 10. This proves very little about the linearity of the system.

Answer 5

I think the answer depends on what you want to do with the spring constant you extract, but generally I think the second approach is more useful.

Presumably you want a simple model of a spring that yields a force (via Hooke's law) with a small error, given the displacement. Unless you care about a specific range of displacements, this is precisely what the second approach gives you, minimizing the error in the least squares sense.

The first approach still minimizes some error, but one that may not necessarily be as useful. It minimizes the error in $F/x$ in the least squares sense, which as Dr. Momo points out emphasizes errors in $F$ for small displacements. This means that the spring constant extracted this way is a worse model for large displacements.

If there is a specific range of displacements you care more about, you could also use a weighting function $w(x)$ that prioritizes the more important displacements, yielding $$k = \frac{\sum w^2(x_i) F_i x_i}{\sum w^2(x_i) x_i^2}$$ which minimizes the error in $w(x)F$ in the least squares sense. $w(x) = 1/x$ for your first approach and $w(x) = 1$ for the second.

Answer 6

The difference is in the underlying hypothesis and what one is trying to verify:

If we know for sure that the spring is linear (i.e., $F=kx$ for any $x$), then the first method is just as good as the second. It is also the method to use when making multiple measurements for the same point (even if the spring is nonlinear, but spring constant is evaluated for each point separately).

On the other hand, if the issue of the linearity is itself to be checked, then the second method is the way to go. One could demonstrate this by plotting $F(x)$ for a non-linear restoring force: it is then clear why the first method would fail miserably. The second however is still applicable in the linear region (i.e., by trancating the dataset, to keep only small $x$).

Update
The question is really about statistics, rather than physics (and it could be a good diea to transfer it to Cross Validated).

Point estimate
In the first case we are performing a point estimate of the spring constant, $k$, provided the data $\{k_1, k_2,...,k_N\}$, using as estimator $$ \hat{k}=\frac{1}{N}\sum_{i=1}^N k_i, $$ whose error can be characterized by the standard deviation $$ \sigma_k=\sqrt{\frac{1}{N-1}\sum_{i=1}^N(\hat{k}-k_i)^2}. $$ That $k_i=F_i/X_i$ are not measured directly is a feature of particular setup - one could imagine a measurement device that measures directly $k$ rather than pairs $(x_i, F_i)$. Any such device would have an erorr, which could be accounted for in addition to the estimation error.

Such estimate could be done by repeatedly extending the spring to the same length - the procedure thus works, in principle, for non-linear springs. The focus here is the value of $k$, not the Hooke's law.

Regression analysis
In the second case one is performing regression analysis, verifying the Hooke's law (i.e., the linear relationship between the force and the extension of the spring). One than could use, e.g., the optimal least squares approach, minimizing $$ f(k|\{x_i,F_i\}) = \frac{1}{N}\sum_{i=1}^N(kx_i-F_i)^2 $$ in respect to $k$. The standard deviation coudl then be calculated as the inverse of the Fisher information, and we could perform a statistical test to decide whether the spring is linear or not with the given level of confidence (for a physicist it rpetty much means comparing the estimated value of $k$ with its standard deviation, but there are more formal statistical procedures).

Remark
Optimal least squares is known to be very sensitive to outliers - one stray point may significantly affect the estimate. One therefore ofte resorts to robust regression, such as least absolute deviations, Theil-Sen estimator, or Huber regression. These require some numerics (apart from Theil-Sen), which is available in many statistical packages, e.g., in R or Python.

Answer 7

Your title is a bit inaccurate. Both methods you propose fall under a broad meaning of "average". We could also use other methods, such as $k = \frac{\sum_i F_i}{\sum_i x_i}$, $k = \left(\frac 1 N\sum_i \frac{x_i}{F_i}\right)^{-1}$, or $k = \frac{\sum_i F_i^2}{\sum_i x_iF_i}$. All of these are "averages". When someone says that they "took the average" from some dataset, there's actually a lot of choice as to just what the "average" is, and it's important to be aware of the possibility of quite different answers.

Each of these formulas weight different datapoints differently. With $k = \frac 1 N\sum_i \frac{F_i}{x_i}$, small errors at datapoints with small $x$ values are magnified, while $k = \frac{\sum_i F_i}{\sum_i x_i}$ can give comparatively more weight to large values. What formula to use depends on what you're assuming about how the errors are distributed.