# Data Fitting Exercises

# 2. Linear Least-Squares

### How to define a line of &ldquo;best fit&rdquo;

In the previous exercise, you used the `linregress` function to calculate a line of best fit to your experimental data. But what is this function actually doing, and how do we define the quality of fit for a particular straight line?

In this exercise, you will work through the maths behind linear regression, and then learn how to implement this using functions in Python as a more general optimisation problem.

### The maths

Linear regression addresses the problem of how to find a line of &ldquo;best fit&rdquo; for some data that is appriximately described by a straight line. For example. we might have the data shown in the figure below:



and we want to find the straight line that best describes the relationship between $x$ and $y$:



How do we find the &ldquo;best&rdquo; line for this particular data set?

The equation of a straight line is

\begin{equation}
y=mx+c.
\end{equation}

This is a mathematical function that defines a set of $y$ values from a set of $x$ values, for any pair or parameters, $m$ and $c$.

More generally, we can say that $y$ is a function of some model parameters, $P$, and the input $x$:

\begin{equation}
y = f(P,x).
\end{equation}

In our example of a straight line, the parameter set $P$ contains the slope and intercept of the line:

\begin{equation}
P=\left\{m,c\right\},
\end{equation}

and the function $f$ is

\begin{equation}
f(P,x) = P_0 x + P_1.
\end{equation}

Now $P_0$ is the slope, and $P_1$ is the intercept.

In principle, we could choose any combination of $P_0$ and $P_1$ for our model parameters. Every combination gives a different straight line, which might be a better, or worse, fit to our experimental data. 

To decide which set of model parameters best describes our real data, we need a way to score the agreement between the model $f(P,x)$ and the data. A common way to quantify the error between the model and the real data is to calculate the **sum of squared errors**. For every $x$ value in our input data set, we can calculate the $y_\mathrm{predicted}$ value predicted by the model. The difference between this predicted value and the real $y$ value for this data point is the **error** for that data point. 



To quantify the overall error of this particular model, we can add up the *squares* of all the error terms:

\begin{equation}
\mathrm{error} = \sum_i\left[y_i - f(P,x_i)\right]^2
\end{equation}

Adding *squares* means that positive and negative errors contribute equally to the total error.

We can now determine mathematically whether one model is a better &ldquo;fit&rdquo; to the data than another. A better quality of fit will give us a smaller error. The **best** fit to the data (for this particular model function) is given by the set of parameters that **minimises** the error. The procedure of finding a set of model parameters that minimises the sum of squared errors is often called **least-squares** fitting.


### The code

In Python the model straight-line function can be expressed using the following function:
>```python
def model_function( P, x ):
    m = P[0]
    c = P[1]
    y = m * x + c
    return y
```

If you compare this to the mathematical function definition above, you should see the same structure in both the mathematical description and the Python representation.  

The function `model_function()` also looks similar to the `line()` function you used in Exercise 1.  

```python
def line( m, c, x ):
    y = m * x + c
    return y
```

The only difference is that in `model_function()` we pass in the model **parameters** as a **list**. This is then unpacked inside the function, instead of passing in $m$ and $c$ separately.

<div class="alert alert-success">
In the code cell below, define the functions <span style='font-family:monospace'>model_function()</span> and <span style='font-family:monospace'>line()</span>.  
Using <span style='font-family:monospace'>x = 2</span>, check that both these functions return the same $y$ values for pairs of parameters, $m$ and $c$. For the <span style='font-family:monospace'>model_function()</span> function, remember that the parameters need to be passed in as a **list**.
</div>


```python

```

We can also write a Python function that calculates the error between the predictions of our model function and the real $y$ values:

>```python
def error_function( P, x, y ):
    y_predicted = model_function( P, x )
    error_terms = y - y_predicted
    total_error = np.sum( error_terms**2 )
    return total_error
```

* The first line calls the `model_function()` function with the list of model parameters, $P$, and the input $x$ values. This returns a set of *predicted* $y$ values for this particular model.

* The second line calculates the error terms, i.e. the differences between the predicted and actual $y$ values.

* The third line calculates the sum of squared errors.

* The fourth line returns the error.

<div class="alert alert-success">
By copying the appropriate code from Exercise 1, or writing it from scratch, read in the experimental data from <span style='font-family:monospace'>'data/equilibrium_constant.dat'</span> and plot it as points.  <br/><br/>
Using the <span style='font-family:monospace'>model_function()</span> function, generate a series of lines with different parameters, and plot these against the experimental data.   <br/><br/>
Define the function <span style='font-family:monospace'>error_function()</span>. For each of your guessed model parameter sets, calculate the error. You should see that models that appear closer to the experimental data give you lower errors.
</div>


