如果一个函数
f
(
x
)
f(x)
f(x) 在整个坐标轴上的积分有界：
∫
−
∞
∞
f
(
x
)
d
x
=
M
<
∞
int^{infty}_{-infty} f(x) dx = M < infty
∫−∞∞f(x)dx=M<∞
那么就可以用它来定义一个分布
f
^
hat{f}
f^：
f
^
(
x
)
≡
f
(
x
)
M
⇒
∫
−
∞
∞
f
^
(
x
)
d
x
=
∫
−
∞
∞
f
(
x
)
d
x
M
=
1
begin{aligned} & hat{f}(x) equiv frac{f(x)}{M} \ Rightarrow & int^{infty}_{-infty} hat{f}(x) dx =frac{int^{infty}_{-infty} f(x) dx}{M} = 1 end{aligned}
⇒f^(x)≡Mf(x)∫−∞∞f^(x)dx=M∫−∞∞f(x)dx=1
如何从这个以
f
^
hat{f}
f^ 为概率密度函数的分布中采样呢？
很多时候由于
f
(
x
)
f(x)
f(x)的复杂性，我们甚至不能计算出上述的归一化因子
M
M
M，那么如何采样，可以使得样本的分布正比于
f
(
x
)
f(x)
f(x) 呢？
举个例子，
f
(
x
)
=
e
−
∣
x
−
4
∣
+
1
(
1
+
x
2
)
2
+
log
(
1
+
∣
x
∣
)
(
1
+
∣
x
+
2
∣
)
2
f(x) = e^{-|x-4|} + frac{1}{(1+x^2)^2} + frac{log(1+|x|)}{(1+|x+2|)^2}
f(x)=e−∣x−4∣+(1+x2)21+(1+∣x+2∣)2log(1+∣x∣)
形状有3个尖峰：
import numpy as np
import matplotlib.pyplot as plt
def f(x):
return np.exp(-np.abs(x-4)) + 1/np.square(1 + np.square(x))+ np.log(1+np.abs(x))/(1+np.abs(x+2))**2
x = np.linspace(-100,100,1000)
y = f(x)
plt.plot(x,y)
plt.xlim([-20,20])
大家感兴趣可以试试把积分即归一化常数
M
M
M 算出来
为了方便之后方案的对比验证，我这里直接用 scipy 的数值积分函数，计算出该分布的期望为：0.1575
def expectedValue(f):
from scipy import integrate
M, error = integrate.quad(f, -np.inf, np.inf)
assert(error < 1e-6)
expected, error = integrate.quad(lambda x: x*f(x), -np.inf, np.inf)
assert(error < 1e-6)
expected /= M
return expected
拒绝采样
f
(
x
)
<
M
g
(
x
)
,
∀
x
f(x)< M g(x), forall x
f(x)
