我提出了以下解决方案，它应该比 @cardinal、@whuber 和 @stephan-kolassa 迄今为止的其他解决方案更简单、更高效和/或计算成本更低。

它涉及以下简单步骤：

1) 绘制两个标准均匀样本：

$$u_1 \sim Unif(0,1)\\ u_2 \sim Unif(0,1).$$

应用以下剪切变换 （右下三角中的点反映到左上三角，它们将是“un-反映”在 2b): $\min\{u_1,u_2\}, \max\{u_1,u_2\}$

$$\begin{bmatrix} x\\y \end{bmatrix} = \begin{bmatrix} 1\\1 \end{bmatrix} + \begin{bmatrix} \frac{\sqrt{2}}{2} & -1\\ \frac{\sqrt{2}}{2} - 1 & 0\\ \end{bmatrix} \, \begin{bmatrix} \min\{u_1,u_2\}\\ \max\{u_1,u_2\}\\ \end{bmatrix}.$$

2b)如果和。$x$$y$$u_1 > u_2$

3）如果在单位圆内（接受率应该在72%左右），即：

$$x^2 + y^2 < 1.$$

该算法背后的直觉如图所示。

步骤 2a 和 2b 可以合并为一个步骤：

2) 应用剪切变换并交换

$$x = 1 + \frac{\sqrt{2}}{2} \min(u_1, u_2) - u_2\\ y = 1 + \frac{\sqrt{2}}{2} \min(u_1, u_2) - u_1$$

以下代码实现了上述算法（并使用@whuber 的代码对其进行了测试）。

n.sim <- 1e6
x.time <- system.time({
    # Draw two standard uniform samples
    u_1 <- runif(n.sim)
    u_2 <- runif(n.sim)
    # Apply shear transformation and swap
    tmp <- 1 + sqrt(2)/2 * pmin(u_1, u_2)
    x <- tmp - u_2
    y <- tmp - u_1
    # Reject if inside circle
    accept <- x^2 + y^2 > 1
    x <- x[accept]
    y <- y[accept]
    n <- length(x)
})
message(signif(x.time[3] * 1e6/n, 2), " seconds per million points.")
#
# Plot the result to confirm.
#
plot(c(0,1), c(0,1), type="n", bty="n", asp=1, xlab="x", ylab="y")
rect(-1, -1, 1, 1, col="White", border="#00000040")
m <- sample.int(n, min(n, 1e4))
points(x[m],y[m], pch=19, cex=1/2, col="#0000e010")

一些快速测试会产生以下结果。

算法https://stats.stackexchange.com/a/258349。最佳 3：每百万分 0.33 秒。

这个算法。最佳 3：每百万分 0.18 秒。

每秒200万点会做吗？

分布是对称的：我们只需要计算出整个圆的八分之一的分布，然后将其复制到其他八分圆周围。在极坐标中，随机位置在值的累积分布之间的面积给出延伸到的圆弧。因此它与$(r,\theta)$$\Theta$$(X,Y)$$\theta$$(0,0), (1,0), (1,\tan\theta)$$(1,0)$$(\cos\theta,\sin\theta)$

$$F_\Theta(\theta) = \Pr(\Theta \le \theta) \propto \frac{1}{2}\tan(\theta) - \frac{\theta}{2},$$

它的密度从何而来

$$f_\Theta(\theta) = \frac{d}{d\theta} F_\Theta(\theta) \propto \tan^2(\theta).$$

我们可以使用拒绝方法（效率为）从这个密度中采样。$8/\pi-2 \approx 54.6479\%$

径向坐标和之间成正比。这可以通过 CDF 的简单反演进行采样。$R$$rdr$$r=1$$r=\sec\theta$

如果我们生成独立样本，则转换回笛卡尔坐标对这个八分圆进行采样。因为样本是独立的，所以随机交换坐标会根据需要从第一象限产生一个独立的随机样本。（随机交换只需要生成一个二项式变量来确定要交换多少个实现。）$(r_i,\theta_i)$$(x_i,y_i)$

的每个这样的实现平均需要一个统一变量（对于）加上乘以两个统一变量（对于）和少量（快速）计算。这是每个点的变量（当然，有两个坐标）。完整的详细信息在下面的代码示例中。这个数字绘制了超过 50 万个点中的 10,000 个。$(X,Y)$$R$$1/(8\pi-2)$$\Theta$$4/(\pi-4) \approx 4.66$

这是R生成此模拟并对其进行计时的代码。

n.sim <- 1e6
x.time <- system.time({
  # Generate trial angles `theta`
  theta <- sqrt(runif(n.sim)) * pi/4
  # Rejection step.
  theta <- theta[runif(n.sim) * 4 * theta <= pi * tan(theta)^2]
  # Generate radial coordinates `r`.
  n <- length(theta)
  r <- sqrt(1 + runif(n) * tan(theta)^2)
  # Convert to Cartesian coordinates.
  # (The products will generate a full circle)
  x <- r * cos(theta) #* c(1,1,-1,-1)
  y <- r * sin(theta) #* c(1,-1,1,-1)
  # Swap approximately half the coordinates.
  k <- rbinom(1, n, 1/2)
  if (k > 0) {
    z <- y[1:k]
    y[1:k] <- x[1:k]
    x[1:k] <- z
  }
})
message(signif(x.time[3] * 1e6/n, 2), " seconds per million points.")
#
# Plot the result to confirm.
#
plot(c(0,1), c(0,1), type="n", bty="n", asp=1, xlab="x", ylab="y")
rect(-1, -1, 1, 1, col="White", border="#00000040")
m <- sample.int(n, min(n, 1e4))
points(x[m],y[m], pch=19, cex=1/2, col="#0000e010")

好吧，可以更高效地完成，但我当然希望你不是在寻找更快的.

这个想法是首先对值上方的垂直蓝色切片的长度成比例：$x$$x$

$$f(x) = 1-\sqrt{1-x^2}.$$

Wolfram 帮助您整合：

$$\int_0^x f(y)dy = -\frac{1}{2}x\sqrt{1-x^2}+x-\frac{1}{2}\arcsin x.$$

所以累积分布函数将是这个表达式，缩放到积分到 1（即除以）。$F$$\int_0^1 f(y)dy$

现在，要生成您的值，请选择一个随机数，均匀分布在和之间。然后找到使得。也就是说，我们需要反转 CDF（逆变换采样）。这可以做到，但并不容易。也不快。$x$$t$$0$$1$$x$$F(x)=t$

最后，给定和之间均匀分布的随机。$x$$y$$\sqrt{1-x^2}$$1$

下面是R代码。请注意，我在值的网格中预先评估 CDF，即使这样，这也需要几分钟。$x$

如果您投入一些思考，您可能可以大大加快 CDF 反转的速度。再一次，思考很痛苦。我个人会选择拒绝抽样，它更快且更不容易出错，除非我有很好的理由不这样做。

epsilon <- 1e-6
xx <- seq(0,1,by=epsilon)
x.cdf <- function(x) x-(x*sqrt(1-x^2)+asin(x))/2
xx.cdf <- x.cdf(xx)/x.cdf(1)

nn <- 1e4
rr <- matrix(nrow=nn,ncol=2)
set.seed(1)
pb <- winProgressBar(max=nn)
for ( ii in 1:nn ) {
    setWinProgressBar(pb,ii,paste(ii,"of",nn))
    x <- max(xx[xx.cdf<runif(1)])
    y <- runif(1,sqrt(1-x^2),1)
    rr[ii,] <- c(x,y)
}
close(pb)

plot(rr,pch=19,cex=.3,xlab="",ylab="")