Variance of $\frac{1}{X+1}$ when $X$ has Poisson distribution I'm trying to find the variance of $\frac{1}{X+1}$  where $X \sim \mbox{Poisson}(\lambda)$.  I found the previous answers here on how to get the expected value of $\frac{1}{X+1}$, and they were helpful and easy to follow. (original: Expectation of $\frac{1}{x+1}$ of Poisson distribution Duplicate: Expected value $\frac{1}{x+1}$ of Poisson distribution )
However, I was unable to generalize this to finding the variance of the same random variable. I tried using the identity $Var(X) = E[X^2] - E[X]^2$ but was unable to derive a neat closed form of $E[\frac{1}{(1+X)^2}]$. While at times it felt like I was getting close, I think I may have hit a dead end. Any help is greatly appreciated.
 A: Let the random variable (r.v.)  $X$ on a probability space $(\Omega,\mathcal F,\Bbb P)$ correspond to a Poisson r.v. with parameter $\lambda$, so $\Bbb P(X=n) =\exp(-\lambda)\lambda^n /n!$. Then:
$$
\begin{aligned}
\Bbb E\left[\frac 1{(X+1)^2}\right]
&=
e^{-\lambda}
\sum_{k\ge 0}\frac1{(k+1)^2}\cdot\frac{\lambda^k}{k!}
\\
&\qquad\text{ we use now $n=k+1$}
\\
&=
e^{-\lambda}
\cdot\frac1\lambda
\sum_{n\ge 1}\frac1{n}\cdot\frac{\lambda^n}{n!}
\\
&=
e^{-\lambda}
\cdot\frac1\lambda
\int_0^\lambda
\frac{e^s-1}s\; ds\ .
\end{aligned}
$$
and this is a complicated (hypergeometric) function (${}_2F_2(1,1;2,2;\lambda)$ ).

Computer check / answer:
sage: var('a,k');
sage: sum( 1/(k+1)^2/factorial(k)*a^k, k, 0, oo )
hypergeometric((1, 1), (2, 2), a)

A: Thank you @dan_fulea. 
It took me a while to parse your answer, so I wanted to expand on it to ensure I'm not mistaken and to clarify for any others who have this question.
As you said:
$$ E\left[\frac 1{(X+1)^2}\right] = \frac{e^{-\lambda}}{\lambda} \int_0^\lambda
\frac{e^s-1}{s} ds$$
After some trial and error I realized what you meant was 
$$ {}_2F_2(1,1;2,2;\lambda) = \frac{1}{\lambda} \int_0^\lambda
\frac{e^s-1}{s} ds $$
So this means that 
$$ Var\left(\frac{1}{1+X}\right) = {}_2F_2(1,1;2,2;\lambda)*e^{-\lambda} - E\left[\frac{1}{1+X}\right]^2 $$
Which combined with the answer to earlier user's question (linked in original question) gives us:
$$ Var\left(\frac{1}{1+X}\right) = {}_2F_2(1,1;2,2;\lambda)*e^{-\lambda} - \frac{(1-e^{-\lambda})^2}{\lambda^2} $$ 
Just to be safe, I confirmed this answer by simulating around 20,000 Poisson draws for a few different lambdas and in each case it closely matched the measured variance. 
Thanks again! 
For those interested here's link to associated Wolfram page with y being used as the $\lambda$ in this case https://www.wolframalpha.com/input/?i=HypergeometricPFQ%5B%7B1,1%7D,+%7B2,2%7D,+y%5D*e%5E-y+++-+(1-e%5E-y)%5E2%2F(y%5E2)
