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

原创已于 2024-12-01 21:08:55 修改 · 1k 阅读

20 ·

CC 4.0 BY-SA版权

文章标签：

#scipy #beta #beta分布 #贝塔分布 #numpy

于 2024-04-12 22:25:40 首次发布

机器学习同时被 2 个专栏收录

121 篇文章

订阅专栏

数学

38 篇文章

订阅专栏

文章讨论了Beta分布的不同形态，包括在参数α和β取不同值时的表现，特别关注了0、1附近的数值问题。作者通过实例展示了如何使用scipy.stats.beta.pdf函数处理这些边界情况，并指出在编程中应限制x值在[ε,1-ε]之间以避免无限大值。此外，还通过numpy.clip函数验证了数值稳定性。

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

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

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

invalid value encountered in divide

此处给出各种形状（ $\alpha,\beta$ 组合）下，变量 x 在各种取值时， $\Beta(x;\alpha,\beta)$ 的值（尤其是变量 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 的值限定在 $[\epsilon, 1 - \epsilon]$ 之间，其中 $\epsilon$ = 1e-7。

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

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 却不能。

（2024.12.2）本文所讲的数值问题在 [5] 早已发现（cvpr’24 deadline 早于本文首发时间），而在方法中只是一笔带过。

References

(CVPR 2023) BiCro: Noisy Correspondence Rectification for Multi-modality Data via Bi-directional Cross-modal Similarity Consistency - paper, code
EM算法估计beta混合模型参数
贝塔分布
Beta Distribution | MIT Mathlets
(cvpr’24) Learning to Rematch Mismatched Pairs for Robust Cross-Modal Retrieval - paper, code