Calculating a MLE for the combinatorial probability distribution: $\sum_k p^k(1-p)^{N-j-k}\binom{N-j}{k} \cdot(1-p)^{i-k}p^{j-i+k}\binom{j}{i-k}$ I have a relatively complicated discrete probability distribution:
$$\begin{aligned}
P(i;j) &= \sum_k p^k(1-p)^{N-j-k}\binom{N-j}{k}\cdot(1-p)^{i-k}p^{j-i+k}\binom{j}{i-k} \\
&= \frac{p^j}{p^i} \frac{(1-p)^i}{(1-p)^j} \sum_{k=0}^i p^{2k}(1-p)^{N-2k} \binom{N-j}{k} \binom{j}{i-k},
\end{aligned}$$
and I need to compute the maximum for the likelihood $\mathcal{L}(\vec{x}; j)$ for some $\vec{x}=\left<x_1,\ldots,x_n \right>, 0\leq x_t\leq N$ a sequence of naturals $i$ less than $N.$ To be clear: I'm fixing $p$ and I want to find the MLE for $j$. 
As usual: 
$$\begin{aligned}
\mathcal{L}(\vec{x}; j) = \prod_{x_t}P(x_t;j)=P(0;j)^{a_0}P(1;j)^{a_1} \cdots P(N;j)^{a_N},
\end{aligned}$$
for $a_0 + a_1 + \cdots + a_N = n$.
Unfortunately, I cannot see a simple way to simplify $P(i;j)$, so computing this maximum likelihood is even worse. The usual trick of taking log doesn't seem to work because of the sum on the front. So:


*

*Any advice on simplifying my expression for $P(i;j)$? (note: it's almost a Vandermonde).

*Any tricks for computing the Maximum likelihood?


I'm ideally looking for analytic solutions. 
Two small notes:


*

*$N$ and $p$ are fixed in the above, $N$ is a fixed value between 4 and 32 and $p$ is a probability (in my case it's $0.05$), in case that's helpful.

*$P(i;j)$ is coming from a convolution, which is probably obvious.

 A: Update
From simulating data from the multinomial distribution it appears that the nearest integer to $\bar{x}$ is the maximum likelihood estimate of $j$ a very high percentage of the time (as in maybe 99% of the time but not 100% of the time).  (I won't have time for a while to give all of the details of that inference but maybe in a few days I'll have time to document that appropriately.)
End of update
If $j$ is the parameter of interest, then I don't see that there is much simplification for the log of the likelihood.  In addition, because the legitimate values of $j$ range from $0$ to $N$, I think you'll just need to calculate the log likelihood for each of those values of $j$ and pick the (hopefully unique) value that maximizes the log of the likelihood.  (Setting the  derivative of the log likelihood with respect to $j$ to zero and solving for $j$ won't get you an integer value and there are not that many values of $j$ to consider in the first place.)
Here is some Mathematica code to do so:
First define $\log(P(i;j))$:
logP[i_, j_, p_, NN_] := (j - i) Log[p] + (i - j) Log[1 - p] +
  Log[Sum[p^(2 k) (1 - p)^(NN - 2 k) Binomial[NN - j, k] Binomial[j, i - k], {k, 0, i}]]

Now define the log of the likelihood for all legitimate values of $j$ given $p$ and the counts $a$:
logL[p_, a_] := Table[{j, N[a.Table[logP[i, j, p, Length[a] - 1], {i, 0, Length[a] - 1}]]},
  {j, 0, Length[a] - 1}]

As an example consider $p=1/20$ and $a=(1,2,4,5)$:
logLikelihood = logL[1/20, {1, 2, 4, 5}]
(* {{0, -68.8659}, {1, -43.224}, {2, -26.9353}, {3, -27.6437}} *)
maxLogLikelihood = Select[logLikelihood, #[[2]] == Max[logLikelihood[[All, 2]]] &]
(* {{2, -26.9353}} *)

ListPlot[{logLikelihood, maxLogLikelihood}, PlotRangeClipping -> False,
  PlotStyle -> {{Blue, PointSize[0.02]}, {Red, PointSize[0.02]}}] 


Addition:
Sometimes simplification is in the eye of the beholder.  The probability $P(i;j)$ can be written as
$$\binom{j}{i} p^{j-i} (1-p)^{i-j+N} \, _2F_1\left(-i,j-N;-i+j+1;\frac{p^2}{(p-1)^2}\right)$$
where $_2F_1$ is the hypergeometric function.
A: It looks like the distribution of i is approximately normal for almost all N, j I tried.
See example

And after some tests, it looks like the maximum is always essentially at
$$i_{max} = p N + (1-2 p) j$$
See  the comparison between this equation and the experimentally determined max for $1<N<70$, $1<j<N$ and $p \sim Uniform(0,1)$

Essentially the true max is never further than one from the formula.
I am sure there is some clever math that proves it because the formula is so simple... Maybe seeing the answer, somebody here can figure out why the maximum is there.
EDIT: After writing the post, I figured that I wasn't quite answering the right question, as you cared about MLE not just maximum of P(i|j,N). But certainly you can try to adopt a Normal approximation, using the equation for the max as a center of the Gaussian. One only has to determine what is the width of the Gaussian.
