Data analysis and fitting

Why it matters

A measurement without an uncertainty is not a result. Physicists fit models to data every day: a decay constant from counts, a spring constant from extensions, g from a pendulum. The fit gives the best parameters, but the questions that matter are how sure are we and does the model actually fit. Weighted least squares and chi-squared answer both.

Uncertainties first

Each data point y_i has an uncertainty sigma_i, its standard deviation. Points with small sigma should count more. Two rules for propagating uncertainties you will use constantly (for independent errors):

sum or difference:  sigma_f^2 = sigma_a^2 + sigma_b^2
product or ratio:   (sigma_f / f)^2 = (sigma_a / a)^2 + (sigma_b / b)^2

Weighted least squares for a line

To fit y = m x + b, choose m and b that minimise

chi^2 = sum of ((y_i - (m x_i + b)) / sigma_i)^2

With weights w_i = 1 / sigma_i^2, the answer has a closed form:

S = sum w     Sx = sum w x     Sy = sum w y
Sxx = sum w x^2     Sxy = sum w x y     D = S*Sxx - Sx^2

m = (S*Sxy - Sx*Sy) / D          sigma_m = sqrt(S / D)
b = (Sxx*Sy - Sx*Sxy) / D        sigma_b = sqrt(Sxx / D)

Judging the fit: chi-squared

Each term in chi^2 asks: how many sigmas away is this point from the line? If the model is right and the sigma values are honest, each term is about 1 on average. With N points and 2 fitted parameters there are N - 2 degrees of freedom, and

reduced chi^2 = chi^2 / (N - 2)   ~ 1 for a good fit

• Much larger than 1: the model is wrong, or the uncertainties are underestimated.
• Much smaller than 1: the uncertainties are probably overestimated.
• Always look at the residuals too: a pattern (a curve, a trend) means the model misses something even if chi^2 looks fine.

Worked example: a cart on a track

Position of a cart every second, each with sigma = 0.03 m:

t (s): 0     1     2     3     4     5
x (m): 0.02  0.51  1.03  1.48  2.05  2.49

The weighted fit gives

m = 0.498 +/- 0.007 m/s     (the speed)
b = 0.019 +/- 0.022 m       (consistent with starting at 0)
chi^2 = 3.60 for 4 degrees of freedom -> reduced chi^2 = 0.90

The reduced chi-squared near 1 says a constant speed describes the data at this precision. Report it as v = 0.498 +/- 0.007 m/s: round the uncertainty to one or two significant figures and the value to match.

With NumPy and SciPy

import numpy as np

t = np.array([0, 1, 2, 3, 4, 5.0])
x = np.array([0.02, 0.51, 1.03, 1.48, 2.05, 2.49])
sigma = np.full(6, 0.03)

coef, cov = np.polyfit(t, x, 1, w=1 / sigma, cov="unscaled")
print(coef, np.sqrt(np.diag(cov)))   # [0.498 0.019] [0.0072 0.0217]

Note that np.polyfit takes weights 1/sigma, not 1/sigma^2, and cov="unscaled" returns errors based on your sigmas. For nonlinear models, scipy.optimize.curve_fit(f, x, y, sigma=sigma, absolute_sigma=True) does the same job; absolute_sigma=True tells it your sigmas are real uncertainties rather than relative weights.

Common mistakes

• Unweighted fits for data with different uncertainties, which lets noisy points pull the line.
• Passing 1/sigma^2 to polyfit (it squares the weights itself).
• Quoting too many digits, such as 0.4977143 +/- 0.0071714.
• Ignoring a reduced chi-squared of 5 because the line "looks fine".
• Forgetting the degrees of freedom: divide by N - p, where p is the number of fitted parameters.

Recap

• Weighted least squares minimises chi-squared with weights 1/sigma^2.
• Parameter errors come from the same sums: sigma_m = sqrt(S/D), sigma_b = sqrt(Sxx/D).
• Reduced chi-squared near 1 plus patternless residuals means a believable fit.

Try it

Open the Do tab to write a weighted line fit with parameter uncertainties and chi-squared, checked against NumPy.