# Integrating Over Parameters, Tractability, MLE, MAP, Bayesian, and Probability Theory

## Probability Theory Basics

**Rule I**: $\displaystyle p(a)=\int p(a,b)\,db$

**Rule II**: $\displaystyle p(a\mid b)=\frac{p(a,b)}{p(b)}$

**Bayes Rule**: $\displaystyle p(a\mid b)=\frac{p(b\mid a)\cdot p(a)}{p(b)}$

## Maximum Likelihood Estimation (MLE)

Method for estimating the parameter(s) of a distribution (?).

Given is a set of points $\mathcal{D}=\{y_n\}_{n=1}^N$.
We _choose to model_ the points with a normal distribution $\mathcal{N}(\mu,\sigma)$ so we can assess the probability of a new point $y$ as $p(y;\mu,\sigma)\sim\mathcal{N}(\mu,\sigma)$.
The goal is to be able to compute $p(y;\mu,\sigma)$ for any $y$.
To be able to do that we need to estimate $\mu$ and $\sigma$ and that is what we use MLE for.

The data is all information we have, so we want $p(\mathcal{D};\mu,\sigma)$ to be high.
We assume the data is i.i.d. given the parameters $\mu$ and $\sigma$.
We can therefore write

$$p(\mathcal{D};\mu,\sigma)=\prod_{n=1}^Np(y_n;\mu,\sigma)$$

What we actually look for are the parameters.
We are interested in the particular set of parameters for which $p(\mathcal{D};\mu,\sigma)$ is maximized.
Let $\mu^\star$ and $\sigma^\star$ denote this maximizer.
We search for
\begin{align}
\mu^\star=&\arg\max_\mu p(\mathcal{D};\mu,\sigma)\\
=&\arg\max_\mu\prod_{n=1}^Np(y_n;\mu,\sigma)\\
=&\arg\max_\mu\log\prod_{n=1}^Np(y_n;\mu,\sigma)\\
=&\arg\max_\mu\sum_{n=1}^N\log p(y_n;\mu,\sigma)
\end{align}

We can apply a $\log$ to the product because it will not change the maximizer.

The same formulation can be written down for $\sigma^\star$ to find the optimal parameter here.

In order to _actually_ get a value for $\mu^\star$ we replace $p(y_n;\mu,\sigma)$ with the definition of our model.
Above it is a normal distribution, so it would be TODO.

Depending on the model there may be a closed-form solution for computing maximizer (i.e., the parameterization for which the data is the most likely) or one can restort to gradient descent to find it.

TODO: apply this to linear regression for example. How can we get $p(y;w)$ there?

## Maximum Aposteriori (MAP)

Method for estimating the parameter(s) of a distribution (?).
Just like MLE, except the probabilistic approach is a bit different.

Previously, with MLE, we modeled $p(\mathcal{D};\mu,\sigma)$, that is, the probability of the data given the parameters.
We then looked for the parameters for which this probability is the highest.
Now we model the probability $p(\mu\mid\mathcal{D})$, that is, the probability of the parameters given the data.
Note that $\sigma$ is also a parameter, it's left out for readability.

Again, analogous to MLE, we seek for the argument that maximizes the probability.
Mathematically speaking, that is

$$\mu^\star=\arg\max_\mu p(\mu\mid\mathcal{D})\,.$$

Following Bayes rule, 

$$p(\mu\mid\mathcal{D})=\frac{p(\mathcal{D}\mid\mu)\cdot p(\mu)}{p(\mathcal{D})}\,.$$

We can plug the rewritten version of $p(\mu\mid\mathcal{D})$ into the optimization objective from above:

\begin{align}
\mu^\star=&\arg\max_\mu p(\mu\mid\mathcal{D})\\
=&\arg\max_\mu\frac{p(\mathcal{D}\mid\mu)\cdot p(\mu)}{p(\mathcal{D})}\\
=&\arg\max_\mu p(\mathcal{D}\mid\mu)\cdot p(\mu)
\end{align}

Note that the denominator $p(\mathcal{D})$ can be dropped because it does not depend on $\mu$.

Just like we did with MLE, we first expand $p(\mathcal{D}\mid\mu)$ into a product over all data samples available to us, and second, apply a $\log$ to ease the optimization, without affecting the position of the maximum (the best parameter configuration).

