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

本文主要是介绍scipy beta分布与numpy clip的数值问题，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

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

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

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

invalid value encountered in divide

此处给出各种形状（$\alpha ,\beta$ 组合）下，变量 x 在各种取值时，$\mathrm{B}(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 01 - eps1, 1 - eps2, 1, 1 + eps2, 1 + eps1, # near 12, # >>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-ϵ]$ 之间，其中 $ϵ$ = 1e-7。

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

import numpy as npzero = 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/1assert (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

这篇关于scipy beta分布与numpy clip的数值问题的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！

---