Introduction

This chapter cover the basics concepts. Only hot takes and notes on slightly more complex concepts will be reported.

Polynomial Curve Fitting Example

Hot takes from the "Polynomial Curve Fitting" example:

• The squared error loss function is squared with respect to the polynomial weights. This means that its derivative is linear with respect to the polynomial weights and can be solved in a closed form. By finding a unique solution that minimizes the loss, we can determine the optimal polynomial curve fit.
• A higher degree \(M\) corresponds to increased flexibility in the polynomial, but it also makes the model more susceptible to overfitting. Overfitting occurs when the model captures noise or random fluctuations in the training data, leading to poor generalization to new data.
• To mitigate overfitting, regularization can be employed. Regularization limits the magnitude of the parameters by incorporating a penalty term in the loss function. When the penalty term is quadratic, the regression is known as ridge regression.
• The \(\lambda\) term in the regularization controls the complexity of the model. A higher \(\lambda\) value results in a more constrained model with reduced flexibility, helping to prevent overfitting.
• Another approach to reducing overfitting is by increasing the size of the training set. By including more diverse examples in the training data, the model can learn more general patterns and avoid over-relying on specific instances.

Probability theory

In page 18, on the paragraph about probability densities, author states that under a non-linear change of variable \(g\), a probability density transforms differently from a simple function, due to the Jacobian factor. To understand it better, let's report the answer provided in math.stackexchange.com.

Conceptually the idea is the following: if \(f\) is a probability density function, it satisfies certain properties, like \(f \ge 0\) and

\[ \int_{-\infty}^{\infty} f(x)dx = 1 \]

If we look at a transformation \(f(g(y))\), firstly the second property might be not true anymore. Hence, \(f(g(\cdot))\) is likely no probability density function. Secondly, what we roughly try to describe with \(f(g(\cdot))\) is the probability distribution of a random variable \(Y\) that is given such that \(g(Y)\) follows the distribution represented by \(f\). Say \(g\) is invertible and sufficiently smooth. The distribution of \(Y\) is given by

\[ P(Y \in A) = P(g(Y) \in g(A)) = \int_{g(A)} f(x)dx \]

Which is not very useful in practice, as this integral is on sets of the form \(g(A)\). According to the integration by substitution formula, we can compute a probability density function ℎ such that

\[ P(Y \in A) = \int_A h(y)dy \]

Here,

\[ h = f(g(\cdot)) \det J_g(\cdot) \]

where \(J_g\) is the Jacobian of \(g\).

A consequence of this observation is that the maximum of a probability density is dependent on the choice of variable.

Bayesian curve fitting

Let's revise the polynomial curve fitting example to also estimate uncertainty. We suppose that each prediction is normally distributed with the mean centered at the prediction, and the variance \(\beta^{-1}\) estimated along with the polynomial weights \(\bar w\): \[ p(t \mid x, \bar w, \beta) = \mathcal N (t \mid y(x, \bar w), \beta^{-1}) \] Supposing that the \(N\) points from the dataset \(\langle X, T \rangle\) are drawn independently, the likelihood function of a prediction is given by: \[ p(T \mid X, \bar w, \beta) = \prod_{n=1}^N \mathcal N (t_n \mid y(x_n, \bar w), \beta^{-1}) \] It is simpler and numerically convenient to work with the log-likelihood: \[ \ln p(T \mid X, \bar w, \beta) = -\frac{\beta}2 \sum_{n=1}^N (y(x_n, \bar w) - t_n)^2 + \frac{N}2\ln \beta - \frac{N}2\ln(2\pi) \] We should find \((\bar w, \beta)\) (\(\beta\) is the precision, inverse of variance) which maximize the likelihood. If we consider first the parameters \(\bar w\), then the maximization problem can drop the last 2 terms since they do not depend on \(\bar w\), which is the same as minimizing the sum of squares error loss.

Once found \(\bar w_{ML}\), we can find \(\beta\) as \[ \frac{1}{\beta_{ML}} = \frac1N \sum_{n=1}^N (y(x_n, \bar w) - t_n)^2 \] By using \(\bar w\) and \(\beta\), instead of providing a single prediction, we can provide a full distribution \(p(t \mid x, \bar w_{ML}, \beta_{ML})\) over the values of \(t\) for each \(x\).

Let \(p(\bar w \mid \alpha)\) be the prior distribution over the weights \(\bar w\), for simplicity: \[ p(\bar w \mid \alpha) = \mathcal N(\bar w \mid \bar 0, \alpha^{-1} I) \] Where the hyperparameter \(\alpha\) is the precision of the distribution. Using the Bayes theorem, the posterior distribution \(p(\bar w \mid X, T, \alpha, \beta)\) is proportional to the product of the likelihood and the prior: \[ p(\bar w \mid X, T, \alpha, \beta) \propto p(T \mid X, \bar w, \beta) p(\bar w \mid \alpha) \] By now we can find the most probable weights \(\bar w\) by maximizing the posterior distribution. This approach is called MAP (MAximum Posterior). We find that maximizing the posterior defined previously is the same as minimizing \[ \frac{\beta}{2}\sum_{n=1}^N (y(x_n, \bar w) - t_n)^2 + \frac\alpha2\bar w^T \bar w \] So maximizing the posterior is the same as minimizing the regularized sum of squares error loss function, where the regularization parameter is \(\lambda = \frac\alpha\beta\).

For a fully Bayesian treatment, we should evaluate \(p(t \mid x, X, T)\). This requires to integrate over all the possible \(\bar w\). By using the product and sum rules, we can write: \[ p(t \mid x, X, T) = \int p(t \mid x, \bar w, \beta) p(\bar w \mid X, T, \alpha, \beta) d\bar w \] Which can be solved analytically, finding that: \[ p(t \mid x, X, T) = \mathcal N (t \mid m(x), s^2(x)) \] where \[ m(x) = \beta \phi(x)^T S \sum_{n=1}^N \phi(x_n) t_n \] and \[ s^2(x) = \beta^{-1} + \phi(x)^T S \phi(x) \] and \[ S^{-1} = \alpha I + \beta \sum_{n=1}^N \phi(x_n) \phi(x_n)^T \] where \(\phi(x) = \langle 1, x, x^2, \dots, x^M\rangle\).

Decision theory

Probability theory is useful to quantify uncertainty, while decision theory allows to make optimal decisions in situations involving uncertainty.

Medical diagnosis example

Consider a medical diagnosis problem in which we have taken an X-ray image of a patient. Our algorithm has to predict if the patient has cancer or not. The input vector $\bar{x}$ is the set of pixel intensities, while the output variable $t$ is binary ($C_1=0$ no cancer, $C_2=1$ cancer).

The general inference problem involves determining the joint distribution $p(\bar{x},t)$. Given it, we then must decide to give treatments to the patient or not, and we would like this choice to be optimal. This is called the decision step.

If our goal is to make as few misclassifications as possible, then it is sufficient to study the posterior probability $p(C_k\mid \bar{x})$. Assigning each observation to a class means dividing the input space in different decision regions $R_k$ ($k=1,2$) where the boundaries are called decision boundaries / surfaces. Regions are not constrained to be contiguous but can be comprised of different disjoint sub-regions.

Given our decision regions, the probability of making a mistake is quantified by: \[ \begin{split} p(\text{mistake}) &= p(\bar x \in R_1, C_2) + p(\bar x \in R_2, C_1) \\ &= \int_{R_1} p(\bar x, C_2) d\bar x + \int_{R_2} p(\bar x, C_1) d \bar x \\ &= \int_{R_1} p(C_2 \mid \bar x)p(\bar x) d\bar x + \int_{R_2} p(C_1\mid \bar x)p(\bar x) d \bar x \end{split} \] Since the prior $p(\bar{x})$ is common in both terms, we can say that the minimum probability of making a mistake is given if each observation $\bar{x}$ is assigned to the class $C_k$ for which the posterior probability $p(C_k\mid \bar{x})$ is the largest.