\begin{align}
\mu^\star=&\arg\max_\mu p(\mathcal{D}\mid\mu)\cdot p(\mu)\\
=&\arg\max_\mu\prod_{n=1}^Np(y_n\mid\mu)\cdot p(\mu)\\
=&\arg\max_\mu\log\left(\prod_{n=1}^N p(y_n\mid\mu)\cdot p(\mu)\right)\\
=&\arg\max_\mu\sum_{n=1}^N\log p(y_n\mid\mu)+\log p(\mu)\\
\end{align}

This is very similar to MLE, with the only difference that there is a new component in the optimization problem that we need to consider.
It is the probability of the parameters, $p(\mu)$, also known as the **prior**.

This is where we have to _explicitly_ make an assumption, for example, by saying our parameter is normally distributed. This is in contrast to MLE where this assumption was made implicitly. (?)

TODO: plug in normal distribution definition and show how it resolves to MSE loss with L2 regularization.

## Bayesian Modeling

Instead of just estimating a single value per parameter (like above just a scalar for the mean of our normal distribution), we now **model the parameters as random variables** themselves.

What was previously $\mu$ is now $\mu$ but stands for a random variable, which is defined by a parameterized distribution we choose.

**Terminology**

* Parameters: $\mu$
* Data: $\mathcal{D}$
* Prior: $p(\mu)$
* Posterior: generally _after having seen the data_, i.e., $p(\cdot\mid\mathcal{D})$, for example the updated parameters $p(\mu\mid\mathcal{D})$
* Likelihood: model of the data given the parameters, i.e., $p(y\mid\mu,\mathcal{D})=p(y\mid\mu)$
* Posterior predictive distribution: $p(y\mid\mathcal{D})$
* Evidence: probability of the data $p(\mathcal{D})$
* Marginal: simplified form in which a random variable was integrated away, e.g., $p(y)$ where the $\mu$ was integrated over
* Marginal likelihood: ?

**Posterior** is always after having seen the data.

What we are interested in is modeling the probability of a new (test) data point $y$ given all the information we have, namely the data $\mathcal{D}$. It is called the **posterior predictive distribution**.

\begin{align}
p(y\mid\mathcal{D})=&\int p(y\mid\mu,\mathcal{D})\cdot p(\mu\mid\mathcal{D})\,d\mu\\
=&\int p(y\mid\mu)\cdot p(\mu\mid\mathcal{D})\,d\mu
\end{align}

We constructed the integral using Rule I.
The choice of $\mu$ is an arbitrary one, one could have integrated over anything, but it is a _helpful_ one. In the second equality we remove $\mathcal{D}$ because the data is i.i.d. given the parameter. Having the parameter plus the data does not provide additional information.

For example: We know the mean of a normal distribution. Gathering additional samples from it will not provide us with additional information.

There are not two components in the integral.
The first is the **likelihood** $p(y\mid\mu,\mathcal{D})$ and something we can compute, because it is our model. It is *easy* because it was our model choice, how to model $y$s.

For example: $p(y\mid\mu)=\mathcal{N}(\mu,\sigma=1)$; $\sigma=1$ for brevity.

As for the second part, the so called **posterior distribution of the parameters** $p(\mu\mid\mathcal{D})$, we can again apply Bayes rule as done in MAP above:

$$p(\mu\mid\mathcal{D})=\frac{p(\mathcal{D}\mid\mu)\cdot p(\mu)}{p(\mathcal{D})}\,.$$

This time around we cannot just discard $p(\mathcal{D})$, because we are not looking for the $\arg\max_\mu$ but the actual probability $p(\mu\mid\mathcal{D})$.

There are two difficulties: One is the numerator, one the denominator.

**Numerator**. Because the data is i.i.d. given the parameters we can break it up into a product of the likelihoods of the data samples. Likelihood, because the part that we condition on is unknown, not the one in the front (where it would be a probability).

$$p(\mathcal{D}\mid\mu)\cdot p(\mu)=\prod_{i=1}^Np(y_i\mid\mu)\cdot p(\mu)$$

The two distributions which we multiply with each other can easily (once multiplied) into a new distribution which is hard to handle / unknown. In some cases, though, for example, for two normal distributions we know what the product is (also normal).

More generally, if the two distributions (likelihood and prior) are **conjugate** to each other, the resulting posterior has the same as the distribution as the prior. Examples are: 
* normal and normal -> normal,
* bernoulli and beta -> beta,
* multinoulli and dirichlet --> dirichlet,
* ...

Example.

