---
jupyter:
  jupytext:
    notebook_metadata_filter: all,-language_info
    split_at_heading: true
    text_representation:
      extension: .Rmd
      format_name: rmarkdown
      format_version: '1.2'
      jupytext_version: 1.14.1
  kernelspec:
    display_name: Python 3 (ipykernel)
    language: python
    name: python3
---

# Sum of squares, root mean square

The [mean and slopes page](mean_and_slopes) used the Sum of Squared Error (SSE)
as the measure of how well a particular slope fits the data.

Here we will think about another, derived measure, called the Root Mean Square Error.

First, let us go back to the original problem.

## Hemoglobin and Packed Cell Volume, again

```{python}
# Run this cell.
import numpy as np
import matplotlib.pyplot as plt
# Make plots look a little bit more fancy
plt.style.use('fivethirtyeight')
# Print to 4 decimal places, show tiny values as 0
np.set_printoptions(precision=4, suppress=True)
import pandas as pd
pd.set_option('mode.copy_on_write', True)
```

We go back to the [data on chronic kidney
disease](../data/chronic_kidney_disease).

Download the data to your computer via this link: {download}`ckd_clean.csv
<../data/ckd_clean.csv>`.

We were interested to find the "best" slope for a line that relates the blood
hemoglobin measures in each patient, with their corresponding Packed Cell
Volume (PCV) blood measure:

```{python}
# Run this cell
ckd = pd.read_csv('ckd_clean.csv')
hgb = np.array(ckd['Hemoglobin'])
pcv = np.array(ckd['Packed Cell Volume'])
# Plot HGB on the x axis, PCV on the y axis
plt.plot(hgb, pcv, 'o')
plt.xlabel('Hemoglobin concentration')
plt.ylabel('Packed cell volume')
```

We are assuming our line goes through `hgb == 0` and `pcv == 0`, on the basis
that the red blood cells measured by PCV contain all the hemoglobin, so if PCV
is 0 then hemoglobin must be 0 and vice versa.  Our *predicted* values for PCV
for each patient are the given by some slope`s` multiplied by the corresponding
hemoglobin value for that patient.

For example, say the slope is 2.25, then the predicted PCV values are:

```{python}
s = 2.25
predicted_pcv = s * hgb
predicted_pcv[:10]
```

The *errors* for these predictions are:

```{python}
errors = pcv - predicted_pcv
errors[:10]
```

## Sum of Squared error

We decided to use the Sum of Squared Error (SSE) as a measure of the quality of
fit.

The sum of squared error is nothing but:

```{python}
sse = np.sum(errors ** 2)
sse
```

Here's a function to return the SSE.  We call this the *cost function* because
it calculates a *cost* in terms of error for a given slope.

```{python}
def calc_sse(slope):
    predicted_pcv = hgb * slope  # 'hgb' comes from the top level
    errors = pcv - predicted_pcv # 'pcv' comes from the top level
    return np.sum(errors ** 2)
```

```{python}
# Recalculate SSE for slope of interest.
calc_sse(2.25)
```

Then we used `minimize` to find the slope that minimizes the SSE:

```{python}
from scipy.optimize import minimize
# 2.25 below is the slope value to start the search.
res_sse = minimize(calc_sse, 2.25)
res_sse
```

Notice the warning - there was some inaccuracy in the calculation.   We  may try slightly different cost functions to avoid that.

## Root mean squared error

One problem with the SSE is that it can get very large, especially if there are
a large number of points.  When the values get large, this can cause
calculation problems.  Very large numbers for SSE can make
it harder to think about the SSE values.

One way of allowing for the number of points, is to use the *mean* square error
(MSE), instead of the SSE.  Here is the MSE cost function.

```{python}
def calc_mse(slope):
    predicted_pcv = hgb * slope  # 'hgb' comes from the top level
    errors = pcv - predicted_pcv # 'pcv' comes from the top level
    return np.mean(errors ** 2)  # Notice mean instead of sum
```

We try `minimize` again, to find the best slope.

```{python}
# 2.25 below is the slope value to start the search.
res_mse = minimize(calc_mse, 2.25)
res_mse
```

Notice there is no warning this time in the optimization message.


Notice that running `minimize` on the MSE found the same slope as it did for
the SSE. The resulting cost function `fun` value is the original SSE value
divided by the number of points. (It's not exactly the same because of the
inaccuracy of the SSE estimate).

```{python}
res_sse.fun / len(hgb)
```

It's also common to use the *square root* of the MSE, or the Root Mean Square
Error.  This has one advantage that it's easier to compare the RMSE to the
individual errors — we don't normally think easily in terms of squared error.
This is the RMSE cost function — another tiny variation.

```{python}
def calc_rmse(slope):
    predicted_pcv = hgb * slope  # 'hgb' comes from the top level
    errors = pcv - predicted_pcv # 'pcv' comes from the top level
    return np.sqrt(np.mean(errors ** 2))  # Notice square root of mean.
```

```{python}
# 2.25 below is the slope value to start the search.
res_rmse = minimize(calc_rmse, 2.25)
res_rmse
```

The cost function `fun` value is just the square root of the MSE `fun` value:

```{python}
np.sqrt(res_mse.fun)
```

Notice again that `minimize` on RMSE found exactly the same slope as for MSE and SSE.   Why?

Reflect a little, and then read on.

(monotonicity)=

## Monotonicity

SSE, MSE, and RMSE give the same parameters from `minimize` because they are
*monotonic* with respect to each other.  Put another way, the SSE and the RMSE
have a monotonic relationship.  Put yet another way,
the RMSE is a [monotonic
function](https://en.wikipedia.org/wiki/Monotonic_function) of the SSE.

Remember, the Root Mean Square Error (RMSE) is just the square root of the Mean
Square Error (MSE), and the MSE is the Sum of Squared Error (SSE) divided by
the number of values in the x array (which is the same as the number of values
in the y array).

For example, imagine we have the following set of SSE values:

```{python}
n_slopes = 10000
some_slopes = np.linspace(2, 4, n_slopes)
SSEs = np.zeros(n_slopes)
for i in np.arange(n_slopes):
    SSEs[i] = calc_sse(some_slopes[i])

plt.scatter(some_slopes, SSEs)
plt.xlabel('Slope')
plt.xlabel('SSE')
plt.title('SSE as a function of slope');
```

We can calculate the matching MSE values by dividing by `n` (the number of
`hgb` or `pcv` values).

```{python}
n = len(hgb)
MSEs = SSEs / n  # Corresponding Mean Square Error
```

When we plot the MSE values as a function of (y-axis) the SSE values (x axis),
we see that there is a straight line relationship, as we would expect, with a
slope of `1 / n`.

```{python}
plt.scatter(SSEs, MSEs)
plt.xlabel('SSE values')
plt.ylabel('Corresponding MSE values')
plt.title('MSE as a function of SSE');
```

This straight line relationship is one example of a *strictly increasing*
monotonic relationship.

Consider any two SSE values, call these `S1` and `S2`.  There are corresponding
MSE values `M1` and `M2` (in fact given by `S1 / n` and `S2 / n`).  The
relationship between SSE and MSE is strictly increasing monotonic because, for
any values `S1` and `S2`:

* when `S1 > S2` then it is also true that`M1 > M2`
* when `S1 < S2` then `M1 < M2`.
* when `S1 == S2`, `M1 == M2`.

You can see that relationship in the plot above, because MSE always goes up, as
SSE goes up.

The same is true of the RMSE values as a function of the SSE values:

```{python}
RMSEs = np.sqrt(MSEs)  # Corresponding Root Mean Square Error
plt.scatter(SSEs, RMSEs)
plt.xlabel('SSE values')
plt.ylabel('Corresponding RMSE values')
plt.title('RMSE as a function of SSE');
```

Notice, again, that the RMSE always goes up as SSE goes up.

Now consider `minimize`.  It is looking at the parameters to minimize the value
for  — say - SSE.  But, you can probably see from the argument above, that if
we found the parameters to minimize SSE, those must also be the parameters to
minimize MSE, and to minimize RMSE, because once we have found the smallest
SSE, that must also be the smallest MSE or RMSE.

Therefore, we can minimize on the RMSE, and we are
guaranteed to get the same parameters, within calculation precision, as we
would if we had minimized on the SSE or the MSE.

In the rest of these textbook pages, we will nearly always use RMSE, because
it's a common metric, and one that is relatively easy to interpret as a typical
error for single point.