For the case of $K$ classes, is easier to maximize the probability of being correct: $$>p(\text{correct})=\sum _{k=1}^{K}p(\bar{x}\in R_k,C_k)=\sum _{k=1}^K{\int }_{R_k}p(\bar{x},C_k)d\bar{x}=\sum _{k=1}^K{\int }_{R_k}p(C_k\mid \bar{x})p(\bar{x})d\bar{x}>$$ The maximum probability is obtained when each observation is assigned to the class with highest posterior.

Minimizing the expected loss

In our example, we have two different misclassification:

• A healthy patient being classified as having cancer (which is bad)
• A patient with cancer being classified as healthy (which is worse due to late treatments)

We can formalize the severity of each misclassification by constructing a loss function (or cost function) which is a single measure of loss incurred in taking any of the available decisions or actions. Let's construct a loss matrix $L$: $$L=\left[\begin{array}{cc}0& 1\\ 1000& 0\end{array}\right]$$ Where $L_{ij}$ indicates the loss severity of classifying an observation of class $i$ with class $j$. The diagonal indicates correct classifications, $L_{21}=1000$ is the loss of classifying a patient with cancer ($C_2$) to be healthy ($C_1$), vice versa for $L_{12}=1$.

The optimal solution is the one which minimizes the average loss function: $$\mathbb{E}[L]=\sum _k\sum _j{\int }_{R_j}{L}_{kj}p(\bar{x},C_k)d\bar{x}$$ Our goal is to choose the regions $R_j$ to minimize the expected loss, which can be formulated with priors instead of the joint probability. Again, the decision rule which minimize the loss is the one that assigns each observation to the class with highest posterior.

If the posterior is too low or comparable with the other posteriors, our classification has higher uncertainty of being true. We can use a threshold $\theta$ such that if $p(C_k\mid \bar{x})\le \theta$ then we avoid making a decision. This is known as the reject option.

Inference and decision

The classification is now broken in two stages:

• Inference stage, where we learn the model $p(C_k\mid \bar{x})$
• Decision stage, where we use posteriors to make optimal assignments

We can identify also 3 approaches to solve the decision stage:

1. Generative models. Model the joint distribution $p(\bar{x},C_k)$, obtain posteriors, make decisions based on posteriors. The joint distribution is useful to simulate samples from the modeled population and for outliers / novelty detection.
2. Discriminative models. Determine only the posterior class probabilities and assign each observation to the class with highest posterior accordingly.
3. Discriminative functions. Model a function $f(\bar{x})=C_k$ directly, where probabilities play no role.

The benefits of computing the posterior probabilities (avoiding approach n.3) are:

• We can modify the loss matrix without re-training the model (minimizing risk)
• We can use the reject option
• When balancing the dataset, we can compensate for class priors (A)
• We can combine models (B)

(A) Suppose we have 1000 samples from class $C_1$ and 100 samples from class $C_2$, so the real priors are $p(C_1)=0.91$ and $p(C_2)=0.09$. Suppose we want to balance our datasets to 100 samples for each class. We know the real priors and hence we can replace them when using the Bayes theorem.

(B) Suppose we have X-ray data ${\bar{x}}_I$ and blood test data ${\bar{x}}_B$, we can develop two models (one for each) and assume that the distributions of the inputs are independent given the class (conditional independence) $$>p({\bar{x}}_I,{\bar{x}}_B\mid C_k)=p({\bar{x}}_I\mid C_k)p({\bar{x}}_B\mid C_k)>$$ The posterior is then $$>\begin{array}{rl}>P(C_k\mid {\bar{x}}_I{\bar{x}}_B)& \propto p({\bar{x}}_I,{\bar{x}}_B\mid C_k)p(C_k)\\ >& \propto p({\bar{x}}_I\mid C_k)p({\bar{x}}_B\mid C_k)p(C_k)\\ >& \propto \frac{p(C_k\mid {\bar{x}}_I)p(C_k\mid {\bar{x}}_B)}{p(C_k)}>\end{array}>$$ Hence the combination of 2+ models is trivial.

Loss functions for regression

The decision stage for regression problems consists of choosing a specific estimate $y(x)$ of the value $t$ for each input $x$. In doing so, we incur a loss $L(t,y(x))$. The expected loss is then given by:

$$\mathbb{E}[L]=\int \int L(t,y(x))p(x,t)dxdt$$

A common choice of loss is the squared loss

$$\mathbb{E}[L]=\int \int \{y(x)-t{\}}^{2}p(x,t)dxdt$$

Supposing that $y$ is a flexible function, this can be solved by using the calculus of variations, obtaining the following regression function

$$y(x)=\int t\cdot p(t\mid x)dt={\mathbb{E}}_t[t\mid x]$$

Notes on Information theory section

KL divergence and Jensen Inequality

Let $p$ and $q$ be the exact and an approximate distributions of the random variable $x$. The KL divergence defines the additional amount of information (in nats for $\ln$ bits for ${\log }_2$) required to transmit $x$ assuming its distribution is $q$ (approximated) instead of $p$ (exact). The mathematical definition is:

$$KL(p||q)=H(p||q)-H(p)=-\int p(x)\ln \frac{q(x)}{p(x)}dx$$

From the Jensen Inequality, we know that for a convex function $f$

$$f(\mathbb{E}[x])\le \mathbb{E}(f(x))$$

Since $-\ln$ is a strictly convex function, we can apply the Jensen Inequality to the KL divergence:

$$\begin{array}{rl}KL(p||q)& =-\int p(x)\ln \frac{q(x)}{p(x)}dx\\ & \ge -\ln \int p(x)\frac{q(x)}{p(x)}dx\\ & \ge -\ln \int q(x)dx=-\ln 1=0\end{array}$$

Since the function $-ln$ is strictly convex, the equality $KL(p||q)=0$ holds if and only if $q(x)=p(x)$ for all values of $x$.

Minimizing KL = Maximizing likelihood

Suppose we want to approximate an unknown distribution $p(x)$ using a parametric distribution $q(x\mid \theta )$. Suppose we have $N$ training points drawn from the distribution $p(x)$, we can approximate the KL divergence as:

$$KL(p||q)\approx \sum _{n=1}^N\left[-\ln q(x_n\mid \theta )+\ln p(x_n)\right]$$

To find the best values of $\theta$ to approximate $p$, we can minimize the approximated KL divergence:

$$\begin{array}{rl}\underset{\theta }{\min }KL(p||q)& \approx \underset{\theta }{\min }\sum _{n=1}^N\left[-\ln q(x_n\mid \theta )+\ln p(x_n)\right]\\ & =\underset{\theta }{\min }-\sum _{n=1}^N\ln q(x_n\mid \theta )\\ & =\underset{\theta }{\max }\sum _{n=1}^N\ln q(x_n\mid \theta )\\ & =\underset{\theta }{\max }\sum _{n=1}^{N}q(x_n\mid \theta )\end{array}$$

Therefore, minimizing the KL divergence is the same as maximizing the likelihood.

Mutual information

Let $x,y$ be two random variables. We don't know if they are indipendent or not, but if they are then

$$p(x,y)=p(x)p(y)$$

