有很多方法可以模拟这个二元随机向量。可能最有效的方法是导出一个变量的边际分布和另一个变量的条件分布，然后使用这些分布单独模拟变量。

另一种效率较低但不需要导出边际分布和条件分布的替代方法是使用拒绝抽样。在这种情况下，最简单的方法是在单位正方形上使用均匀的生成分布。我们有一个二元连续随机向量 ，在支持上具有有界密度。因此，我们可以生成然后以接受概率接受生成的值：$(X,Y)$$f$$\mathcal{S} \subset [0,1]^2$$X_*,Y_* \sim \text{IID U}[0,1]$

$$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_*).$$

将此方法编程到R. 在下面的代码中，我们创建了一个函数，该函数SIMULATE接受一个输入n并生成一个矩阵，其中包含所讨论的二元随机向量的这么多输出。（矩阵有两列代表两个变量；每一行是随机向量的一个模拟值。）

#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; }

我们可以使用这个函数来模拟来自这个分布的任意数量的二元输出。下面我们模拟输出。$n=10^3$

#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