• α, β < 0: Invalid region
• α = β = 1: Constant distribution, B(x; 1, 1) ≡ 1
• α, β > 1: Bell-shaped (unimodal)
• 0 < α < 1 ≤ β: L-shaped
• 0 < β < 1 ≤ α: J-shaped
• 0 < α, β < 1: U-shaped

The last three categories approach positive infinity at the boundaries (0 or 1), which can cause numerical problems in programming:

invalid value encountered in divide

The following demonstrates PDF values for various parameter combinations, particularly focusing on behavior near the boundaries:

import scipy.stats as stats
import numpy as np

# Boundary thresholds
delta_small = 1e-7
delta_tiny = 1e-8

# Test points spanning critical regions
sample_points = np.array([
    -1,  # much less than 0
    -delta_small, -delta_tiny, 0, delta_tiny, delta_small,  # near 0
    1 - delta_small, 1 - delta_tiny, 1, 1 + delta_tiny, 1 + delta_small,  # near 1
    2,  # much greater than 1
], dtype=np.float32)

print("Sample points:", sample_points)

print("\n--- Invalid region: alpha, beta < 0 ---")
print("alpha < 0:", stats.beta.pdf(sample_points, -0.5, 1))
print("beta < 0:", stats.beta.pdf(sample_points, 1, -0.5))

print("\n--- U-shaped: 0 < alpha, beta < 1 ---")
print(stats.beta.pdf(sample_points, 0.5, 0.5))

print("\n--- L-shaped: 0 < alpha < 1 <= beta ---")
print(stats.beta.pdf(sample_points, 0.5, 1))

print("\n--- J-shaped: 0 < beta < 1 <= alpha ---")
print(stats.beta.pdf(sample_points, 1, 0.5))

print("\n--- Constant: alpha = beta = 1 ---")
print(stats.beta.pdf(sample_points, 1, 1))

print("\n--- Bell-shaped (unimodal): 1 < alpha, beta ---")
print(stats.beta.pdf(sample_points, 2, 2))

Sample output:

[-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]

To avoid infinite values, one solution is to clamp input values to a safe range before computing the PDF. A boundary threshold of ε = 1e-7 works well in practice.

The following verifise numerical stability of numpy.clip for different epsilon values:

import numpy as np

# Generate mixed data with zeros and ones
padding_zeros = np.zeros([500], dtype=np.float32)
padding_ones = np.ones([500], dtype=np.float32)
mixed_data = np.concatenate([padding_zeros, padding_ones], axis=0)

# Test numpy.clip stability across different thresholds
for threshold in (1e-7, 1e-8):
    print(f"Testing threshold: {threshold}")
    for iteration in range(100):
        clipped = np.clip(mixed_data.copy(), threshold, 1 - threshold)
        # After clipping, no values should remain at exactly 0 or 1
        assert (0 != clipped).all() and (1 != clipped).all()
    print(f"  Passed: threshold {threshold} is stable")

Results show that 1e-7 provides stable clipping behavior, while 1e-8 leads to numerical instability where boundary values may not be properly clipped.