# Music Machine Learning - Probability distributions

### Author: Philippe Esling (esling@ircam.fr)

In this course we will cover
1. A [quick recap](#recap) on simple probability concepts
2. An introduction to [probability distributions](#distributions)
3. An explanation on how to [sample](#sampling) from distributions ourselves

<a id="recap"> </a>

## Quick recap on probability

The field of probability aims to model random or uncertain events through random variables $X$ denoting a quantity that is uncertain, which take values $\omega$ in a sample space $\omega \in \Omega$. If we observe several occurrences of this variable $\{\mathbf{x}_{i}\}_{i=1}$, we might try to model the _probability distribution_ $p(\mathbf{x})$ of that variable. We recall here that the probability of an event $a$ is a real number $p(a)$, with $0 \leq p(a) \leq 1$, knowing that $p(\Omega)=1$ and $p(\{\})=0$. The probability of two events occuring simultaneously is defined as $p\left(a \cup b \right)$. Therefore, the probability of one event **or** the other occuring is defined as
$
\begin{equation}
p\left(a \cap b \right) = p(a) + p(b) - p\left(a \cup b \right)
\end{equation}
$

The *conditional probability* of an event $a$ occuring *given* another event $b$ is denoted $p \left(a \mid b \right)$ and is defined as
$$
\begin{equation}
p \left(a \mid b \right) = \frac{p \left(a , b \right)}{p \left(b \right)}
\end{equation}
$$

This can be understood as the probability of event $a$ to occur if we restrict the world of possibilities to event $b$. The *chain rule* defines the probabilities of a set of events to co-occur simultaneously
$$
\begin{equation}
p \left(x_{1},...,x_{n} \right) = \prod_{i=n}^{1}{p \left(x_{i}\mid x_{i-1},..., x_{1} \right)}
\end{equation}
$$

Finally, we say that two events are independent if $p(a\mid b) = p(a)$.

To understand these concepts graphically, we will rely on `PyTorch` and specifically the `distributions` package.


```python
import torch.distributions as distribution
import torch.distributions.transforms as transform
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from helper_plot import hdr_plot_style
hdr_plot_style()
```

<a id="distributions"> </a>
### Probability distributions

Let $X$ be a random variable associated with a *probability distribution function* that assigns probabilities to the different outcomes $X$ can take in the sample space. We can divide random variables into three different types:

- **$X$ is discrete**: Discrete random variables may only assume values on a specified list. 
- **$X$ is continuous**: Continuous random variable can take on arbitrarily exact values. 
- **$X$ is mixed**: Mixed random variables assign probabilities to both discrete and continuous random variables, (i.e. a combination of the above two categories). 

**Expected Value**
The expected value $\mathbb{E}\left[X\right]$ for a given probability distribution can be described as _"the mean expected value for many repeated samples from that distribution."_ As the number of repeated observation goes to infinity, the difference between the average outcome and the expected value becomes arbitrarily small. Formally, it is defined for discrete variables as

$$
\mathbb{E}\left[X\right] = \sum\limits_{i}\bx_i p(\bx=\bx_i)
$$

#### Discrete distributions

If $X$ is discrete, then its distribution is called a *probability mass function* (pmf), which measures the probability $X$ takes on the value $x_{i}$, denoted $P(X=x_{i})$. Let $\mathbf{x}$ be a discrete random variable with range $R_{X}=\{x_1,\cdots,x_n\}$ (finite or countably infinite). The function

$$
p_{X}(x_{i})=p(X=x_{i}), \forall i\in\{1,\cdots,n\}
$$

is called the probability mass function (PMF) of $X$.

Hence, the PMF defines the probabilities of all possible values for a random variable. The above notation allows to express that the PMF is defined for the random variable $X$, so that $p_{X}(1)$ gives the probability that $X=1$. For discrete random variables, the PMF is also called the \textit{probability distribution}. The PMF is a probability measure, therefore it satisfies all the corresponding properties
- $0 \leq p_{X}(x_i) < 1, \forall x_i$
- $\sum_{x_i\in R_{X}} p_{X}(x_i) = 1$
- $\forall A \subset R_{X}, p(X \in A)=\sum_{x_a \in A}p_{X}(x_a)$

A very simple example of discrete distribution is the `Bernoulli` distribution. With this distribution, _we can model a coin flip_, if it has equal probability. More formally, a Bernoulli distribution is defined as

$$ 
Bernoulli(x)= p^x (1-p)^{(1-x)} 
$$

with $p$ controlling the probability of the two classes. Hence, a fair coin should have $p=0.5$, and if we throw the coin a very large number of times, we hope to see on average an equal amount of _heads_ and _tails_.


```python
bernoulli = distribution.Bernoulli(0.5)
samples = bernoulli.sample((10000,))
plt.figure(figsize=(10,8))
sns.distplot(samples)
plt.title("Samples from a Bernoulli (coin toss)")
plt.show()
```

However, we can also _sample_ from the distribution to have individual values of a single throw. In that case, we obtain a series of separate events that _follow_ the distribution 


```python
vals = ['heads', 'tails']
samples = bernoulli.sample((10, ))
for s in samples:
    print('Coin is tossed on ' + vals[int(s)])
```

Now, we can mess up our probability to model an unfair (loaded) coin, as shown in the following example (where we use a cheated coin that should give us a lot more of _heads_ than _tails_) 


```python
bernoulli = distribution.Bernoulli(0.3)
samples = bernoulli.sample((10000,))
plt.figure(figsize=(10, 8)); sns.distplot(samples)
plt.title("Samples from an unfair coin"); plt.show()
vals = ['heads', 'tails']
samples = bernoulli.sample((10, ))
for s in samples:
    print('Coin is tossed on ' + vals[int(s)])
```

#### Poisson distribution

Let's introduce one of the (many) useful probability mass functions. We say $Z$ is *Poisson*-distributed if:

$$
\begin{equation}
P\left(Z = k\right) =\frac{ \lambda^k e^{-\lambda} }{k!}, \; \; k \in \mathbb{N^{+}} 
\end{equation}
$$

$\lambda \in \mathbb{R}$ is a parameter of the distribution that controls its shape (usually termed the *intensity* of the Poisson distribution). By increasing $\lambda$, we add more probability to larger values. One can describe $\lambda$ as the *intensity* of the Poisson distribution. If a random variable $Z$ has a Poisson mass distribution, we denote it by

$$
\begin{equation}
Z \sim \text{Poi}(\lambda) 
\end{equation}
$$

One useful property of the Poisson distribution is that its expected value is equal to its parameter, i.e.:

$$E\large[ \;Z\; | \; \lambda \;\large] = \lambda $$

We will plot the probability mass distribution for different $\lambda$ values.

### Continuous distributions

The same ideas apply to _continuous_ random variables, which can model for instance the height of human beings. If we try to guess the height of someone that we do not know, there is a higher probability that this person will be around 1m70, instead of 20cm or 3m. For the rest of this course, we will use the shorthand notation $p(\mathbf{x})$ for the distribution $p(\mathbf{x}=x_{i})$,  which expresses for a real-valued random variable $\mathbf{x}$, evaluated at $x_{i}$, the probability that $\mathbf{x}$ takes the value $x_i$.

One notorious example of such distributions is the Gaussian (or Normal) distribution, which is defined as 
\begin{equation}
    p(x)=\mathcal{N}(\mu,\sigma)=\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}
\end{equation}

Similarly as before, we can observe the behavior of this distribution with the following code (in our height example)


```python
normal = distribution.Normal(150., 30.)
samples = normal.sample((10000, ))
sns.distplot(samples)
plt.title("Samples from a standard Normal")
plt.show()
```

If we have access to this complete probability distribution (its exact parameterization and function), we can generate samples (in this case "new humans") that follow the correct distribution. You can experiment with this in the following code space.


```python

######################
# YOUR CODE GOES HERE
######################
    
```

An example of continuous random variable is a random variable with *exponential density*

$$
\begin{equation}
f_Z(z | \lambda) = \lambda e^{-\lambda z }, \;\; z\ge 0
\end{equation}
$$

When a random variable $Z$ has an exponential distribution with parameter $\lambda$, we say *$Z$ is exponential*

$$Z \sim \text{Exp}(\lambda)$$

Given a specific $\lambda$, the expected value of an exponential random variable is equal to the inverse of $\lambda$, that is

$$E[\; Z \;|\; \lambda \;] = \frac{1}{\lambda}$$ 

<a id="distribs"></a>
### PyTorch distributions

Here, we rely on the [PyTorch distributions module](https://pytorch.org/docs/stable/_modules/torch/distributions/), which is defined in `torch.distributions`. Most notably, we are going to rely both on the `Distribution` and `Transform` objects.


```python
# Imports for plotting
import numpy as np
import matplotlib.pyplot as plt
# Define grids of points (for later plots)
x = np.linspace(-4, 4, 1000)
z = np.array(np.meshgrid(x, x)).transpose(1, 2, 0)
z = np.reshape(z, [z.shape[0] * z.shape[1], -1])
```

Inside this toolbox, we can already find some of the major probability distributions that we are used to deal with



```python
p = distrib.Normal(loc=0, scale=1)
p = distrib.Bernoulli(probs=torch.tensor([0.5]))
p = distrib.Beta(concentration1=torch.tensor([0.5]), concentration0=torch.tensor([0.5]))
p = distrib.Gamma(concentration=torch.tensor([1.0]), rate=torch.tensor([1.0]))
p = distrib.Pareto(alpha=torch.tensor([1.0]), scale=torch.tensor([1.0]))
```

The interesting aspect of these `Distribution` objects is that we can both obtain some samples from it through the `sample` (or `sample_n`) function, but we can also obtain the analytical density at any given point through the `log_prob` function


```python
# Based on a normal
n = distrib.Normal(0, 1)
# Obtain some samples
samples = n.sample((1000, ))
# Evaluate true density at given points
density = torch.exp(n.log_prob(torch.Tensor(x))).numpy()
# Plot both samples and density
fig, (ax1, ax2) = plt.subplots(1, 2, sharex=True, figsize=(15, 4))
ax1.hist(samples, 50, alpha=0.8);
ax1.set_title('Empirical samples', fontsize=18);
ax2.plot(x, density); ax2.fill_between(x, density, 0, alpha=0.5)
ax2.set_title('True density', fontsize=18);
```

Here you can experiment with different distributions, and try to compare how they behave depending on their parameters and also on how much _samples_ you draw from these.


```python

######################
# YOUR CODE GOES HERE
######################
    
```

<a id="sampling"></a>
## Sampling from distributions

The advantage of using probability distributions is that we can *sample* from these to obtain examples that follow the distribution. For instance, if we perform sampling repeatedly (up to infinity) from a Gaussian PDF, the different values will be distributed following the exact Gaussian distribution. However, although we know the PDF, we need to compute the *Cumulative Distribution Function* (CDF), and then its inverse to obtain the sampling function. Therefore, if we denote the PDF as $f_{X}(x)$, we need to compute the CDF 

$$
\begin{equation}
F_{X}\left(x\right)=\intop_{\infty}^{x}f_{X}\left(t\right)dt
\end{equation}
$$

Then, the *inverse sampling method* consists in solving and applying the inverse CDF $F_{X}^{-1}\left(x\right)$. Here, we recall that the Exponential probability is defined with the following function.

$$
p_{\lambda}(x) = \lambda e^{-\lambda x}
$$

with $\lambda$ defining the _rate_ parameter. Therefore, to be able to define our own `sample` method, we need to solve for

$$
F^{-1}_Y(x) = \left(\int_0^x \lambda e^{-\lambda y} dy\right)^{-1}
$$

***

**Exercise**

1. Code the exponential probability density functions 
2. Perform the inverse transform method on the exponential distribution PDF
3. Code the `sample_exponential` function to sample from the exponential
4. (optional) Perform the same operations for sampling from the Beta distribution

***


```python
from scipy.stats import expon
nb_samples = 500
nb_bins = 50

def sample_exponential(mu, n, m):
    
    ######################
    # YOUR CODE GOES HERE
    ######################
    
    ######################
    # Solution

    samples = np.random.uniform(0, 1, size=(m,n))
    samples = - np.log(1 - samples) 
    samples *= mu
    
    ######################
    
    return samples

# Exponential distribution
mu = 2
samples = np.random.exponential(mu, nb_samples)
samples_ex = sample_exponential(mu, nb_samples, 1)
# Compute the PDF
X = np.linspace(0, np.max(samples), int(np.max(samples)) * 100)
y1 = expon.pdf(X) #* (nb_samples / nb_bins) * int(np.max(samples) * 1.5)
# Display both
plt.figure(figsize=(12, 8))
plt.subplot(2,1,1)
plt.hist(samples, 50, label='Samples')
plt.plot(X,y1,ls='--',c='r',linewidth=2, label='Exponential PDF')
plt.legend(loc=1)
plt.subplot(2,1,2)
plt.hist(samples_ex, 50, label='Samples')
plt.plot(X,y1,ls='--',c='r',linewidth=2, label='Our exponential sampples')
plt.legend(loc=1)
```


```python
def sampleBeta(a, b, M, N):
    
    ######################
    # YOUR CODE GOES HERE
    ######################
    
    ######################
    # Solution

    samples = np.random.beta(a, b, size=(M, N))
    
    ######################
    
    return samples

from scipy.stats import beta
# Beta distribution
a = 0.6
b = 1.5
samples = np.random.beta(a, b, nbSamples)
samplesBeta = sampleBeta(a, b, nbSamples, 1)
# Compute the PDF
X = np.linspace(0, 1, 100)
y1 = beta.pdf(X, a, b) * (nbSamples / nbBins) * (np.max(samples) * 1.5)
# Display both
plt.figure(figsize=(12, 8))
plt.subplot(2,1,1)
plt.hist(samples, 50, label='Samples')
plt.plot(X,y1,ls='--',c='r',linewidth=2, label='Beta PDF')
plt.legend(loc=1)
plt.subplot(2,1,2)
plt.hist(samplesBeta, 50, label='Samples')
plt.plot(X,y1,ls='--',c='r',linewidth=2, label='Beta PDF')
plt.legend(loc=1)
```

### Comparing distributions (KL divergence)
$
\newcommand{\R}{\mathbb{R}}
\newcommand{\bb}[1]{\mathbf{#1}}
\newcommand{\bx}{\bb{x}}
\newcommand{\by}{\bb{y}}
\newcommand{\bz}{\bb{z}}
\newcommand{\KL}[2]{\mathcal{D}_{\text{KL}}\left[#1 \| #2\right]}$
Originally defined in the field of information theory, the _Kullback-Leibler (KL) divergence_ (usually noted $\KL{p(\bx)}{q(\bx)}$) is a dissimilarity measure between two probability distributions $p(\bx)$ and $q(\bx)$. In the view of information theory, it can be understood as the cost in number of bits necessary for coding samples from $p(\bx)$ by using a code optimized for $q(\bx)$ rather than the code optimized for $p(\bx)$. In the view of probability theory, it represents the amount of information lost when we use $q(\bx)$ to approximate the true distribution $p(\bx)$. %that explicit the cost incurred if events were generated by $p(\bx)$ but charged under $q(\bx)$


Given two probability distributions $p(\bx)$ and $q(\bx)$, the Kullback-Leibler divergence of $q(\bx)$ _from_ $p(\bx)$ is defined to be
\begin{equation}
    \KL{p(\bx)}{q(\bx)}=\int_{\R} p(\bx) \log \frac{p(\bx)}{q(\bx)}d\bx
\end{equation}

Note that this dissimilarity measure is _asymmetric_, therefore, we have
\begin{equation}
    \KL{p(\bx)}{q(\bx)}\neq \KL{q(\bx)}{p(\bx)}
\end{equation}
This asymmetry also describes an interesting behavior of the KL divergence, depending on the order to which it is evaluated. The KL divergence can either be a _mode-seeking_ or _mode-coverage_ measure (we will come back to these notions in the _approximate inference_ course)


```python

```

8cd934b6febdda2186cdda4118dea1246ccf472a
ipynb
Jupyter Notebook
04b_distributions.ipynb
piptouque/atiam_ml
9da637eae179237d30a15dd9ce3e95a2a956c385
MIT
04b_distributions.ipynb
piptouque/atiam_ml
9da637eae179237d30a15dd9ce3e95a2a956c385
MIT
04b_distributions.ipynb
piptouque/atiam_ml
9da637eae179237d30a15dd9ce3e95a2a956c385
MIT
Qwen/Qwen-72B
1. YES
2. YES