Find the regression line of $X$ on $Y$ 
Suppose that $$f (x, y) = xe^{−x(y+1)}$$ where  $0 ≤ x < ∞$, $0 ≤ y < ∞$

Now how to find the regression line.
My approach ,is to find marginal density functions $f_{X} (x)$ and $f_{Y} (y)$ then using that i will findout the standard deviation $\sigma_{x} \sigma_{y}$ and $\sigma_{xy}$
I got
$$f_Y (y)=\int_0^{\infty} xe^{-x(y+1)}dx$$ $$=-\frac x {y+1}e^{-x(y+1)}|_0^{\infty} +\frac 1 {y+1} \int_0^{\infty}e^{-x(y+1)}dx$$ $$=\frac 1{(y+1)^{2}}.$$
$f_X (x)=\int_0^{\infty} xe^{-x(y+1)}dy=e^{-x}$
Now the mean around origini is coming out to be $1$ but i am unable to find $\sigma^2_y$
By defination $\sigma^2_y=\int_{0}^{\infty}\frac{(y-1)^2}{(y+1)^2} dy$ which is divergent.
Please help me out ,i think i am making some calculation mistake.
 A: First, let's do a quick recap. As StubbornAtom commented, the regression line of $X$ on $Y$ is formally defined as the conditional expectation
$$\text{E}[X | Y = y] = \int x f_{X|Y}(x|y) dx$$
For $(X, Y)$ being bivariate normally distributed, one can further show that the regression line of $X$ on $Y$ is given by
$$\text{E}[X | Y = y] = \mu_x + \rho_{xy}(\sigma_x/\sigma_y)(y - \mu_y),$$
where $\rho_{xy} = \sigma_{xy}/(\sigma_x\sigma_y)$ is the correlation between $X$ and $Y$.
In this problem, we will directly calculate $\text{E}[X | Y = y]$ to obtain the regression line of $X$ on $Y$. As you have already (correctly) computed, the marginal density of $Y$ is
$$f_Y(y) = \int_0^{\infty} xe^{-x(y+1)}dx = \frac{1}{(y+1)^2},$$
such that the conditional density of $X$ given $Y$ is
$$f_{X|Y}(x|y) = \frac{f_{X,Y}(x,y)}{f_Y(y)} = x(y+1)^2 e^{−x(y+1)}.$$
The next step is to calculate
$$\text{E}[X | Y = y] = \int x f_{X|Y}(x|y) dx = \int_0^\infty x^2(y+1)^2 e^{−x(y+1)} dx.$$
This integral can be obtained by using Integration by parts twice. Here we borrow some useful facts from Exponential distribution to make the computation easier. Recall that if $W \sim \text{Exp}(\lambda)$ with $f_W(w) = \lambda e^{-\lambda w}$ for $w \geq 0$. Then we have $\text{E}[W] = 1/\lambda$ and $\text{Var}[W] = 1/\lambda^2$, which implies
$$\text{E}[W^2] = \int_0^\infty w^2 \lambda e^{-\lambda w} dw = (\text{E}[W])^2 + \text{Var}[W] = \left(\frac{1}{\lambda}\right)^2 + \frac{1}{\lambda^2} = \frac{2}{\lambda^2}.$$
Back to this problem, we then have
\begin{align*}
\text{E}[X | Y = y] &= \int_0^\infty x^2(y+1)^2 e^{−x(y+1)} dx\\
&= (y+1) \int_0^\infty x^2(y+1) e^{−x(y+1)} dx\\
&= (y+1) \cdot \frac{2}{(y+1)^2} \quad (\text{let } \lambda = y+1)\\
&= \frac{2}{y+1},
\end{align*}
which is the regression line of $X$ on $Y$.
Remark: When $(X, Y)$ is not bivariate normal (like in this problem), the regression line of $X$ on $Y$ depends on the joint distribution of $(X, Y)$, and may not necessarily be a straight line.
