从自定义分布生成随机数

Question

我正在尝试从自定义分布中生成随机数，我已经找到了这个问题：Simulate from an (arbitrary) continuous probability distribution，但不幸的是它并不能帮助我，因为那里提出的方法需要一个分布函数的公式。我的分布是多个均匀分布的组合，基本上分布函数看起来像一个直方图。例如：

f(x) = { 
    0     for  x < 1
    0.5   for  1 <= x < 2
    0.25  for  2 <= x < 4
    0     for  4 <= x
}

- Simon Eismann

Answer 1

您只需要使用反向CDF方法即可：

samplef <- function (n) {
  x <- runif(n)
  ifelse(x < 0.5, 2 * x + 1, 4 * x)
  }

计算累积分布函数来验证：

F(x) = 0                 x < 1
       0.5 * x - 0.5     1 < x < 2
       0.25 * x          2 < x < 4
       1                 x > 4

因此，它的反函数是：

invF(x) = 2 * x + 1      0 < x < 0.5
          4 * x          0.5 < x < 1

- Zheyuan Li

---

谢谢，这很好用。您能解释一下我为什么需要翻转 cdf 吗？ - Simon Eismann

哦，现在我明白了， x <- runif(n) ifelse(x < 0.5, 2 * x + 1, 4 * x) 这是两行代码（我总是忘记在R中不需要分号）。 首先生成一个[0,1]的均匀分布，然后进行映射。再次感谢。 - Simon Eismann

Answer 2

你可以将各种有效的抽样方法与连续均匀分布相结合，用于离散分布。也就是说，从变量的整数部分Y=[X]模拟，该整数部分具有等于每个区间内概率的离散分布（例如通过表格法-别名方法），然后只需添加随机均匀分布[0,1$]，即 X=Y+U。

在你的例子中，Y取值为1、2、3，概率分别为0.5、0.25和0.25（这等价于以相等的概率抽取1、1、2、3），然后再加上一个随机均匀分布。

如果你的“直方图”非常大，这可能是一种非常快速的方法。

在R中，你可以通过简单的方式（虽然效率不高）实现这一点。

sample(c(1,1,2,3))+runif(1)

or

sample(c(1,1,2,3),n,replace=TRUE)+runif(n)

更普遍地说，您可以在sample中使用概率权重参数。如果您需要比这更快的速度（对于某些应用程序可能需要），特别是在具有大量直方图和非常大的样本大小的情况下，您可以使用链接中提到的方法来加速离散部分，并在较低级别语言（例如C）中编程该函数的工作核心部分。尽管如此，即使只是使用上述代码与相当“大”的直方图 - 数十到数百个箱子 - 这种方法似乎甚至可以在我的相当普通的笔记本电脑上在不到一秒钟的时间内生成一百万个随机值，因此对于许多应用程序，这将是可以接受的。

- Glen_b

---

谢谢，如果所有的“箱子”宽度相同，您的方法似乎非常直观。我不太理解您的代码行sample(c(1,1,2,3),n,replace=TRUE)+runif(n)。 - Simon Eismann

顺便提一下，你可以像这样在sample函数中使用概率： sample(c(1, 2, 3), size=3000000, replace=TRUE, prob=c(0.5, 0.25, 0.25))，这样你就不需要使用两个1的解决方法了 :-) - Simon Eismann

(1,1, ... 部分是因为1个箱子的出现频率是其他箱子的两倍；如果样本实现得好，这种方式应该比更一般的概率加权方法更快。...2,3) 部分是将2-4号箱子分割成相同宽度的部分。同样，这是为了提高速度。如果箱子的高度和宽度不都是有理数（尽管问题中没有提到），你可能需要稍微慢一些但更一般的方法，可以通过适当使用 sample 和 runif 函数来实现。 - Glen_b

谢谢您的解释，我明白了 :-) - Simon Eismann

Answer 3

如果您从经验观察定义的自定义分布中获取随机数，则可以使用fishmethods::remp()包和函数直接从经验分布中获取。

https://rdrr.io/cran/fishmethods/man/remp.html

- Roasty247