To check if they are indipendent we can calculate the KL divergence between $p(x,y)$ and $p(x)p(y)$. Given the properties of the KL divergence, if $KL=0$ then $p(x,y)=p(x)p(y)$. This quantity is called mutual information between the variables $x,y$

$$I[x,y]=KL(p(x,y)||p(x)p(y))=-\int \int p(x,y)\ln \frac{p(x)p(y)}{p(x,y)}dxdy$$

It is related to the conditional entropy through:

$$I[x,y]=H[x]-H[x\mid y]=H[y]-H[y\mid x]$$

The mutual information $I[x,y]$ represents the reduction in uncertainty about $x$ as a consequence of the new observation $y$ (or vice versa).

Chapter 1 Exercises

Solution to exercise 1.1

Given

\[ E(\bar w) = \frac12 \sum_{n=1}^N \left( x_n^M w_M + \dots + x_n w_1 + w_0 - t_n \right)^2 \]

Compute the partial derivative

\[ \frac{\partial E}{\partial w_i} = \sum_{n=1}^N x_n^i \left( x_n^M w_M + \dots + x_n w_1 + w_0 - t_n \right) \]

Set the partial derivatives to zero to get $w$ that minimizes \(E\)

\[ \frac{\partial E}{\partial w_i} = \sum_{n=1}^N x_n^i \left( x_n^M w_M + \dots + x_n w_1 + w_0 - t_n \right) = 0 \]

\[ = \sum_{n=1}^N x_n^{M+i} w_M + \dots + x_n^{1+i} w_1 + x_n^{0 + i}w_0 - x_n^{i} t_n = 0 \]

\[ = \sum_{n=1}^N x_n^{M+i} w_M + \dots + x_n^{1+i} w_1 + x_n^{0 + i}w_0 = \sum_{n=1}^N x_n^{i} t_n \]

\[ = w_M \sum_{n=1}^N x_n^{M+i} + \dots + w_j \sum_{n=1}^N x_n^{j+i} + \dots + w_0 \sum_{n=1}^N x_n^{0 + i} = \sum_{n=1}^N x_n^{i} t_n \]

\[ = A_{iM} w_M + \dots + A_{ij} w_j + \dots + A_{i0} w_0 = T_i \]

\[ = \sum_{j=1}^M A_{ij}w_j = T_i \]

Solution to exercise 1.2

Given

\[ E(\bar w) = \frac12 \sum_{n=1}^N \left( x_n^M w_M + \dots + x_n w_1 + w_0 - t_n \right)^2 + \frac\lambda2 || \bar w ||^2 \]

Compute the partial derivative

\[ \frac{\partial E}{\partial w_i} = \left[ \sum_{n=1}^N x_n^i \left( x_n^M w_M + \dots + x_n w_1 + w_0 - t_n \right) \right] + \lambda w_i \]

Set the partial derivatives to zero to get $w$ that minimizes \(E\)

\[ \frac{\partial E}{\partial w_i} = \sum_{n=1}^N x_n^i \left( x_n^M w_M + \dots + x_n w_1 + w_0 - t_n \right) + \lambda w_i = 0 \]

\[ = \sum_{n=1}^N x_n^{M+i} w_M + \dots + x_n^{1+i} w_1 + x_n^{0 + i}w_0 - x_n^{i} t_n + \lambda w_i = 0 \]

\[ = \sum_{n=1}^N x_n^{M+i} w_M + \dots + x_n^{1+i} w_1 + x_n^{0 + i}w_0 + \lambda w_i = \sum_{n=1}^N x_n^{i} t_n \]

\[ = w_M \sum_{n=1}^N x_n^{M+i} + \dots + w_i (\lambda + \sum_{n=1}^N x_n^{2i}) + \dots + w_0 \sum_{n=1}^N x_n^{0 + i} = \sum_{n=1}^N x_n^{i} t_n \]

\[ = \hat A_{iM} w_M + \dots + \hat A_{ij} w_j + \dots + \hat A_{i0} w_0 = T_i \]

\[ = \sum_{j=1}^M \hat A_{ij}w_j = T_i \]

where

\[ \hat A_{ij} = \begin{cases} \lambda + \sum_{n=1}^N x_n^{2i} & \text{if } i = j \\ \sum_{n=1}^N x_n^{j+i} & \text{otherwise} \end{cases} \]

or simply

\[ \hat A = \lambda I_M + A \]

Solution to exercise 1.3

The probability of selecting an apple is

\[ \begin{split} P(a) &= P(a, r) + P(a, b) + P(a, g) \\ &= P(a \mid r)P(r) + P(a \mid b)P(b) + P(a \mid g)P(g) \\ &= 3/10 \times 0.2 + \frac12 \times 0.2 + 3/10 \times 0.6 = 0.34 \end{split} \]

Following the Bayes theorem, the probability of the selected box being green given that the selected fruit is an orange is

\[ P(g \mid o) = \frac{P(o \mid g) P(g)}{P(o)} = \frac{0.3 \times 0.6}{0.36} = 0.5 \]

where

\[ \begin{split} P(o) &= P(o, r) + P(o, b) + P(o, g) \\ &= P(o \mid r)P(r) + P(o \mid b)P(b) + P(o \mid g)P(g) \\ &= 4/10 \times 0.2 + \frac12 \times 0.2 + 3/10 \times 0.6 = 0.36 \end{split} \]

Solution to exercise 1.4

Revised solution from Bishop's solution manual.

Let $g$ be a non-linear change of variable $x=g(y)$, for probability density functions we know that $$p_y(y)=p_x(g(y))\cdot |g^{\prime }(y)|$$ Let $\hat{x},\hat{y}$ be the maximum of $p_x,p_y$ densities respectively. Let $s=\text{sign}(g^{\prime }(y))\in \{-1,1\}$ and re-write: $$p_y(y)=p_x(g(y))\cdot s{g}^{\prime }(y)$$ Differentiate both sides: $$p^{\prime }(y)=s{p}_x^{\prime }(g(y))[g^{\prime }(y){]}^2+s{p}_x(g(y))g^{\prime \prime }(y)$$ Suppose that $\hat{x}=g(\hat{y})$, then $$\begin{array}{rl}{p}^{\prime }(\hat{y})& =s{p}_x^{\prime }(g(\hat{y}))[g^{\prime }(\hat{y}){]}^2+s{p}_x(g(\hat{y}))g^{\prime \prime }(\hat{y})\\ & =s{p}_x^{\prime }(\hat{x})[g^{\prime }(\hat{y}){]}^2+s{p}_x(\hat{x})g^{\prime \prime }(\hat{y})\\ & =s\cdot 0\cdot [g^{\prime }(\hat{y}){]}^2+s{p}_x(\hat{x})g^{\prime \prime }(\hat{y})\\ & =s{p}_x(\hat{x})g^{\prime \prime }(\hat{y})=0\end{array}$$ Where:

1. $s\in \{-1,1\}$ cannot be zero
2. $p_x(\hat{x})$ is the maximum probability, thus cannot be zero

This means $\frac{{\partial }^{2}g(\hat{y})}{\partial y^2}$ has to be 0, which depends on $g$, hence the relation $\hat{x}=g(\hat{y})$ may not hold. If $g$ is linear, then the second derivative of $g$ is 0 and the relation $\hat{x}=g(\hat{y})$ is valid.

Solution to exercise 1.5

