$\newcommand{\ind}[1]{\left[#1\right]}$

# Model

- Data set
$$ {\cal D} = \{ x_1, \dots  x_N \} $$
- Model with parameter $\theta$
$$ p(\cal D | \theta) $$




## Maximum Likelihood

- Maximum Likelihood (ML)
$$ \theta^{\text{ML}} = \arg\max_{\theta} \log p({\cal D} | \theta) $$
- Predictive distribution
$$ p(x_{N+1} |  {\cal D} ) \approx  p(x_{N+1} |  \theta^{\text{ML}})  $$

## Maximum Aposteriori

- Prior
$$ p(\theta) $$

- Maximum a-posteriori (MAP) : Regularised Maximum Likelihood
$$
\theta^{\text{MAP}} = \arg\max_{\theta} \log p({\cal D} | \theta) p(\theta)
$$

- Predictive distribution
$$ p(x_{N+1} |  {\cal D} ) \approx  p(x_{N+1} |  \theta^{\text{MAP}})  $$

## Bayesian Learning

- We treat parameters on the same footing as all other variables
- We integrate over unknown parameters rather than using point estimates (remember the many-dice example)
 - Self-regularisation, avoids overfitting
 - Natural setup for online adaptation
 - Model selection


- Predictive distribution
\begin{eqnarray}
p(x_{N+1} ,  {\cal D} ) &=& \int d\theta \;\; p(x_{N+1} |  \theta) p( {\cal D}| \theta) p(\theta)  \\
 &=& \int d\theta \;\; p(x_{N+1}|  \theta) p( {\cal D}, \theta)   \\
 &=& \int d\theta \;\; p(x_{N+1}|  \theta) p(  \theta| {\cal D}) p({\cal D})   \\
 &=&  p({\cal D}) \int d\theta \;\; p(x_{N+1}|  \theta) p(  \theta| {\cal D})  \\
p(x_{N+1} |  {\cal D} ) &=& \int d\theta \;\; p(x_{N+1} |  \theta) p(\theta | {\cal D}) 
\end{eqnarray}

The interpretation is that past data provides an 'update' to the recent prior to be used for the current prediction.

- Bayesian learning is just inference ...




# Learning the parameter of a (possibly fake) coin

Suppose we have a coin, flipped several times independently. A vague question one can ask is if one can predict the outcome of the next flip.

It depends. If we already know that the coin is fair, there is nothing that we can learn from past data and indeed the future flips are independent of the previous flips. However, if we don't know the probability of the coin, we could estimate the parameter from past data to create a better prediction. Mathematically, the model is identical to  



Here, $\theta$ is the parameter of the coin.

## Maximum Likelihood Estimation

We observe the outcome of $N$ coin flips $\{x^{(n)}\}_{n=1\dots N}$ where $x^{(n)} \in \left\{0,1\right\}$. The model is a Bernoulli distribution with parameter $\pi = (\pi_0, \pi_1)$. We have $\pi_0 = 1 - \pi_1$ where $0 \leq \pi_1 \leq 1$. 

\begin{eqnarray}
x^{(n)} & \sim & p(x|\pi) = (1-\pi_1)^{1-x^{(n)} } \pi_1^{x^{(n)} }
\end{eqnarray}

The loglikelihood is 

\begin{eqnarray}
{\cal L}(\pi_1) & = & \sum_{n=1}^N (1- x^{(n)}) \log (1 - \pi_1) + \sum_{n=1}^N x^{(n)} \log (\pi_1)  \\
& = & \log (1 - \pi_1) \sum_{n=1}^N (1- x^{(n)})  + \log (\pi_1) \sum_{n=1}^N x^{(n)}  
\end{eqnarray}

We define the number of $0$'s  
\begin{eqnarray}
c_0 = \sum_{n=1}^N (1- x^{(n)})
\end{eqnarray}
and $1$'s as
\begin{eqnarray}
c_1 = \sum_{n=1}^N  x^{(n)}
\end{eqnarray}

\begin{eqnarray}
{\cal L}(\pi_1) & = & \log (1 - \pi_1) c_0  + \log (\pi_1) c_1  
\end{eqnarray}

We compute the gradient 
\begin{eqnarray}
\frac{\partial}{\partial \pi_1} {\cal L}(\pi_1) & = & - \frac{c_0}{1 - \pi_1} + \frac{c_1}{\pi_1}  = 0 
\end{eqnarray}

The solution is quite predictable
\begin{eqnarray}
\pi_1 & = &\frac{c_1}{c_0 + c_1}  = \frac{c_1}{N}  
\end{eqnarray}

## Maximum A-posteriori estimation

We need a prior over the probability parameter. One choice is the beta distribution

\begin{eqnarray}
p(\pi_1) & = & \mathcal{B}(\pi_1; \alpha, \beta) =  \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta) } \pi_1^{\alpha-1} (1-\pi_1)^{\beta-1}
\end{eqnarray}

The log joint ditribution of data is
\begin{eqnarray}
\log p(X, \pi_1) & = & \log p(\pi_1) + \sum_{n=1}^N \log p(x^{(n)}|\pi_1) \\
& = & \log \Gamma(\alpha + \beta) -\log \Gamma(\alpha) - \log \Gamma(\beta) \\
& & + (\alpha-1) \log \pi_1 + (\beta-1) \log(1-\pi_1) \\
& & + c_1 \log (\pi_1)  + c_0 \log (1 - \pi_1) \\
& = & \log \Gamma(\alpha + \beta) -\log \Gamma(\alpha) - \log \Gamma(\beta) \\
& & + (\alpha + c_1 -1) \log \pi_1 + (\beta + c_0 -1) \log(1-\pi_1) 
\end{eqnarray}

The gradient is 

\begin{eqnarray}
\frac{\partial}{\partial \pi_1} \log p(X, \pi_1)  & = & - \frac{\beta + c_0 -1}{1 - \pi_1} + \frac{\alpha + c_1 -1}{\pi_1}  = 0 
\end{eqnarray}

We can solve for the parameter.
\begin{eqnarray}
\pi_1 (\beta + c_0 -1) & = & (1 - \pi_1) (\alpha + c_1 -1) \\ 
\pi_1 \beta + \pi_1 c_0 - \pi_1 & = & \alpha  + c_1  - 1 - \pi_1 \alpha - \pi_1 c_1 + \pi_1 \\ 
\pi_1  & = & \frac{\alpha - 1  + c_1}{\alpha + \beta  - 2 + c_0 + c_1}    \\ 
\end{eqnarray}

When the prior is flat, i.e., when $\alpha = \beta = 1$, MAP and ML solutions coincide.

## Full Bayesian inference

We infer the posterior

\begin{eqnarray}
p(\pi_1| X) & = & \frac{p(\pi_1, X)}{p(X)} 
\end{eqnarray}


The log joint density is 
\begin{eqnarray}
\log p(X, \pi_1) & = & \log \Gamma(\alpha + \beta) -\log \Gamma(\alpha) - \log \Gamma(\beta) \\
& & + (\alpha + c_1 -1) \log \pi_1 + (\beta + c_0 -1) \log(1-\pi_1) 
\end{eqnarray}

At this stage, we may try to evaluate the integral 
$$
p(X)  =  \int d\pi_1 p(X, \pi_1) 
$$

Rather than trying to evaluate this integral directly, a simple approach is known as 'completing the square': we add an substract terms to obtain an expression that corresponds to a known, normalized density. This typically involves adding and substracting an expression that will make us identify a normalized density. 

