拒绝抽样没有给出正确的分布

Question

我正在编写代码以根据柯西分布执行拒绝抽样。

x = np.linspace(-4, 4, 100)

dist = scipy.stats.cauchy
global_max = dist.pdf(0)

这很简单。 scipy.stats中柯西分布的默认参数是(0, 1)。在下面的代码中，我生成了一百万个随机点，并且只接受 y 坐标位于曲线下方的点，即小于相应 x 坐标处的 Cauchy PDF 值的点。

num_samples=1000000

sample_y = np.random.uniform(0, global_max, size=num_samples)
sample_x = np.random.uniform(-4, 4, size=num_samples)

X = sample_x[sample_y <= dist.pdf(sample_x)]
params = scipy.stats.cauchy.fit(X)

最后我将实际分布与估计分布进行比较：

print('Estimated parameters =', params)

plt.hist(X, bins=50, density=True, alpha=0.3)
plt.plot( x, scipy.stats.cauchy( params[0], params[1] ).pdf(x), color='red' )
plt.plot( x, dist.pdf(x), color='green' );

输出:

Actual parameters    = (0, 1)
Estimated parameters = (-0.0030743926369336217, 0.7362620502669363)

我不明白为什么会这样。方差明显不同。我在这里错过了什么？

Answer 1

柯西分布以“肥尾”着称，但您将样本限制在 [-4, +4] 范围内。太窄了。

尝试更大的范围。例如

num_samples=100000000
max_x = 2000

sample_y = np.random.uniform(0, global_max, size=num_samples)
sample_x = np.random.uniform(-max_x, max_x, size=num_samples)

拟合参数为(-0.008691511218505928, 0.9918572787654598)

X = sample_x[sample_y <= dist.pdf(sample_x)]

params = sps.cauchy.fit(X)
_X = X[np.logical_and(X < 4, X>-4)]

plt.hist(_X, bins=50, density=True, alpha=0.3)
plt.plot( x, sps.cauchy( params[0], params[1] ).pdf(x), color='red' )
plt.plot( x, dist.pdf(x), color='green' );

如果您不喜欢它的外观，请与参考 Cauchy-distributed PRNG 进行比较。

ref_x = np.random.standard_cauchy(size=int(num_samples/max_x))
plt.hist(ref_x[np.logical_and(ref_x < 4, ref_x>-4)], bins=50, density=True, alpha=0.3)
plt.plot( x, dist.pdf(x), color='green' );

这些照片看起来很像，不是吗？

拒绝抽样没有给出正确的分布

Rejection sampling not giving the correct distribution

python

numpy

simulation

probability-distribution