scipy beta分布与numpy clip的数值问题

时间：2024-04-13 11:09:22

[1] 用到混合 Beta 分布，估计参数的方法见 [2]。由 [3] 可见 Beta 分布在其参数 α , β \alpha,\beta α,β 在不同取值范围时存在几种形态：

  • α , β < 0 \alpha,\beta < 0 α,β<0：不合法；
  • α = β = 1 \alpha=\beta=1 α=β=1：常数， B ( x ; 1 , 1 ) ≡ 1 \Beta(x;1,1)\equiv 1 B(x;1,1)≡1；
  • α , β > 1 \alpha, \beta > 1 α,β>1：钟形（bell shape），即单峰（unimodal）；
  • 0 < α < 1 ≤ β 0<\alpha<1\leq\beta 0<α<1≤β：L 形；
  • 0 < β < 1 ≤ α 0<\beta<1\leq\alpha 0<β<1≤α：J 形；
  • 0 < α , β < 1 0<\alpha,\beta<1 0<α,β<1：U 形。

其中后三种在 0、1 处会取到正无穷，可能在编程时引起问题，如：

invalid value encountered in divide

此处给出各种形状（ α , β \alpha,\beta α,β 组合）下，变量 x 在各种取值时， B ( x ; α , β ) \Beta(x;\alpha,\beta) B(x;α,β) 的值（尤其是变量 x 在 0、1 附近时）作为参考：

  • 调包：scipy.stats.beta.pdf 
import scipy.stats as stats
import numpy as np

# 临界 epsilon
eps1 = 1e-7
eps2 = 1e-8

# 变量
x = np.array([
    -1, # <<0
    - eps1, - eps2, 0, eps2, eps1, # near 0
    1 - eps1, 1 - eps2, 1, 1 + eps2, 1 + eps1, # near 1
    2, # >>1
], dtype=np.float32)
print(x)

print("\tinvalid: alpha, beta < 0")
print("alpha < 0:", stats.beta.pdf(x, -0.5, 1))
print("beta < 0:", stats.beta.pdf(x, 1, -0.5))

print("\tU-shape: 0 < alpha, beta < 1")
print(stats.beta.pdf(x, 0.5, 0.5))

print("\tL-shape: 0 < alpha < 1 <= beta")
print(stats.beta.pdf(x, 0.5, 1))

print("\tJ-shape: 0 < beta < 1 <= alpha")
print(stats.beta.pdf(x, 1, 0.5))

print("\tconstant: alpha = beta = 1")
print(stats.beta.pdf(x, 1, 1))

print("\tbell-shape (unimodal): 1 < alpha, beta")
print(stats.beta.pdf(x, 2, 2))

输出：

[-1, -1e-7, -1e-8, 0, 1e-8, 1e-7, 0.99999988, 1.0000000e+00, 1, 1.0000000e+00, 1.0000001, 2]

        invalid: alpha, beta < 0
alpha < 0: [nan nan nan nan nan nan nan nan nan nan nan nan]
beta < 0: [nan nan nan nan nan nan nan nan nan nan nan nan]

        U-shape: 0 < alpha, beta < 1
[0, 0, 0, inf, 5.1460, 4.0876, 4.0164, inf, inf, inf, 0, 0]

        L-shape: 0 < alpha < 1 <= beta
[0, 0, 0, inf, 6.2062, 4.9298, 0.1997, 0, 0, 0, 0, 0]

        J-shape: 0 < beta < 1 <= alpha
[0, 0, 0, 0, 0.1558, 0.1962, 4.8439, inf, inf, inf, 0, 0]

        constant: alpha = beta = 1
[0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0]

        bell-shape (unimodal): 1 < alpha, beta
[0, 0, 0, 0, 5.99999990e-08, 5.99999947e-07, 7.15255652e-07, 0, 0, 0, 0, 0]

考察其中出现 inf 的位置，可以考虑在调用 scipy.stats.beta.pdf 时将 x 的值限定在 [ ϵ , 1 − ϵ ] [\epsilon, 1 - \epsilon] [ϵ,1−ϵ] 之间，其中 ϵ \epsilon ϵ = 1e-7。

除了上面的测试，此 ϵ \epsilon ϵ 还能如此验证：用 numpy.clip 重复实验，将 0/1 截断到 [ ϵ , 1 − ϵ ] [\epsilon, 1 - \epsilon] [ϵ,1−ϵ] 之间，看从哪个精度开始数值开始不稳定。代码：

import numpy as np

zero = np.zeros([500], dtype=np.float32)
one = np.ones([500], dtype=np.float32)
# 有 0 有 1 的数据
x = np.concatenate([zero, one], axis=0)

# 测试 numpy.clip 对各 epsilon 的稳定性
for eps in (1e-7, 1e-8):
    print(eps)
    for _ in range(100):
        y = np.clip(x.copy(), eps, 1 - eps) # deep copy, then clip
        # 若成功截断，则不应再有 0/1
        assert (0 != y).all() and (1 != y).all()

实验表明，1e-7 能让 numpy.clip 稳定截断，而 1e-8 却不能。

References

  1. (CVPR 2023) BiCro: Noisy Correspondence Rectification for Multi-modality Data via Bi-directional Cross-modal Similarity Consistency - paper, code
  2. EM算法估计beta混合模型参数
  3. 贝塔分布
  4. Beta Distribution | MIT Mathlets

相关文章