Simulate a Continuous Joint pdf in R Using Known Distributions I am trying to simulate a continuous joint pdf as follows, so that I can verify my solutions to a problem set. The joint pdf is $$f(x,y) = 24x(1-y) \text{ for } 0 \leq x \leq y \leq 1$$ and $0$ otherwise. I have found the marginal pdf of $y$ to be $12y^2(1-y)$, a nice beta distribution. My idea was to simulate the $x$ values using the conditional CDF of $F_{X|Y}(X|Y)$ which I have found to be $\frac{x^2}{y^2}$, but I am having trouble using a known distribution for these x-values. Any ideas? 
I have simulated the y values as follows: 
y = rbeta(10000,3,2)
 A: There are many ways you can simulate this bivariate random vector.  Probably the most efficient way is to derive the marginal distribution of one variable, and the conditional distribution of the other, and then simulate the variables individually using these distributions.
An alternative method, which is less efficient, but does not require you to derive the marginal and conditional distributions, is to use rejection sampling.  The simples method in this case is to use a uniform generating distribution over the unit square.  We have a bivariate continuous random vector $(X,Y)$ with a bounded density $f$ over the support $\mathcal{S} \subset [0,1]^2$.  Thus, we can generate $X_*,Y_* \sim \text{IID U}[0,1]$ and then accept the generated value with acceptance probability:
$$A(x_*,y_*) \equiv \frac{f(x_*,y_*)}{\sup_{(x,y) \in \mathcal{S}} f(x,y)} = \frac{24 x_* (1-y_*) \cdot \mathbb{I}(x_* \leqslant y_*)}{24} = x_* (1-y_*) \cdot \mathbb{I}(x_* \leqslant y_*).$$
It is fairly simple to program this method into R.  In the code below we create a function SIMULATE that takes an input n and produces a matrix with this many outputs of the bivariate random vector in question.  (The matrix has two columns for the two variables; each row is one simulated value of the random vector.)
#Create function to simulate vectors from specified distribution
SIMULATE <- function(n) {

  #Set output matrix
  OUT <- matrix(NA, nrow = n, ncol = 2);
  colnames(OUT) <- c('X','Y');

  #Undertake rejection sampling
  for (i in 1:n) {

    ACCEPT <- FALSE;
    while (!ACCEPT) {

      #Simulate proposed values
      X <- runif(1);
      Y <- runif(1);  

      #Determine acceptance
      AA <- X*(1-Y)*(X <= Y);
      ACCEPT <- (runif(1) <= AA); } 

    OUT[i,] <- c(X,Y); }

OUT; }

We can use this function to simulate any number of bivariate outputs from this distribution.  Below we simulate $n=10^3$ outputs of the random vector.
#Set the seed
set.seed(75375211);

#Generate simulations and show the first few values
SIMULATIONS <- SIMULATE(1000);
head(SIMULATIONS);

              X         Y
[1,] 0.4124875 0.4681140
[2,] 0.1345465 0.1565690
[3,] 0.4703997 0.4810464
[4,] 0.6532923 0.8114625
[5,] 0.5971606 0.6286653
[6,] 0.6476007 0.8088133

