Maximal Entropy

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma
import seaborn as sns

sns.set_theme()

Laplace’s Method

\(L(\mu) \approx L(\mu_*) + \frac{\partial}{\partial \mu} L(\mu) \Bigr|_{\mu = \mu_*} (\mu - \mu_*)+ \frac{1}{2} \frac{\partial^2}{\partial \mu^2} L(\mu) \Bigr|_{\mu = \mu_*} (\mu - \mu_*)^2 + ...\) We only take the terms up to the second derivative as a good enough approximation. Because then $\frac{\partial}{\partial \mu} L(\mu) \Bigr|{\mu = \mu} = 0$, because at $\mu_$ the posterior has a maxima, we get

Then setting $L’’(\mu_*) = -\frac{1}{\sigma^2}$ and taking the exponent we get a gaussian:

Expanding then around the logarithm we get

Again the first term vanishes because of the maxima at $ \mathbf{\alpha}^*$. The second derivate gives the so called hessian matrix.

If we set $B_{ij} = -H_{ij}$ we get the approximated posterior by a multivariate gaussian

The covariance matrix is then given by the inverse of the Hessian. If now integrate all parameters but one, the marginal distribution again is a gaussian with variance $\sigma^2i = B{ij}^{-1}$:

Lagrange

Assume we have $N$ genes in total, where each gene can have $m=0,1,2,…,$ mRNA copies in a given cell. On average, a gene has $\bar{m}$ mRNAs per cell, i.e. there are $M=\bar{m}N$ mRNAs per cell. We now would like to know, what is the probability $ \mathbb{P}(m | I)$ that one particular gene has $m$ mRNAs in the cell.

1. Standard way

For each gene $g$ we specify the number $m_g$ of mRNAs that it has in the cell, i.e. the mutually exclusive and exhaustive possibilities are vectors $ \mathbf{m} = (m_1, m_2, … , m_N)$. We can then incorporate our prior information that among these potentially possible states, our information $I$ specifies that only those states have nonzero probability for which we have:

This we have a uniform prior which restricts our space to the space which has the previous given property

and for a single gene we have

For large N these calculations become difficult for large N, which is why we look at an approximate. Instead of looking at th number of mRNAs $m_g$ that each gene $g$ has, we look at the vector with the number of genes that have precisely $m$ mRNA copies $ \mathbf{n} = (n_0, n_1, n_2, …, n_M)$, which counts how many genes there are with 0 mRNA ($n_0$), with 1 mRNA ($n_1$) etc. For a given vector $ \mathbf{n}$ there are many vectors $ \mathbf{m}$. Let then $W(\mathbf{n})$ denote the number of vectors $ \mathbf{m}$ that all havce the same count vector $ \mathbf{n}$. The numbers $W(\mathbf{n})$ are given by the multinomial coefficients.

\(W(\mathbf{n}) = \frac{N!}{n_0 ! n_1 ! ... n_M!}\) With our constraint $M( \mathbf{n}) = \sum_{m=1}^{\infty} mn_m = M = N \bar{m}$.Then our probability distribution becomes

Using now the stirling approximation for factorials $n! \approx n^n e^{-n}$:

With $f_m = \frac{n_m}{N}$, we get the entropy function $H[\mathbf{f}] = - \sum_m f_m \log f_m$. Thus we get that $W(\mathbf{f}) = e^{N H [\mathbf{f} ]} $ and $ \mathbb{P}(\mathbf{f} | I) \propto e^{NH [\mathbf{f} ]} \delta (m(\mathbf{f} | - \bar{m}))$ with $m(\mathbf{f} ) = \sum_{m} mf_m$
For $N \rightarrow \infty$ we gat the probability distribution is dominated by the entropy of the vector $ \mathbf{f}_* $ which has the highest entropy.

