Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
206 changes: 154 additions & 52 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,67 +6,123 @@ Use:
nox -s test
```
While coding, use `nox -s lint`, and `nox -s typing` to check your code.
Autoformatting help is available via `nox -s format`.
Autoformatting help is available via `nox -s` format.
Feel free to read more about nox at https://nox.thea.codes/en/stable/ .

Today's goal is to get more familiar with the three important concepts of linear regression, fitting a polynomial of higher order and regularization.
You will be given some data in the form of pairs of (a, b)-values.
So for each a-value there is one b-value.
The main idea is to find an easy function (i.e. a polynomial) that best explains our data.
That means, that if we plug in an a-value, the result should be close enough to the corresponding b-value.

### Part 1: Proof of concept
Use `b = pandas.read_csv('./data/noisy_signal.tab')` to load a noisy signal.
The first part will be concerned with modelling this signal using polynomials.
A general polynomial $f$ of order $n$ is given by:

#### Regression:
Linear regression is usually a good first step. Start by implementing the function
`set_up_point_matrix` from the `src/regularization.py` module.
The function should produce polynomial-coordinate matrices $\mathbf{A}_n$ of the form:
$$ b = f(a) = c_1 + c_2 \cdot a^1 + c_3 \cdot a^2 + ... + c_n \cdot a^{n-1} $$

The $c$-values are the coefficients of the polynomial - these numbers are to be estimated!
The $a-$ and $b-$values are already given.
If we plug all the given values into the general form of the polynomial we get a system of linear equations which can be reformulated as a matrix multiplication:

$$
\mathbf{A}_n =
\begin{pmatrix}
1 & a_1^1 & a_1^2 & \dots & a_1^{n-1} \\\\
1 & a_2^1 & a_2^2 & \dots & a_2^{n-1} \\\\
1 & a_3^1 & a_3^2 & \dots & a_3^{n-1} \\\\
\vdots & \vdots & \vdots & \ddots & \vdots \\\\
1 & a_m^1 & a_m^2 & \dots & a_m^{n-1} \\\\
\end{pmatrix} \cdot
\begin{pmatrix}
c_1 \\\\
\vdots \\\\
c_n \\\\
\end{pmatrix} = \begin{pmatrix}
b_1 \\\\
b_2 \\\\
b_3 \\\\
\vdots \\\\
b_m \\\\
\end{pmatrix}
$$

With n=2,
Or in short:

$$\mathbf{A}_n\mathbf{c} = \mathbf{b}$$

The optimal $\mathbf{c}$ is given by

$$\mathbf{A}_2^{\dagger}\mathbf{b} = \mathbf{x}$$
$$\mathbf{A}_n^{\dagger}\mathbf{b} = \mathbf{c} ,$$

where $\mathbf{A}_n^{\dagger}$ is the Pseudo-Inverse.

For linear regression, the polynomial will be of order n=2 and will look like this:

$$ f(a) = c_1 + c_2 \cdot a $$

This will just be a straight line and there are only two coefficients $c_1$ and $c_2$ that have to be estimated.
Very often, this line is too simple to explain the data sufficiently.
That is why we want to fit polynomials of higher order, so $n > 2$.
Unfortunately, the more complex the model gets (i.e. the higher the order of the polynomial gets), the more noise will be tracked. Here we can make use of regularization techniques.
In the first and following part, you will be given artificial data and in the second part you will make use of real data!

### Part 1: Proof of concept
The line `b = pandas.read_csv('./data/noisy_signal.tab')` is used to load a noisy signal.
The line `x_axis = np.linspace(0, 1, num=len(b_noise))` will provide you with corresponding x-values (these are the a-values from above and the lecture).
This is some artificial data that serves as a means to try out the concepts you have learned about in the lecture.
The first part will be concerned with modeling this signal using polynomials.

will produce the coefficients for a straight line. Evaluate your first-degree polynomial via $ax+b$.
Plot the result using `matplotlib.pyplot`'s `plot` function.
#### ⊙ Task 1.1: Regression
Linear regression is usually a good first step.

1. Start by implementing the function `set_up_point_matrix` from the `src/regularization.py` module.
The function should produce polynomial-coordinate matrices $\mathbf{A}_n$ of the form:
$$
\mathbf{A}_n =
\begin{pmatrix}
1 & a_1^1 & a_1^2 & \dots & a_1^{n-1} \\\\
1 & a_2^1 & a_2^2 & \dots & a_2^{n-1} \\\\
1 & a_3^1 & a_3^2 & \dots & a_3^{n-1} \\\\
\vdots & \vdots & \vdots & \ddots & \vdots \\\\
1 & a_m^1 & a_m^2 & \dots & a_m^{n-1} \\\\
\end{pmatrix}
$$

2. Go to the main-function and use the function you just implemented to create the point-matrix A for n=2.

#### Fitting a Polynomial to a function:
The straight line above is insufficient to model the data. Using your
implementation of `set_up_point_matrix` set $n=300$ (to set up a square matrix) and fit the polynomial
by computing
3. Now,
$$\mathbf{A}_2^{\dagger}\mathbf{b} = \mathbf{c} = \begin{pmatrix}
c_1 \\\\
c_2 \\\\
\end{pmatrix} $$
will produce the coefficients for a straight line.

$$\mathbf{A}^{\dagger}\mathbf{b} = \mathbf{x}_{\text{fit}}.$$
4. Evaluate your first-degree polynomial via $c_1 + c_2 \cdot x$ and plot the result as well as the original data using `matplotlib.pyplot`'s `plot` function.

Having estimated the coefficients

$$\mathbf{A} \mathbf{x}_{\text{fit}}$$

computes the function values. Plot the original points and the function values using matplotlib.
#### ⊙Task 1.2: Fitting a Polynomial to a function
The straight line above is insufficient to model the data.
So perform the very same steps as above, but change the degree of the polynomial to n=300 (to set up a square matrix since we have 300 data-points):
1. Set up the point matrix by setting n=300.
2. Estimate the coefficients via $$\mathbf{A}^{\dagger}\mathbf{b} = \mathbf{x}_{\text{fit}}.$$
3. Having estimated the coefficients $$\mathbf{A} \mathbf{x}_{\text{fit}}$$ computes the function values.
Plot the original points and the function values using matplotlib.
What do you see?



#### Regularization:
Unfortunately, the fit is not ideal. The polynomial tracks the noise.

#### ⊙Task 1.3: Regularization
Unfortunately, the fit is not ideal. The polynomial is too complex and tracks the noise.
The singular value decomposition (SVD) can help!
Recall that the SVD turns a matrix

$$\mathbf{A} \in \mathbb{R}^{m,n}$$

into the form:

$$ \mathbf{A} = \mathbf{U} \Sigma \mathbf{V}^T
$$
$$\mathbf{A} = \mathbf{U} \Sigma \mathbf{V}^T$$

In the SVD-Form computing, the inverse is simple. Swap $U$ and $V$ and replace every of the m singular values with it's inverse
In the SVD-Form, computing the pseudoinverse is simple! Swap U and V and replace every of the m singular values with it's inverse

$$1/\sigma_i .$$

Expand All @@ -75,26 +131,36 @@ This results in the matrix
\Sigma^\dagger = \begin{pmatrix}
\sigma_1^{-1} & & & \\\\
& \ddots & \\\\
& & \sigma_n^{-1} \\\\ \hline
& & \sigma_m^{-1} \\\\ \hline
& 0 &
\end{pmatrix}^T
\end{pmatrix}
```

