EM algorithm on discrete random variable Im trying to find the iteration formula for the maximum likelihood estimator of the parameter $\theta$ of the following discrete random variable, through the EM algorithm:
$$P(Y = y ) = \sum_{i=0}^{\infty} \frac{\theta^{i + y}(1+i)^y}{i!y!}e^{-2\theta-\theta i}, y = 0, 1, 2 \ldots$$
I know that we can use this as the augmented variable of the EM algorithm:
$$f(y,\theta) = \frac{\theta^{i + y}(1+i)^y}{i!y!}e^{-2\theta-\theta i}$$.
But I don't know how to proceed with the algorithm from this point, how do I find the function that I need to maximize? How do I construct the closed iterative formula?

I'm going to generalise your problem a little bit, to deal with a broader class of distributions.  Specifically, I will add an additional parameter $\alpha$ to the problem; your question follows as a special case with $\alpha=1$.

The EM algorithm is used when you want to compute an MLE in a situation where we have a latent variable in the analysis.  The basic idea of the algorithm is that it allows us to maximise the likelihood function for marginal observations by maximising the underlying likelihood function for joint observations that include the latent variable (which is very useful).  So the first thing we want to do to use the EM algorithm is to break down the marginal density to determine the underlying hierarchical model using the latent variable.  To do this, we observe that the marginal density can be written as:
$$\begin{align}
\mathbb{P}(Y=y|\theta) 
&= \sum_{z=0}^\infty \frac{\theta^{z + y}(1+z)^y}{z!y!}e^{-(\alpha+1) \theta-\theta z} \\[6pt]
&= \sum_{z=0}^\infty \frac{[(1+z)\theta]^y e^{-(\alpha+z)\theta}}{y!} \cdot \frac{\theta^z e^{-\theta}}{z!} \\[6pt]
&= \sum_{z=0}^\infty \text{Pois}(y|(\alpha+z)\theta) \cdot \text{Pois}(z|\theta), \\[6pt]
\end{align}$$
which means that it is generated from the following hierarchical model using the latent variable $Z$:
$$\begin{align}
Y|Z,\theta &\sim \text{Pois}((\alpha+Z)\theta), \\[6pt]
Z|\theta &\sim \text{Pois}(\theta). \\[6pt]
\end{align}$$
Your particular question is the special case of this model where we have $\alpha=1$.

The EM algorithm: Now that we understand the latent variable in the model, we can proceed to implement the EM algorithm.  Suppose that we have IID data points $\mathbf{y} = (y_1,...,y_n)$ and $\mathbf{z} = (z_1,...,z_n)$ generated from the above model and let $\ell_{\mathbf{y}}$ and $\ell_{\mathbf{y},\mathbf{z}}$ be the marginal and joint log-likelihood functions respectively.  The latter can be written as:
$$\begin{align}
\ell_{\mathbf{y},\mathbf{z}}(\theta) 
&= \sum_{i=1}^n \log \bigg( \frac{\theta^{z_i + y_i}(1+z_i)^{y_i}}{z_i!y_i!}e^{-(\alpha+1)\theta-\theta z_i} \bigg) \\[6pt]
&= \text{const} + \sum_{i=1}^n [(z_i + y_i) \log(\theta) - (z_i+\alpha+1) \theta ] \\[6pt]
&= \text{const} + n(\bar{z}_n + \bar{y}_n) \log(\theta) - n(\bar{z}_n+\alpha+1) \theta. \\[6pt]
\end{align}$$
For the expectation-step in the EM algorithm, we form the function:
$$\begin{align}
Q(\theta|\theta^{(t)}) 
&\equiv \mathbb{E}(\ell_{\mathbf{y},\mathbf{Z}}(\theta) | \theta^{(t)}) \\[12pt]
&= \text{const} + n(\mathbb{E}(\bar{Z}_n|\theta^{(t)}) + \bar{y}_n) \log(\theta) - n(\mathbb{E}(\bar{Z}_n|\theta^{(t)})+\alpha+1) \theta \\[8pt]
&= \text{const} + n \bigg[ (\theta^{(t)} + \bar{y}_n) \log(\theta) - (\theta^{(t)}+\alpha+1) \theta \bigg]. \\[6pt]
\end{align}$$
This function has first and second derivatives:
$$\begin{align}
\frac{dQ}{d \theta}(\theta|\theta^{(t)}) 
&= n \bigg[ \frac{\theta^{(t)} + \bar{y}_n}{\theta} - (\theta^{(t)}+\alpha+1) \bigg], \\[6pt]
\frac{d^2 Q}{d \theta^2}(\theta|\theta^{(t)}) 
&= -n \cdot \frac{\theta^{(t)} + \bar{y}_n}{\theta^2}, \\[6pt]
\end{align}$$
which establishes that it is a strictly concave function with a maximising point at its unique critical point.  For the maximisation-step in the EM algorithm, we then have:
$$\begin{align}
\theta^{(t+1)} 
= \underset{\theta}{\text{arg max}} \ Q(\theta|\theta^{(t)}) 
= \frac{\theta^{(t)} + \bar{y}_n}{\theta^{(t)} + \alpha+1}. \\[6pt]
\end{align}$$
Taking $t \rightarrow \infty$ gives a simple quadratic form for the MLE, which then yields the explicit solution:
$$\hat{\theta}_\text{MLE} = \frac{1}{2} \bigg[ \sqrt{\alpha^2 + 4\bar{y}_n} - \alpha \bigg].$$