To now find the vector $ \mathbf{f}* $ which maximizes the entropy, we use the lagrange multipliers. Denote $ \mathbf{\nu} $ as vector in the dimension $M + 1$ suscht that if we move away from $ \mathbf{f}* $ in the direction of $ \mathbf{\nu} $ then the average $m (\mathbf{f} )$ is unchanged. Our optimum is characterized by $ \mathbf{\nu} \nabla H[ \mathbf{f} ] = 0 \ \forall \mathbf{\nu}, \ \mathbf{\nu} \nabla m (\mathbf{f} ) = 0$

Because by the multivariate taylor expansion we get :

Thus then if $ \mathbf{\nu} \nabla m( \mathbf{f} ) = 0 $ it follows that $m( \mathbf{f} + \epsilon \mathbf{\nu}) = m ( \mathbf{f} ) $, meaning that it is constant in the direction of $ \mathbf{\nu} $.

Then the same for our entropy function

\(H[ \mathbf{f} + \epsilon \mathbf{\nu} ] = H[ \mathbf{f} ] + \epsilon \mathbf{\nu} \nabla H[\mathbf{f} ]\) Because we want our optimum, $ \mathbf{\nu} \nabla H[ \mathbf{f} ] = 0 \ \forall \mathbf{\nu} $ with $ \mathbf{\nu} \nabla m( \mathbf{f} ) = 0 $. If we definde $X[\mathbf{f}] = H[\mathbf{f}] + \lambda m(\mathbf{f})$ then $\nabla X[\mathbf{f] = 0 \Leftrightarrow \mathbf{\nu} \nabla H[\mathbf{f}] = 0 $ when $ \mathbf{\nu} \nabla m(\mathbf{f}) = 0 $.
When $\nabla X[\mathbf{f}] = 0 \Leftrightarrow \nabla H[\mathbf{f}] = - \lambda \nabla m(\mathbf{f})$. Thus at $ \mathbf{f}_* $ we have that for any direction $ \mathbf{\nu} $ Here we get

If we add another constraint that $C(\mathbf{f}) = \sum_m f_m = 1 $ we get $X(\mathbf{f} ) = H[\mathbf{f}] + \lambda_1 m(\mathbf{f}) + \lambda_2 C(\mathbf{f})$ which then gives us $f_m = e^{\lambda_1 m + \lambda_2 - 1}$

This can be generalized to

Looking back, because we get that $f_m = e^{\lambda_1 m + \lambda_2 - 1}$ we set our $\lambda_2$ by demanding that $\sum_m f_m = 1$

for $\lambda_1$ we demand that $\sum_m m f_m = \bar{m} = \sum_m m e^{\lambda_1 m} (1 - e^{\lambda_1})$. For this we introduce the partition function

Which then gives $\langle m \rangle = \frac{d \log Z(\lambda)}{d \lambda} = \frac{e^{\lambda}}{1 - e^{\lambda}} = \bar{m} \Rightarrow \lambda = \log \left( \frac{\bar{m}}{\bar{m} + 1} \right) $. Thus in the end we get our probability of finding a gene having m mRNAs is $ \mathbb{P}(m | I) = f_m = \frac{1}{\bar{m} + 1} \left( \frac{\bar{m}}{\bar{m} + 1} \right)^m $

Information theory

Axioms:

1. There exist a funtion $H[\mathbb{P} ]$ that assigns a real number to each probability distribution $ \mathbb{P} $ which quantifie the ignorance / uncertainty / not-knowing associated with that probability distribution.
2. It is a continous function of it’s arguments
3. Fr the uniform distribution over n variables, the function $h(n) = H \left( \frac{1}{n}, \frac{1}{n}, …, \frac{1}{n} \right)$ should increase with n.
4. It should be consistent in that, if it can be calculated in multiple ways, it always gives the same result.

Assuming a uniform distribution, we split our n possibilities into $G$ groups with the first $g_1$ hypotheses in the first group, the $g_2$ in the second etc. Then we get

\(h(n) = H \left[ \frac{g_1}{n}, \frac{g_2}{n},..., \frac{g_G}{n},\right] + \sum_{i=1}^G \frac{g_i}{n} h(g_i)\) The first term here signifies the uncertainty of which group we are in, and then the sum signifies the uncertainty of which of the $g_i$ are we in the group $g_i$ . If we then set all $ \frac{g_i}{n} = \frac{g}{n} = \frac{1}{G}$ we get $h(gG) = h(G) + h(g)$. The solution to this equation is then $h(n) = k \log n$, where $k \in \mathbb{R}$. Substituting this back we get

Setting then $p_i = \frac{g_i}{n}$ we get our again our entropy function $H[p_1, p_2, …, p_n] = - \sum_{i=1}^n p_i \log p_i$.

Guessing game

For this we divide all possible values of $k$ into two subsets of equal probability $K$ and $K’$, i.e.

\(\sum_{k \in K} \mathbb{P}(k | I) = \sum_{k' \in K'} \mathbb{P}(k' | I) = \frac{1}{2}\) We then ask the oracle if $k$ lies in $K$ or $K’$. By this we are maximizing the information we receive at each step. On average then the number of questions is given by $H = - \sum_k \mathbb{P}(k | I ) \log_2 \mathbb{P}(k | I )$

Then given the list of optimal questions, the sender of the message can already foresee what questions an optimal receiver would ask, and send the list of yes/no answers as a binary bit string. This way the value $k$ can be transmitter as a bit string that has on average $H$ bits.

Conditional entropy

For a two dimensional joint distribution we get the entropy

Imagine now that we have sampled from the joint distribution $ \mathbb{P}(x, y)$ and are told the value of $y$, and then calculate the entropy of $x$ given $y$. What is the average of this entropy?

This is the conditional entropy of $x$ given $y$ and it tells us how much uncertainty there is on average left about $x$ if we’re told about the value of $y$. We can further rewrite our conditional entropy to

\(\text{-} \sum_{x, y} \mathbb{P}(x,y) \log \mathbb{P}(x|y) = - \sum_{x, y} \mathbb{P}(x,y) \log \frac{\mathbb{P}(x,y)}{\mathbb{P}(y)} = - \sum_{x, y} \mathbb{P}(x,y) \log \mathbb{P}(x,y) - \sum_{x, y} \mathbb{P}(x,y) \log \mathbb{P}(y) = - \sum_{x, y} \mathbb{P}(x,y) \log \mathbb{P}(x,y) - \sum_{y} \mathbb{P}(y) \log \mathbb{P}(y) = H(x, y) - H(y)\) Then if we rewrite this for both the conditional $x$ and $y$ we get the so called mutual information $I(x, y) = H(x) + H(y) - H(x, y)$ or $H(y) - H(y|x) = H(x) - H(y|x) = I(x, y)$. The mutual information quantifies how much information $x$ and $y$ contain about each other. The mutual information is zero if and only if knowing $x$ does not reduce the number of questions about $y$ and vice versa, $x$ and $y$ are then independent and we have

Thus only if $\mathbb{P}(x, y) = \mathbb{P}(x)\mathbb{P}(y) $ is there no mutual information. Mutual information is thus often a measure of dependence between variables.

Maximal entropy

Often we can encode constraint certain observables in our prior information i.e. functions $f_k(x)$ are constrained to have certain values.

Given this information, our prior probability should be the one that satisfies these constraints and maximizes the entropy

The maximum entropy distribution $ \mathbb{P}(x)$ satisfies

If we introduce the partition function $Z(\lambda_1, \lambda_2, …, \lambda_n) = \sum_x \exp \left( - \sum_{k=1}^n \lambda_k f_k(x) \right)$, then our solution can be written as

Where the lagrange multipliers are set by solving

The second derivatives of the logarithm of the partition function give covariances between the constraints

But because the first derivatives gave the constraints, these second derivatives also give the derivatives of the mean of one constraint in respect to the lagrange multiplier of the other

Thus we get the relationship

Sufficient statistics

In progress

Jane Widget problem

Imagine a company makes red, yellow and green widgets. Only 200 widgets of a single color can be made each day. How do we decide which color widget to make to maximize the probability of being able to fill all orders coming in that day.

1. Stock

Imagine we only know the current stock of the widgets. We have 100 red, 150 yellow and 50 green. To find out which color we should make, we look at the universe of possibilities. Each color has a count and in total we have a triplet of integers which signify the amount ordered of each color.

We assume nothing for the different colors, i.e. each triplet has the same probability $ \mathbb{P}(n_r, n_y, n_g) = const$. We then want to maximize the probability to fill all orders for that day given the amount for each color we have in stock $(s_r, s_y, s_g)$

Thus we want all triplets for which it holds that $s_i \geq n_i$. Thus we get the total probability of

For our given stock we get then, when making only on color:

Thus our best guess is to make only green widgets.

1. Average order amopunt

We’re given additional information that on average per day we get orders for 50 red, 100 yellow and 10 green. We can setup this probability be maximizing our probability of orders under constraints of the averages. So we choose the distribution which maximizes our entropy distribution, thus the distribution with the least amount of ‘surprise’. Given $ \mathbb{P}(n_r, n_y, n_g)$ our lagrange problem becomes

This is solved by differentiating in respect to $ \mathbb{P}$ and setting it zero.

By normalising with the function

Note that through dividing of our $ \mathbb{P}$ by our normalising constant, $e^{-\lambda - 1}$ cancels itself out, thus we get:

Thus we get our final distribution, which also shows independence of the variables.

To find our lambdas, we need that $- \frac{d \log Z}{d_i} = \langle n_i \rangle $. For our red’s we get that

Which then gives that

Our final distribution is then

Because they are independant we can look at each color individually and see the probability of being able to fill all orders

This is the same for all the other colors and we get

The probabilities are then

Thus making yellow is our best option.

1. Average order for color

We now are given additional information which gives the total number of independent orders per day. On average we have 2/3 reds, 10 yellow and 1/2 green per day. For this we now need to distinguish orders of different sizes. Our universe of possibilities now specifies how many orders $n_k$ for $k$ widgets of each color. The possible set of orders for today can be specified by a vactor $(n_1, n_2, n_3, …)$, where $n_1$ is the amount of orders for 1 widget etc.
The total number of widgets ordered today is then $t = \sum_i i n_i$ and the amount of orderes is given by $n = \sum_i n_i$. We again find the maximal entropy distribution

Which gives us

The normalisation constant is given by

For our lagrange multipliers we get

We must then fit these our lagrange multipliers to the given averages $\bar{n_r} = 2/3$, $\bar{n_y} = 10$, $\bar{n_g} = 1/2$ and $\bar{t_r} = 50$, $\bar{t_y} = 100$ and $\bar{t_g} = 10$. For red we must have then

Numerically solving this for all the colors gives us

We now need for each color the distribution $ \mathbb{P}(t)$ for the total number ordered. Because the the number of orders for each amount are mutually exclusive we can write

Because $t$ is the sum of many independent terms we can approximate the distribution $ \mathbb{P}(t) $ by the central limit theorem with a gaussian distribution $ \mathbb{P}(t) \approx C \exp \left( -\frac{(t - \langle t \rangle)^2}{2 \sigma^2} \right)$.

For the mean we have $\langle t \rangle = \sum_i i \langle n_i \rangle $ and the variance $var(t) = \sum_i i^2 var(n_i) $. For $\langle n_i \rangle$ we have

For $var(n_i)$ we get

We then get $\sigma_r = 85.8$, $\sigma = 46.0$ and $\sigma_g = 19.3$. Then we get the distributions

The probability to then fill all orders is given by

Which gives us then

### Priors

Imagine we use a uniform prior in the variable $x$, which gives us the prior distribution $ \mathbb{P}(x | I) = const$. If we now perform a coordinate transformation $x \mapsto y = x^3$. Our uniform under this coordinate transformation becomes

Thus through a coordinate transformation our initial want to be unbias / ignorant causes our transformed distribution to not reflect this. Thus we basically want to find priors such that under coordinate transformation, their initial structure is preserved.

1. Shifts

Let our transformation be a shift, i.e. $x \mapsto y = x + a$ then for our consistency we want that

The general solution to this is the function $ \mathbb{P}(x) \propto e^{\lambda x}$ with $\lambda \in \mathbb{R}$. Because

For $\lambda = 0$ this becomes the uniform prior

1. Rescaling

Let our transformation be a rescaling, i.e. $x \mapsto y = \alpha x $ then for our consistency we want

The general solution of this is given by $ \mathbb{P}(x) \propto x^{\lambda}$. Because

For $\lambda = 0$ we get our uniform prior. For rescaling and shifts, the only distribution which has both is the uniform distribution.

1. Changing dimensions

Let our transformation be a change of dimensions, i.e. $x \mapsto y = x^{\alpha} $ then for our consistency we want

The only solution to this equation is $ \mathbb{P}(x) dx = \frac{\log^c x}{x} dx$. Because

The only distribution that is invariant under changing dimensions and rescaling is $ \frac{1}{x}$. This distribution is simply uniform in the logarithm of the parameter $ \mathbb{P}(x) dx \propto d \log x$.

Benfords law

Assume we obtain a very large table of numbers. We then would like to ask, what fraction of the numbers would we expect to start with digit $d$, i.e. $x = d(…)$? We assume that the numbers are distributed uniform in the log-space $ \mathbb{P}(x) dx \propto d \log x$. To calculate this we look at the range [0.1, 0.2], [1, 2], [10, 20], etc.

In generel for any digit $d$ we get $ \mathbb{P}(d) \propto \log \left( \frac{d + 1}{d} \right) $. In general our distribution becomes

This benfords law. It gives the frequencies of digits.

Dirichlet Priors

Assume our systems contains $n$ different objects in unknown proportions. Our space of hypotheses is thus the $n$-dimensional space of proportions

If we now demand that our prior is invariant to scaling, i.e $\rho_i \mapsto \rho_i’ = \alpha_i \rho_i$. The prior that is invariant under these transformation then has the general form

If we now demand symmetrie between all objects then $\lambda_i \overset{!}{=} \lambda_j$ for $i \neq j$. Then our general, normalized distribution becomes

lam = [0.1, 0.5, 1, 5, 10]
f = np.arange(0, 1, 0.01)

for lam_i in lam:
    prob = (gamma(2 * lam_i) / gamma(lam_i)**2) * f**(lam_i - 1) * (1 - f)**(lam_i - 1)
    plt.plot(f, prob, label=f'{lam_i}')
plt.legend()
plt.show()

C:\Users\mathi\AppData\Local\Temp\ipykernel_24320\3866864515.py:5: RuntimeWarning: divide by zero encountered in power
  prob = (gamma(2 * lam_i) / gamma(lam_i)**2) * f**(lam_i - 1) * (1 - f)**(lam_i - 1)

from matplotlib import cm

x = np.arange(0, 1, 0.001)

X, Y = np.meshgrid(x, x)
lam = 2

for i, x_i in enumerate(X):
    for j, y_j in enumerate(Y):
        if X[i, j] + Y[i , j] > 1 + 1e-3:
            X[i, j] = 0
            Y[i , j] = 0


Z = (gamma(3 * lam) / gamma(lam)**3) * X**(lam - 1) * Y**(lam - 1) * (1 - X - Y)**(lam - 1)
fig = plt.figure(figsize=plt.figaspect(0.5))

ax = fig.add_subplot(1, 2, 1, projection='3d')

ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False)

ax = fig.add_subplot(1, 2, 2)
countour = ax.contourf(X, Y, Z)
fig.colorbar(countour, ax=ax)
plt.show()