Our prior is $$p(\mu)=\operatorname{Beta}(1,1)=\frac{\mu^{\alpha-1}(1-\mu)^{\beta-1}}{B(\alpha,\beta)}$$

with ${\displaystyle \mathrm {B} (\alpha ,\beta )={\frac {\Gamma (\alpha )\Gamma (\beta )}{\Gamma (\alpha +\beta )}}}$.

The likelihood is $$p(y\mid\mu)=\operatorname{Ber}(\mu)=\mu^{y}(1-\mu)^{1-y}$$

This gives the numerator

$$\mu^{y}(1-\mu)^{1-y}\cdot\frac{\mu^{\alpha-1}(1-\mu)^{\beta-1}}{B(\alpha,\beta)}\,.$$

For a data point $y=1$ it is

\begin{align}
&\mu\cdot\frac{\mu^{\alpha-1}(1-\mu)^{\beta-1}}{B(\alpha,\beta)}\\
=&\frac{1}{B(\alpha,\beta)}\cdot\mu\cdot\mu^{\alpha-1}(1-\mu)^{\beta-1}\\
=&\frac{1}{B(\alpha,\beta)}\cdot\mu^{\underbrace{(\alpha+1)}_{\text{new }\alpha}-1}(1-\mu)^{\beta-1}\\
\end{align}

It is apparent how the multiplication with Bernoulli did not change the distribution type.
We still have a Beta distribution and its parameter $\alpha$ has changed slightly to turn into $\alpha+1$.
Had our data point been $y=0$, would $\alpha$ have stayed the same whereas $\beta$ would have been incremented.

**Denominator**. $p(\mathcal{D})$ is just a constant. So the part that actually determines the _distribution part_ of $p(\mu\mid\mathcal{D})$ is the numerator.

In the following, $\mathcal{D}=\{y\}$, i.e., there is only a single data point. Taking the numerator for the example from above and plugging it back into the fraction we get

\begin{align}
p(\mu\mid\mathcal{D})=&\frac{\mu^{y}(1-\mu)^{1-y}\cdot\frac{\mu^{\alpha-1}(1-\mu)^{\beta-1}}{B(\alpha,\beta)}}{p(\mathcal{D})}\\
=&\underbrace{\frac{1}{p(\mathcal{D})\cdot B(\alpha,\beta)}}_{\text{numbers}}\cdot\underbrace{\mu^{y}(1-\mu)^{1-y}\cdot\mu^{\alpha-1}(1-\mu)^{\beta-1}}_{\text{depends on }\mu}\\
=&\underbrace{\frac{1}{p(\mathcal{D})\cdot B(\alpha,\beta)}}_{\text{numbers}}\cdot\underbrace{\mu^{(\alpha+y)-1}(1-\mu)^{(\beta+1-y)-1}}_{\text{depends on }\mu}
\end{align}

The additional information comes from the fact that the distribution has to integrate to 1. That helps us determine the unknown component $p(\mathcal{D})$.

\begin{align}
\int p(\mu\mid\mathcal{D})\,d\mu=&1\\
\int \frac{1}{p(\mathcal{D})\cdot B(\alpha,\beta)}\cdot\mu^{(\alpha+y)-1}(1-\mu)^{(\beta+1-y)-1}\,d\mu=&1\\
\frac{1}{p(\mathcal{D})\cdot B(\alpha,\beta)}\cdot\int\mu^{(\alpha+y)-1}(1-\mu)^{(\beta+1-y)-1}\,d\mu=&1\\
\int\mu^{(\alpha+y)-1}(1-\mu)^{(\beta+1-y)-1}\,d\mu=&\underbrace{p(\mathcal{D})}_\text{unknown}\cdot B(\alpha,\beta)\\
\int\mu^{(\alpha+y)-1}(1-\mu)^{(\beta+1-y)-1}\,d\mu=&B(\alpha+y,\beta+1-y)
\end{align}

The last equality holds because we know that the left side is an unnormalized $\operatorname{Beta}(\alpha+y,\beta+1-y)$ distribution PDF. So we define the right-hand side to be its corresponding normalization constant, which is $B(\alpha+y,\beta+1-y)$. And $p(\mathcal{D})$ is just a part of it.

$$p(\mathcal{D})=\frac{B(\alpha+y,\beta+1-y)}{B(\alpha,\beta)}$$

---

Now we jump back to the posterior predictive distribution, where we aim to determine the probability of a test point $y^*$ given the data. The above paragraphs helped determine the $p(\mu\mid\mathcal{D})$ part of it, and the likelihood is given by $p(y\mid\mu)=\operatorname{Ber}(\mu)=\mu^{y}(1-\mu)^{1-y}$

