```python
import scipy.io as sio
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.sparse import csc_matrix
data = sio.loadmat('heightWeight.mat')
hwdata = data['heightWeightData']
y = hwdata[:,0]
x = hwdata[:,1:]
maleNdx = np.where(y==1)
femaleNdx = np.where(y==2)
```

### Gaussian Models
The probability density function for an Multivariate Normal in D dimensions is defined as follows,
$$\mathcal{N}(\mathbf{x}|\mathbf{\mu},\mathbf{\Sigma})\triangleq \frac{1}{(2\pi)^{\frac{D}{2}}|\mathbf{\Sigma}|^{\frac{1}{2}}}\exp\bigg[-\frac{1}{2}(\mathbf{x}-\mathbf{\mu})^\top\mathbf{\Sigma}^{-1}(\mathbf{x}-\mathbf{\mu})\bigg]$$

The expression inside the exponent is the **Mahalanobis distance** between a data vector $\mathbf{x}$ and the mean vector $\mathbf{\mu}$. To gain a better understanding of this quantity by performing an **eigendecomposition** of $\mathbf{\Sigma}=\mathbf{U\Lambda U^\top}$, where $\mathbf{U}$ is an orthonormal matrix of eigenvectors satisfying $\mathbf{U^\top U=I}$, and $\mathbf{\Lambda}$ is a diagonal matrix of eigenvalues. That is, 

$$\mathbf{\Sigma}^{-1}=\mathbf{U^{-\top}\Lambda^{-1}U^{-1}}=\mathbf{U\Sigma^{-1}U^\top}=\sum^D_{i=1}\frac{1}{\lambda_i}\mathbf{u_iu_i^\top}$$
where $\mathbf{u_i}$ is the $i^{th}$ column of $\mathbf{U}$, containing the $i^{th}$ eigenvector. Hence, the Mahalanobis distance can be rewritten as,

\begin{align}
(\mathbf{x}-\mathbf{\mu})^\top\mathbf{\Sigma}^{-1}(\mathbf{x}-\mathbf{\mu})&= (\mathbf{x}-\mathbf{\mu})^\top\sum^D_{i=1}\frac{1}{\lambda_i}\mathbf{u_iu_i^\top}(\mathbf{x}-\mathbf{\mu}) \\
&= \sum^D_{i=1}\frac{1}{\lambda}(\mathbf{x}-\mathbf{\mu})^\top \mathbf{u_i}\mathbf{u_i}^\top(\mathbf{x}-\mathbf{\mu}) = \sum^D_{i=1}\frac{y_i^2}{\lambda_i}
\end{align}
where $y_i\triangleq \mathbf{u_i}^\top(\mathbf{x}-\mathbf{\mu})$. This matches the elliptical equation in 2d, 

$$\frac{y_1^2}{\lambda_1}+\frac{y_2^2}{\lambda_2}=1$$



#### MLE for an MVN
The parameters are given by,

\begin{align}
\mathbf{\mu} &= \frac{1}{N}\sum^N_{i=1}\mathbf{x}_i \triangleq \bar{\mathbf{x}} \\
\mathbf{\Sigma} &= \frac{1}{N}\sum^N_{i=1}(\mathbf{x}_i-\bar{\mathbf{x}})(\mathbf{x}_i-\bar{\mathbf{x}})^\top = \frac{1}{N}\sum^N_{i=1}\mathbf{x_i}\mathbf{x_i}^\top - \bar{\mathbf{x}}\bar{\mathbf{x}}^\top
\end{align}


```python
mx = x[maleNdx,:]
mm = np.mean(mx,axis=1); mv = np.cov(np.transpose(mx.squeeze()))
fx = x[femaleNdx,:]
fm = np.mean(fx,axis=1); fv = np.cov(fx.squeeze().T)
```


```python
def gaussPlot2d(mu,Sigma):
    n = 100; mu = mu.T; k = np.sqrt(6)
    d,u = np.linalg.eig(Sigma)
    t = np.linspace(0,2*math.pi,n)
    xy = np.vstack((np.cos(t),np.sin(t))) 
    w = np.matmul((k*np.matmul(u,np.sqrt(np.diag(d)))),xy)
    z = np.tile(mu,(1,n))+w
    plt.plot(z[0,:],z[1,:])
```


```python
plt.scatter(x[maleNdx,0],x[maleNdx,1],c='b')
gaussPlot2d(mm,mv)
plt.scatter(x[femaleNdx,0],x[femaleNdx,1],c='r')
gaussPlot2d(fm,fv)
```

### Proof
The log-likelihood is 

$$\mathcal{l}(\mathbf{\mu},\mathbf{\Sigma}) = \log p(\mathcal{D}|\mathbf{\mu},\mathbf{\Sigma}) = \frac{N}{2}\log|\mathbf{\Lambda}|-\frac{1}{2}\sum^N_{i=1}(\mathbf{x_i}-\mathbf{\mu})^\top \mathbf{\Lambda}(\mathbf{x_i}-\mathbf{\mu})$$
where $\mathbf{\Lambda}=\mathbf{\Sigma}^{-1}$ is the precision matrix.

Using the substitution $\mathbf{y_i}=\mathbf{x_i}-\mathbf{\mu}$, 

\begin{align}
\frac{\partial}{\partial \mathbf{\mu}}(\mathbf{x_i}-\mathbf{\mu})^\top\mathbf{\Sigma}^{-1}(\mathbf{x_i}-\mathbf{\mu})=\frac{\partial}{\partial\mathbf{\mu}}\mathbf{y_i}^\top\mathbf{\Sigma}^{-1}\mathbf{y_i}
\end{align}


