[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 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 − ϵ ] [\epsilon, 1 - \epsilon] [ϵ,1−ϵ] 之间，其中 ϵ \epsilon ϵ = 1e-7。
除了上面的测试，此 ϵ \epsilon ϵ 还能如此验证：用 numpy.clip 重复实验，将 0/1 截断到 [ ϵ , 1 − ϵ ] [\epsilon, 1 - \epsilon] [ϵ,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
- (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
)