\begin{align}
p(y^*\mid\mathcal{D})=&\int p(y^*\mid\mu)\cdot p(\mu\mid\mathcal{D})\,d\mu=\mathbb{E}_{p(\mu\mid\mathcal{D})}\left[p(y^*\mid\mu)\right]\\
=&\int p(y^*\mid\mu)\cdot\frac{1}{p(\mathcal{D})\cdot B(\alpha,\beta)}\cdot\mu^{(\alpha+y)-1}(1-\mu)^{(\beta+1-y)-1}\,d\mu\\
=&\int p(y^*\mid\mu)\cdot\frac{1}{B(\alpha+y,\beta+1-y)}\cdot\mu^{(\alpha+y)-1}(1-\mu)^{(\beta+1-y)-1}\,d\mu\\
=&\int \mu^{y^*}(1-\mu)^{1-y^*}\cdot\frac{1}{B(\alpha+y,\beta+1-y)}\cdot\mu^{(\alpha+y)-1}(1-\mu)^{(\beta+1-y)-1}\,d\mu\\
=&\frac{1}{B(\alpha+y,\beta+1-y)}\cdot\int \mu^{y^*}(1-\mu)^{1-y^*}\cdot\mu^{(\alpha+y)-1}(1-\mu)^{(\beta+1-y)-1}\,d\mu\\
=&\frac{1}{B(\alpha+y,\beta+1-y)}\cdot\int \mu^{y^*}(1-\mu)^{1-y^*}\cdot\mu^{(\alpha+y)-1}(1-\mu)^{(\beta+1-y)-1}\,d\mu\\
=&\frac{1}{B(\alpha+y,\beta+1-y)}\cdot\int \mu^{(\alpha+y+y^*)-1}(1-\mu)^{(\beta+1-y+1-y^*)-1}\,d\mu\\
=&\frac{1}{B(\alpha+y,\beta+1-y)}\cdot B(\alpha+y+y^*,\beta+1-y+1-y^*)\\
=&\frac{B(\alpha+y+y^*,\beta+1-y+1-y^*)}{B(\alpha+y,\beta+1-y)}\\
\end{align}

\begin{align}
p(y^*\mid\mathcal{D})=&\begin{cases}
\underbrace{\frac{B(\alpha+y,\beta-y+2)}{B(\alpha+y,\beta+1-y)}}_{1-\mu?} & \text{if }y^*=0\\
\underbrace{\frac{B(\alpha+y+1,\beta+1-y)}{B(\alpha+y,\beta+1-y)}}_{\mu?} & \text{if }y^*=1
\end{cases}\\
=&\operatorname{Ber}\left(\frac{B(\alpha+y+1,\beta+1-y)}{B(\alpha+y,\beta+1-y)}\right)
\end{align}

Note: Here $y$ is the data that we learned from (a single point).

Here we would have almost given up. What we hope here is (as **successfully** validated in the following) that 

\begin{align}
\frac{B(\alpha+y,\beta-y+2)}{B(\alpha+y,\beta-y+1)}+\frac{B(\alpha+y+1,\beta-y+1)}{B(\alpha+y,\beta-y+1)}\overset{?}{=}&1\\\\
\frac{B(\alpha+y,\beta-y+2)+B(\alpha+y+1,\beta-y+1)}{B(\alpha+y,\beta-y+1)}\overset{?}{=}&1\\
\end{align}

With the definition of $B$ 

$${\displaystyle \mathrm {B} (\alpha ,\beta )={\frac {\Gamma (\alpha )\Gamma (\beta )}{\Gamma (\alpha +\beta )}}}$$

we get the fraction

