实现密度函数

Question

我正在浏览我的书，它说“为这个密度函数编写一个采样算法”

y=x^2+(2/3)*x+1/3; 0 < ???? < 1

或者我可以使用蒙特卡洛？ 任何帮助将不胜感激！

Answer 1

我假设你有一个函数 y(x)，它取一个 [0,1] 之间的值并返回 y 的值。你只需要提供一个随机的 x 值并返回对应的 y 值。

def getSample():
  #get uniform random number
  x = numpy.random.random()

  #sample my custom function
  return y(x)

Answer 2

我假设您的意思是要生成具有由密度 y(x) 指定的分布的随机 x 值。

通常需要通过对密度进行积分来导出累积分布函数，并使用inverse transform sampling 生成x 值。在您的情况下，CDF 是一个三阶多项式，它不会产生简单的立方根解，因此您必须使用数值求解器来找到逆。是时候考虑替代方案了。

另一种选择是使用acceptance/rejection method。检查导数后，很明显你的密度是凸的，所以很容易通过从f(0) 到f(1) 画一条直线来创建边界函数b(x)。这产生b(x) = 1/3 + 5x/3。此边界函数的面积为 7/6，而您的 f(x) 的面积为 1，因为它是有效密度。因此，在b(x) 下统一生成的 6/7 点也将落在f(x) 下，并且拒绝方案中只有七分之一的尝试会失败。这是f(x) 和b(x) 的图：

由于b(x) 是线性的，因此很容易生成x 值，将其用作在按6/7 缩放后的分布，使其成为有效的分布函数。算法，用伪代码表示，然后变成：

function generate():
  while TRUE:
    x <- (sqrt(1 + 35 * U(0,1)) - 1) / 5     # inverse CDF transform of b(x)
    if U(0, b(x)) <= f(x):
      return x
  end while
end function

其中U(a,b) 表示生成一个均匀分布在a 和b 之间的值，f(x) 是你的密度，b(x) 是上面描述的边界函数。

我实现了上述算法以生成 100,000 个候选值，其中 14,199 个（约 1/7）如预期的那样被拒绝。最终结果显示在以下直方图中，您可以将其与上图中的f(x) 进行比较。