Convolution of random variables: unimodality of the likelihood function Let $X_1, X_2,...X_k$ be random independent variables, each $X_i$ drawn from a Geometric distribution $\mathcal{G}(p_i)$, and let its convolution, or sum, be $Y = \sum_{i=1}^k X_i$.
The likelihood function of each $X_i$ with respect to $p_i$ is unimodal (the $r$ parameter is known). But what can we say about the likelihood of the sum (with respect to  $p_1,...,p_k$)?
$$
\mathcal{L}(\boldsymbol{p} = p_1,...p_k; Y)
$$
The context of my problem is the following: I have an EM algorithm to estimate the parameters $\boldsymbol{p} = p_1,p_2,...,p_k$ (given a set of observations $Y_1,...,Y_n$), and I would like to know whether I can say that this is the unique global maximum or whether it might be a local optimum.
Are there some distributions, or family of distributions, for which we have theoretical results?*
I'm looking either for a direct answer or a pointer to some references.
*Appart from the trivial cases such as sums of i.i.d Poissons, sums of Geometrics with same $p$, etc which give known distributions.
 A: The short answer is no, it is not a unique global maximum, but there is a longer answer as well.
Consider the sum of two Geometric variates with parameters $p_1 = 0.1,\space p_2 = 0.9$.  If you have access to the individual variates, all is well, but if you only have access to the sum, it should be intuitively clear that you will not be able to distinguish the two cases $p_1 =0.1,\space p_2 = 0.9$ and $p_1 = 0.9,\space p_2 = 0.1$.  More generally, you will not be able to distinguish $p_1 = a,\space p_2 = b$ from $p_1 = b,\space p_2 = a$, because the distribution of the sum is the same either way.  Because of this, depending on your starting values, your algorithm could converge to either (for example) $\hat{p}_1 = 0.134, \space \hat{p}_2 = 0.802$ or $\hat{p}_1 = 0.802, \space \hat{p}_2 = 0.134$, and the associated values of the likelihood function will be the same.  This is because the distribution of the sum is independent of the arrangement of the distributions of the random variates across the indices of the random variates being summed.  This will always be the case if the components are independent.
In another context - estimation of mixture models - this is related to the "label-switching" problem.
If, on the other hand, you don't care about $p_1,\space p_2,\dots,p_k$ as they relate to the specific unobserved random variables, then things are different.  You just want a collection of $k$ estimates, and don't care which is associated with index 1, index 2, etc.  You can achieve this by generating your estimates, then ordering them from smallest to largest, for example (there are plenty of other ways of doing it), so $p_1^*$ is your estimate of the smallest $p_i$, etc.  In this case, the ordered estimates do form a unique global maximum of the likelihood function - assuming at least $k$ data points - subject to the constraint.  But note that you achieve this by, in effect, imposing a constraint on the estimation procedure that forces the ordering $p_1 \leq p_2 \leq \dots \leq p_k$ on the estimates themselves.  
A: jbowman provided an answer regarding the label-switching problem. But putting that aside it made me wonder about the proof for the unique global maximum given the constraint of $p_1 < p_2 < etc$. How can we be sure that the set of p-values for the maximum $\mathcal{L}$ is unique?
I could get an answer in the case of a single measurement $Y$. But still wonder about multiple measurements of $Y$.

Consider a particular set $(X_1, X_2, X_3 ...)$ to obtain the sum $Y=\sum X_i$ 
Then for each individual variable $X_i$ the maximum likelihood of the associated $p_i$ is $\hat{p_i} = \tfrac{1}{X_i}$ with value $\mathcal{L}(\hat{p_i} \vert X_i) = (1-\tfrac{1}{X_i})^{X_i-1}\tfrac{1}{X_i}$
And for the total set of variables $X_i$ the likelihood (if for each $X_i$ we choose the maximum $p_i$) is $\mathcal{L}(\hat{p_1}, \hat{p_2}, \hat{p_3}, etc. \vert   X_1,X_2,X_3,etc.) = \prod \mathcal{L}(\hat{p_i} \vert X_i) = (1-\tfrac{1}{X_1})^{X_1-1}\tfrac{1}{X_1} (1-\tfrac{1}{X_2})^{X_2-1}\tfrac{1}{X_2} (1-\tfrac{1}{X_3})^{X_3-1}\tfrac{1}{X_3} etc$
For a given set $X$ we can analyze two separate terms in the $\mathcal{L}$-function. And consider whether the likelihood increases if we would choose two different terms with the same sum. This relates to evaluating the maximum of the related terms parameterized by $m$ and $x$ (with $m$ the midpoint of the two variables $X_i$ and $x$ the difference of the variables from the midpoint):
$\frac{1}{m+x}\frac{1}{m-x}(1-\frac{1}{m+x})^{m+x-1}(1-\frac{1}{m-x})^{m-x-1}$
whose derivative for $x$ is equal to
$\frac{(1-\frac{1}{m-x})^{(m-x)}(1-\frac{1}{m+x})^{(m+x)}\left( log(1-\frac{1}{m+x}) log(1-\frac{1}{m-x}) \right) }{(m-x-1)(m+x-1)}$
from which we can conclude that the maximum is present $x=0$ and that the closer the variables in the set $X$, or the more uniform, the higher the likelihood. 
And thus there is a unique set of $p_i$ (the most homogeneous set $X$ of variables $X_i$ such that we can't select any pair from $X$ for which we could decrease the value of $x$ according to the previous discussed pattern) for the maximum likelihood given a single $Y$ and a given number size of the set of variables $X$.