$$\begin{array}{rl}var[X]& =\int p(x)[f(x)-\mathbb{E}[f(x)]{]}^{2}dx\\ & =\int p(x)[f(x{)}^2-2\mathbb{E}[f(x)]f(x)+\mathbb{E}[f(x){]}^2]dx\\ & =\int (p(x)f(x{)}^2-2p(x)\mathbb{E}[f(x)]f(x)+p(x)\mathbb{E}[f(x){]}^2)dx\\ & =\int p(x)f(x{)}^{2}dx-2\int p(x)\mathbb{E}[f(x)]f(x)dx+\int p(x)\mathbb{E}[f(x){]}^{2}dx\\ & =\mathbb{E}[f(x{)}^2]-2\mathbb{E}[f(x)]\int p(x)f(x)dx+\mathbb{E}[f(x){]}^2\int p(x)dx\\ & =\mathbb{E}[f(x{)}^2]-2\mathbb{E}[f(x){]}^2+\mathbb{E}[f(x){]}^2\\ & =\mathbb{E}[f(x{)}^2]-\mathbb{E}[f(x){]}^2\end{array}$$

Solution to exercise 1.6

From 1.41

$$cov(x,y)={\mathbb{E}}_{x,y}[xy]-\mathbb{E}[x]\mathbb{E}[y]$$

But if x and y are indipendent, then

$$\begin{array}{rl}{\mathbb{E}}_{x,y}[xy]& =\int \int p(x,y)xydxdy\\ & =\int \int xp(x)yp(y)dxdy\\ & =\mathbb{E}[x]\int yp(y)dy\\ & =\mathbb{E}[x]\mathbb{E}[y]\end{array}$$

Therefore $cov(x,y)=0$.

Solution to exercise 1.32

Let $x\sim p_x(x)$ and let $y=Ax$ be a linear change of variable. In that case, the jacobian factor is the determinant $|A|$ and we can write

$$p_y(y)=p_x(y)|A|$$

So we can write

$$\begin{array}{rl}H(y)& =-\int p_y(y)\ln p_y(y)dy\\ & =-\int p_y(y)\ln [p_x(y)|A|]dy\\ & =-\int p_y(y)\left[\ln p_x(y)+\ln |A|\right]dy\\ & =[-\int p_y(y)\ln p_x(y)dy]-\ln |A|\int p_y(y)dy\\ & =[-\int p_y(y)\ln p_x(y)dy]-\ln |A|\end{array}$$

(last steps: solve the integral on the left-hand side using the substitution $A^{-1}y=x$ remembering that A is non-singular)

Probability distributions

This chapter will focus on the problem of density estimation, which consists in finding / estimating the probability distribution $p(x)$ from $N$ independent and identically distributed datapoints $x_1,x_2,\dots ,x_N$ drawn from $p(x)$. There are two main ways of doing that: the first way is to use parametric density estimation, where you choose one known parametric distribution (e.g., Gaussian) and try to get the right parameters that fit the data. This method assumes that the parametric distribution we use it's suitable for the data, which is not always the case. Another way of doing that is using non-parametric density estimation techniques (e.g., histograms, nn, kernels).

Bernoulli experiment

Suppose we have a data set $D=\{x_1,\dots ,x_N\}$ of i.i.d. observed values of $x\sim Bern(x\mid \mu )$. We can estimate the $\mu$ parameter from the sample in a frequentist way, by maximizing the likelihood (or the log-likelihood):

$$\begin{array}{rl}\ln p(D\mid \mu )& =\ln \prod _{n=1}^{N}p(x_i\mid \mu )\\ & =\sum _{n=1}^N\ln [{\mu }^{x_n}(1-\mu {)}^{(1-x_n)}]\\ & =\sum _{n=1}^N[x_n\ln \mu +(1-x_n)\ln (1-\mu )]\end{array}$$

To find $\mu$, let's set the log-likelihood derivative w.r.t. $\mu$ to 0:

$$\begin{array}{rl}\frac{\partial }{\partial \mu }\ln p(D\mid \mu )& =0\\ \sum _{n=1}^N\left(\frac{x_n}{\mu }+\frac{1-x_n}{1-\mu }\right)& =0\\ \sum _{n=1}^N\frac{x_n-\mu }{\mu +{\mu }^2}& =0\end{array}$$

Since $\frac{x_n-\mu }{\mu +{\mu }^2}=0\iff x_n=\mu$ then:

$$\begin{array}{rl}\sum _{n=1}^N{x}_n-\mu & =0\\ \sum _{n=1}^N{x}_n-\sum _{n=1}^N\mu & =0\\ \sum _{n=1}^N{x}_n& =N\mu \\ \frac{1}{N}\sum _{n=1}^N{x}_n& =\mu \end{array}$$

$\mu$ is estimated from the sample mean. In this case, the sample mean is an example of sufficient statistic for the model, i.e. calculating other statistics from the sample will not add more information than that.

Binomial distribution

Sticking with the coin flips example, the binomial distribution models the probability of obtaining $m$ heads out of $N$ total coin flips:

$$Bin(m\mid N,\mu )=\left(\genfrac{}{}{0ex}{}{N}{m}\right){\mu }^m(1-\mu {)}^{N-m}$$

Where $\left(\genfrac{}{}{0ex}{}{N}{m}\right)$ represents all the possible ways of obtaining $m$ heads out of $N$ coin flips. The mean and variance of a binomial variabile can be estimated by knowning that for i.i.d events the mean of the sum is the sum of the mean, and the variance of the sum is the sum of the variances. Because $m=x_1+\cdots +x_N$ then:

$$\mathbb{E}[m]=N\mu var[m]=N\mu (1-\mu )$$

Beta distribution

Please read the estimating_parameters_using_a_bayesian_approach notebook. Some quick notes here:

$$\beta (\mu )=\frac{\Gamma (a+b)}{\Gamma (a)+\Gamma (b)}{\mu }^{(a-1)}(1-\mu {)}^{(b-1)}$$

and

$$\mathbb{E}[\mu ]=\frac{a}{a+b}var[\mu ]=\frac{ab}{(a+b{)}^2(a+b+1)}$$

Multinomial variables

It's a generalization of the Bernoulli distribution where a random variable has $K$ possible values instead of being binary. We can represent the variable as a $K$-dimensional binary vector $x=⟨x_1,x_2,\dots ,x_K⟩$ where only one component can be asserted:

$$\sum _{k=1}^K{x}_k=1$$

The probability of each component to be asserted is regulated by a probability vector $\mu =⟨{\mu }_1,{\mu }_2,\dots ,{\mu }_K⟩$, so that basically $x_k\sim Bern({\mu }_k)$. Since the $\mu$ vector represents a probability distribution, then:

$$\sum _{k=1}^K{\mu }_k=1$$

The multinomial distribution of $x$ is given by:

$$p(x\mid \mu )=\prod _{k=1}^K{\mu }_k^{x_k}$$

And the expected values is $\mathbb{E}[x\mid \mu ]=\mu$. Let's consider a dataset of N independente observations, then the likelihood function is:

$$p(D\mid \mu )=\prod _{n=1}^N\prod _{k=1}^K{\mu }_k^{x_{nk}}=\prod _{k=1}^K{\mu }_k^{({\sum }_{n=1}^N{x}_{nk})}=\prod _{k=1}^K{\mu }_k^{m_k}$$

where $m_k={\sum }_{n=1}^N{x}_{nk}$.

If we want to find $\mu$ from $D$ by maximizing the (log) likelihood, we have to constrain that to be a probability distribution and therefore we can use the Lagrangian multiplier $\lambda$

$$\sum _{k=1}^K{m}_k\ln {\mu }_k+\lambda \left(\sum _{k=1}^K{\mu }_k-1\right)$$

Setting the derivative w.r.t. ${\mu }_k$ to zero we get ${\mu }_k=-m_k/\lambda$. We can solve for the Lagrangian multiplier $\lambda$ by replacing this result in the equation and then we get that $\lambda =-N$ and ${\mu }_k^{ML}=m_k/N$.