With $\frac{\partial}{\partial \mathbf{a}}(\mathbf{a}^\top\mathbf{A}\mathbf{a})=(\mathbf{A}+\mathbf{A}^\top)\mathbf{a}$, and chain rule $\frac{\partial}{\partial \mathbf{\mu}}=\frac{\partial}{\partial\mathbf{y_i}}\frac{\partial\mathbf{y_i}}{\partial\mathbf{\mu}}$,

$$\frac{\partial}{\partial\mathbf{\mu}}\mathbf{y_i}^\top\mathbf{\Sigma^{-1}}\mathbf{y_i} = \frac{\partial \mathbf{y_i}}{\partial\mathbf{\mu}}\frac{\partial}{\partial \mathbf{y_i}}\mathbf{y_i}^\top\mathbf{\Sigma}^{-1}\mathbf{y_i}=-(\mathbf{\Sigma}^{-1}+\mathbf{\Sigma}^{-\top})\mathbf{y_i}=-2\mathbf{\Sigma}^{-1}\mathbf{y_i}$$

Hence

\begin{align}
\frac{\partial}{\partial\mathbf{\mu}}\mathcal{l}(\mathbf{\mu},\mathbf{\Sigma})&=-\frac{1}{2}\sum^N_{i=1}-2\mathbf{\Sigma}^{-1}(\mathbf{x_i}-\mathbf{\mu})=\mathbf{\Sigma}^{-1}\sum^N_{i=1}(\mathbf{x_i}-\mathbf{\mu})=0 \\
\hat{\mathbf{\mu}}\sum^N_{i=1} &= N\hat{\mathbf{\mu}} = \sum^N_{i=1}\mathbf{x_i} \Longrightarrow\hat{\mu} = \frac{1}{N}\sum^N_{i=1}\mathbf{x_i}=\mathbf{\bar{x}}
\end{align}

Given the trace-trick $\mathbf{x^\top Ax}=Tr(\mathbf{x^\top Ax})=Tr(\mathbf{xx^\top A})=Tr(\mathbf{Axx^\top})$,

\begin{align}
\mathcal{l}(\mathbf{\Lambda})&=\frac{N}{2}\log|\mathbf{\Lambda}|-\frac{1}{2}\sum_i Tr[(\mathbf{x_i}-\mathbf{\mu})(\mathbf{x_i}-\mathbf{\mu})^\top\mathbf{\Lambda}] \\
&= \frac{N}{2}\log|\mathbf{\Lambda}|-\frac{1}{2}Tr|\mathbf{S_\mu}\mathbf{\Lambda}|
\end{align}
where $$\mathbf{S_\mu} \triangleq \sum^N_{i=1}(\mathbf{x_i}-\mathbf{\mu})(\mathbf{x_i}-\mathbf{\mu})^\top$$ is the scatter matrix centered on $\mathbf{\mu}$.

Given $\frac{\partial}{\partial\mathbf{A}}\log|\mathbf{A}| = \mathbf{A}^{-\top}$ and $\frac{\partial}{\partial\mathbf{A}}Tr(\mathbf{BA})=\mathbf{B}^\top$, taking derivatives of this expression with respect to $\mathbf{\Lambda}$ yields,

\begin{align}
\frac{\partial\mathcal{l}(\mathbf{\Lambda})}{\partial\mathbf{\Lambda}}&=\frac{N}{2}\frac{\partial}{\partial\mathbf{\Lambda}}\log|\mathbf{\Lambda}|-\frac{1}{2}\frac{\partial}{\partial\mathbf{\Lambda}}Tr|\mathbf{S_\mu\Lambda}| = \frac{N}{2}\mathbf{\Lambda}^{-\top}-\frac{1}{2}\mathbf{S_\mu^\top} = 0 \\
\mathbf{\Lambda}^{-\top} &= \mathbf{\Lambda}^{-1} = \mathbf{\Sigma}=\frac{1}{N}\mathbf{S_\mu}
\end{align}

so
$$\mathbf{\hat{\Sigma}}=\frac{1}{N}\sum^N_{i=1}(\mathbf{x_i}-\mathbf{\mu})(\mathbf{x_i}-\mathbf{\mu})^\top$$

efb6bb320fcfd5ed2bf47905fceeb6899865341e
ipynb
Jupyter Notebook
Machine Learning A Probabilistic Perspective/4GaussianModels/F4.2/.ipynb_checkpoints/4.2gaussHeightWeight-checkpoint.ipynb
zcemycl/ProbabilisticPerspectiveMachineLearning
8291bc6cb935c5b5f9a88f7b436e6e42716c21ae
MIT
2019-11-20T10:20:29.000Z
2021-11-09T11:15:23.000Z
Machine Learning A Probabilistic Perspective/4GaussianModels/F4.2/.ipynb_checkpoints/4.2gaussHeightWeight-checkpoint.ipynb
zcemycl/ProbabilisticPerspectiveMachineLearning
8291bc6cb935c5b5f9a88f7b436e6e42716c21ae
MIT
Machine Learning A Probabilistic Perspective/4GaussianModels/F4.2/.ipynb_checkpoints/4.2gaussHeightWeight-checkpoint.ipynb
zcemycl/ProbabilisticPerspectiveMachineLearning
8291bc6cb935c5b5f9a88f7b436e6e42716c21ae
MIT
2020-05-27T03:56:38.000Z
2021-05-02T13:15:42.000Z
Qwen/Qwen-72B
1. YES
2. YES