Linear Model Fitting with scipy.optimize.minimize()

Having seen that \(\chi^2\) varies as we change our model parameters, and that it is smallest when the parameters are close to the true values, we now have a strategy for finding the best-fit parameters: we need to search through parameter space to find the values that minimise χ².

This is an optimisation problem, and we can solve it using the same tool you encountered in Week 5 for geometry optimisation: scipy.optimize.minimize(). Just as you used minimize() to find the bond length that minimised potential energy, you can now use it to find the model parameters that minimise \(\chi^2\).

In this section, you will implement the complete model fitting workflow:

1. Define a model function

2. Define an error function that calculates \(\chi^2\) for given parameters

3. Use minimize() to find the best-fit parameters

4. Visualise the results

Note

This section assumes you are continuing directly from the previous page and have the x and y_data arrays already defined. If you need to regenerate the data, refer back to the previous page.

Defining the Model Function

Before we can fit a model to our data, we need to define the mathematical form of our model. For a linear model, we want a function that calculates predicted \(y\) values given input \(x\) values and parameters.

Exercise: Write a Linear Model Function

Write a function called linear_model() that implements the equation:

\[y = mx + c\]

Your function should:

• Take three arguments: x (the independent variable), m (the slope), and c (the intercept)

• Return the predicted y value(s)

• Work with both single values and numpy arrays for x

def linear_model(x, m, c):
    """Calculate y values for a linear model y = mx + c.
    
    Args:
        x (float): Independent variable value(s).
        m (float): Slope parameter.
        c (float): Intercept parameter.
        
    Returns:
        float: Predicted y value(s).
    """
    # Your code here

Test your function by calculating and plotting y_model = linear_model(x, 2.0, 3.0) alongside your data points to verify it works correctly.

Defining the Error Function

To use scipy.optimize.minimize(), we need an error function that:

• Takes the parameters we want to optimise as its first argument

• Takes the observed data as additional arguments

• Returns a single number representing the quality of fit (which we want to minimise)

For least-squares fitting, this error function should calculate \(\chi^2\).

Exercise: Write an Error Function

Write a function called error_function() that calculates the sum-of-squares error.

Your function should:

• Take params as the first argument (a list or array containing [m, c])

• Take x_data and y_data as additional arguments

• Extract the slope and intercept from params

• Use your linear_model() function to calculate predicted \(y\) values, given the input model slope and intercept

• Calculate and return \(\chi^2\) as the sum of squared residuals

def error_function(params, x_data, y_data):
    """Calculate sum-of-squares error for linear model.
    
    Args:
        params (list): List/array containing [slope, intercept].
        x_data (list): Observed x values.
        y_data (list): Observed y values.
        
    Returns:
        Chi-squared error value.
    """
    m, c = params
    # Your code here

Test your function by calculating the error for the true parameters [2.0, 3.0] and for some incorrect parameters like [1.5, 4.0]. The error should be smaller for the true parameters.

Finding the Best-Fit Parameters

Now we have everything needed to perform the optimisation. We will use scipy.optimize.minimize() to find the parameter values that minimise our error function.

Exercise: Optimise the Parameters

Use scipy.optimize.minimize() to find the best-fit parameters.

1. Import minimize:

   from scipy.optimize import minimize
   

2. Choose initial guesses for the parameters. Try initial_guess = [1.0, 1.0] (deliberately quite different from the true values).

3. Call minimize, passing:

   • Your error function

   • Your initial guess

   • The observed data as additional arguments using args=(x, y_data)

   result = minimize(error_function, initial_guess, args=(x, y_data))
   

4. Extract the optimised parameters from result.x and print them:

   m_fit, c_fit = result.x
   print(f"Best-fit parameters: m = {m_fit:.4f}, c = {c_fit:.4f}")
   

5. Print the final \(\chi^2\) value from result.fun.

6. Check whether the optimisation was successful by printing result.success.

We have just found “best-fit” parameters by minimising \(\chi^2\). It’s worth reflecting that this definition of “best fit” represents a choice we made. We could have defined our error function differently—for example, minimising the sum of absolute deviations or the maximum deviation—and obtained different “best-fit” parameters. Least-squares fitting (minimising \(\chi^2\)) is a common and often well-motivated choice; it is particularly well-suited to data where the measurement errors are normally distributed, as they are here. However, recognising that “best fit” is defined by our choice of error function, rather than being an inherent property of the data, is important.

Questions to consider:

• How close are the fitted parameters to the true values (\(m\) = 2.0, \(c\) = 3.0)?

• Why are they not exactly the same?

• How many iterations did the optimiser take? (Check result.nit)

Visualising the Fitted Model

The final step in any model fitting workflow is to visualise how well the fitted model matches the data.

Exercise: Plot the Data and Fitted Model

Create a plot that shows:

1. The measured data as scatter points with label “Data”

2. The true model line (using \(m\) = 2.0, \(c\) = 3.0) as a dashed grey line with label “True model”

3. The fitted model line (using your optimised parameters) as a solid line with label “Fitted model”

4. Add a legend, axis labels, and a title.

5. Optionally, add a text box showing the fitted parameters using:

   plt.text(0.5, 20, f'Fitted: y = {m_fit:.3f}x + {c_fit:.3f}', fontsize=10)
   

The fitted line should lie very close to (or on top of) the true model line, confirming that our optimisation successfully recovered the underlying relationship from the noisy data.

Testing Different Initial Guesses

A good optimisation algorithm should converge to the same minimum regardless of where it starts. Let us verify this for our linear fitting problem.

Exercise: Explore Different Initial Guesses

Repeat the optimisation with several different initial guesses:

1. Create a list of initial guesses to test:

   initial_guesses = [
       [1.0, 1.0],
       [0.0, 0.0],
       [5.0, 10.0],
       [-1.0, 8.0]
   ]
   

2. For each initial guess:

   • Run the optimisation

   • Print the initial guess, the optimised parameters, and the final \(\chi^2\)

3. Format your output clearly to compare the results.

Questions to consider:

• Do all starting points converge to approximately the same parameters?

• Do all starting points give approximately the same final \(\chi^2\)?

• Does the number of iterations vary depending on the starting point?

• Why is linear fitting generally more robust to initial guesses than some other optimisation problems?

Summary

You have now completed the full workflow for model fitting:

1. Define your model – A mathematical function that describes the relationship you expect

2. Define an error metric – A function that quantifies how well parameters fit the data (typically \(\chi^2\))

3. Optimise the parameters – Use minimize() to find parameters that minimise the error

4. Visualise the results – Plot the data and fitted model to assess quality

This same pattern applies to any model fitting problem, whether the model is linear or non-linear, and whether you have two parameters or twenty. The exercises that follow will give you practice applying this workflow to real chemical systems, where the models and data have physical significance.

Note

For linear regression specifically, there exist analytical solutions and specialized functions like scipy.stats.linregress() that are more efficient than using minimize(). However, the approach you have learned here is completely general and works for any model – including non-linear models where no analytical solution exists. This generality makes it a powerful tool for fitting diverse chemical models to experimental data.