\begin{eqnarray}
\log p(X, \pi_1) & = & \log \Gamma(\alpha + \beta) -\log \Gamma(\alpha) - \log \Gamma(\beta) \\
& & - \log \Gamma(\alpha + \beta + c_0 + c_1) + \log \Gamma(\alpha + c_1) + \log \Gamma(\beta + c_0) \\
& & + \log \Gamma(\alpha + \beta + c_0 + c_1) - \log \Gamma(\alpha + c_1) - \log \Gamma(\beta + c_0) \\
& & + (\alpha + c_1 -1) \log \pi_1 + (\beta + c_0 -1) \log(1-\pi_1) \\
& = & \log \Gamma(\alpha + \beta) -\log \Gamma(\alpha) - \log \Gamma(\beta) \\
& & - \log \Gamma(\alpha + \beta + c_0 + c_1) + \log \Gamma(\alpha + c_1) + \log \Gamma(\beta + c_0) \\
& & + \log \mathcal{B}(\alpha + c_1, \beta + c_0) \\
& = & \log p(X) + \log p(\pi_1| X)
\end{eqnarray}

From the resulting expression, taking the exponent on both sides we see that 
\begin{eqnarray}
p(\pi_1| X) & = & \mathcal{B}(\alpha + c_1, \beta + c_0) \\
p(X) & = & \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \frac{\Gamma(\alpha + c_1)\Gamma(\beta + c_0)}{\Gamma(\alpha + \beta + c_0 + c_1)}
\end{eqnarray}

Predictive distribution: Let $a=\alpha + c_1$ and $b=\beta+c_0$

\begin{eqnarray}
\int d\pi_1 p(x|\pi_1) p(\pi_1| X) & = & \int d\pi_1 \mathcal{BE}(x; \pi_1) \mathcal{B}(\pi_1; a, b) \\
 & = & \int d\pi_1 \pi_1^{\ind{x=1}}(1-\pi_1)^{\ind{x=0}}\mathcal{B}(a, b) \\
\end{eqnarray}

\begin{eqnarray}
p(x) & = & \frac{\Gamma(a + b)}{\Gamma(a)\Gamma(b)} \frac{\Gamma(a + \ind{x=1})\Gamma(b + \ind{x=0})}{\Gamma(b + a + 1)}
\end{eqnarray}

\begin{eqnarray}
p(x=1) & = & \frac{\Gamma(a + b)}{\Gamma(a)\Gamma(b)} \frac{\Gamma(a + 1)\Gamma(b)}{\Gamma(b + a + 1)} \\
& = & \frac{\Gamma(a + b)}{\Gamma(a)\Gamma(b)} \frac{a\Gamma(a)\Gamma(b)}{(b+a)\Gamma(b + a)} \\
& = & \frac{a}{a+b}
\end{eqnarray}


## Alternative Derivation
Alternatively, we may directly write
\begin{eqnarray}
p(X, \pi_1) & = & \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)}   \pi_1^{(\alpha + c_1 -1)} (1-\pi_1)^{(\beta + c_0 -1)} 
\end{eqnarray}

\begin{eqnarray}
p(X) &=& \int d\pi_1 p(X, \pi_1)  =  \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \int d\pi_1  \pi_1^{(\alpha + c_1 -1)} (1-\pi_1)^{(\beta + c_0 -1)} 
\end{eqnarray}


From the definition of the beta distribution, we can arrive at the 'formula' for the integral 
\begin{eqnarray}
1 &=& \int d\pi \mathcal{B}(\pi; a, b) \\
& = & \int d\pi \frac{\Gamma(a + b)}{\Gamma(a)\Gamma(b)} \pi^{(a -1)} (1-\pi)^{(b -1)} \\
\frac{\Gamma(a)\Gamma(b)}{\Gamma(a + b)} & = & \int d\pi \pi^{(a -1)} (1-\pi)^{(b -1)}
\end{eqnarray}
Just substitute $a = \alpha + c_1$ and $b = \beta + c_0$

## An Approximation
For large $x$, we have the following approximation
\begin{eqnarray}
\log \Gamma(x + a) - \log \Gamma(x) & \approx & a \log(x) \\
\Gamma(x + a)  & \approx & \Gamma(x) x^a \\
\end{eqnarray}

When $c_0$ and $c_1$ are large, we obtain:

\begin{eqnarray}
p(X) & \approx & \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \frac{\Gamma(c_1)\Gamma(c_0)c_0^{\beta}c_1^{\alpha}}{\Gamma(c_0 + c_1)(c_0+c_1)^{\alpha + \beta}}
\end{eqnarray}

Let $\hat{\pi}_1 = c_1/(c_0+c_1)$ and $N = c_0 + c_1$, we have
\begin{eqnarray}
p(X) & \approx &  \frac{\Gamma(c_1)\Gamma(c_0)}{\Gamma(c_0 + c_1)} (1-\hat{\pi}_1) \hat{\pi}_1 \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} (1-\hat{\pi}_1)^{\beta-1}\hat{\pi}_1^{\alpha-1}
\end{eqnarray}

### Illustration: Bayesian update of a Beta Distribution 


```python
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import gammaln

def log_beta_pdf(x, a, b):
    return - gammaln(a) - gammaln(b) + gammaln(a+b) + np.log(x)*(a-1) + np.log(1-x)*(b-1) 

x = np.arange(0.01,1,0.01)

a = 1
b = 1
c_0 = 1
c_1 = 1
N = c_0 + c_1

pi_ML = c_1/N

plt.figure(figsize=(8,4))
plt.plot(x, np.exp(log_beta_pdf(x, a, b)), 'b')
plt.plot(x, np.exp(log_beta_pdf(x, a+c_1, b+c_0)), 'r')
yl = plt.gca().get_ylim()
plt.plot([pi_ML, pi_ML], yl , 'k:')
plt.legend(["Prior $\cal B$ $(a={}, b={})$".format(a,b), "Posterior $\cal B$ $(a={}, b={})$".format(a+c_1, b+c_0)], loc="best")
plt.show()
```

### Illustration: Learning from a sequence of coin flips


```python
%matplotlib inline
import matplotlib.pyplot as plt
from scipy.special import gammaln
import numpy as np
import scipy.special as sps

#savefigs = 

def log_beta_pdf(x, a, b):
    
    y = - gammaln(a) - gammaln(b) + gammaln(a+b) + np.log(x)*(a-1) + np.log(1-x)*(b-1) 
    idx = np.where(x == 0)
    if a==1:
        y[idx] = - gammaln(a) - gammaln(b) + gammaln(a+b)
    elif a<1:
        y[idx] = np.Inf
    else:
        y[idx] = - np.Inf
    
    idx = np.where(x == 1)
    if b==1:
        y[idx] = - gammaln(a) - gammaln(b) + gammaln(a+b)
    elif b<1:
        y[idx] = np.Inf
    else:
        y[idx] = - np.Inf
    return y

a = 1
b = 1

xx = [1,1,1,1,1,0,1,1, 1,0,1,0,1,1,1,1,0, 1,1,1,1,1,1,1]

p = np.arange(0.0,1.01,0.01)

c = [0,0]
N = 0
plt.figure(figsize=(5,3))
plt.plot(p, np.exp(log_beta_pdf(p, a+c[1], b+c[0])), 'r')
plt.xticks([0,1])
plt.yticks([])
plt.xlim([-0.1,1.1])
plt.ylim([0,6])
pi_ML =  (a+c[1])/(a+b+N)
yl = plt.gca().get_ylim()
plt.plot([pi_ML, pi_ML], yl , 'k:')

plt.xlabel('$\lambda$')
#plt.savefig('/Users/cemgil/Dropbox/tex/cam/talks/cmpe547/beta{n}.eps'.format(n=N), bbox_inches='tight')
plt.show()

for x in xx:
    c[x] += 1
    N += 1
    plt.figure(figsize=(5,3))
    plt.plot(p, np.exp(log_beta_pdf(p, a+c[1], b+c[0])), 'r')

    
    pi_ML =  (a+c[1])/(a+b+N)
    pi_mode = (a+c[1]-1)/float(a+b+N-2)
    tmp = str(int(a+c[1]-1))+'/'+str(a+b+N-2)
    plt.xticks([0,pi_mode],('0',tmp))
    #plt.xticks([0,pi_mode])
    plt.yticks([])
    plt.xlim([-0.1,1.1])
    plt.ylim([0,6])
    
    yl = plt.gca().get_ylim()
    plt.plot([pi_ML, pi_ML], yl , 'k:')
    plt.plot([pi_mode, pi_mode], yl , 'b:')
    
    plt.xlabel('$\lambda$')
    #plt.savefig('/Users/cemgil/Dropbox/tex/cam/talks/cmpe547/beta{n}.eps'.format(n=N), bbox_inches='tight')
    plt.show()

```

