Skip to content
Merged
Show file tree
Hide file tree
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
106 changes: 81 additions & 25 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,28 +1,34 @@
# Sample Statistics and Gaussians Exercise

### Exercise 1: Mean and Variance
This exercise first explores sample statistics like mean and variance.
Today's exercise first explores sample statistics like mean and variance.
It continues with the definition of the gaussian distributions.
It ends with applying gaussian mixture models (gmm) and the expectation-maximization algorithm used to optimize gmm.

### ⊙ Task 1: Mean and Variance

- To get started, take a look at the `src/sample_mean_corr_rhein.py` file.

Implement functions to compute sample mean and the standard deviation.
Use

$$ \hat{\mu} = \frac{1}{n} \sum_{i=1}^n x_i , $$
1. Use

to calculate the mean and
$$ \hat{\mu} = \frac{1}{n} \sum_{i=1}^n x_i , $$

$$ \hat{\sigma} = \sqrt{\frac{1}{n-1} \sum_{i=1}^n (x_i - \hat{\mu})^2} $$
to calculate the mean.

to compute the standard deviation. $x \in \mathbb{R}$ denotes individual sample elements, and $n \in \mathbb{N}$ the size of the sample.
2. Use

Return to the Rhine data-set. Load the data from `./data/pegel.tab`. Compute the water level mean and standard deviation before and after the year 2000.
$$ \hat{\sigma} = \sqrt{\frac{1}{n-1} \sum_{i=1}^n (x_i - \hat{\mu})^2} $$

to compute the standard deviation.
$x_i \in \mathbb{R}$ for $i \in \{1, ... , n\}$ denotes individual sample elements, and $n \in \mathbb{N}$ the size of the sample.
Don't use the pre-build functions np.mean() or np.std() to solve these tasks.

Return to the Rhine data-set and go to the `__main__`-function.
The data from `./data/pegel.tab` is already loaded and processed and is ready to be used.
3. Compute the water level mean and standard deviation before the year 2000.
4. And now compute the water level mean and standard deviation after the year 2000.

### Exercise 2: Autocorrelation
### ⊙ Task 2: Autocorrelation
We now want to use autocorrelation to analyse the discrete time signal of the rhine level measurements. Implement the `auto_corr` function in `src/sample_mean_corr_rhein.py`. It should implement the engineering version without the normalization and return the autocorrelation

$$ R_{xx} = (c_{-N+1},\ldots,c_{1}, c_0, c_{1}, \ldots, c_{N-1}) $$
Expand All @@ -31,36 +37,81 @@ with

$$ c_{k} = \sum_{t=1}^{N-|k|} n_t n_{t + |k|}$$

with $n$ the normalized version of your signal of length $N$. The time shift $k$ moves from $-(N-1)$ to $N-1$. Therefore, the resulting array has a length of $2N-1$. For example the autocorrelation of an input signal $x=(2,3,-1)$ is $R_{xx}=(c_{-2}, c_{-1}, c_0, c_1, c_2)=(-2, 3, 14, 3, -2)$ and is symmetrical. Make sure that you normalize the signal *before* giving it to your `auto_corr` function. Once you have checked your implementation using `nox -s test`, you can use `np.correlate` for efficiency. Plot the autocorrelation for the rhine level measurements since 2000.
Normalize your data via
with $n$ the normalized version of your signal of length $N$. The time shift $k$ moves from $-(N-1)$ to $N-1$. Therefore, the resulting array has a length of $2N-1$.

For example the autocorrelation of an input signal $x=(2,3,-1)$ is $R_{xx}=(c_{-2}, c_{-1}, c_0, c_1, c_2)=(-2, 3, 14, 3, -2)$ and is symmetrical.

>>> In the following table you can see an illustrative depicition on how the $c_k$'s are calculated.
The header contains the input signal x padded with 0's on its sides.
For autocorrelation, we compute correlation between $x$ and $x$ itself.
So in visual terms, we slide $x$ from left to right across itself.
At each step we compute one $c_k$ by first multiplying the numbers that are aligned with the input signal in the header. Then, these products will be summed up. The result is written in the respective cell of the last column.

| 0 | 0 | 2 | 3 | -1 | 0 | 0 | $c_k$ |
| -------- | ------- | ------- | ------- | ------- | ------- | ------- | ------- |
| 2 | 3 | -1 | | | | | 0 + 0 - 2 = -2
| | 2 | 3 | -1 | | | | 0 + 6 - 3 = 3
| | | 2 | 3 | -1 | | | 4 + 9 + 1 = 14
| | | | 2 | 3 | -1 | | 6 - 3 + 0 = 3
| | | | | 2 | 3 | -1 | -2 + 0 + 0 = -2
>>>As you can see, when reading from top to bottom we get the correct solution $R_{xx}=(-2, 3, 14, 3, -2)$.

