Is there a relationship between Maximum Likelihood Estimation and the Maximum Entropy Principle? I know that both techniques can be used to estimate distribution from the data, but I didn't see anything in common between the two and I haven't found anything yet for the internet that relates the two.
 A: Yes there is a connection. Namely, fitting maximum likelihood in exponential (Gibbs) family models is equivalent to doing maximum entropy with some constraints. Formally speaking, there is a primal-dual relationship between max entropy and max likelihood - see for example Section 3.6 in Wainwright and Jordan (2008).
To see why this is the case less formally (i.e. without going into convex analysis and dual representations) let's first look at the maximum entropy principle and its solution. The principle tells us to choose the distribution $p^*$ which is "consistent with the data" and has maximum entropy, i.e.
\begin{align}
& p^* = \arg\max_{p \in \mathcal{P}} H(p) \\
 & \text{subject to} \; \mathbb{E}_p \phi_k(x) = \hat{\mu}_k \quad \forall k  = {1, \dots, K} \quad (\dagger)
\end{align}
where $H(.)$ is Shannon entropy, $\mathcal{P}$ is the set of all possible distributions and the constraints encode the aforementioned consistency with data.
The optimal solution $p^*$ is of the form (e.g. see Section 3.1 in Wainwright and Jordan (2008) or Wikipedia or derive by introducing Lagrange multipliers)
$$
p_\theta(x) \propto \exp\left\{ \sum_{k=1}^K \theta_k \phi_k(x) \right\}  \quad \quad (1),
$$
where $\theta_k$ are chosen to satisfy the constraints.
Now suppose you do maximum likelihood, assuming a $k$-parameter exponential family for our data, i.e.
$$
p_\theta(x)= \exp\left\{\sum_{k=1}^K \theta_k \phi_k(x) - A(\theta) \right\}  \quad \quad (2),
$$
where the parameters $\{\theta_k\}$ are unknown, the functions $\{\phi_k\}$ are given (these are called sufficient statistics), and $A(\theta) = \log\left(\int \exp\left\{\sum_{k=1}^K \theta_k \phi_k(x)\right\}dx\right)$ is the log-normalising constant (aka log-partition function).
Note that Equation (2) is of the same form as Equation (1). Hence, if we can show that the maximum likelihood solution satisfies the moment matching constraints $(\dagger)$ then we are done.
Start by writing the log-likelihood of $n$ datapoints, setting $\theta=\{\theta_k\}$ and $\phi(.) = \{\phi_k(.)\}$ for cleaner notation:
$$
\ell(\theta; x_{1:n}) = \sum_{i=1}^n \log p_\theta(x_i) = \sum_{i=1}^n \theta^T\phi(x_i) - nA(\theta),
$$
then differentiating and setting to zero gives us the maximum likelihood solution:
$$
\nabla_\theta A(\theta) = \frac{1}{n}\sum\phi(x_i).
$$
Finally note that the LHS $\nabla_\theta A(\theta) = \mathbb{E}_{p(\theta)}\phi(x)$, while the RHS is exactly its empirical  average and so we recover the constraints $(\dagger)$.
