我正在编写代码来执行柯西分布的拒绝采样。
x = np.linspace(-4, 4, 100)
dist = scipy.stats.cauchy
global_max = dist.pdf(0)
这很简单。 scipy.stats
中柯西分布的默认参数为 (0, 1)
。在下面的代码中,我生成了一百万个随机点,并且只接受 y 坐标位于曲线下方的点,即小于相应 x 坐标处的柯西 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)
我无法理解为什么会发生这种情况。方差显着不同。我在这里缺少什么?
最佳答案
柯西分布以“肥尾”而闻名,但您将样本限制在 [-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' );
如果您不喜欢它的外观,请与引用柯西分布式 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' );
这些照片看起来很相似,不是吗?
关于python - 拒绝抽样未给出正确的分布,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/67864985/