So here are your tasks:
1. Implement the `auto_corr` function as described above.
>>> The function expects $x$ to be normalized. That means, that no normalization is done inside the `auto_corr`-function. Instead you normalize the input signal before calling `auto_corr`.
2. Check your implementation using `nox -s test`.
If the test passes you can use `np.correlate` for efficiency in the following exercises!

$$ n_{t} = \frac{x_{t} - \hat{\mu}}{\hat{\sigma}} ,$$
Now go back to the `__main__`-function and consider the Rhine data-set after 2000.
3. Normalize the data of the Rhine level measurements since 2000 via

for all $t$ measurements until the signal length N. Before running the autocorrelation computation. Compare the autocorrelation to a random signal from `np.random.randn` by plotting both results with `plt.plot`.
$$ n_{t} = \frac{x_{t} - \hat{\mu}}{\hat{\sigma}} ,$$

for all $t$ measurements until the signal length N.
4. Compute and plot the autocorrelation for the Rhine level measurements since 2000.

### Exercise 3: Distributions
Now we want to compare this autocorrelation to the one of a random signal.
5. Create a random signal from `np.random.randn` that has the same shape as the Rhine level measurements
6. Normalize the random signal.
7. Compute the autocorrelation of the normalized, random signal.
8. Plot both autocorrelations of the Rhine level measurements and the random signal using `plt.plot` and compare the results.


- Consider the `src/plot_gaussian.py` module. Implement the `gaussian_pdf` function.
### ⊙ Task 3: Distributions

In one dimension gaussian probability density function is defined as
1. Consider the `src/plot_gaussian.py` module. Implement the `gaussian_pdf` function.

$$\phi_1(x | \mu, \sigma) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp({-\frac{1}{2}(\frac{x - \mu}{\sigma})^2}) .$$
In one dimension gaussian probability density function is defined as

$\pi \in \mathbb{R}$ denotes Pi, $\mu \in \mathbb{R}$ the mean and $\sigma \in \mathbb{R}$ the standard deviation for a random variable $X$. $e^x$ denotes the exponential function. $X$ having a gaussian pdf is described as gaussion or normal distribution $\mathcal{N}$. Explore the behavior of $\mathcal{N}(\mu, \sigma)$ for different values of $\mu$ and $\sigma$.
$$\phi_1(x | \mu, \sigma) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp({-\frac{1}{2}(\frac{x - \mu}{\sigma})^2}) .$$

$\pi \in \mathbb{R}$ denotes Pi, $\mu \in \mathbb{R}$ the mean and $\sigma \in \mathbb{R}$ the standard deviation for a random variable $X$.
$e^x$ denotes the exponential function.
A random variable $X$ having a gaussian pdf is described as gaussion or normal distribution $\mathcal{N}$.
>>> Remark: In the notation $\phi_1(x | \mu, \sigma)$, x is the variable that is plugged into the function and $\mu$ and $\sigma$ are parameters which are needed to define the function and that are determined beforehand.
2. Explore the behavior of $\mathcal{N}(\mu, \sigma)$ for different values of $\mu$ and $\sigma$.
The Code for plotting the pdf's is already given.

- Consider the `src/mixture_concpets.py` module.
Implement a two-dimensional gaussian pdf following,

$$ \phi_2(\mathbf{x} | \mu_g, \Sigma_g) = \frac{1}{\sqrt{(2\pi)^2 \| \Sigma_g \|}} \exp({-\frac{1}{2}(\mathbf{x}-\mu_g)^T \Sigma_g^{-1}(\mathbf{x}-\mu_g)}).$$

$\mu_g \in \mathbb{R}^2$ denotes the two dimensional mean vector, $\Sigma_g \in \mathbb{R}^{2\times2}$ the covariance matrix, $^{-1}$ the matrix inverse, $T$ the transpose and $g \in \mathbb{N}$ the number of the distrubtion, which will be important later.
Plot a 2d-bell curve with $\mu_1 = [-1.5, 2]$ and $\Sigma_1 = [[1, 0], [0, 1]]$ using the `plt.imshow` function. `np.linspace` and `np.meshgrid` will help you.
3. Consider the `src/mixture_concepts.py` module.
Go to the `twod_gaussian_pdf`-function and implement a two-dimensional gaussian pdf following,

$$ \phi_2(\mathbf{x} | \mu_g, \Sigma_g) = \frac{1}{\sqrt{(2\pi)^2 \| \Sigma_g \|}} \exp({-\frac{1}{2}(\mathbf{x}-\mu_g)^T \Sigma_g^{-1}(\mathbf{x}-\mu_g)}).$$

### Exercise 4: Gaussian mixture models (optional)
$\mu_g \in \mathbb{R}^2$ denotes the two dimensional mean vector, $\Sigma_g \in \mathbb{R}^{2\times2}$ the covariance matrix, $^{-1}$ the matrix inverse, $T$ the transpose, $\| \|$ the determinant and $g \in \mathbb{N}$ the number of the distrubtion, which will be important later.

