def METHOD_NAME(
    x: np.ndarray,
    a: Union[np.ndarray, float],
    b: Union[np.ndarray, float],
    loc: Union[np.ndarray, float] = 0,
    scale: Union[np.ndarray, float] = 1,
) -> np.ndarray:
    x = (x - loc) / scale
    x, a, b = np.atleast_1d(x, a, b)
    out = _norm_logpdf(x) - _log_gauss_mass(a, b)
    x, a, b = np.broadcast_arrays(x, a, b)
    out[(x < a) | (b < x)] = -np.inf
    out[a == b] = math.nan
    return out