Proof of an identity involving a double sum of binomial factors.

Question

In the course of analyzing a certain Markov model described in Solving another non-trivial recurrence relation we discovered a following identity. Let $\Lambda \ge 0$ and $\lambda^M \ge 0$ and $0 \le q \le 1$ and let $n \in {\mathbb N}$ and $n\ge 1$. Then we have: \begin{eqnarray} &&q \sum\limits_{p=0}^n \sum\limits_{\eta=0}^p \binom{\lambda^M+n}{\eta} \binom{-\lambda^M}{p-\eta} \cdot \left(\lambda^M+1+\Lambda+p-\eta\right)^{(n-1-p)} \cdot \left( \frac{q \Lambda}{1-q} -(p-\eta)+1\right)^{(p-1)} \cdot \left( 1_{n=p} (\Lambda+n+\lambda^M-\eta) + \Lambda 1_{n\neq p}\right) =\\ && \left( \frac{\Lambda}{1-q} + \lambda^M+1 \right)^{(n-1)} \quad (1) \end{eqnarray} Below is the Mathematica code that verifies that identity for $n \le 5$. We have:

In[346]:= (*Simplifying distribution for r\[Equal]1*)
n = RandomInteger[{1, 7}]; lM =.; L =.; q =.; r =.;
eX1 = Table[
   Sum[Binomial[lM + n, eta] Binomial[-lM, p - eta]  Pochhammer[
      lM + 1 + L + p - eta, 
      n - 1 - p] Pochhammer[(q L)/(1 - q) - (p - eta) + 1, p - 1] If[
      n == p, -eta + lM + L + n, L] q, {p, 0, n}, {eta, 0, p}], {n, 1,
     5}];
Collect[Take[eX1, 5], L, FullSimplify]
eX11 = Factor[eX1];
eX2 = Table[Pochhammer[L/(1 - q) + lM + 1, n - 1] , {n, 1, 5}];
Simplify[eX11 - eX2]

Out[348]= {1, 
 1 + lM + L/(1 - q), (1 + lM) (2 + lM) + (L (3 + 2 lM))/(1 - q) + 
  L^2/(-1 + q)^2, (1 + lM) (2 + lM) (3 + lM) - L^3/(-1 + q)^3 + (
  3 L^2 (2 + lM))/(-1 + q)^2 + (
  L (-11 - 3 lM (4 + lM)))/(-1 + 
   q), (1 + lM) (2 + lM) (3 + lM) (4 + lM) + L^4/(-1 + q)^4 - (
  2 L^3 (5 + 2 lM))/(-1 + q)^3 + (
  L^2 (35 + 6 lM (5 + lM)))/(-1 + q)^2 - (
  2 L (5 + 2 lM) (5 + lM (5 + lM)))/(-1 + q)}

Out[351]= {0, 0, 0, 0, 0}

Now, in the attempt to proving the identity in question for any value of $n$ I have done the following. Firstly we notice that both sides of the equation are polynomials of order $(n-1)$ in the variable $\Lambda$ and as such we can compare coefficients at the leading power, i.e. at $\Lambda^{n-1}$. We have: \begin{eqnarray} lhs &=& q \sum\limits_{p=0}^n \sum\limits_{\eta=0}^p \binom{\lambda^M+n}{\eta} \binom{-\lambda^M}{p-\eta} \cdot \Lambda^{n-1} \left( \frac{q}{1-q} \right)^{p-1} \\ &=& (1-q) \cdot \left(1+\frac{q}{1-q} \right)^n \cdot \Lambda^{n-1}= \frac{\Lambda^{n-1}}{(1-q)^{n-1}} = rhs \end{eqnarray} Here in the second line we firstly did the sum over $\eta$ using the Chu-Vandermonde identity and then the sum over $p$ using the binomial expansion.

Being encouraged by this success now we move on to comparing the coefficient at $\Lambda^0$ in both sides. Because of the fourth term under the sum in the left hand side in $(1)$ there is only very few terms that are non zero. Firstly it is the term ($\eta=p$ and $p=n$ and $p\ge 1$) then there is the term ($\eta=0$ and $p=n$ and $p\ge 1$) and finally there is the term ($\eta=0$ and $p=0$). We have: \begin{eqnarray} \left.\mbox{term in sum}\right|_{\eta=p \; \& \; p=n \; \& \; p \ge 1} &=& q \cdot \binom{\lambda^M+n}{n} \cdot 1 \cdot (n-1)! \\ \left.\mbox{term in sum}\right|_{\eta=0 \; \& \; p=n \; \& \; p \ge 1} &=& q \cdot \binom{-\lambda^M}{n} \cdot (-1)^{n-1} \cdot (n-1)! \\ \left.\mbox{term in sum}\right|_{\eta=0 \; \& \; p=0 } &=& (1-q) \cdot (\lambda^M+1)^{(n-1)} \end{eqnarray} Therefore we have: \begin{eqnarray} &&lhs = q \cdot (n-1)!\left[ \binom{\lambda^M+n}{n}- \binom{\lambda^M+n-1}{n} \right] + (1-q) \cdot (1+\lambda^M)^{(n-1)} =\\ (\lambda^M+1)^{(n-1)} = rhs \end{eqnarray} where in the second line we used the Pascal triangle identity.

Now my question would be the following. How do we complete the proof? So far we have only compared two extreme coefficients of the polynomial in question . Can we prove that all coefficients match? Or maybe there is some other easier way of proving the identity in question.

Now let us take a look at the coefficient at $\Lambda^1$ in both sides of the equation. The only non-zero contributions to that coeffiecient are listed below: \begin{eqnarray} \left.\mbox{term in sum}\right|_{p=n \; \& \; \eta=0,\cdots,p} &=& \frac{q^2 (\lambda^{M}+1)^{(n-1)} (\psi ^{(0)}(\lambda^{M}+n)-\psi ^{(0)}(\lambda^{M}+1))}{1-q} \\ \left.\mbox{term in sum}\right|_{1\le p \le n-1 \; \& \; \eta=0} &=& q (\lambda^{M}+1)^{(n-1)} \left(-H_{n-1}+\psi ^{(0)}(\lambda^{M}+n)-\psi ^{(0)}(\lambda^{M}+1)\right) \\ \left.\mbox{term in sum}\right|_{1\le p \le n-1 \; \& \; \eta=p} &=& q (\lambda^{M}+1)^{(n-1)} \left(H_{n-1}+\psi ^{(0)}(\lambda^{M}+n)-\psi ^{(0)}(\lambda^{M}+1)\right) \\ \left.\mbox{term in sum}\right|_{p=0 \; \& \; \eta=0} &=& (1-q) (\lambda^{M}+1)^{(n-1)} (\psi ^{(0)}(\lambda^{M}+n)-\psi ^{(0)}(\lambda^{M}+1)) \end{eqnarray} Therefore summing up those terms we get : \begin{equation} lhs = \frac{1}{1-q} \cdot (\lambda^{M}+1)^{(n-1)} (\psi ^{(0)}(\lambda^{M}+n)-\psi ^{(0)}(\lambda^{M}+1)) = rhs \end{equation}