I need to show that the following equality holds for any $i,j,m,n \in \Bbb N$ and $0 \le p \le 1$ ($p$ is a probability). Could you please help me? I think I should use hyper-geometric function but I could not do it.
$$ \sum _{i=0}^m\sum _{j=0}^n \binom{m}{i}\binom{n}{j}\frac{(i+j)!(m+n-i-j)!}{(m+n+1)!}p^i(1-p)^{m-i}(j)=\frac{n+nmp}{(m+2)}$$
Finally, I solved it and wanted to share with you! Let $$LHS=\sum _{i=0}^m\binom{m}{i}p^i(1-p)^{m-i}\sum _{j=0}^n \binom{n}{j}\frac{(i+j)!(m+n-i-j)!}{(m+n+1)!} j$$ and let $$B=\sum _{j=0}^n \binom{n}{j}\frac{(i+j)!(m+n-i-j)!}{(m+n+1)!} j=\sum _{j=0}^n\frac{n(n-1)!}{(n-j)!(j-1)!}\frac{(i+j)!(m+n-i-j)!}{(m+n+1)!} .$$ We can simplify $B$ by using the Vandermonde's identity for $k=m+n-i, r=m-i, s=i+1$ and we have $$\sum _{j=0}^k \begin{pmatrix}m+n-i-j \\m-i\end{pmatrix}\begin{pmatrix}i+j\\i+1\end{pmatrix}=\begin{pmatrix}m+n+1\\m+2\end{pmatrix}$$ where the index $j$ runs to $n$ since for $j>n$ the first binomial coefficient in the summation vanishes. Now it is just a matter of writing the binomial coefficients in terms of factorials and recombining (P.S. I copied this sentence from Lozenges' solution to my previous question): $$\sum _{j=0}^n \frac{(m+n-i-j)!}{(m-i)!(n-j)!}\frac{(i+j)!}{(i+1)!(j-1)!}=\frac{(m+n+1)!}{(m+2)!(n-1)!}$$ $$\sum _{j=0}^n \frac{(n-1)!}{(n-j)!(j-1)!}\frac{(i+j)!(m+n-i-j)!}{(m+n+1)!}=\frac{(i+1)!(m-i)!}{(m+2)!}$$ $$\sum _{j=0}^n \begin{pmatrix}n-1 \\j-1 \end{pmatrix} \frac{(i+j)!(m+n-i-j)!}{(m+n+1)!}=\frac{(i+1!)(m-i)!}{(m+2)!}=\frac{1}{(m+2)\begin{pmatrix}m+1 \\i+1 \end{pmatrix}} .$$ Hence, $$B=\frac{n}{(m+2)\begin{pmatrix}m+1 \\i+1 \end{pmatrix}} .$$ Now, $$LHS=\sum _{i=0}^m\binom{m}{i}p^i(1-p)^{m-i}B=\sum _{i=0}^m\binom{m}{i}p^i(1-p)^{m-i}\frac{n}{(m+2)\begin{pmatrix}m+1 \\i+1 \end{pmatrix}}$$ $$LHS=\frac{n}{(m+2)}[\sum _{i=0}^m\binom{m}{i}p^i(1-p)^{m-i}(i+1)]=\frac{n}{(m+2)}[1+\sum _{i=0}^m\binom{m}{i}p^i(1-p)^{m-i}(i)] $$ $$LHS=\frac{n}{(m+2)}[1+\sum _{i=0}^m\frac{m!}{(m-i)!(i)!}p^i(1-p)^{m-i}(i)]=\frac{n}{(m+2)}[1+\sum _{i=0}^m\frac{m!}{(m-i)!(i-1)!}p^i(1-p)^{m-i}]$$ $$LHS=\frac{n}{(m+2)}[1+m\sum _{i=0}^m\binom{m-1}{i-1}p.p^{i-1}(1-p)^{m-i}]=\frac{n}{(m+2)}[1+mp\sum _{i=0}^m\binom{m-1}{i-1}p^{i-1}(1-p)^{m-i}] .$$ Since function in the summation is the pdf of binomial distribution with parameters $n-1$ and $p$, the summation equals $1$. Hence, $$LHS=\frac{n(1+mp)}{(m+2)}$$ which is the desired result.