实现密度函数

Question

我正在翻阅我的书，它指出“为此密度函数编写一个采样算法”

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

或者我可以使用 Monte Carlo？如有任何帮助，我们将不胜感激！

Answer 1

我假设您有一个函数 y(x)，它的值介于 [0,1] 和 return 之间，即 y 的值。您只需要提供一个随机值 x 和 return 相应的 y 值。

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

  #sample my custom function
  return y(x)

Answer 2

我假设你的意思是你想要生成具有密度指定分布的随机 x 值 y(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) 下，并且 7 次尝试中只有 1 次会在拒绝方案中失败。这是 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) 进行比较。

实现密度函数

Implement density function

algorithm

function

distribution

montecarlo