How to generate random variables from a bivariate known PDF in R?

Question

I have this bivariate probability density function in a DX x DY rectangular region:

Link to my pdf

I am using R. How can I generate random (x,y) points within the rectangle following this pdf distribution?

I have read many answers regarding "inverse transform sampling", but I don't have a univariate pdf. I have seen this method, but looks tedious and extremely difficult to implement.

Thank you in advance!

Answer 1

This seems to be more of a statistical than a programming question. Anyway, rejection sampling should work in your case. You can read up the details at wikipedia. The following function performs rejection sampling (n...number of samples; pdf...probability density function; maxval...maximum (or bigger) value your pdf can yield; xlim, ylim... boundary box wherein density function yields values greater zero):

reject.sample.2d <- function(n,pdf,maxval,xlim,ylim)
{
  smpl <- data.frame(x=numeric(n),y=numeric(n))
  i <- 0
  while (i<n)
  {
    xval <- runif(1,xlim[1],xlim[2])
    yval <- runif(1,ylim[1],ylim[2])
    if (runif(1)<pdf(xval,yval)/maxval)
    {
      i <- i+1
      smpl[i,] <- c(xval,yval)
    }
  }
  return(smpl)
}

For instance, it can be used for the 2d normal distribution

mydens <- function(x,y)
{
  dnorm(x)*dnorm(y)
}

in that way (the maximum value of the here defined 2d normal distribution is at x=0,y=0 and a bit smaller than 0.16):

res <- reject.sample.2d(5000,mydens,0.16,c(-5,5),c(-5,5))

Some checking:

> sd(res[["x"]])
[1] 1.015413
> sd(res[["y"]])
[1] 0.9981738

> shapiro.test(res[["x"]])

    Shapiro-Wilk normality test

data:  res[["x"]]
W = 0.9995, p-value = 0.1603

> shapiro.test(res[["y"]])

    Shapiro-Wilk normality test

data:  res[["y"]]
W = 0.9997, p-value = 0.8304

p.s. You can also use inverse transform sampling. Calculate the marginal distribution along i.e. the x-axis. Use this distribution for inverse transform sampling to get an x-value. For this x-value use the conditional probability distribution (probability for y given the obtained x) and use inverse transform sampling again, now for the y-value.