Simulation for the random vector (X,Y) with density f(x,y,a) I need to generate data from a random vector with joint density $(X,Y)$ with density: 
$$
\frac{x^\alpha(x+y+2)}{x+y+1}e^{-x-y}~~~~~,\alpha,x,y> 0
$$
Do you have any hints on how to start?
 A: The accept-reject simulation algorithm requires a density function $g(x,y)$ such that $f(x,y)/g(x,y)$ is bounded, i.e., that there exists $M>0$ such that$$\dfrac{f(x,y)}{g(x,y)}\le M$$everywhere (with the convention that $\frac{0}{0}=0$). If we look at your density
\begin{align*}f(x,y)&=\frac{x^\alpha(x+y+2)}{x+y+1}e^{-x-y}\\&=\frac{x^\alpha(x+y+1+1)}{x+y+1}e^{-x-y}\\&=x^\alpha e^{-x-y}+\frac{x^\alpha}{x+y+1}e^{-x-y}\end{align*}
Hence your target density can be written as a mixture of two positive functions, thus of two densities. The first term in the sum is the product of a Gamma $\text{Ga}(\alpha+1,1)$ density on $x$ and of an Exponential $\text{E}(1)$ density on $y$, missing the normalising constant $\Gamma(\alpha+1)$. The second term is a more unusual joint density that requires a specific simulation algorithm. For instance,
\begin{align*}\frac{x^{\alpha-1}\,x}{x+y+1}e^{-x-y}&=\frac{x^\alpha}{x+y+1}e^{-x-y}
\\&\le x^{\alpha-1} e^{-x-y}\end{align*}
which exhibits a product of a Gamma $\text{Ga}(\alpha,1)$ density on $x$ and of an Exponential $\text{E}(1)$ density on $y$, missing the normalising constant $\Gamma(\alpha)$. In conclusion, we have that
$$f(x,y)\le \Gamma(\alpha+1)\dfrac{x^\alpha}{\Gamma(\alpha+1)} e^{-x-y}+\Gamma(\alpha)\frac{x^{\alpha-1}}{\Gamma(\alpha)}e^{-x-y}=Mg(x,y)$$
Generating from $g$ is straightforward:

  
*
  
*choose between $i=1$ and $i=0$ with probabilities $\alpha/1+\alpha$ and $1/1+\alpha$;
  
*generate $x\sim \text{Ga}(\alpha+i,1)$ and $y\sim \text{E}(1)$
  

and therefore generating from $f$ is as straightforward:

  
*
  
*generate $(x,y)\sim g(x,y)$
  
*generate $u\sim\text{U}(0,1)$ a uniform variate
  
*accept $(x,y)$ if$$u\le \dfrac{f(x,y)}{Mg(x,y)}$$else goto 1.
  

A: Welcome to cross validated! This is an interesting first question!
There are three classical ways to achieve this.
The first one is to try to factor the densities into two factors, one only including $x$ and another only including $y$. Then the $X$ and the $Y$ random variables would be independent and things would be easy. I do not see a way to factor this so I guess this option is out in this case.
The second option is inverse transform sampling. That requires you to factor your density as:
$$
f_{XY}(x,y) = f_Y(y|x)\cdot f_X(x)
$$
where $f_{XY}(x,y)$ is your joint density. Then you sample from the marginal density of $X$ (namely $f_X(x)$) and use that value to sample from the conditional density of $Y$ given $X$ based on the value you sampled from the density of $X$. This requires you to calculate the inverse of the corresponding CDFs, which would be a bit tricky in this case. I found the marginal density for $X$, but it involved the exponential integral. Calculating the inverse of the CDF seems to be nontrivial (to me at least).
The third option is the usage of rejection sampling. That is probably your best bet in this case. 
Here is a solution from stackoverflow with some R code. The user that answers uses a multivariate normal as an example and the instrumental distribution is a bivariate uniform. You will just need to replace your pdf in the mydens function. In this case instrumental distribution has bounded support, so your sample will be within that support. This will not give you an exact solution, but will probably be good enough if the support for the uniform is large enough. 
You may want to pick another instrumental distribution in your case. It is easiest to find something that is as similar to your distribution as possible.
Hope this helps!
EDIT: I have no idea what your parameter $\alpha$ is, so for the purpose of this example I will assume that it is 2. I need a specific value to make some of these plots. First let's examine how your density looks:
mydens <- function(x,y)
{
  alpha <- 2
  val <- ((x^alpha)*(x+y+2)/(x+y+1))*exp(-x-y)
  return(val)
}

x <- seq(0,7,len=50)
y <- seq(0,7,len=50)

f<-outer(x,y,mydens)
persp(x,y,f,theta=140)

This gives me the following plot:

Visually I think this looks as if I can assume that the majority of the probability mass is contained in the region $[0,10]\times[0,10]$ for $\alpha=2$. Now let's use reject.sample.2d function from the other answer:
samp <- reject.sample.2d(1000,mydens,10,c(0,10),c(0,10))
plot(samp[,1],samp[,2])

This gives the following plot:

This looks visually already very much like the other one. Let's look at an image of the contour of the density and a kernel estimate of the density with 10000 samples:


These look very similar. Here the instrumental distribution is a uniform, you will definitely want to find something better than that!
