MLE for Poisson-binomial distribution I am looking for the maximum likelihood estimator (MLE) for the Poisson-binomial distribution.  I understand the derivation of the MLE for a Poisson distribution and a binomial distribution, but I am unable to derive the MLE equation for the Poisson-binomial distribution.
 A: Consider a Poisson-binomial distribution with probability parameter $\boldsymbol{\theta} \equiv (\theta_1, ..., \theta_m) \in [0,1]^m$.  The probability density function (PDF) for this distribution is given by:
$$p(k | \boldsymbol{\theta}) = \sum_{\mathcal{A} \in \mathfrak{P}_k(m)} \Bigg( \prod_{j \in \mathcal{A}} \theta_j \Bigg) \Bigg( \prod_{j \notin \mathcal{A}} (1-\theta_j) \Bigg),$$
where $\mathfrak{P}_k(m)$ is the class of all subsets of $k$ labels from $\{ 1, 2, ..., m \}$.  The sampling density is clearly invariant to permutations of the probability vector $\boldsymbol{\theta}$, and so this parameter is unidentifiable in the distribution.  The minimal sufficient parameter for this distribution is the vector of ordered probability values $\boldsymbol{\theta}_*= (\theta_{(1)} \geqslant ... \geqslant \theta_{(m)})$.  This means that the MLE for $\boldsymbol{\theta}$ cannot be estimated uniquely, and for any MLE vector, any permutation of that vector will also be an MLE.

Setting aside this identifiability issue, we can obtain equations for the MLE.  Letting $\boldsymbol{\theta}_{-a} \in [0,1]^{m-1}$ be the probability parameter excluding the element $\theta_a$, we have the recursive formula:
$$p(k | \boldsymbol{\theta}) = \theta_a p(k-1 | \boldsymbol{\theta}_{-a}) + (1-\theta_a) p(k | \boldsymbol{\theta}_{-a}).$$
We therefore have the useful preliminary result:
$$\frac{\partial p}{\partial \theta_a} (k | \boldsymbol{\theta}) = p(k-1 | \boldsymbol{\theta}_{-a}) - p(k | \boldsymbol{\theta}_{-a}).$$
Now, given the observed vector $\boldsymbol{k} \equiv (k_1, ..., k_n)$ taken from IID draws from the Poisson-binomial distribution, we have log-likelihood function $l_{\boldsymbol{k}} (\boldsymbol{\theta}) = \sum_{i=1}^n \ln p(k_i | \boldsymbol{\theta})$ which gives us the score:
$$\frac{\partial l_\boldsymbol{k}}{\partial \theta_a} (\boldsymbol{\theta}) = \sum_{i=1}^n \frac{\partial}{\partial \theta_a} \ln p(k_i | \boldsymbol{\theta}) = \sum_{i=1}^n \frac{p(k_i-1 | \boldsymbol{\theta}_{-a}) - p(k | \boldsymbol{\theta}_{-a})}{p(k_i | \boldsymbol{\theta})}.$$
The MLE occurs at any point $\hat{\boldsymbol{\theta}}$ satisfying:
$$\sum_{i=1}^n \frac{p(k_i-1 | \hat{\boldsymbol{\theta}}_{-a}) - p(k | \hat{\boldsymbol{\theta}}_{-a})}{p(k_i | \hat{\boldsymbol{\theta}})} = 0 \text{ } \text{ } \text{ } \text{ for all } a=1,2,...,m.$$
We have already noted that the MLE is invariant to permutations, and therefore we only require a "representative" of the class of permutations.  We can obtain this by imposing the additional order constraint $\theta_{1} \geqslant ... \geqslant \theta_{m}$.  For $n \geqslant m$ this set of $m$ score equations --plus the ordering constraint-- should yield a unique solution for $\hat{\boldsymbol{\theta}}$.  This solution represents the class of $m!$ permutations of this vector that are all MLEs.