- As you can see, the x-parameter of the function is a grid of shape (grid_height, grid_width, 2).
That means, that you get not only one but grid_height*grid_width many 2-dimensional values that should be evaluated.
It's up to you how you want to approach this task.
But Broadcasting might be an elegant way to evaluate all these values at the same time (https://numpy.org/doc/stable/user/basics.broadcasting.html). In this case, you might take a look at `np.swapaxes` to deal with the transponation.
- At the very end of this document, we included a diagram that depicts how the shapes should develop when using Broadcasting. This is purely optional and just for when you need some guidance regarding the relevant shapes.

4. Plot a 2d-bell curve with $\mu_1 = [-1.5, 2]$ and $\Sigma_1 = [[1, 0], [0, 1]]$ using the `plt.imshow` function. `np.linspace` and `np.meshgrid` will help you.



### ✪ Task 4: Gaussian mixture models (optional)

We can use bell-curve sums for classification! A Gaussian mixture model has the density

Expand Down Expand Up @@ -113,3 +164,8 @@ Train a gmm to find the diabetic patients.


- Standard packages like sci-kit-learn implement GMMs, too. Take a minute to read https://scikit-learn.org/stable/modules/mixture.html .

---
If you need some inspiration for broadcasting in Task 3.3:
![broadcasting_hint](./figures/broadcasting.png)

4 changes: 3 additions & 1 deletion src/mixture_concepts.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ def twod_gaussian_pdf(x: np.ndarray, mu: np.ndarray, sigma: np.ndarray) -> np.nd
Returns:
np.ndarray: The two dimensional gaussian distribution.
"""
# TODO: Implement me.
# TODO: 3.3 Implement me.
return np.zeros_like(x)


Expand Down Expand Up @@ -66,6 +66,8 @@ def fit_gmm(points: np.ndarray, init_params_list: List) -> List:

if __name__ == "__main__":
np.random.seed(42)
#--------------------------------------------------------------------------------------#

dist1 = np.random.normal(loc=(2, 2), scale=(1.0, 1.0), size=(100, 2))
dist2 = np.random.normal(loc=(-2, -2), scale=(1.0, 1.0), size=(100, 2))

Expand Down
5 changes: 4 additions & 1 deletion src/plot_gaussian.py
Original file line number Diff line number Diff line change
Expand Up @@ -45,12 +45,15 @@ def forward_euler(x: np.ndarray, fun: Callable, int_0: float = 0.0) -> np.ndarra


if __name__ == "__main__":
params = (
# 3.2 Explore the behavior of tha gaussian pdf by trying out different parameter values
params = ( # each parameter-tuple contains a mean and a std: (mu, std)
(0.0, np.sqrt(0.2)),
(0.0, np.sqrt(1.0)),
(0.0, np.sqrt(5.0)),
(-2.0, np.sqrt(0.5)),
)

# Plot the gaussian pdfs for the different parameters specified above
for param in params:
x = np.linspace(-5.0, 5.0, 500)
plt.plot(
Expand Down
19 changes: 14 additions & 5 deletions src/sample_mean_corr_rhein.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,13 +10,13 @@

def my_mean(data_sample) -> float:
"""Implement a function to find the mean of an input List."""
# TODO: Implement me.
# TODO: 1.1 Implement me.
return 0.


def my_std(data_sample) -> float:
"""Implement a function to find the standard deviation of a sample in a List."""
# TODO: Implement me.
# TODO: 1.2 Implement me.
return 0.


Expand All @@ -29,7 +29,8 @@ def auto_corr(x: np.ndarray) -> np.ndarray:
Returns:
np.ndarray: Autocorrelation of input signal of shape (signal_length*2 - 1,).
"""
# TODO: Implement me.
# TODO: 2.1 Implement me.
# TODO: 2.2 Check your implementation via nox -s test.
return np.zeros_like(x)


Expand Down Expand Up @@ -57,6 +58,14 @@ def auto_corr(x: np.ndarray) -> np.ndarray:
]


# TODO: compute the mean and standard deviation before and after 2000.
# TODO: 1.3 Compute the mean and standard deviation before 2000.
# TODO: 1.4 Compute the mean and standatd deviation after 2000.

# TODO: Compare the autocorrelation functions of the rhine data and of a random signal.
#----------------------------------------------------------------------------------------------#
# TODO: 2.3 Normalize the data of the Rhine level measurements since 2000.
# TODO: 2.4 Compute and plot the autocorrelation.

# TODO: 2.5 Create a random signal.
# TODO: 2.6 Normalize the random signal.
# TODO: 2.7 Compute the autocorrelation.
# TODO: 2.8. Plot both autocorrelations and compare the results.
Loading