\begin{align}
\frac{\frac{\Gamma(\alpha+y)\Gamma(\beta-y+2)}{\Gamma(\alpha+y+\beta-y+2)}+\frac{\Gamma(\alpha+y+1)\Gamma(\beta-y+1)}{\Gamma(\alpha+y+1+\beta-y+1)}}{\frac{\Gamma(\alpha+y)\Gamma(\beta-y+1)}{\Gamma(\alpha+y+\beta-y+1)}}\overset{?}{=}&1\\\\
\frac{\frac{\Gamma(\alpha+y)\Gamma(\beta-y+2)}{\Gamma(\alpha+\beta+2)}+\frac{\Gamma(\alpha+y+1)\Gamma(\beta-y+1)}{\Gamma(\alpha+\beta+2)}}{\frac{\Gamma(\alpha+y)\Gamma(\beta-y+1)}{\Gamma(\alpha+\beta+1)}}\overset{?}{=}&1\\\\
\frac{\frac{\Gamma(\alpha+y)\Gamma(\beta-y+2)+\Gamma(\alpha+y+1)\Gamma(\beta-y+1)}{\Gamma(\alpha+\beta+2)}}{\frac{\Gamma(\alpha+y)\Gamma(\beta-y+1)}{\Gamma(\alpha+\beta+1)}}\overset{?}{=}&1\\
\end{align}

With the definition of ${\displaystyle \Gamma (n)=(n-1)!}$ we continue to pursue our goals at 00:43 am...

$$\Gamma (n+1)=\Gamma(n)\cdot n=n!$$

\begin{align}
\frac{\frac{\Gamma(\alpha+y)\Gamma(\beta-y+2)+\Gamma(\alpha+y+1)\Gamma(\beta-y+1)}{\Gamma(\alpha+\beta+2)}}{\frac{\Gamma(\alpha+y)\Gamma(\beta-y+1)}{\Gamma(\alpha+\beta+1)}}\overset{?}{=}&1\\\\
\frac{
    (\Gamma(\alpha+y)
    \Gamma(\beta-y+2)+\Gamma(\alpha+y+1)\Gamma(\beta-y+1))\cdot\Gamma(\alpha+\beta+1)}
    {\Gamma(\alpha+\beta+2)\cdot\Gamma(\alpha+y)\cdot\Gamma(\beta-y+1)}\overset{?}{=}&1\\\\
\frac{
    (\Gamma(\alpha+y)
    \Gamma(\beta-y+2)+\Gamma(\alpha+y+1)\Gamma(\beta-y+1))}
    {(\alpha+\beta+1)\cdot\Gamma(\alpha+y)\cdot\Gamma(\beta-y+1)}\overset{?}{=}&1\\\\
\frac{
    (
    \Gamma(\beta-y+2)+(\alpha+y)\Gamma(\beta-y+1))}
    {(\alpha+\beta+1)\cdot\Gamma(\beta-y+1)}\overset{?}{=}&1\\\\
\frac{\beta-y+1+\alpha+y}{\alpha+\beta+1}\overset{?}{=}&1\\\\
\frac{\beta+1+\alpha}{\alpha+\beta+1}\overset{!}{=}&1\\\\
\end{align}

---

What we did above is a special case. It's an example for a conjugate prior (bernoulli likelihood and beta prior) in which a closed-form solution was available. In cases where there is no closed-form solution, one can solve

$$p(y^*\mid\mathcal{D})=\int p(y^*\mid\mu)\cdot p(\mu\mid\mathcal{D})\,d\mu=\mathbb{E}_{p(\mu\mid\mathcal{D})}\left[p(y^*\mid\mu)\right]$$

by sampling from $p(\mu\mid\mathcal{D})$ to approximate the expectation (**Monte Carlo**). So here we would approximate $p(y^*\mid\mu)$.

Alternatively, one can (probably?) approximate $p(\mathcal{D}\mid\mu)$ in

$$p(\mu\mid\mathcal{D})=\frac{p(\mathcal{D}\mid\mu)\cdot p(\mu)}{p(\mathcal{D})}\,.$$

**Variational inference** (probably?) fits into this picture where one says we define $q$ to approximate $p(y^*\mid\mathcal{D})$ by minimizing the KL between the two.

c4c5b7c4e552176521cf1b424b2208133ee93a67
ipynb
Jupyter Notebook
probability-theory/integrating-over-parameters.ipynb
Simsso/Machine-Learning-Tinkering
0a024aab0bb1ac5fbdd2f77380ab36d192278701
MIT
2020-02-19T16:30:04.000Z
2020-03-13T02:22:29.000Z
probability-theory/integrating-over-parameters.ipynb
Simsso/Machine-Learning-Tinkering
0a024aab0bb1ac5fbdd2f77380ab36d192278701
MIT
2021-02-17T14:00:33.000Z
2022-02-10T03:05:35.000Z
probability-theory/integrating-over-parameters.ipynb
Simsso/Machine-Learning-Tinkering
0a024aab0bb1ac5fbdd2f77380ab36d192278701
MIT
Qwen/Qwen-72B
1. YES
2. YES