We can also consider the distribution of the quantities $m_1,\dots ,m_k$ (Multinomial distribution) conditioned on the parameter $\mu$ and on the number $N$ of observations:

$$Mult(m_1,dots,m_K\mid \mu ,N)=\left(\genfrac{}{}{0ex}{}{N}{m_1{m}_2\dots m_K}\right)\prod _{k=1}^K{\mu }_k^{m_k}$$

where

$$\left(\genfrac{}{}{0ex}{}{N}{m_1{m}_2\dots m_K}\right)=\frac{N!}{m_1!m_2!\dots m_K!}\sum _{k=1}^K{m}_k=N$$

Short description of Lagrangian Multiplier utility taken from Quora: You are trying to maximize or minimize some function $f$ (distance to treasure), while keeping some other function $g$ fixed at a certain value $c$ (stay on the path). At this point, the gradient $\nabla f$ (the compass needle) must be parallel to the gradient $\nabla g$ (the arrows on the signs), but the two vectors will not generally have the same length. The test for whether or not they’re parallel is $\nabla f=\lambda \nabla g$, where $\lambda$ is whatever multiplier is needed to have them match; it will still only be able to be equal if they’re parallel (you can resize the compass needle however you want to make it match the sign arrow, but you have to be at a spot with the right direction).

Dirichlet distribution

While the beta distribution is a prior of the Bernoulli parameter $\mu$, the Dirichlet distribution is a prior of the Multinomial probability vector $\bar{\mu }$. The definition is:

$$Dir(\mu \mid \alpha )=\frac{\Gamma ({\alpha }_0)}{\Gamma ({\alpha }_1)\dots \Gamma ({\alpha }_K)}\prod _{k=1}^K{\mu }_k^{{\alpha }_k-1}$$

Where ${\alpha }_0={\sum }_{k=1}^K{\alpha }_k$. Since the $\bar{\mu }$ parameters are bounded to ${\sum }_k{\mu }_k=1$, then the distribution is confined to a simplex in the $K-1$ space.

By multiplying the likelihood function (which is the multinomial distribution) by the prior (which is a Dirichlet distribution) we get something that is proportional to the posterior $p(\mu \mid D,\alpha )$. Assuming a conjugate prior, the posterior has the same form and hence we can derive the normalization constant by comparison with the dirichlet distribution definition. The posterior is defined as:

$$p(\mu \mid D,\alpha )=\frac{\Gamma ({\alpha }_0+N)}{\Gamma ({\alpha }_1+m_1)\dots \Gamma ({\alpha }_K+m_k)}\prod _{k=1}^K{\mu }_k^{{\alpha }_k+m_k-1}$$

Gaussian distribution

Univariate Gaussian distribution:

$$\mathcal{N}(x\mid \mu ,{\sigma }^2)=\frac{1}{(2\pi {\sigma }^2{)}^{1/2}}\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)$$

Where $\mu$ and ${\sigma }^2$ are the mean and variance of the population.

Multivariate Gaussian distribution:

$$\mathcal{N}(\bar{x}\mid \bar{\mu },\Sigma )=\frac{1}{(2\pi {)}^{D/2}|\Sigma {|}^{1/2}}\exp \left(-\frac{(\bar{x}-\bar{\mu }{)}^T{\Sigma }^{-1}(\bar{x}-\bar{\mu })}{2}\right)$$

Where $D$ is the dimensionality of $\bar{x}$, $\mu$ is the mean vector and $\Sigma$ is the covariance matrix.

Central Limit Theorem. Subject to certain mild conditions, the sum of a set of random variables, which is of course itself a random variable, has a distribution that becomes increasingly Gaussian as the number of terms in the sum increases.

Observation n.1 - The covariance matrix $\Sigma$ is always positive semi-definite. This means that the eigenvalues are non-negative.

Observation n.2 - To be well-defined, a Gaussian must have a positive definite covariance matrix, which means that all the eigenvalues are strictly positive. If the covariance matrix has one or more null eigenvalues (positive semi-definite), then the distribution is singular and is confined to a subspace of lower dimensionality.

Observation n.3 - Given a 2D Gaussian distribution, we find that elements $x$ with constant density are distributed as ellipses, where the axis of the ellipse is given by the eigenvectors, and the length of the axis is proportional to the square root of the corresponding eigenvalue. The ellipse defined by the axis having a length equal to the square root of the eigenvalues, we find all the elements with a density of $\exp (-1/2)$.

This is a hint on my preferred interpretation of eigenvectors and eigenvalues calculated from data: eigenvectors represent that capture most of the variability, and the corresponding eigenvalues are an indicator of the variability on that axis (see PCA).

Observation n.4 - A multivariate Gaussian can be decomposed as a product of multiple univariate Gaussians.

Observation n.5 - The book provides a formal proof to find that $\mathbb{E}[\bar{x}]=\bar{\mu }$ and $\text{var}[\bar{x}]=\Sigma$.

Observation n.6 - A general symmetric covariance matrix $\Sigma$ has $D(D+1)/2$ independent parameters, and there are another $D$ independent parameters in $\mu$, giving $D(D+3)/2$ parameters in total. This means that the number of parameters grows quadratically with the dimension $D$.

One way to reduce the number of parameters is to use restricted forms of the covariance matrix. For example, with a diagonal covariance matrix (figure b) we a linear dependency between parameters and dimensionality, but we explain the variability only from the features axis. Another possibility is to use the isotropic covariance (Figure c) ${\sigma }^{2}I$ where we have only one parameter, but we discard the variability along multiple axis. We have a trade-off between model complexity and flexibility that must be adressed based on the application.

Limitations of the Gaussian - Gaussian is unimodal (only 1 maximum) and thus is not good at representing multimodal data.

Conditional & Marginal Gaussian Distributions