```python

```

You should have seen how changing the parameters for your model lets you generate different lines, that each do better or worse jobs of describing the experimental data. In each case, the quality of fit between your model and the data is given by the error, calculated with your error function.

The final step in the problem is to find the **least-squares** solution. This corresponds to the set of model parameters that **minimise** the result from your error function.

Finding the minimum or maximum of a function is a common mathematical problem, and there are a large number of algorithms for doing this numerically, using computers.

The `scipy` module contains functions for doing exactly this, in `scipy.optimize`. For this exercise, you will import the `minimize()` function, and use this to find the parameters that minimise the output of your error function. You can import the `minimize()` function as follows:

>```python
from scipy.optimize import minimize
```

The `minimize()` function can be given a [large number of options](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html), but the simplest usage requires three arguments:
>```python
minimize( function_to_minimise, initial_guess, other_arguments )
```

- The first argument, `function_to_minimise`, is the name of the function to minimise. For this exercise, this is `error_function`.
- The second argument, `initial_guess`, is a starting guess for the parameters that you want to fit. e.g. you might guess that `P=[1.0,1.0]`.
- The third argument is a list of any other arguments that need to be passed into the function you are minimising. In this exercise, your function `error_function` takes three arguments: `error_function( P, x, y )`: The first is the set of model parameters, which you want to change during the fitting. The second and third are the **observed** $x$ and $y$ values (here these are $1/T$ and $\ln{K}$) for each data point. These do not change during the fit, and get passed into `minimize()` as a list in the `other_arguments` position.

e.g. if you have stored the $1/T$ ($x$ values) and $\ln{K}$ ($y$ values) for the experimental data in numpy arrays called, respectively, `inverse_temperature` and `ln_K` then you could write

>```python
P_initial = [ 1.0, 1.0 ]
other_args = inverse_temperature, ln_K
minimize( error_function, P_initial, other_args )
```

The `other_arg = inverse_temperature, ln_K` command stores the $x$ and $y$ data values in a data type called a **tuple**. From the perspective of this exercise, a tuple is like a list, in that it can store multiple values (here you store two `numpy` arrays).

<div class="alert alert-success">
In the following code cell, import the <span style='font-family:monospace'>minimize</span> function, and use it to minimise your error function, with the experimental data you loaded previously.
</div>


```python

```

Your output should look like:
```
      fun: 0.026018593683863903
 hess_inv: array([[  1.53892987e+06,  -4.95784961e+03],
       [ -4.95784961e+03,   1.60346929e+01]])
      jac: array([  1.53901055e-07,  -1.45286322e-07])
  message: 'Optimization terminated successfully.'
     nfev: 60
      nit: 4
     njev: 15
   status: 0
  success: True
        x: array([ 6917.390252  ,   -21.05540855])
```       
You will see that a lot of information about the fit is produced. For this exercise, the important parts of the output are the `message`, which tells you whether the optimisation was successful, and the result labelled `x`, which gives the input parameters that minimise your error function. The value of the error function with these optimised model parameters is given by the `fun` result.

You can access specific parts of the output by appending the variable names to the end of your `minimize()` call, e.g. to get the optimal parameters you can write

>```python
minimize( error_function, P_initial, other_args ).x
```

<div class="alert alert-success">
Repeat your optimisation, and store the optimised model parameters in a variable <span style='font-family:monospace'>P_opt</span>. Check these values against the fitted slope and intercept you found in Exercise 1.
</div>

21e4f0482649d0a6c83292947cdc588e651fc1b3
ipynb
Jupyter Notebook
Data_Fitting_Exercise_S2_2.ipynb
pythoninchemistry/chem_data_analysis_jupyter
4af545f1a8acdded28d96508bb5adc8929da92cc
CC-BY-4.0
2019-05-05T00:21:55.000Z
2021-09-16T14:15:15.000Z
Data_Fitting_Exercise_S2_2.ipynb
pythoninchemistry/chem_data_analysis_jupyter
4af545f1a8acdded28d96508bb5adc8929da92cc
CC-BY-4.0
Data_Fitting_Exercise_S2_2.ipynb
pythoninchemistry/chem_data_analysis_jupyter
4af545f1a8acdded28d96508bb5adc8929da92cc
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES