python - 逆 (CDF) 变换采样的错误分布

Question

我正在尝试使用逆 CDF 方法来模拟几何分布，但是我得到的结果略有错误，我不确定原因。

更具体地说，形状因子 p = 0.8 的几何分布应具有以下特征:

mean: 1.25 
variance: 0.31

但是，运行下面的代码，我得到:

mean: 0.6224363901913519
var: 0.391813011265263
[Finished in 0.3s]

正如您所看到的，与预期相比，我得到的平均值截然不同。

np.log(uniform[i])/np.log(1-p) is the result of solving the equation: F(X) = R for X in terms of R, F(X) = CDF of geometric distribution = 1 - (1 - p)^k.

R 是区间 (0,1) 上的均匀分布。

因此解决它会产生以下结果:

X = ln(1-R)/ln(1-p)

但是，由于 1-R 和 R 均均匀分布在 (0,1) 上，因此我们可以进行以下简化:

X = ln(R)/ln(1-p)

上述方程是正确的，应该会产生几何分布样本。

import numpy as np

n = 10000
p = 0.8
geo_dist = np.zeros(n,dtype = np.float64)
uniform = np.random.uniform(0, 1, n)
for i in range(n):
    geo_dist[i] = np.log(uniform[i])/np.log(1-p)
print("mean: " +str(geo_dist.mean()))
print("var: " +str(geo_dist.var())) 

我尝试通过使用 np.float64 来提高计算精度，拼命尝试修复本来应该是微不足道的脚本，但无济于事。

我还尝试使用 scipy Uniform.rvs() 而不是 np.uniform 生成均匀分布，但问题仍然存在。

如果 p = 0.5:

expected mean: 2
expected variance : 2

但是，我编写的代码具有以下结果:

mean: 1.4440009653569306
var: 2.0421079966161093
[Finished in 0.3s]

有人知道为什么这不起作用吗？ 谢谢。

Answer 1

您实际上是在连续采样 exponential distribution lambda 等于 -1/ln(1-p)

好的，这是正确采样的代码，上限应用于指数输出

import numpy as np

N = 100000
p = 0.8

q = np.random.random(N)
g = np.ceil(np.log(1.0 - q)/np.log(1.0-p))

print(np.mean(g))
print(np.var(g))

打印内容

1.25055
0.3146946975

请注意:

1. 您最好使用 NumPy 向量化功能，而无需显式循环

2. 从 U(0,1) 采样的 R 的替换 (1-R) -> R 对于 NumPy RNG 不起作用 - 它会返回值在半封闭范围 [0...1) 中，这意味着您可能会时不时地得到 log(0) 和 FP 异常。