Linear Curve Fitting with numpy.polyfit

When engineers and scientists collect data from experiments or sensors, the measurements are almost never perfectly clean. There is always some degree of noise, instrument error, or natural variability. Curve fitting is the process of finding a mathematical model — a smooth, continuous function — that best represents the underlying trend in a set of noisy, discrete data points. Rather than connecting the dots (which would simply reproduce the noise), a fitted curve generalizes from the data, capturing the dominant behavior while ignoring random fluctuations.

The standard mathematical strategy for finding the "best" curve is to minimize the sum of squared residuals (SSR). A residual is the vertical distance between an observed data point and the value predicted by the model at the same x-coordinate. Squaring each residual before summing ensures that positive and negative errors do not cancel each other out, and it penalizes large deviations more heavily than small ones. The set of model parameters that produces the smallest possible SSR is called the least-squares solution. This approach, known as ordinary least squares (OLS), has a closed-form analytical solution for polynomial models, which is exactly what NumPy exploits internally.

A critical concept in curve fitting is the bias-variance trade-off. A model that is too simple (high bias) will systematically miss the trend in the data — this is called underfitting. A model that is too complex (high variance) will chase every noise spike and fail to generalize to new measurements — this is called overfitting. A good fit sits in the middle: complex enough to capture the real signal, but simple enough to ignore noise. Choosing the right model complexity is often the most important decision in curve fitting, and it requires judgment, visualization, and sometimes statistical tests.

NumPy provides numpy.polyfit as a straightforward, powerful tool for polynomial regression. Its call signature is:

coefficients = numpy.polyfit(x, y, deg)

Here, x is a 1-D array of independent variable values, y is a 1-D array of corresponding dependent variable measurements, and deg is a non-negative integer specifying the degree of the polynomial to fit. The function returns a 1-D array of deg + 1 coefficients in descending power order. For example, a degree-2 fit returns [a, b, c] representing the polynomial a·x² + b·x + c.

For a degree-1 (linear) fit, polyfit returns two coefficients [m, b] representing the slope and intercept of the best-fit line y = m·x + b. This is the most common use case in engineering when a linear relationship is expected between two variables. For example, suppose you measure the extension of a spring at different applied forces and want to verify Hooke's Law:

import numpy as np
import matplotlib.pyplot as plt

force   = np.array([0.5, 1.0, 1.5, 2.0, 2.5, 3.0])   # Newtons
stretch = np.array([0.12, 0.25, 0.35, 0.51, 0.60, 0.75])  # meters (with noise)

coeffs = np.polyfit(force, stretch, 1)
print(coeffs)   # e.g., [0.2504, 0.0017]
# Interpretation: slope ≈ 0.25 m/N  (spring constant k = 1/0.25 = 4 N/m)
#                 intercept ≈ 0.002 m (small offset, nearly zero as expected)

Once you have the coefficients, you need to evaluate the polynomial at new x values to draw the fitted curve or make predictions. NumPy provides two convenient tools for this:

• numpy.polyval(coeffs, x_new) — evaluates the polynomial defined by coeffs at each point in x_new. It is a simple function call, ideal for a quick evaluation.
• numpy.poly1d(coeffs) — creates a callable polynomial object. Once created, you call it like a regular Python function. It also has useful attributes like .roots and supports arithmetic between polynomials. This is often more convenient in larger programs.

# Using poly1d
p = np.poly1d(coeffs)        # create callable object
x_smooth = np.linspace(0.5, 3.0, 200)
y_fitted = p(x_smooth)       # evaluate at 200 evenly spaced points

# Equivalently with polyval
y_fitted = np.polyval(coeffs, x_smooth)

It is important to understand that increasing the polynomial degree always reduces training error. A polynomial of degree n-1 can pass exactly through n data points, giving zero residuals but also perfectly encoding all the noise. This is why a higher-degree fit should not be automatically preferred.

Consider the problem of selecting the right polynomial degree. If your data follows a clearly linear trend, using deg=1 is appropriate. If there is a known curvature — for example, the trajectory of a projectile under gravity is quadratic — then deg=2 may be correct. In general, you should choose the lowest degree that captures the meaningful trend.

