Least Squares: Fitting a Straight Line with Python

Given a set of points, how do we find the straight line that best represents them? Least squares turns this problem into a computable objective. We derive the formula using five observed points, then use Python to compute the line and plot the residuals.

Writing "best fit" as an objective

Let the $i$-th observed point be $(x_i,y_i)$, and let the fitted line be $\widehat y=ax+b$. At the same $x_i$, the vertical difference between the observed value and the fitted value $e_i=y_i-(ax_i+b)$ is called the residual.

Least squares chooses the slope and intercept that minimize $S(a,b)$. After squaring, positive and negative residuals do not cancel out. Larger residuals are also penalized more heavily, so outliers can noticeably pull the fitted line.

Why the fitted line passes through the sample center

Let $\bar x$ and $\bar y$ be the sample means. Split the error into a centered part and a constant part, then use the fact that the sum of the centered variables is zero to eliminate the cross term.

When $a$ is fixed, the first term is unchanged, and the second term reaches zero at $b=\bar y-a\bar x$. Therefore the fitted line must pass through $(\bar x,\bar y)$. Substituting this intercept back in, taking the derivative with respect to $a$ and setting the derivative to zero gives

As long as the $x_i$ are not all the same, the denominator is positive, and what remains is an upward-opening quadratic, so we get a unique global optimum. If all the $x_i$ are the same, the data can only determine the best predicted value at that x-coordinate, $\bar y$, and the slope and intercept cannot be determined separately.

Fitting with Python and plotting the residuals

Click the run button on the code block. The scatter points are the observed values, the solid line is the fitted line, and the dashed segments mark the vertical residuals of each point. The code directly computes the formula above and prints the sum of squared residuals.

import matplotlib.pyplot as plt

x = [0, 1, 2, 3, 4]
y = [1.0, 2.1, 2.7, 4.2, 5.1]

n = len(x)
x_mean = sum(x) / n
y_mean = sum(y) / n
sxx = sum((xi - x_mean) ** 2 for xi in x)
if sxx == 0:
    raise ValueError("At least two distinct x values are required.")

sxy = sum((xi - x_mean) * (yi - y_mean) for xi, yi in zip(x, y))
slope = sxy / sxx
intercept = y_mean - slope * x_mean
fitted = [slope * xi + intercept for xi in x]
residuals = [yi - fit for yi, fit in zip(y, fitted)]
sse = sum(residual ** 2 for residual in residuals)

print(f"Slope: {slope:.3f}")
print(f"Intercept: {intercept:.3f}")
print(f"Sum of squared residuals: {sse:.4f}")

fig, ax = plt.subplots(figsize=(6.4, 3.8), layout="constrained")
ax.plot(x, fitted, color="#1d4ed8", label="Least-squares line")
ax.vlines(x, fitted, y, colors="#b45309", linestyles="--",
          linewidth=2, label="Residuals")
ax.scatter(x, y, color="#0f172a", zorder=3, label="Observations")
ax.set(xlabel="x", ylabel="y", title="Least-squares fit")
ax.grid(alpha=0.18)
ax.legend(loc="upper left")

For this dataset we get $\widehat y=1.030x+0.960$, and the sum of squared residuals is $0.1390$. You can change the last observed value from $5.1$ to $8.0$ and run it again to see how the point on the right pulls the slope up. Goodness of fit only describes how well the line fits this dataset. When using it to predict values outside the observed range, you also need to judge whether there is reason for the linear relationship to continue to hold.