# Finding if a Coin is Fair or Fake 

We consider the folowing problem: Given a sequence of coin tosses $X = \{x^{(n)}\}_{n=1\dots N}$, determine if the coin is fair or fake.

This can be cast as a model selection problem:

\begin{eqnarray}
\pi_1|m & \sim & \left\{ \begin{array}{cc} \delta(\pi_1 - 0.5) & m = 0\\ \mathcal{B}(\pi; a, b) & m = 1 \end{array} \right.
\end{eqnarray}
For $n = 1\dots N$
\begin{eqnarray}
x^{(n)}| \pi_1 & \sim & \mathcal{BE}(x; \pi_1)
\end{eqnarray}

This model defines the following:
- The indicator $m$, that denotes if the coin is fake,
- What a fake coin is: a fake coin is one that has an arbitrary probability $\pi_1$ between $0$ and $1$. 
- What a fair coin is: a fair coin has $\pi_1 = 0.5$

We need to calculate the marginal likelihoods for $m=0$ and $m=1$
\begin{eqnarray}
p(X| m) & = & \int d\pi_1 p(X | \pi_1) p(\pi_1|m)
\end{eqnarray}

###### Not Fake
\begin{eqnarray}
p(X| m) & = & \int d\pi_1 p(X| \pi_1) \delta(\pi_1 - 0.5) \\
& = & \prod_{n=1}^N \left(\frac{1}{2}\right)^{x^{(n)}} \left(\frac{1}{2}\right)^{1-x^{(n)}} = \frac{1}{2^N}
\end{eqnarray}

###### Fake

\begin{eqnarray}
p(X| m) & = & \int d\pi_1 p(\pi_1; a, b) \prod_{n=1}^{N} p(x^{(n)}| \pi_1)  \\
& = &  \int d\pi_1 \left(\prod_{n=1}^N \left(1-\pi_1\right)^{1-x^{(n)}} \pi_1^{x^{(n)}} \right) \mathcal{B}(\pi; a, b) \\
& = & \frac{\Gamma(a + b)}{\Gamma(a)\Gamma(b)} \int d\pi_1 \left(1-\pi_1\right)^{c_0+a-1} \pi_1^{c_1+b-1} \\
& = & \frac{\Gamma(a + b)}{\Gamma(a)\Gamma(b)} \frac{\Gamma(c_0+a)\Gamma(c_1+b)}{\Gamma(c_0 + c_1 +a + b)}
\end{eqnarray}

The log-odds is the ratio of marginal likelihoods

$$
l(X) = \log\left( \frac{p(X|m = \text{Fair})}{p(X|m = \text{Fake})} \right)
$$

If $l(X)>0$, we may conclude that the coin is fair and biased when $l<0$.


```python
import numpy as np
import scipy.special as sps

def log_odds(c_0, c_1, a, b):
    # Total number of tosses
    N = c_0 + c_1
    
    M_fair = N*np.log(0.5)
    M_fake = sps.gammaln(a+b) - sps.gammaln(a) - sps.gammaln(b) +  sps.gammaln(c_0+a) + sps.gammaln(c_1+b) - sps.gammaln(N+a + b) 
    return M_fair - M_fake

# Number of Zeros observed
c_0 = 6
# Number of Ones
c_1 = 1

# Prior
a = 1
b = 1


print('log_odds = ', log_odds(c_0, c_1, a, b) )

```

    log_odds =  -0.826678573184



```python
a = 1
b = 1
N = 10

l = np.zeros(N+1)

for c in range(0,N+1):
    l[c] = log_odds(N-c, c, a, b)

plt.plot(range(0,N+1), l, 'o')
plt.plot(range(0,N+1), np.zeros(N+1), 'k:')
ax = plt.gca()
ax.set_xlabel('Number of ones $c_1$')
ax.set_ylabel('log-odds $l(X)$')
plt.show()
```

We can visualize the region where we would decide that the coin is fake by plotting the points where the log-odds is negative.


```python
a = 1
b = 1

for N in range(1, 25):

    l = np.zeros(N+1)

    for c in range(0,N+1):
        l[c] = log_odds(N-c, c, a, b)
    
    
    idx = np.where( np.array(l)<0 )
    p = np.arange(0,N+1)/N
    plt.plot(N*np.ones_like(p), p, '.k',markersize=4)    
    plt.plot(N*np.ones_like(p[idx]), p[idx], '.r',markersize=20)
    

ax = plt.gca()
ax.set_ylim((0,1))
ax.set_xlabel('$N$')
ax.set_ylabel('$c_1/N$')
plt.show()
```

# Estimation of a Categorical distribution

## Maximum Likelihood Estimation

We observe a dataset $\{x^{(n)}\}_{n=1\dots N}$. The model for a single observation is a categorical distribution with parameter $\pi = (\pi_1, \dots, \pi_S)$ where 

\begin{eqnarray}
x^{(n)} & \sim & p(x|\pi) = \prod_{s=1}^{S} \pi_s^{\ind{s = x^{(n)}}}
\end{eqnarray}
where $\sum_s \pi_s  = 1$.

The loglikelihood of the entire dataset is

\begin{eqnarray}
{\cal L}(\pi_1,\dots,\pi_S) & = & \sum_{n=1}^N\sum_{s=1}^S \ind{s = x^{(n)}} \log \pi_s
\end{eqnarray}
This is a constrained optimisation problem.
Form the Lagrangian
\begin{eqnarray}
\Lambda(\pi, \lambda) & = & \sum_{n=1}^N\sum_{s'=1}^S \ind{s' = x^{(n)}} \log \pi_{s'}  + \lambda \left( 1 - \sum_{s'} \pi_{s'} \right ) \\
\frac{\partial \Lambda(\pi, \lambda)}{\partial \pi_s} & = & \sum_{n=1}^N \ind{s = x^{(n)}} \frac{1}{\pi_s} - \lambda = 0 \\
\pi_s & = & \frac{\sum_{n=1}^N \ind{s = x^{(n)}}}{\lambda}
\end{eqnarray}

We solve for $\lambda$
\begin{eqnarray}
1 & = & \sum_s \pi_s = \frac{\sum_{s=1}^S \sum_{n=1}^N \ind{s = x^{(n)}}}{\lambda} \\
\lambda & = & \sum_{s=1}^S \sum_{n=1}^N \ind{s = x^{(n)}} =  \sum_{n=1}^N 1 = N
\end{eqnarray}

Hence
\begin{eqnarray}
\pi_s & = & \frac{\sum_{n=1}^N \ind{s = x^{(n)}}}{N}
\end{eqnarray}


```python
# %load template_equations.py
from IPython.display import display, Math, Latex, HTML
import notes_utilities as nut
from importlib import reload
reload(nut)
Latex('$\DeclareMathOperator{\trace}{Tr}$')

L = nut.pdf2latex_dirichlet(x=r'\pi', a=r'a',N=r'S', i='s')

display(HTML(nut.eqs2html_table(L)))
```


<table width=100%><tr><td align="center">\begin{eqnarray}\mathcal{D}(\pi_{1:S}; a_{1:S} )\end{eqnarray}</td><td>\mathcal{D}(\pi_{1:S}; a_{1:S} )</td></tr><tr><td align="center">\begin{eqnarray}\frac{\Gamma(\sum_{s} a_{s})}{\prod_{s} \Gamma(a_{s})} \prod_{{s}=1}^{S} {\pi}_{s}^{a_{s} - 1} \end{eqnarray}</td><td>\frac{\Gamma(\sum_{s} a_{s})}{\prod_{s} \Gamma(a_{s})} \prod_{{s}=1}^{S} {\pi}_{s}^{a_{s} - 1} </td></tr><tr><td align="center">\begin{eqnarray}\exp\left(\log{\Gamma(\sum_{s} a_{s})} - {\sum_{s} \log \Gamma(a_{s})} + \sum_{{s}=1}^{S} (a_{s} - 1) \log{\pi}_{s} \right)\end{eqnarray}</td><td>\exp\left(\log{\Gamma(\sum_{s} a_{s})} - {\sum_{s} \log \Gamma(a_{s})} + \sum_{{s}=1}^{S} (a_{s} - 1) \log{\pi}_{s} \right)</td></tr><tr><td align="center">\begin{eqnarray}\log{\Gamma(\sum_{s} a_{s})} - {\sum_{s} \log \Gamma(a_{s})} + \sum_{{s}=1}^{S} (a_{s} - 1) \log{\pi}_{s} \end{eqnarray}</td><td>\log{\Gamma(\sum_{s} a_{s})} - {\sum_{s} \log \Gamma(a_{s})} + \sum_{{s}=1}^{S} (a_{s} - 1) \log{\pi}_{s} </td></tr></table>


#### Maximum A-Posteriori Estimation

$$
\pi \sim \mathcal{D}(\pi_{1:S}; a_{1:S} )
$$
where $\sum_s \pi_s  = 1$. For $n = 1\dots N$
\begin{eqnarray}
x^{(n)} & \sim & p(x|\pi) = \prod_{s=1}^{S} \pi_s^{\ind{s = x^{(n)}}}
\end{eqnarray}
$X = \{x^{(1)},\dots,x^{(N)} \}$

The posterior is 
\begin{eqnarray}
\log p(\pi_{1:S}| X) & =^+ & \log p(\pi_{1:S}, X) \\
& = & +\log{\Gamma(\sum_{s} a_{s})} - {\sum_{s} \log \Gamma(a_{s})} + \sum_{{s}=1}^{S} (a_{s} - 1) \log{\pi}_{s}  + \sum_{s=1}^S\sum_{n=1}^N \ind{s = x^{(n)}} \log \pi_s \\
 & =^+ &   \sum_{s=1}^S \left(a_s - 1 +  \sum_{n=1}^N \ind{s = x^{(n)}}\right) \log \pi_s 
\end{eqnarray}
Finding the parameter vector $\pi_{1:S}$ that maximizes the posterior density is a constrained optimisation problem. After omitting constant terms that do not depend on $\pi$, we form the Lagrangian
\begin{eqnarray}
\Lambda(\pi, \lambda) & = & \sum_{s=1}^S \left(a_s - 1 +  \sum_{n=1}^N \ind{s = x^{(n)}}\right) \log \pi_s   + \lambda \left( 1 - \sum_{s'} \pi_{s'} \right ) \\
\frac{\partial \Lambda(\pi, \lambda)}{\partial \pi_s} & = & \left(a_s - 1 + \sum_{n=1}^N \ind{s = x^{(n)}}\right) \frac{1}{\pi_s} - \lambda = 0 \\
\pi_s & = & \frac{a_s - 1 + \sum_{n=1}^N \ind{s = x^{(n)}}}{\lambda}
\end{eqnarray}


We solve for $\lambda$
\begin{eqnarray}
1 & = & \sum_s \pi_s = \frac{- S + \sum_{s=1}^S \left( a_s + \sum_{n=1}^N \ind{s = x^{(n)}  \right) }}{\lambda} \\
\lambda & = &  N - S + \sum_{s=1}^S a_s
\end{eqnarray}

Setting the count of observations equal to $s$ as $C_s \equiv \sum_{n=1}^N \ind{s = x^{(n)}}$, we obtain 

Hence
\begin{eqnarray}
\pi_s & = & \frac{C_s + a_s - 1}{N + \sum_{s=1}^S a_s - S}
\end{eqnarray}



#### Full Bayesian Inference

Setting the count of observations equal to $s$ as $C_s \equiv \sum_{n=1}^N \ind{s = x^{(n)}}$, we obtain 


The posterior is 
\begin{eqnarray}
\log p(\pi_{1:S}, X) & = & \log{\Gamma(\sum_{s} a_{s})} - {\sum_{s} \log \Gamma(a_{s})} +  \sum_{s=1}^S \left(\left(a_s + \sum_{n=1}^N \ind{s = x^{(n)}}\right) - 1  \right) \log \pi_s \\
& =  & \log{\Gamma(\sum_{s} a_{s})} - {\sum_{s} \log \Gamma(a_{s})}\\
& &  + \sum_{s=1}^S (a_s + C_s - 1) \log \pi_s \\
 & & + \log{\Gamma(\sum_{s} (a_{s} + C_s) )} - {\sum_{s} \log \Gamma(a_{s} + C_s)} \\
 & & - \log{\Gamma(\sum_{s} (a_{s} + C_s) )} + {\sum_{s} \log \Gamma(a_{s} + C_s)} \\
 & = & \log \mathcal{D}(\pi_{1:S}, \alpha_{1:S} ) + \log p(X) 
\end{eqnarray}

\begin{eqnarray}
\log p(X) & = & \log{\Gamma(\sum_{s} a_{s})} - {\sum_{s} \log \Gamma(a_{s})} - \log{\Gamma(\sum_{s} (a_{s} + C_s) )} + {\sum_{s} \log \Gamma(a_{s} + C_s)}
\end{eqnarray}




```python
from IPython.display import display, Math, Latex, HTML
import html_utils as htm
wd = '65px'
L = [[htm.TableCell('', width=wd), htm.TableCell('$x_2=1$', width=wd), htm.TableCell('$x_2=j$', width=wd), htm.TableCell('$x_2=S_2$', width='80px')],
     [r'$x_1=1$',r'$C(1,1)$',r'',r'$C(1,S_2)$'],
     [r'$x_1=i$',r'',r'$C(i,j)$',r''],
     [r'$x_1=S_1$',r'$C(S_1,1)$',r'',r'$C(S_1,S_2)$']]

t = htm.make_htmlTable(L)
#display(HTML(str(t)))
print(str(t))
```

## Are $x_1$ and $x_2$ independent?
Suppose we observe a dataset of $(x_1, x_2)$ pairs where
$x_1 \in \{1,\dots,S_1\}$ and
$x_2 \in \{1,\dots,S_2\}$.


| $x_1^{(1)}$ | $x_2^{(1)}$ | 
| --- | --- |
| $\vdots$ | $\vdots$ | 
| $x_1^{(n)}$ | $x_2^{(n)}$ | 
| $\vdots$ | $\vdots$ | 
| $x_1^{(N)}$ | $x_2^{(N)}$ | 


We are given the counts of observations where $x_1 = i$ while $x_2 = j$. These counts can be stored as an array, that is known as a contingency table
$
C(i,j) = \sum_{n=1}^N \ind{x_1^{(n)} = i}\ind{x_2^{(n)} = j}
$

<TABLE cellpadding="4" style="border: 1px solid #000000; border-collapse: collapse;" border="1">
 <TR>
  <TD width="65px">&nbsp;</TD>
  <TD width="65px">$x_2=1$</TD>
  <TD width="65px">$x_2=j$</TD>
  <TD width="80px">$x_2=S_2$</TD>
 </TR>
 <TR>
  <TD>$x_1=1$</TD>
  <TD>$C(1,1)$</TD>
  <TD>&nbsp;</TD>
  <TD>$C(1,S_2)$</TD>
 </TR>
 <TR>
  <TD>$x_1=i$</TD>
  <TD>&nbsp;</TD>
  <TD>$C(i,j)$</TD>
  <TD>&nbsp;</TD>
 </TR>
 <TR>
  <TD>$x_1=S_1$</TD>
  <TD>$C(S_1,1)$</TD>
  <TD>&nbsp;</TD>
  <TD>$C(S_1,S_2)$</TD>
 </TR>
</TABLE>

The goal is deciding if the random variables are independent or dependent.

||$x_2$||
|-|-|-|-|
|$x_1$ |$3$|$5$|$9$|
| |$7$|$9$|$17$|

### Independent model $M_1$

\begin{equation}
p(x_1, x_2) = p(x_1) p(x_2) 
\end{equation}

\begin{eqnarray}
\pi_1 & \sim & \mathcal{D}(\pi_1; a_1) \\
\pi_2 & \sim & \mathcal{D}(\pi_2; a_2) \\
x_1^{(n)} & \sim & \mathcal{C}(x_1; \pi_1) \\
x_2^{(n)} & \sim & \mathcal{C}(x_2; \pi_2) 
\end{eqnarray}

We let 

* $C_1(i) = \sum_j C(i,j)$ 
* $C_2(j) = \sum_i C(i,j)$
* $A_1 = \sum_i a_{1}(i)$ 
* $A_2 = \sum_j a_{2}(j)$

The marginal likelihood can be found as 
\begin{eqnarray}
\log p(X|M_1) & = & \log{\Gamma(\sum_{i} a_{1}(i))} - {\sum_{i} \log \Gamma(a_{1}(i)} - \log{\Gamma(\sum_{i} (a_{1}(i) + \sum_{j} C(i,j)) )} + {\sum_{i} \log \Gamma(a_{1}(i) + \sum_{j} C(i,j))} \\
& & +\log{\Gamma(\sum_{j} a_{2}(j)} - {\sum_{j} \log \Gamma(a_{2}(j))} - \log{\Gamma(\sum_{j} (a_{2}(j) + \sum_{i} C(i,j)) )} + {\sum_{j} \log \Gamma(a_{2}(j) + \sum_{i} C(i,j))} \\
& = & \log{\Gamma(A_1)} - \sum_{i} \log \Gamma(a_{1}(i)) - \log{\Gamma(A_1+ N)} + {\sum_{i} \log \Gamma(a_{1}(i) + C_1(i))} \\
& & + \log{\Gamma(A_2)} - {\sum_{j} \log \Gamma(a_{2}(j))} - \log{\Gamma(A_2 + N )} + {\sum_{j} \log \Gamma(a_{2}(j) +  C_2(j))} \\
\end{eqnarray}

### Dependent model $M_2$
\begin{equation}
p(x_1, x_2)
\end{equation}

$\pi_{1,2}$ is a $S_1 \times S_2$ matrix where the joint distribution of entries is Dirichlet $\mathcal{D}(\pi_{1,2}; a_{1,2})$ with $S_1 \times S_2$ parameter matrix $a_{1,2}$. Then, the probability that $p(x_1 = i, x_2 = j|\pi_{1,2}) = \pi_{1,2}(i,j)$.

\begin{eqnarray}
\pi_{1,2} & \sim & \mathcal{D}(\pi_{1,2}; a_{1,2}) \\
(x_1, x_2)^{(n)} & \sim & \mathcal{C}((x_1,x_2); \pi_{1,2}) \\
\end{eqnarray}

\begin{eqnarray}
\log p(X|M_2) & = & \log{\Gamma(A_{1,2})} - {\sum_{i,j} \log \Gamma(a_{1,2}(i,j))} - \log{\Gamma(A_{1,2}+ N)} + {\sum_{i,j} \log \Gamma(a_{1,2}(i,j) + C(i,j))} 
\end{eqnarray}



### Dependent model $M_3$
\begin{equation}
p(x_1, x_2) = p(x_1) p(x_2|x_1) 
\end{equation}

\begin{eqnarray}
\pi_1 & \sim & \mathcal{D}(\pi_1; a_1) \\
\pi_{2,1} & \sim & \mathcal{D}(\pi_2; a_2) \\
\vdots \\
\pi_{2,S_1} & \sim & \mathcal{D}(\pi_2; a_2) \\
x_1^{(n)} & \sim & \mathcal{C}(x_1; \pi_1) \\
x_2^{(n)} & \sim & \mathcal{C}(x_2; \pi_{2}(x_1^{(n)},:)) 
\end{eqnarray}

\begin{eqnarray}
\log p(x_1^{(1:N)}|\pi_1) & = & \sum_n \sum_i \sum_j \ind{x_1^{(n)} = i} \ind{x_2^{(n)} = j} \log \pi_{1}(i)  =  \sum_i \sum_j C(i,j) \log \pi_{1}(i)
\end{eqnarray}
\begin{eqnarray}
\log p(x_2^{(1:N)}|\pi_2, x_1^{(1:N)} ) & = & \sum_n \sum_i \sum_j \ind{x_1^{(n)} = i} \ind{x_2^{(n)} = j} \log \pi_{2}(i,j) =  \sum_i \sum_j C(i,j) \log \pi_{2}(i,j)
\end{eqnarray}



\begin{eqnarray}
\log p(\pi_1) & = & \log{\Gamma(\sum_{i} a_{1}(i))} - {\sum_{i} \log \Gamma(a_{1}(i))} + \sum_{{i}=1}^{S_1} (a_{1}(i) - 1) \log{\pi_1}(i) 
\end{eqnarray}
\begin{eqnarray}
\log p(\pi_2) & = & \sum_i \left(\log{\Gamma(\sum_{j} a_{2}(i,j))} - {\sum_{j} \log \Gamma(a_{2}(i,j))} + \sum_{{j}=1}^{S_2} (a_{2}(i,j) - 1) \log{\pi_2}(i,j) \right)
\end{eqnarray}

The joint distribution is
\begin{eqnarray}
\log p(X, \pi| M_2)&= & \log{\Gamma(\sum_{i} a_{1}(i))} - {\sum_{i} \log \Gamma(a_{1}(i))} + \sum_{{i}=1}^{S_1} (a_{1}(i) + C_1(i) - 1) \log{\pi_1}(i)  \\
& & + \sum_i \left(\log{\Gamma(\sum_{j} a_{2}(i,j))} - {\sum_{j} \log \Gamma(a_{2}(i,j))} + \sum_{{j}=1}^{S_2} (a_{2}(i,j) + C(i,j) - 1) \log{\pi_2}(i,j) \right) 
\end{eqnarray}

We will assume $a_2(i,j) = a_2(i',j)$ for all $i$ and $i'$.
\begin{eqnarray}
\log p(X| M_2)&= & \log{\Gamma(\sum_{i} a_{1}(i))} - {\sum_{i} \log \Gamma(a_{1}(i))} - \log{\Gamma(\sum_{{i}} a_{1}(i) + C_1(i) )} + \sum_{{i}} \log \Gamma(a_{1}(i) + C_1(i))  \\
& & + \sum_i \left( \log\Gamma(\sum_{j} a_{2}(i,j)) - \sum_{j} \log \Gamma(a_{2}(i,j)) - \log\Gamma( \sum_{j} a_{2}(i,j) + C(i,j)) + \sum_j \log\Gamma( a_{2}(i,j) + C(i,j) ) \right) \\
& = & \log{\Gamma(A_1)} - {\sum_{i} \log \Gamma(a_{1}(i))} - \log{\Gamma(A_1+ N)} + {\sum_{i} \log \Gamma(a_{1}(i) + C_1(i))} \\
& & + \sum_i \left( \log{\Gamma(A_2)} - {\sum_{j} \log \Gamma(a_{2}(j))} - \log{\Gamma(A_2 + C_1(i) )} + {\sum_{j} \log \Gamma(a_{2}(j) + C(i,j))} \right) \\
\end{eqnarray}



### Dependent model $M_3b$
The derivation is similar and corresponds to the factorization:

\begin{equation}
p(x_1, x_2) = p(x_2) p(x_1|x_2) 
\end{equation}





```python
import numpy as np
from notes_utilities import randgen, log_sum_exp, normalize_exp, normalize

import scipy as sc
from scipy.special import gammaln

#C = np.array([[3,5,9],[7,9,17]])
C = 1*np.array([[1,1,3],[1,1,7]])

#C = np.array([[0,1,1],[1,0,2]])

C_i = np.sum(C, axis=1)
C_j = np.sum(C, axis=0)

N = np.sum(C)
S_1 = C.shape[0]
S_2 = C.shape[1]

#M1 Parameter
M1 = {'a_1': np.ones(S_1), 'a_2': np.ones(S_2), 'A_1': None, 'A_2': None}
M1['A_1'] = np.sum(M1['a_1'])
M1['A_2'] = np.sum(M1['a_2'])
#p(x_1) p(x_2)
log_marglik_M1 = gammaln(M1['A_1']) - np.sum(gammaln(M1['a_1'])) - gammaln(M1['A_1'] + N) + np.sum(gammaln(M1['a_1'] + C_i)) \
            + gammaln(M1['A_2']) - np.sum(gammaln(M1['a_2'])) - gammaln(M1['A_2'] + N) + np.sum(gammaln(M1['a_2'] + C_j))

# p(x_1, x_2)
M2 = {'a_12': np.ones((S_1,S_2)), 'A_12':None}
M2['A_12'] = np.sum(M2['a_12'])
log_marglik_M2 = gammaln(M2['A_12']) - np.sum(gammaln(M2['a_12'])) - gammaln(M2['A_12'] + N) + np.sum(gammaln(M2['a_12'] + C)) 

    
    
M3 = {'a_1': S_2*np.ones(S_1), 'a_2': np.ones(S_2), 'A_1': None, 'A_2': None}
M3['A_1'] = np.sum(M3['a_1'])
M3['A_2'] = np.sum(M3['a_2'])

#p(x_1) p(x_2|x_1)
log_marglik_M3 = gammaln(M3['A_1']) - np.sum(gammaln(M3['a_1'])) - gammaln(M3['A_1'] + N) + np.sum(gammaln(M3['a_1'] + C_i)) 
for i in range(S_1):
    log_marglik_M3 +=  gammaln(M3['A_2']) - np.sum(gammaln(M3['a_2'])) - gammaln(M3['A_2'] + C_i[i]) + np.sum(gammaln(M3['a_2'] + C[i,:]))

# Beware the prior parameters
M3b = {'a_1': np.ones(S_1), 'a_2': S_1*np.ones(S_2), 'A_1': None, 'A_2': None}
M3b['A_1'] = np.sum(M3b['a_1'])
M3b['A_2'] = np.sum(M3b['a_2'])

#p(x_2) p(x_1|x_2)
log_marglik_M3b = gammaln(M3b['A_2']) - np.sum(gammaln(M3b['a_2'])) - gammaln(M3b['A_2'] + N) + np.sum(gammaln(M3b['a_2'] + C_j)) 
for j in range(S_2):
    log_marglik_M3b +=  gammaln(M3b['A_1']) - np.sum(gammaln(M3b['a_1'])) - gammaln(M3b['A_1'] + C_j[j]) + np.sum(gammaln(M3b['a_1'] + C[:,j]))


print('M1:', log_marglik_M1)
print('M2:', log_marglik_M2)
print('M3:', log_marglik_M3)
print('M3b:', log_marglik_M3b)

print('Log Odds, M1-M2')
print(log_marglik_M1 - log_marglik_M2)
print(normalize_exp([log_marglik_M1, log_marglik_M2]))

print('Log Odds, M1-M3')
print(log_marglik_M1 - log_marglik_M3)
print(normalize_exp([log_marglik_M1, log_marglik_M3]))

print('Log Odds, M1-M3b')
print(log_marglik_M1 - log_marglik_M3b)
print(normalize_exp([log_marglik_M1, log_marglik_M3b]))

```

    M1: -23.7979581523
    M2: -24.2354716141
    M3: -24.2354716141
    M3b: -24.2354716141
    Log Odds, M1-M2
    0.437513461821
    [ 0.60766638  0.39233362]
    Log Odds, M1-M3
    0.437513461821
    [ 0.60766638  0.39233362]
    Log Odds, M1-M3b
    0.437513461821
    [ 0.60766638  0.39233362]


Conceptually $M_2$, $M_3$ and $M_3b$ should have the same marginal likelihood score, as the dependence should not depend on how we parametrize the conditional probability tables. However, this is dependent on the choice of the prior parameters. 

How should the prior parameters of $M_2$, $M_3$ and $M_3b$ be chosen such that we get the same evidence score?

The models 
$M_2$, $M_3$ and $M_3b$ are all equivlent, if the prior parameters are chosen appropriately. For $M_2$ and $M_3$, we need to take $a_1(i) = \sum_j a_{1,2}(i,j)$.

For example, if in $M_2$, the prior parameters $a_{1,2}$ are chosen as
\begin{eqnarray}
a_{1,2} & = & \left(\begin{array}{ccc} 1 & 1 & 1\\ 1 & 1 & 1 \end{array} \right)
\end{eqnarray}

we need to choose in model $M_3$
\begin{eqnarray}
a_{1} & = & \left(\begin{array}{c} 3 \\ 3  \end{array} \right)
\end{eqnarray}
\begin{eqnarray}
a_{2} & = & \left(\begin{array}{ccc} 1 & 1 & 1\\ 1 & 1 & 1 \end{array} \right)
\end{eqnarray}

and in model $M_3b$
\begin{eqnarray}
a_{2} & = & \left(\begin{array}{ccc} 2 & 2 & 2  \end{array} \right)
\end{eqnarray}
\begin{eqnarray}
a_{1} & = & \left(\begin{array}{ccc} 1 & 1 & 1\\ 1 & 1 & 1 \end{array} \right)
\end{eqnarray}

This is due to fact that the marginals of a Dirichlet distribution are also Dirichlet. In particular,
if a probability vector $x$ and corresponding parameter vector $a$ are partitioned as $x = (x_\iota, x_{-\iota})$ 
and $a = (a_\iota, a_{-\iota})$, the Dirichlet distribution
$$
\mathcal{D}(x_\iota, x_{-\iota}; a_\iota, a_{-\iota})
$$
has marginals 
$$
\mathcal{D}(X_{\iota}, X_{-\iota}; A_{\iota}, A_{-\iota})
$$
where $X_\iota = \sum_{i \in \iota}  x_i$ and $X_{-\iota} = \sum_{i \in -\iota}  x_i$, where $A_\iota = \sum_{i \in \iota}  a_i$ and $A_{-\iota} = \sum_{i \in -\iota}  a_i$. The script below verifies that the marginals are indeed distributed according to this formula.



```python
S_1 = 2
S_2 = 10
M = S_1*S_2

a = np.ones(M)

N = 100
P = np.random.dirichlet(a, size=N)
A = np.zeros((N,S_1))
B = np.zeros((N*S_1,S_2))

for n in range(N):
    temp = P[n,:].reshape((S_1,S_2))
    A[n,:] = np.sum(temp, axis=1)
    for i in range(S_1):
        B[(n*S_1+i),:] = temp[i,:]/A[n,i]

import pylab as plt

plt.hist(A[:,0],bins=20)
plt.gca().set_xlim([0,1])
#plt.plot(B[:,0],B[:,1],'.')
plt.show()

P2 = np.random.dirichlet(S_2*np.ones(S_1), size=N)
plt.hist(P2[:,0],bins=20)
plt.gca().set_xlim([0,1])
plt.show()
```

Are two given histograms drawn from the same distribution?

$[3,5,12,4]$

$[8, 14, 31, 14]$

## Visualizing the Dirichlet Distribution

[http://blog.bogatron.net/blog/2014/02/02/visualizing-dirichlet-distributions/]


```python
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
from functools import reduce
from scipy.special import gammaln

corners = np.array([[0, 0], [1, 0], [0.5, 0.75**0.5]])
triangle = tri.Triangulation(corners[:, 0], corners[:, 1])

refiner = tri.UniformTriRefiner(triangle)
trimesh = refiner.refine_triangulation(subdiv=4)


# Mid-points of triangle sides opposite of each corner
midpoints = [(corners[(i + 1) % 3] + corners[(i + 2) % 3]) / 2.0 \
             for i in range(3)]

def xy2bc(xy, tol=1.e-3):
    '''Converts 2D Cartesian coordinates to barycentric.'''
    s = [(corners[i] - midpoints[i]).dot(xy - midpoints[i]) / 0.75 \
         for i in range(3)]
    return np.clip(s, tol, 1.0 - tol)

class Dirichlet(object):
    def __init__(self, alpha):
        self._alpha = np.array(alpha)
        self._coef = gammaln(np.sum(self._alpha)) - np.sum(gammaln(self._alpha))
    def log_pdf(self, x):
        return self._coef + np.sum(np.log(x)*(self._alpha - 1)) 
    def pdf(self, x):
        '''Returns pdf value for `x`.'''
        return np.exp(self.log_pdf(x))
        
def draw_pdf_contours(dist, nlevels=200, subdiv=8, **kwargs):

    refiner = tri.UniformTriRefiner(triangle)
    trimesh = refiner.refine_triangulation(subdiv=subdiv)
    pvals = [dist.pdf(xy2bc(xy)) for xy in zip(trimesh.x, trimesh.y)]

    plt.tricontourf(trimesh, pvals, nlevels, **kwargs)
    plt.axis('equal')
    plt.xlim(0, 1)
    plt.ylim(0, 0.75**0.5)
    plt.axis('off')
```


```python
draw_pdf_contours(Dirichlet([1, 3, 1]))
```


```python
draw_pdf_contours(Dirichlet([1.99, 3.99, 10.99]))
```

# 6-faced die with repeated labels

Consider a die where the numbers on each face are labeled, possibly with repetitions, from the set $1\dots 6$. 
A 'normal' die has labels $1,2,3,4,5,6$ but we allow other labelings, for example as $1,1,3,5,5,5$ or $1,1,1,1,1,6$.

Can we construct a method to find how the die has been labeled from a sequence of outcomes?  



# Does my data have a single cluster or are there two clusters?

We observe a dataset of $N$ points and want to decide if there are one or two clusters. For example, the below dataset, when visualized seems to suggest two clusters; however the separation is not very clear; perhaps a single component might have been also sufficient. How can we derive a procedure that leads to a resonable answer in this ambigious situation?



One principled approach is based on Bayesian model selection. Our approach will be describing two alternative generative models for data. Each generative model will reflect our assumption what it means to have clusters. In other words, we should describe two different procedures: how to generate a dataset with a single cluster, or a datset that has two clusters. Once we have a description of each generative procedure, we may hope to convert our qualitative question (how many clusters?) into a well defined computational procedure. 

Each generative procedure will be a different probability model and we will compute the marginal posterior distribution conditioned on an observed dataset.

The single cluster model will have a single cluster center, denoed as $\mu$. Once $\mu$ is generated, each observation is generated by a Gaussian distribution with variance $R$, centered around $\mu$.

Model $M =1$: Single Cluster
\begin{eqnarray}
\mu & \sim & {\mathcal N}(\mu; 0, P) \\
x_i | \mu & \sim & {\mathcal N}(x; \mu, R)
\end{eqnarray}

The parameter $P$ denotes a natural range for the mean, $R$ denotes the variance of data, the amount of spread around the mean. To start simple, we will assume that these parameters are known; we will able to relax this assumption later easily.

Below, we show an example where each data point is a scalar.


```python
# Parameters
P = 100
R = 10

# Number of datapoints
N = 5

mu = np.random.normal(0, np.sqrt(P))
x  = np.random.normal(mu, np.sqrt(R), size=(N))

plt.figure(figsize=(10,1))
plt.plot(mu, 0, 'r.')
plt.plot(x, np.zeros_like(x), 'x')
ax = plt.gca()
ax.set_xlim(3*np.sqrt(P)*np.array([-1,1]))
ax.set_ylim([-0.1,0.1])
ax.axis('off')
plt.show()
```

Model $M = 2$: Two Clusters
\begin{eqnarray}
\mu_0 & \sim & {\mathcal{N}}(\mu; 0, P) \\
\mu_1 & \sim & {\mathcal{N}}(\mu; 0, P) \\
c_i & \sim & {\mathcal{BE}}(r; 0.5) \\
x_i | \mu_0, \mu_1, c_i & \sim & {\mathcal N}(x; \mu_0, R)^{1-c_i} {\mathcal N}(x; \mu_1, R)^{c_i}
\end{eqnarray}

The parameter $P$ denotes a natural range for both means, $R$ denotes the variance of data in each cluster. The variables $r_i$ are the indicators that show the assignment of each datapoint to one of the clusters.


```python
# Parameters
P = 100
R = 2

# Number of datapoints
N = 10

# Number of clusters
M = 2

mu = np.random.normal(0, np.sqrt(P), size=(M))
c = np.random.binomial(1, 0.5, size=N)
x = np.zeros(N)
for i in range(N):
    x[i]  = np.random.normal(mu[c[i]], np.sqrt(R))

plt.figure(figsize=(10,1))
#plt.plot(mu, np.zeros_like(mu), 'r.')
plt.plot(x, np.zeros_like(x), 'x')
ax = plt.gca()
ax.set_xlim(3*np.sqrt(P)*np.array([-1,1]))
ax.set_ylim([-0.1,0.1])
ax.axis('off')
plt.show()
```

#### Extension
We dont know the component variances


#### Combining into a single model

Capture

Recapture

### Change point
Coin switch

Coal Mining Data
 Single Change Point
 Multiple Change Point

# Bivariate Gaussian model selection


Suppose we are given a dataset $X = \{x_1, x_2, \dots, x_N \}$ where $x_n \in \mathbb{R}^K$ for $n=1 \dots N$ and consider two competing models:

- ### Model $m=1$:
$\newcommand{\diag}{\text{diag}}$

- Observation
\begin{eqnarray}
\text{for}\; n=1\dots N&& \\
x_n| s_{1:K} & \sim &  \mathcal{N}\left(x; 0, \diag\{s_1, \dots, s_K\}\right) = \mathcal{N}\left(x; 0, \left(\begin{array}{ccc} s_1 & 0 & 0\\0 & \ddots & 0 \\ 0 & \dots & s_K \end{array} \right) \right)
\end{eqnarray}

$$
p(x_n| s_{1:K})  = 
\prod_{k=1}^K \mathcal{N}(x_{k,n}; 0, s_k) 
$$

$$
p(X|s_{1:K} ) = \prod_{n=1}^N p(x_n| s_{1:K})  = 
\prod_{n=1}^N \prod_{k=1}^K \mathcal{N}(x_{k,n}; 0, s_k) 
$$


- Prior
\begin{eqnarray}
\text{for}\; k=1\dots K&& \\
s_k & \sim & \mathcal{IG}(s_k; \alpha, \beta) 
\end{eqnarray}



- ### Model $m=2$:

-- Observation

\begin{eqnarray}
x_n \sim \mathcal{N}(x_n; 0, \Sigma)=\left|{ 2\pi \Sigma } \right|^{-1/2} \exp\left(-\frac12 {x_n}^\top {\Sigma}^{-1} {x_n} \right)=\exp\left( -\frac{1}{2}\trace {\Sigma}^{-1} {x_n}{x_n}^\top -\frac{1}{2}\log \left|2{\pi}\Sigma\right|\right)
\end{eqnarray}

$$
{\cal IW}(\Sigma; 2a, 2B)  =  \exp( - (a + (k+1)/2) \log |\Sigma| - \trace B\Sigma^{-1}  - \log\Gamma_k(a) + a\log |B|) \\
$$


```python
# %load template_equations.py
from IPython.display import display, Math, Latex, HTML
import notes_utilities as nut
from importlib import reload
reload(nut)
Latex('$\DeclareMathOperator{\trace}{Tr}$')

#L = nut.pdf2latex_gauss(x=r's', m=r'\mu',v=r'v')
#L = nut.pdf2latex_mvnormal(x=r's', m=r'\mu',v=r'\Sigma')
L = nut.pdf2latex_mvnormal(x=r'x_n', m=0,v=r'\Sigma')
#L = nut.pdf2latex_gamma(x=r'x', a=r'a',b=r'b')
#L = nut.pdf2latex_invgamma(x=r'x', a=r'a',b=r'b')
#L = nut.pdf2latex_beta(x=r'\pi', a=r'\alpha',b=r'\beta')

eq = L[0]+'='+L[1]+'='+L[2]
display(Math(eq))
display(Latex(eq))
display(HTML(nut.eqs2html_table(L)))
```


$$\mathcal{N}(x_n; 0, \Sigma)=\left|{ 2\pi \Sigma } \right|^{-1/2} \exp\left(-\frac12 {x_n}^\top {\Sigma}^{-1} {x_n} \right)=\exp\left(  -\frac{1}{2}\trace {\Sigma}^{-1} {x_n}{x_n}^\top  -\frac{1}{2}\log \left|2{\pi}\Sigma\right|\right)$$



\mathcal{N}(x_n; 0, \Sigma)=\left|{ 2\pi \Sigma } \right|^{-1/2} \exp\left(-\frac12 {x_n}^\top {\Sigma}^{-1} {x_n} \right)=\exp\left(  -\frac{1}{2}\trace {\Sigma}^{-1} {x_n}{x_n}^\top  -\frac{1}{2}\log \left|2{\pi}\Sigma\right|\right)



<table width=100%><tr><td align="center">\begin{eqnarray}\mathcal{N}(x_n; 0, \Sigma)\end{eqnarray}</td><td>\mathcal{N}(x_n; 0, \Sigma)</td></tr><tr><td align="center">\begin{eqnarray}\left|{ 2\pi \Sigma } \right|^{-1/2} \exp\left(-\frac12 {x_n}^\top {\Sigma}^{-1} {x_n} \right)\end{eqnarray}</td><td>\left|{ 2\pi \Sigma } \right|^{-1/2} \exp\left(-\frac12 {x_n}^\top {\Sigma}^{-1} {x_n} \right)</td></tr><tr><td align="center">\begin{eqnarray}\exp\left(  -\frac{1}{2}\trace {\Sigma}^{-1} {x_n}{x_n}^\top  -\frac{1}{2}\log \left|2{\pi}\Sigma\right|\right)\end{eqnarray}</td><td>\exp\left(  -\frac{1}{2}\trace {\Sigma}^{-1} {x_n}{x_n}^\top  -\frac{1}{2}\log \left|2{\pi}\Sigma\right|\right)</td></tr><tr><td align="center">\begin{eqnarray}  -\frac{1}{2}\trace {\Sigma}^{-1} {x_n}{x_n}^\top  -\frac{1}{2}\log \left|2{\pi}\Sigma\right|\end{eqnarray}</td><td>  -\frac{1}{2}\trace {\Sigma}^{-1} {x_n}{x_n}^\top  -\frac{1}{2}\log \left|2{\pi}\Sigma\right|</td></tr></table>


Computing the marginal likelihood

* #### Model 1

\begin{eqnarray}
p(X| m=1) & = & \int d{s_{1:K}} p(X|s_{1:K}) p(s_{1:K}) = 
\end{eqnarray}



```python
# %load template_equations.py
from IPython.display import display, Math, Latex, HTML
import notes_utilities as nut
from importlib import reload
reload(nut)
Latex('$\DeclareMathOperator{\trace}{Tr}$')

#L = nut.pdf2latex_gauss(x=r's', m=r'\mu',v=r'v')
L = nut.pdf2latex_mvnormal(x=r'x_t', m=r'(Ax_{t-1})',v=r'Q')
#L = nut.pdf2latex_mvnormal(x=r's', m=0,v=r'I')
#L = nut.pdf2latex_gamma(x=r'x', a=r'a',b=r'b')
#L = nut.pdf2latex_invgamma(x=r'x', a=r'a',b=r'b')
#L = nut.pdf2latex_beta(x=r'\pi', a=r'\alpha',b=r'\beta')

eq = L[0]+'='+L[1]+'='+L[2]
display(Math(eq))

L = nut.pdf2latex_mvnormal(x=r'y_t', m=r'(Cx_{t})',v=r'R')
eq = L[0]+'='+L[1]+'='+L[2]
display(Math(eq))

```


$$\mathcal{N}(x_t; (Ax_{t-1}), Q)=\left|{ 2\pi Q } \right|^{-1/2} \exp\left(-\frac12 ({x_t} - {(Ax_{t-1})} )^\top {Q}^{-1} ({x_t} - {(Ax_{t-1})} ) \right)=\exp\left(  -\frac{1}{2}\trace {Q}^{-1} {x_t}{x_t}^\top + \trace {Q}^{-1} {x_t}{(Ax_{t-1})}^\top -\frac{1}{2}\trace {Q}^{-1} {(Ax_{t-1})}{(Ax_{t-1})}^\top  -\frac{1}{2}\log \left|2{\pi}Q\right|\right)$$



$$\mathcal{N}(y_t; (Cx_{t}), R)=\left|{ 2\pi R } \right|^{-1/2} \exp\left(-\frac12 ({y_t} - {(Cx_{t})} )^\top {R}^{-1} ({y_t} - {(Cx_{t})} ) \right)=\exp\left(  -\frac{1}{2}\trace {R}^{-1} {y_t}{y_t}^\top + \trace {R}^{-1} {y_t}{(Cx_{t})}^\top -\frac{1}{2}\trace {R}^{-1} {(Cx_{t})}{(Cx_{t})}^\top  -\frac{1}{2}\log \left|2{\pi}R\right|\right)$$



```python
%connect_info
```

    {
      "key": "5fe2a052-2599-465d-96fa-793a58b02ea2",
      "transport": "tcp",
      "signature_scheme": "hmac-sha256",
      "shell_port": 65197,
      "stdin_port": 65199,
      "ip": "127.0.0.1",
      "hb_port": 65201,
      "control_port": 65200,
      "iopub_port": 65198
    }
    
    Paste the above JSON into a file, and connect with:
        $> ipython <app> --existing <file>
    or, if you are local, you can connect with just:
        $> ipython <app> --existing kernel-457030b1-3d25-4222-9094-06fb4bbdbefe.json 
    or even just:
        $> ipython <app> --existing 
    if this is the most recent IPython session you have started.


3d45602fb55ac0abffcb32ff2fb2b2683fe82a4e
ipynb
Jupyter Notebook
ModelingExamples.ipynb
bkoyuncu/notes
0e660f46b7d17fdfddc2cad1bb60dcf847f5d1e4
MIT
ModelingExamples.ipynb
bkoyuncu/notes
0e660f46b7d17fdfddc2cad1bb60dcf847f5d1e4
MIT
ModelingExamples.ipynb
bkoyuncu/notes
0e660f46b7d17fdfddc2cad1bb60dcf847f5d1e4
MIT
2019-09-11T11:46:36.000Z
2019-09-11T11:46:36.000Z
Qwen/Qwen-72B
1. YES
2. YES