Understanding Residuals and Quality of Fit

Understanding Residuals and Quality of Fit#

Before applying model fitting to real experimental data, we need to understand what we mean by the “quality of fit” and how to quantify it. In this exercise, you will work with synthetic data where we know the true underlying relationship. This allows us to see clearly how residuals and the sum-of-squares error metric behave as we try different model parameters.

Generating Synthetic Data#

Let us start by creating some artificial data from a known linear relationship with added random noise. This simulates what happens in a real experiment: there is an underlying relationship between variables, but our measurements include some random experimental error.

Exercise: Generate Noisy Linear Data#

We will generate data from the linear relationship:

\[y = 2.0x + 3.0\]

and then add some random noise to simulate experimental error.

Create an array x containing 20 evenly spaced values between 0 and 10 using np.linspace().
Calculate the “true” y values using the equation above. Store these in an array called y_true.
Add random noise to create our “measured” data. Use the following code:

   np.random.seed(42)
   noise = np.random.normal(0, 1.5, size=20)
   y_data = y_true + noise

The function np.random.normal(0, 1.5, size=20) generates 20 random numbers drawn from a normal (Gaussian) distribution with mean 0 and standard deviation 1.5. These random values simulate the experimental errors that occur in real measurements.

Create a plot showing:
- The true line (plot x vs y_true as a solid line)
- The noisy measured data (plot x vs y_data as scattered points)
- Appropriate axis labels and a legend

This plot should show that whilst our measured data approximately follow the true linear relationship, individual points are scattered around the line due to random measurement error.

Calculating Residuals#

The residual for each data point is the difference between the observed y value and the y value predicted by our model. If we have a perfect model with perfect parameters, and our data had no noise, all residuals would be zero. In practice, we want to find parameters that make the residuals as small as possible.

For a linear model \(y = mx + c\), the residual for the \(i\)-th data point is:

\[r_i = y_i^\mathrm{data} - y_i^\mathrm{model} = y_i^\mathrm{data} - (mx_i + c)\]

Exercise: Calculate Residuals for the True Parameters#

Let us calculate the residuals when we use the true parameters \(m = 2.0\) and \(c = 3.0\).

Calculate the predicted y values using the true parameters: y_model = 2.0 * x + 3.0
Calculate the residuals: residuals = y_data - y_model
Print the first 5 residuals to see what they look like.
Create a plot showing:
- The model line (plot x vs y_model)
- The measured data points (plot x vs y_data as scatter)
- Vertical lines showing the residuals for each point using:
```
for i in range(len(x)):
    plt.plot([x[i], x[i]], [y_data[i], y_model[i]], 'r-', alpha=0.5)
```

Questions to consider:

Are the residuals positive or negative? What does the sign tell you?
Are the residuals roughly the same size as the noise we added (standard deviation 1.5)?
Do the residuals appear randomly distributed around zero, or is there a pattern?

Quantifying Quality of Fit with \(\chi^2\)#

Whilst individual residuals tell us how well our model fits each point, we need a single number to quantify the overall quality of fit. The sum-of-squares error, \(\chi^2\), achieves this by adding up the squared residuals:

\[\chi^2 = \sum_i r_i^2 = \sum_i \left[y_i^\mathrm{data} - y_i^\mathrm{model}\right]^2\]

Squaring the residuals ensures that positive and negative deviations both contribute positively to the error, and also gives more weight to larger deviations.

Exercise: Calculate \(\chi^2\) for the True Parameters#

Using the residuals you calculated above, compute \(\chi^2\) as the sum of squared residuals. You can use either:
- chi_squared = np.sum(residuals**2)
- chi_squared = np.sum((y_data - y_model)**2)
Print the value of \(\chi^2\). This represents how well the true model parameters fit our noisy data.

Note: Since we are using the true parameters that generated the data, this \(\chi^2\) value reflects only the random noise we added, not any systematic mismatch between model and data.

Exploring How χ² Changes with Parameters#

\(\chi^2\) provides a single number that quantifies the overall quality of fit: smaller χ\(\chi^2\) values indicate that the model predictions are closer to the observed data. This suggests a natural way to think about what we mean by the “best-fit” parameters.

If we want our model to fit the data as closely as possible, we should look for the parameters that make \(\chi^2\) as small as possible. This leads us to define the “best-fit” parameters as those that minimise \(\chi^2\).

Let us verify that \(\chi^2\) does indeed vary as we change our parameters, and observe that it is smallest for parameters close to the true values that generated our data.

Exercise: Calculate \(\chi^2\) for Different Parameter Values#

We will test several different combinations of slope \(m\) and intercept \(c\) to see how \(\chi^2\) varies.

Create a list of parameter pairs to test:

parameters_to_test = [
    (2.0, 3.0),   # True parameters
    (1.5, 3.0),   # Wrong slope
    (2.0, 4.0),   # Wrong intercept  
    (2.5, 2.0),   # Both wrong
    (1.0, 5.0),   # Both very wrong
]

For each parameter pair (m, c) in this list:
- Calculate the predicted y values: y_model = m * x + c
- Calculate \(\chi^2\): chi_squared = np.sum((y_data - y_model)**2)
- Print the parameters and the resulting \(\chi^2\) value

Format your output neatly, for example:

m = 2.0, c = 3.0: χ² = ...
m = 1.5, c = 3.0: χ² = ...