Given a joint Gaussian distribution $\mathcal{N}(x\mid \mu ,\Sigma$ with $\Lambda ={\Sigma }^{-1}$ (precision matrix) and

$$x=\left[\begin{array}{c}{x}_a\\ x_b\end{array}\right]\mu =\left[\begin{array}{c}{\mu }_a\\ {\mu }_b\end{array}\right]$$

$$\Sigma =\left[\begin{array}{cc}{\Sigma }_{aa}& {\Sigma }_{ab}\\ {\Sigma }_{ba}& {\Sigma }_{bb}\end{array}\right]\Lambda =\left[\begin{array}{cc}{\Lambda }_{aa}& {\Lambda }_{ab}\\ {\Lambda }_{ba}& {\Lambda }_{bb}\end{array}\right]$$

where ${\Sigma }_{ab}={\Sigma }_{ba}^T$ and ${\Lambda }_{ab}={\Lambda }_{ba}^T$.

Then we have that the conditional distribution $p(x_a\mid x_b)$ is

$$\begin{array}{rl}p(x_a\mid x_b)& =\mathcal{N}({\mu }_{a\mid b},{\Lambda }_{aa}^{-1})\\ {\mu }_{a\mid b}& ={\mu }_a+{\Sigma }_{ab}{\Sigma }_{bb}^{-1}(x_b-{\mu }_b)\end{array}$$

Note that ${\mu }_{a\mid b}$ is a linear function of $x_b$.

where

$${\Lambda }_{aa}^{-1}={\Sigma }_{aa}-{\Sigma }_{ab}{\Sigma }_{bb}^{-1}{\Sigma }_{ba}$$

And that the marginal distribution $p(x_a)$:

$$p(x_a)=\mathcal{N}(x_a\mid {\mu }_a,{\Sigma }_{aa})$$

All the derivations focus on the quadratic relationship of the exponential factor to the $x$ and are detailed starting from page 85. The point is that the conditional and marginal distribution of a joint Gaussian distribution are again Gaussian distributions.

Given a marginal Gaussian distribution for $x$ and $y$ and a conditional Gaussian distribution for $y$ given $x$ in the form:

$$\begin{array}{r}p(x)=\mathcal{N}(x\mid \mu ,{\Lambda }^{-1})\\ p(y\mid x)=\mathcal{N}(Ax+b,L^{-1})\end{array}$$

Where $Ax+b$ expresses the fact that the mean of the conditional distribution $y$ given $x$ is a linear function of $x$, and $L^{-1}$ is another precision matrix. The marginal distribution of $y$ and the conditional distribution of $x$ given $y$ are given by:

$$\begin{array}{rl}p(y)& =\mathcal{N}(y\mid A\mu +b,L^{-1}+A\Lambda A^T)\\ p(x\mid y)& =\mathcal{N}(x\mid \Sigma [A^{T}L(y-b)+\Lambda \mu ],\Sigma )\end{array}$$

where

$$\Sigma =(\Lambda +A^{T}LA{)}^{-1}$$

Maximum likelihood for the Gaussian

Let $X=(x_1,\dots ,x_N{)}^T$ be a dataset of observation $x_n$ drawn independently from a multivariate Gaussian distribution. We can estimate the parameters by maximizing the log-likelihood:

$$p(X\mid \mu ,\Sigma )=-\frac{ND}{2}\ln (2\pi )-\frac{N}{2}\ln |\Sigma |-\frac{1}{2}\sum _{n=1}^N(x_n-\mu {)}^T{\Sigma }^{-1}(x_n-\mu )$$

By setting the derivative to zero, we can compute:

$${\mu }_{ML}=\frac{1}{N}\sum _{n=1}^N{x}_n{\Sigma }_{ML}=\frac{1}{N}\sum _{n=1}^N(x_n-{\mu }_{ML})(x_n-{\mu }_{ML}{)}^T$$

Where we calculate the two parameters in sequential steps since there's no dependency on $\Sigma$ when maximizing with respect to $\mu$. By taking the expectation we see that:

$$\mathbb{E}[{\mu }_{ML}]=\mu \mathbb{E}[{\Sigma }_{ML}]=\frac{N-1}{N}\Sigma$$

We see that ${\Sigma }_{ML}$ is a biased estimator, and we can correct it by:

$${\tilde{\Sigma }}_{ML}=\frac{1}{N-1}\sum _{n=1}^N(x_n-{\mu }_{ML})(x_n-{\mu }_{ML}{)}^T$$

Now $\mathbb{E}[{\tilde{\Sigma }}_{ML}]=\Sigma$ is a correct extimator for the true covariance.

Student's t distribution

In order to estimate the precision $\tau$ of a Gaussian distribution $\mathcal{N}(x\mid \mu ,{\tau }^{-1})$ by using a Bayesian approach, we can use the Gamma distribution as a prior: $$Gam(\tau \mid a,b)=\frac{b^a}{\Gamma (a)}{\tau }^{a-1}\exp (-b\tau )$$ Then we can marginalize to obtain $$p(x\mid \mu ,a,b)={\int }_0^{\infty }\mathcal{N}(x\mid \mu ,{\tau }^{-1})Gam(\tau \mid a,b)d\tau$$ By performing a number of steps, we can say that this is a Student's t distribution: $$St(x\mid \mu ,\lambda ,\nu )=\frac{\Gamma (\frac{\nu }{2}+\frac{1}{2})}{\frac{\nu }{2}}\sqrt{\frac{\lambda }{\pi \nu }}{\left[1+\frac{\lambda (x-\mu {)}^2}{\nu }\right]}^{-\frac{\nu }{2}-\frac{1}{2}}$$ Where $\lambda =a/b$ (precison) and $\nu =2a$ (degrees of freedom). The precision of the t distribution doesn't correspond to the inverse of the variance!

The parameters can be estimated by Expectation-Maximization, and the result has the property of robustness: outliers does not severely affect the distribution.

The multivariate version of the t distribution is: $$St(x\mid \mu ,\Lambda ,\nu )=\frac{\Gamma (\frac{D}{2}+\frac{\nu }{2})}{\frac{\nu }{2}}\frac{\sqrt{|\Lambda |}}{(\pi \nu {)}^{D/2}}{\left[1+\frac{{\Delta }^2}{\nu }\right]}^{-\frac{D}{2}-\frac{\nu }{2}}$$ and $${\Delta }^2=(x-\mu {)}^T\Lambda (x-\mu )$$ The statistics are

$$\begin{gathered} \mathbb{E}[x]=\mu \text{ if }\nu >1 \\ cov[x]=\frac{\nu }{\nu -2}{\Lambda }^{-1}\text{ if }\nu >2 \\ mode[x]=\mu \end{gathered}$$

Exponential family

A distribution that is part of the exponential family can be represented as:

$$p(x\mid \eta )=h(x)g(\eta )\exp ({\eta }^{T}u(x))$$

Where $\eta$ are the natural parameters of the distribution. The function $g(\eta )$ ensures that the distribution is normalized:

$$g(\eta )\int h(x)\exp ({\eta }^{T}u(x))dx=1$$

The Bernoulli ($\eta =\mu$), Multinomial ($\eta =⟨{\mu }_1,\dots ,{\mu }_K⟩$) and Gaussian distributions $(\eta =⟨\mu ,\sigma ⟩)$ are part of this family, and the PRML book proofs this at page 113.

If we want to estimate $\eta$, we can do that by maximum-likelihood. Let's set the gradient of $p(x\mid \eta )$ w.r.t. $\eta$ to 0:

$$\begin{gathered} \nabla \left[(g(\eta )\int h(x)\exp ({\eta }^{T}u(x))dx\right]=0 \\ \nabla g(\eta )\underset{\text{(b)}}{\underbrace{\int h(x)\exp ({\eta }^{T}u(x))dx}}+g(\eta )\underset{\text{(a)}}{\underbrace{\int h(x)\exp ({\eta }^{T}u(x))u(x)dx}}=0 \end{gathered}$$

(a) To understand the underlined expression, remember that the integration is w.r.t. $x$, while the differentiation is w.r.t. $\nabla$, so we directly differentiate the content as the derivative of sums is the sum of derivatives.

(b) recall that $g(\eta )\int h(x)\exp ({\eta }^{T}u(x))dx=1$, this means that $\int h(x)\exp ({\eta }^{T}u(x))dx=1/g(\eta )$

$$\begin{gathered} -\frac{1}{g(\eta )}\nabla g(\eta )=\int h(x)g(\eta )\exp ({\eta }^{T}u(x))u(x)dx \\ -\frac{1}{g(\eta )}\nabla g(\eta )=\int p(x\mid \eta )u(x)dx \\ \underset{\text{(c)}}{\underbrace{-\frac{1}{g(\eta )}\nabla g(\eta )}}=\mathbb{E}[u(x)] \\ -\nabla \ln g(\eta )=\mathbb{E}[u(x)] \end{gathered}$$

Point (c) is provided by the identity $\nabla \ln g(\eta )=\frac{1}{g(\eta )}\nabla g(\eta )$ (chain rule).

We will use this result later on. Now suppose to have $N$ i.i.d. observations drawn from the exponential distribution $X=\{x_1,\dots ,x_N\}$. The likelihood function is given by:

$$p(X\mid \eta )=(\prod _{n=1}^{N}h(x_n))g(\eta {)}^N\exp \{{\eta }^T\sum _{n=1}^{N}u(x_n)\}$$

Setting $\nabla \ln p(X\mid \eta )=0$ we get:

$$-\nabla \ln g({\eta }_{ML})=\frac{1}{N}\sum _{n=1}^{N}U(x_n)$$

which can be solved to obtain ${\eta }_{ML}$. The solution depends on the data only through ${\sum }_{n}u(x_n)$, which is called the sufficient statistic of the exponential distribution. For $N\to \infty$ then ${\eta }_{ML}\to \eta$,

For each exponential distribution of the previous form, there exists a conjugate prior distribution over the parameters $\eta$ of the following form

$$p(\eta \mid \chi ,\nu )=f(\chi ,\nu )g(\eta {)}^{\nu }\exp \left\{\nu {\eta }^T\chi \right\}$$

Where $f(\chi ,\eta )$ is a normalization coefficient, $g(\eta )$ is the same function presented in the exponential distribution, $\nu$ can be interpreted as a effective number of pseudo-observations in the prior, each of which has a value for the sufficient statistic $u(x)$ given by $\chi$.

Again - why do we need conjugate priors??

A prior $p(\eta )$ which is conjugate to the likelihood $p(X\mid \eta )$ produce a posterior $p(\eta \mid X)$ that has the same functional form as the choosen prior. This allows to derive a closed-form expression for the posterior distribution (otherwise you need to compute the normalization coefficient by integration, YOU DON'T WANT TO DO THAT, RIGHT?)

Noninformative prior

If we have no prior information, we want a prior with minimal influence on the inference. We call such a prior a noninformative prior. The Bayes/Laplace postulate, stated about 200 years ago says the following:

The principle of insufficient reason. When nothing is known about $\theta$ in advance, let the prior $\pi (\theta )$ be a uniform distribution, that is, let all possible outcomes of $\theta$ have the same probability.

One noninformative prior could be the uniform distribution, but there are two problems:

1. If the parameter $\theta$ is unbounded, the prior distribution cannot be correctly normalized because the integral over $\lambda$ diverges. In that case, we have an improper prior. In pratice, improper priors are used if the posterior is proper, i.e. is correctly normalized.
2. The second problem is that if we perform a non-linear change of variable, then the resulting density will not be constant (recall the Jacobian multiplier).

Nonparametric Methods

The distribution we have seen are governed by parameters that are estimated from the data. This is called parametric approach to density modelling.

In this section we talk about nonparametric approaches to density estimation (only simple frequentist methods).

Consider a continuous variable $x$, the simplest way to model that distribution is to partition $N$ observations of $x$ in different bins of width ${\Delta }_i$ (often the same for every bin ${\Delta }_i=\Delta$), and then count the number $n_i$ of observation of $x$ falling in bin $i$. To turn this count into a normalized probability density:

$$p_i=\frac{n_i}{N{\Delta }_i}$$

Problems:

• If we choose $\Delta$ too small, the resulting distribution will be too spiky (i.e. will show structures that are not present in the real distribution)
• If we choose $\Delta$ too big, the resulting distribution may fail to capture the structures of the real distribution

Advantages:

• Good visualization of the distribution
• The dataset can be discarded once the histogram is built
• Good setup if data points are arriving sequentially

Disvantages:

• The estimated density has discontinuities in the bin edges
• Does not scale with dimensionality (curse of dimensionality), the amount of data needed to work in high dimensional spaces is prohibitive.

Good ideas:

• To estimate probability density at a particular location, we should consider the data points that lie within some local neighbourhood of that point.

Kernel Density Estimators

Suppose that we want to estimate the value of $p(x)$, where $x\in {\mathbb{R}}^D$, and the PDF $p$ is totally unknown. We want to use a non-parametric estimation method. We suppose that $x$ lies in an Euclidean space.

Let $R$ be a region containing the point $x$, the probability $P$ of falling into the region $x$ is defined by:

$$P={\int }_{R}p(x)dx$$

Let $x_1,\dots ,x_N$ be $N$ i.i.d. data points drawn from $p(x)$, then the probability that $K$ of $N$ data points fall into region $R$ is a binomial distribution:

$$Bin(K\mid N,P)=\left(\genfrac{}{}{0ex}{}{N}{K}\right){\mu }^K(1-\mu {)}^{N-K}$$

Fromt the properties of the Binomial distribution, we have that:

1. $\mathbb{E}[K/N]=P$
2. $var[K/N]=P(1-P)/N$ becomes smaller with $N\to \infty$

We have to suppose (A) that the region $R$ is large enough such that the $N$ points are sufficient to get a sharply peaked binomial distribution.

In that case, for $N\to \infty$ we have that $K/N\approx P$, and also:

$$K\approx NP$$

Now we suppose (B) that the region $R$ is sufficiently small such that the probability $p(x)$ is constant for $x\in R$. In this case, we have that:

$$P\approx p(x)V$$

where $V$ is the volume of the region $R$.

Observe that the two assumptions (A) and (B) are contradictory. Which is a limitation of the KDE methods. By assuming this, we can use the two results to derive:

$$p(x)=\frac{K}{NV}$$

The KDE methods are based on this result. They usually fix $V$ and then get $K$ from the $N$ observations available:

def kde(points: List, region: Region) -> float:
    """
    Generic KDE estimator
    """
    K = 0
    N = len(points)
    for point in points:
        if point in region:
            K += 1
    return K / (N * region.volume)

Parzen estimator

In the function above, we suppose that we can evaluate if the point lies inside the region, but this depends on the shape of the region we use. The Parzen estimator uses a hypercube centered at the point $x$ where we want to evaluate the PDF $p(x)$.

We now introduce a kernel function $k(u)$, also called Parzen Window, defined as follows:

$$k(u)=\{\begin{array}{ll}1& \text{ if }|u_i|\le 1/2,i=1,\dots ,D\\ 0& \text{ otherwise}\end{array}$$

To know if a point $x_i$ lies inside the hypercube of side $h$ centered on $x$, we need to scale the point coordinates using this formula:

$$k\left(\frac{x-x_i}{h}\right)$$

In this way, we can compute $K$ by:

$$K=\sum _{n=1}^{N}k\left(\frac{x-x_n}{h}\right)$$

and since the volume of an hypercube of $D$ dimensions and of edge size $h$ is $V=h^D$, we can replace $K$ and $V$ in the $p(x)$ equation and get:

$$p(x)=\frac{1}{N}\sum _{n=1}^N\frac{1}{h^D}\cdot k\left(\frac{x-x_n}{h}\right)$$

class ParzenWindow(Region):

    def __init__(self, h: float, origin: Point):
        self.h = h
        self.origin = origin
        self.volume = h ** len(origin)

    def __contains__(self, point):
        return self._k( (point - self.origin) / self.h )

    def _k(self, u):        
        return int(all([ u_n < 0.5 for u_n in u ]))

Gaussian Kernel

The Parzen Estimator has discontinuities (boundaries of the cubes). A smoother kernel function is the Gaussian:

$$p(x)=\frac{1}{N}\sum _{n=1}^N\frac{1}{\sqrt{2\pi }{h}^2}\exp \left\{-\frac{||x-x_n|{|}^2}{2h^2}\right\}$$

Where $h$ now represents the standard deviation of the Gaussian components.

In general, we can use every kernel function $k(u)$ that satisfies:

$$k(u)\ge 0\int k(u)du=1$$

Nearest-neighbour methods

A density distribution could have locations with smaller structures and location with bigger structures. KDE methods use the same $h$ parameter everywhere, and this may not be suitable for every location of the distribution.

Nearest-neighbour methods address this problem: recall the general problem of non-parametric density estimation, let's fix K and use the data to find an appropiate value for V.

We consider a small sphere centered on the point $x$ at which we wish to estimate the density $p(x)$, we allow the radius of the sphere to grow until it contains precisely $K$ data points. $V$ will be the volume of that sphere.

The model produced by kNN is not a true density model since the integral over all spaces diverges.

kNN can also be used as a classifier. An interesting property of kNN with $K=1$ is that in the limit $N\to \infty$, the error rate is never more than twice the minimum achievable error rate of an optimal classifier.

Chapter 2 Exercises

Solution to exercise 2.1

Prop n.1.

$$\sum _{x=0}^{1}p(x\mid \mu )=\mu +(1-\mu )=1$$

Prop n.2

$$\mathbb{E}[x]=\sum _{x=0}^{1}xp(x\mid \mu )=\sum _{x=0}^{1}x\cdot {\mu }^x(1-\mu {)}^{1-x}=1\cdot {\mu }^1(1-\mu {)}^{1-1}=\mu$$

Prop n.3

$$\begin{array}{rl}var[x]& =\sum _{x=0}^{1}p(x)(x-\mu {)}^2\\ & =p(0){\mu }^2+p(1)(1-\mu {)}^2\\ & =(1-\mu ){\mu }^2+\mu (1-\mu {)}^2\\ & =\mu ((1-\mu )\mu +(1-\mu {)}^2)\\ & =\mu (\mu -{\mu }^2+1-2\mu +{\mu }^2)\\ & =\mu (1-\mu )\end{array}$$

Solution to exercise 2.30

We know that

$$\mathbb{E}[z]=R^{-1}\left[\begin{array}{c}\Lambda \mu -A^{T}Lb\\ Lb\end{array}\right]$$

and that

$$cov[z]=R^{-1}=\left[\begin{array}{cc}{\Lambda }^{-1}& {\Lambda }^{-1}{A}^T\\ A{\Lambda }^{-1}& L^{-1}+A{\Lambda }^{-1}{A}^T\end{array}\right]$$

So by replacing $R^{-1}$ we get

$$\mathbb{E}[z]=\left[\begin{array}{cc}{\Lambda }^{-1}& {\Lambda }^{-1}{A}^T\\ A{\Lambda }^{-1}& L^{-1}+A{\Lambda }^{-1}{A}^T\end{array}\right]\left[\begin{array}{c}\Lambda \mu -A^{T}Lb\\ Lb\end{array}\right]$$

which results in

$$\left[\begin{array}{c}\mu -{\Lambda }^{-1}{A}^{T}Lb+{\Lambda }^{-1}{A}^{T}Lb\\ A\mu -A{\Lambda }^{-1}{A}^{T}Lb+b+A{\Lambda }^{-1}{A}^{T}Lb\end{array}\right]=\left[\begin{array}{c}\mu \\ A\mu +b\end{array}\right]$$

Solution to exercise 2.31

There are various approaches to compute the marginal distribution $p(y)$ where $y=x+z$ and $x\sim \mathcal{N}(x\mid {\mu }_x,{\Sigma }_x)$, $z\sim \mathcal{N}(z\mid {\mu }_z,{\Sigma }_z)$.

The first approach came to my mind after a video from 3Blue1Brown, that demonstrates exactly that in this case $p(y)=p(x)\ast p(z)$, where $\ast$ is the convolution operator.

$$p(y)=\int p_x(y-t)p_z(t)dt$$

In this way, we consider every way to obtain $y$ from the sum $x+z$. This solution has been adopted by Tommy Odland in his solutions.

But there is a simpler way to do this. Let's consider the conditional distribution $p(y\mid x)$. Since $x$ is fixed, and $y=x+z$, the only variability is up to $z$. We can define this as

$$p(y\mid x)=\mathcal{N}(y\mid {\mu }_z+x,{\Sigma }_z)$$

We can now compare the obtained results with expressions 2.99 and 2.100:

$$\mu ={\mu }_x{\Lambda }^{-1}={\Sigma }_{x}A=Ib={\mu }_z{L}^{-1}={\Sigma }_{z}x=x$$

And using results from 2.109 and 2.110 we know:

$$\mathbb{E}[y]=A\mu +b=I{\mu }_x+{\mu }_z={\mu }_x+{\mu }_z$$

and

$$cov[y]=L^{-1}+A{\Lambda }^{-1}{A}^T={\Sigma }_z+{\Sigma }_x$$

Linear models for regression

The goal of regression is to predict the value of one or more continuous targets variables $t$ given the value of a D-dimensional vector $x$ of input variables.

By linear models we mean that the model is a linear function of the adjustable parameters. E.g. the polynomial curve fitting algorith builds a linear model. The simplest form of linear regression models are also linear functions of the input variables.

We get a more useful class of functions by taking linear combinations of a fixed set of nonlinear functions of the input variables, known as basis functions. Such models are linear functions of the parameters (which gives simple analytical properties) and yet can be nonlinear with respect to the input variables.

Given a dataset of $N$ observations $\{x_n\}$ where $n=1,\dots ,N$, together with the corresponding target values $\{t_n\}$, the goal is to predict $t$ for a new value of $x$.

• Simple approach: Find an appropiate function $y(x)\approx t$
• General approach: Find the predictive distribution $p(t\mid x)$ to get the uncertainty of a prediction

Linear Basis Function Models

The simplest linear model involves a linear combination of the input variables, also called linear regression:

$$y(x,w)=w_0+w_1{x}_1+w_2{x}_2+\cdots +w_D{x}_D$$

This is:

• A linear function of the parameters (good for optimization)
• A linear function of the inputs (bad for expressiveness)

Extend the concept of linear combination to combine fixed nonlinear functions of the input:

$$y(x,w)=w_0+\sum _{j=1}^{M-1}{w}_j{\varphi }_j(x)$$

where ${\varphi }_j(x)$ are known as basis functions. By having $M-1$ components, the total number of parameters is $M$ (consider the bias).

If we consider ${\varphi }_0(x)=1$, then we can write:

$$y(x,w)=\sum _{j=0}^{M-1}{w}_j{\varphi }_j(x)$$

These linear models are:

• A linear function of the parameters (good for optimization)
• A nonlinear function of the inputs (good for expressiveness)

Polynomials are basis functions of the form ${\varphi }_j(x)=x^j$. The problem with polynomial is that they are global functions: a change in a region of the input space affects all the other regions.

Other choices for the basis functions are:

$${\varphi }_j(x)=\exp \{-\frac{(x-{\mu }_j{)}^2}{2s^2}\}$$

Which is the Gaussian basis function. In this case, ${\mu }_j$ is the location in the input space, and $s$ is the scale. This function doesn't have a probabilistic interpretation.

Another possibility is the Sigmoidal basis function:

$${\varphi }_j(x)=\sigma \left(\frac{x-{\mu }_j}{s}\right)\sigma (a)=\frac{1}{1+\exp (-a)}$$

Where $\sigma$ is the logistic function, but we can also use the tanh function.

We can use also Fourier basis functions such that the the regression function is an expansion of sinusoidal functions at a certain frequency. Combining basis functions localized in both space and frequency leads to a class of functions known as wavelets.

Maximum likelihood and least squares

Let $y$ be a deterministic function such that $y(x)\approx t$, let $ϵ\sim \mathcal{N}(ϵ\mid 0,{\beta }^{-1})$ be a random Gaussian variable with precision $\beta$. We assume that the target variable $t$ is given by: