def METHOD_NAME(
    vec: List[Any], alpha: float = 0.01, trials: int = 3000
) -> List[Tuple[int, float]]:
    """Returns the indices with values that are significant outliers, with their p-values"""
    m = len(vec)
    if m == 0:
        return []
    removed = 0
    results = []
    # pv = multinomial_pvalue(vec, trials)
    # Hack: for now, set pv to alpha because computing exact multinomial p-values is too expensive
    pv = alpha
    # We use the Benjamin-Yekutieli procedure to control false-discovery rate.
    # See https://en.wikipedia.org/wiki/False_discovery_rate#Benjamini%E2%80%93Yekutieli_procedure
    c_m = harmonic_number(m)
    if pv <= alpha:
        while removed < m:
            # While we remain below the threshold, remove (zero-out by
            # setting to NaN) the max and add its index to the list of
            # results with its p-value.
            max_index = argmax(vec)
            # See how unlikely this bin is to have occurred at random,
            # assuming a uniform distribution into bins.
            this_pvalue = one_sided_binomial_test_ge(
                int(np.nansum(vec)), vec[max_index], 1 / (m - removed)
            )
            # print("max_index = ", max_index, "p-value = ", this_pvalue)
            if this_pvalue <= (alpha * (removed + 1) / (m * c_m)):
                results.append((max_index, this_pvalue))
                vec[max_index] = np.nan
                removed += 1
            else:
                break
    return results