Underfitting is visible when the fitted curve is systematically above the data in some regions and below it in others, even though the data clearly curves. For example, fitting a straight line to data that follows a distinct parabolic arc will leave large, structured residuals. Overfitting is visible when a high-degree polynomial wiggles wildly between the data points, especially near the edges of the x-range (a phenomenon called Runge's phenomenon).

# Example illustrating underfitting vs. overfitting
np.random.seed(42)
x_data = np.linspace(0, 2*np.pi, 15)
y_data = np.sin(x_data) + np.random.normal(0, 0.15, size=15)

x_plot = np.linspace(0, 2*np.pi, 300)

for deg in [1, 3, 12]:
    coeffs = np.polyfit(x_data, y_data, deg)
    y_fit  = np.polyval(coeffs, x_plot)
    plt.plot(x_plot, y_fit, label=f'deg={deg}')

plt.scatter(x_data, y_data, color='black', zorder=5, label='Data')
plt.legend()
plt.show()
# deg=1  → clearly underfits (straight line through sinusoidal data)
# deg=3  → captures the general wave shape, reasonable fit
# deg=12 → fits training points almost exactly but oscillates wildly

A visual check — plotting the fitted curve over the data — is a fast first diagnostic. However, a more rigorous tool is the residual plot: a scatter plot of the residuals (y_data − y_fitted at each data point) against the corresponding x values. For a good fit, residuals should look like random noise scattered uniformly around zero with no discernible pattern. If residuals show a systematic curve (a "U" shape or "S" shape), this indicates that the model degree is too low and the fit is underfitting. If residuals for a high-degree polynomial have large magnitudes and wild swings, overfitting is likely.

# Residual plot example
coeffs = np.polyfit(x_data, y_data, 3)
y_pred = np.polyval(coeffs, x_data)
residuals = y_data - y_pred

plt.scatter(x_data, residuals, color='red')
plt.axhline(0, color='black', linestyle='--')
plt.xlabel('x')
plt.ylabel('Residuals')
plt.title('Residual Plot — Degree 3 Fit')
plt.show()

Not all real-world relationships are polynomial. Many engineering systems follow exponential or power-law behavior. Rather than fitting these directly with a nonlinear solver, you can often linearize the model and then apply polyfit with deg=1.

For an exponential model of the form y = a · e^(bx), taking the natural logarithm of both sides gives:

ln(y) = ln(a) + b·x

This is linear in x. You fit ln(y) vs. x with polyfit(..., 1). The returned coefficients are [b, ln(a)]. To recover the original parameters, note that a = e^(intercept).

# Exponential fit via linearization
# Suppose y = 2.5 * exp(0.4 * x) + noise
np.random.seed(0)
x_exp = np.linspace(0, 5, 20)
y_exp = 2.5 * np.exp(0.4 * x_exp) + np.random.normal(0, 0.5, 20)

# Linearize: take natural log of y
log_y = np.log(y_exp)          # requires y > 0 for all points

coeffs_lin = np.polyfit(x_exp, log_y, 1)
b_fit  = coeffs_lin[0]         # slope → b
a_fit  = np.exp(coeffs_lin[1]) # intercept → ln(a), so a = exp(intercept)

print(f"a ≈ {a_fit:.3f},  b ≈ {b_fit:.3f}")
# Should recover values close to a=2.5, b=0.4

For a power-law model of the form y = a · x^b, taking the log of both sides gives:

log(y) = log(a) + b · log(x)

This is linear in log(x). You fit log(y) vs. log(x) with polyfit(..., 1). The slope is b and the intercept is log(a), so a = 10^(intercept) if base-10 logarithms are used, or a = e^(intercept) for natural logs.

# Power-law fit via linearization
# Suppose y = 3.0 * x^1.5 + noise
np.random.seed(1)
x_pow = np.linspace(1, 10, 20)
y_pow = 3.0 * x_pow**1.5 + np.random.normal(0, 1.5, 20)

log_x = np.log(x_pow)
log_y_pow = np.log(y_pow)

coeffs_pow = np.polyfit(log_x, log_y_pow, 1)
b_pow = coeffs_pow[0]
a_pow = np.exp(coeffs_pow[1])

print(f"a ≈ {a_pow:.3f},  b ≈ {b_pow:.3f}")
# Should recover values close to a=3.0, b=1.5

An important caveat about linearization: the least-squares criterion is minimized in the transformed (log) space, not in the original space. This implicitly assumes that the measurement noise is multiplicative (proportional to the magnitude of y) rather than additive (constant absolute noise). If the noise is truly additive and the signal spans a large range of magnitudes, the linearized fit can give somewhat different results than a direct nonlinear fit. In practice, for moderate noise levels and well-behaved data, linearization works well and is far simpler to implement.

To produce a professional, engineering-quality plot of the fitted curve alongside the raw data, you should always follow these steps:

• Use numpy.linspace to generate a dense array of x values spanning the data range. This ensures the fitted curve appears smooth, not as a connected polygon between data points. For example, np.linspace(x_data.min(), x_data.max(), 300) creates 300 evenly spaced evaluation points.
• Create a poly1d object from the coefficients and call it on the linspace array to generate y values for the smooth curve.
• Plot the raw data as a scatter plot (using plt.scatter or plt.plot with marker-only format) and the fitted curve as a line with a distinct color or style, so the two are visually distinguishable.
• Always include a title, axis labels with units, and a legend identifying the data series and the model. Without these elements, a figure is incomplete for any engineering report or publication.

import numpy as np
import matplotlib.pyplot as plt

# --- Data ---
force   = np.array([0.5, 1.0, 1.5, 2.0, 2.5, 3.0])
stretch = np.array([0.12, 0.25, 0.35, 0.51, 0.60, 0.75])

# --- Fit ---
coeffs = np.polyfit(force, stretch, 1)
p      = np.poly1d(coeffs)

# --- Smooth evaluation points ---
x_smooth = np.linspace(force.min(), force.max(), 300)
y_smooth = p(x_smooth)

# --- Plot ---
plt.figure(figsize=(7, 4))
plt.scatter(force, stretch, color='steelblue', zorder=5,
            label='Measured data')
plt.plot(x_smooth, y_smooth, color='tomato', linewidth=2,
         label=f'Linear fit: y = {coeffs[0]:.3f}x + {coeffs[1]:.4f}')
plt.xlabel('Applied Force (N)')
plt.ylabel('Spring Extension (m)')
plt.title('Spring Extension vs. Applied Force — Linear Curve Fit')
plt.legend()
plt.tight_layout()
plt.show()

This example brings together every element of the workflow: collecting data into NumPy arrays, calling polyfit with the appropriate degree, creating a callable poly1d object, generating a smooth evaluation grid with linspace, and producing a clearly labeled figure. The coefficient values printed in the legend legend give the reader quantitative information about the model without requiring them to look at a separate table. By following this complete process — fitting, evaluating, and visualizing with proper labels — you produce results that are both technically rigorous and clearly communicated.