Beta Distribution PDF and Statistics

import math
from typing import Dict

def beta_function(a: float, b: float) -> float:
    return math.gamma(a) * math.gamma(b) / math.gamma(a + b)

def beta_pdf(x: float, alpha: float, beta: float) -> float:
    if x < 0 or x > 1:
        return 0.0
    if alpha <= 0 or beta <= 0:
        raise ValueError("Alpha and beta must be positive")

    if x == 0:
        return float('inf') if alpha < 1 else (1.0 if alpha == 1 else 0.0)
    if x == 1:
        return float('inf') if beta < 1 else (1.0 if beta == 1 else 0.0)

    B = beta_function(alpha, beta)
    return (x ** (alpha - 1) * (1 - x) ** (beta - 1)) / B

def beta_statistics(alpha: float, beta: float) -> Dict[str, float]:
    mean = alpha / (alpha + beta)
    variance = (alpha * beta) / ((alpha + beta) ** 2 * (alpha + beta + 1))

    if alpha > 1 and beta > 1:
        mode = (alpha - 1) / (alpha + beta - 2)
    elif alpha <= 1 and beta <= 1 and alpha != 1 and beta != 1:
        mode = None  # Bimodal at 0 and 1
    elif alpha <= 1:
        mode = 0.0
    else:
        mode = 1.0

    return {
        "mean": round(mean, 4),
        "variance": round(variance, 6),
        "mode": round(mode, 4) if mode is not None else None,
        "std_dev": round(variance ** 0.5, 4)
    }