A solution to the overfitting problem is to filter the singular values.
A solution to the overfitting problem is to filter the singular values.
The idea is, that small singular values often correspond to directions in the data where noise dominates. And now we want to create a filter matrix, that gets rid of singular values if they are too small (i.e. if they fall below a threshold $\epsilon$).
Compute a diagonal for a filter matrix by evaluating:

$$f_i = \sigma_i^2 / (\sigma_i^2 + \epsilon)$$

The idea is to compute a loop over $i$ for all of the m singular values.

Roughly speaking multiplication by $f_i$ will filter a singular value when

$$\sigma_i \lt \epsilon .$$
$$\sigma_i \lt \epsilon ,$$
since in this case, $f_i$ will be close to $0$.
If however
$$\sigma_i \geq \epsilon ,$$
$f_i$ will be closer to $1$ and the respective singular value will not be filtered out.

Apply the regularization by computing:


$$
\mathbf{x}_r= \mathbf{V} \mathbf{F} \mathbf{\Sigma}^\dagger
\mathbf{x}_r= \mathbf{V} \mathbf{F} \begin{pmatrix}
\sigma_1^{-1} & & & \\\\
& \ddots & \\\\
& & \sigma_n^{-1} \\\\ \hline
& 0 &
\end{pmatrix}
\mathbf{U}^T \mathbf{b} = \mathbf{V} \mathbf{F} \mathbf{\Sigma}^\dagger
\mathbf{U}^T \mathbf{b}
$$

Expand All @@ -103,39 +169,75 @@ with

$$\mathbf{V} \in \mathbb{R}^{n,n}, \mathbf{F} \in \mathbb{R}^{n,n}, \Sigma^{\dagger} \in \mathbb{R}^{n,m}, \mathbf{U} \in \mathbb{R}^{m,m} \text{ and } \mathbf{b} \in \mathbb{R}^{m,1}.$$

Setting $n=300$ turns $A$ into a square matrix. In this case, the zero block in the sigma-matrix disappears.
Plot the result for epsilon equal to 0.1, 1e-6, and 1e-12.
Setting n=300 turns A into a square matrix. In this case, the zero block in the sigma-matrix disappears and you don't have to worry about transposing $\sigma$ when computing the pseudoinverse.

#### Model Complexity (Optional):
Another solution to the overfitting problem is reducing the complexity of the model.
To assess the quality of polynomial fit to the data, compute and plot the Mean Squared Error (Mean Squared Error (MSE) measure how close the regression line is to data points) for every degree of polynomial up to 20.
To sum it up, your tasks are:
1. Compute the SVD of A.

Perform the following steps 2. - 4. for epsilon equal to 0.1, 1e-6, and 1e-12.
2. Compute the diagonal for the filter matrix and turn it into a matrix.
3. Estimate the regularized coefficients by applying the formula above.
4. Plot the result.

MSE can be calculated using the following equation, where $N$ is the number of samples, $y_i$ is the original point and $\hat{y_i}$ is the predictied output.

#### ✪Task 1.4: Model Complexity (Optional):
Another solution to the overfitting problem is reducing the complexity of the model.
To assess the quality of polynomial fit to the data, compute and plot the Mean Squared Error (Mean Squared Error measure how close the regression line is to data points) for every degree of polynomial upto 20.
So as before:
1. Set up the point matrix for the current degree from 1 to 20.
2. Estimate the coefficients.
3. Compute the predictions.
4. Calculate the MSE.
MSE can be calculated using the following equation, where $N$ is the number of samples, $y_i$ is the original point and $\hat{y_i}$ is the predicted output.
$$MSE=\frac{1}{N} \sum_{i=1}^{N} (y_i-\hat{y_i})^2$$
5. Plot the MSE-error against the degree.
6. Are the degree of the polynomial and the MSE linked?
From the plot, estimate the optimal degree of polynomial and fit the polynomial with this specific degree.


Are the degree of the polynomial and the MSE linked?

### Part 2: Real data analysis

### Part 2: Rhine water level analysis
Now we are ready to deal with real data! Feel free to use your favorite time series data or work with the Rhine level data we provide.
The file `./data/pegel.tab` contains the Rhine water levels measured in Bonn over the last 100 years.
Data source: https://pegel.bonn.de.

#### Regression:
#### ⊙Task 2.1 Regression
The `src/pegel_bonn.py` file already contains code to pre-load the data for you.
Make the Rhine level measurements your new vector $\mathbf{b}$.

Generate a matrix A with n=2 using the timestamps for the data set and compute
The Rhine level measurements will be your new vector $\mathbf{b}$ from before.
Now we want to do the same as in Part 1 and start with linear regression!
1. Generate a matrix A with n=2 using the timestamps for the data set as your x-values.
2. Compute $$\mathbf{A}^{\dagger}\mathbf{b}$$ to estimate the coefficients of the line.
3. Evaluate your polynomial and plot the result.
4. Compute the zero-intercept with the y-axis. When do the regression line and the x-axis intersect?
Or in other words: On which day will the Rhine water level be at 0 cm?
> **Hint:** Plug in $y=0$ into the equation of your line with the estimated coefficients and solve for the date $x$.

$$\mathbf{A}^{\dagger}\mathbf{b}.$$

Plot the result. Compute the zero. When do the regression line and the x-axis intersect?

#### Fitting a higher order Polynomial:
#### ⊙ Task 2.2: Fitting a higher-order Polynomial

Re-using the code you wrote for the proof of concept task, fit a polynomial of degree 20 to the data. Before plotting have a closer look at `datetime_stamps` and its values and scale the axis appropriately.
Plot the result.

#### Regularization:
1. Scale the x-axis.
2. Set up the point matrix for the scaled x-axis with n=20.
3. Compute the coefficients.
4. Evaluate the polynomial.
5. Plot the result.




#### ⊙Task 2.3: Regularization
Focus on the data from the year 2000 onward and filter the singular values.
Matrix A is not square in this case. Consequently, a zero block must appear in your singular value matrix.
Plot filtered eigen-polynomials using epsilon equal to 0.1, 1e-3, 1e-9.
We will use again a degree of 20.
Matrix A is not square in this case, because the degree is smaller than the number of datapoints! Consequently, a zero block must appear in your singular value matrix and when computing the Pseudoinverse from the SVD, $\sigma$ has to be transposed!
Like in Part 1:
1. Compute the SVD of the point matrix from the previous task.

Perform the following steps 2. - 4. for epsilon equal to 0.1, 1e-3, and 1e-9.
2. Compute the filter matrix.
3. Estimate the regularized coefficients by applying the formula from before.
> **Hint:** Remember the zero-block! You need degree-many rows and number-of-datapoints-many columns!
4. Evaluate the regularized polynomial and plot the results.

Loading