def fdr(pvals, alpha=0.05, method='fdr_bh'):
    """P-values FDR correction with Benjamini/Hochberg and
    Benjamini/Yekutieli procedure.

    This covers Benjamini/Hochberg for independent or positively correlated and
    Benjamini/Yekutieli for general or negatively correlated tests.

    Parameters
    ----------
    pvals : array_like
        Array of p-values of the individual tests.
    alpha : float
        Error rate (= alpha level).
    method : str
        FDR correction methods ::

        'fdr_bh' : Benjamini/Hochberg for independent / posit correlated tests
        'fdr_by' : Benjamini/Yekutieli for negatively correlated tests

    Returns
    -------
    reject : array, bool
        True if a hypothesis is rejected, False if not
    pval_corrected : array
        P-values adjusted for multiple hypothesis testing using the BH or BY
        correction.

    See also
    --------
    bonf : Bonferroni correction
    holm : Holm-Bonferroni correction

    Notes
    -----
    From Wikipedia:

    The **Benjamini–Hochberg** procedure (BH step-up procedure) controls the
    false discovery rate (FDR) at level :math:`\\alpha`. It works as follows:

    1. For a given :math:`\\alpha`, find the largest :math:`k` such that
    :math:`P_{(k)}\\leq \\frac {k}{m}\\alpha.`

    2. Reject the null hypothesis (i.e., declare discoveries) for all
    :math:`H_{(i)}` for :math:`i = 1, \\ldots, k`.

    The BH procedure is valid when the m tests are independent, and also in
    various scenarios of dependence, but is not universally valid.

    The **Benjamini–Yekutieli** procedure (BY) controls the FDR under arbitrary
    dependence assumptions. This refinement modifies the threshold and finds
    the largest :math:`k` such that:

    .. math::
        P_{(k)} \\leq \\frac{k}{m \\cdot c(m)} \\alpha

    References
    ----------
    - Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery
      rate: a practical and powerful approach to multiple testing. Journal of
      the Royal Statistical Society Series B, 57, 289–300.

    - Benjamini, Y., and Yekutieli, D. (2001). The control of the false
      discovery rate in multiple testing under dependency. Annals of
      Statistics, 29, 1165–1188.

    - https://en.wikipedia.org/wiki/False_discovery_rate

    Examples
    --------
    FDR correction of an array of p-values

    >>> from pingouin import fdr
    >>> pvals = [.50, .003, .32, .054, .0003]
    >>> reject, pvals_corr = fdr(pvals, alpha=.05)
    >>> print(reject, pvals_corr)
    [False  True False False  True] [0.5    0.0075 0.4    0.09   0.0015]
    """
    assert method.lower() in ['fdr_bh', 'fdr_by']
    # Convert to array and save original shape
    pvals = np.asarray(pvals)
    shape_init = pvals.shape
    pvals = pvals.ravel()
    num_nan = np.isnan(pvals).sum()

    # Sort the (flattened) p-values
    pvals_sortind = np.argsort(pvals)
    pvals_sorted = pvals[pvals_sortind]
    sortrevind = pvals_sortind.argsort()
    ntests = pvals.size - num_nan

    # Empirical CDF factor
    ecdffactor = np.arange(1, ntests + 1) / float(ntests)

    if method.lower() == 'fdr_by':
        cm = np.sum(1. / np.arange(1, ntests + 1))
        ecdffactor /= cm

    # Now we adjust the p-values
    pvals_corr = np.diag(pvals_sorted / ecdffactor[..., None])
    pvals_corr = np.minimum.accumulate(pvals_corr[::-1])[::-1]
    pvals_corr = np.clip(pvals_corr, None, 1)

    # And revert to the original shape and order
    pvals_corr = np.append(pvals_corr, np.full(num_nan, np.nan))
    pvals_corrected = pvals_corr[sortrevind].reshape(shape_init)
    with np.errstate(invalid='ignore'):
        reject = np.less(pvals_corrected, alpha)
    # reject = reject[sortrevind].reshape(shape_init)
    return reject, pvals_corrected
P-values FDR correction with Benjamini/Hochberg and
    Benjamini/Yekutieli procedure.

    This covers Benjamini/Hochberg for independent or positively correlated and
    Benjamini/Yekutieli for general or negatively correlated tests.

    Parameters
    ----------
    pvals : array_like
        Array of p-values of the individual tests.
    alpha : float
        Error rate (= alpha level).
    method : str
        FDR correction methods ::

        'fdr_bh' : Benjamini/Hochberg for independent / posit correlated tests
        'fdr_by' : Benjamini/Yekutieli for negatively correlated tests

    Returns
    -------
    reject : array, bool
        True if a hypothesis is rejected, False if not
    pval_corrected : array
        P-values adjusted for multiple hypothesis testing using the BH or BY
        correction.

    See also
    --------
    bonf : Bonferroni correction
    holm : Holm-Bonferroni correction

    Notes
    -----
    From Wikipedia:

    The **Benjamini–Hochberg** procedure (BH step-up procedure) controls the
    false discovery rate (FDR) at level :math:`\\alpha`. It works as follows:

    1. For a given :math:`\\alpha`, find the largest :math:`k` such that
    :math:`P_{(k)}\\leq \\frac {k}{m}\\alpha.`

    2. Reject the null hypothesis (i.e., declare discoveries) for all
    :math:`H_{(i)}` for :math:`i = 1, \\ldots, k`.

    The BH procedure is valid when the m tests are independent, and also in
    various scenarios of dependence, but is not universally valid.

    The **Benjamini–Yekutieli** procedure (BY) controls the FDR under arbitrary
    dependence assumptions. This refinement modifies the threshold and finds
    the largest :math:`k` such that:

    .. math::
        P_{(k)} \\leq \\frac{k}{m \\cdot c(m)} \\alpha

    References
    ----------
    - Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery
      rate: a practical and powerful approach to multiple testing. Journal of
      the Royal Statistical Society Series B, 57, 289–300.

    - Benjamini, Y., and Yekutieli, D. (2001). The control of the false
      discovery rate in multiple testing under dependency. Annals of
      Statistics, 29, 1165–1188.

    - https://en.wikipedia.org/wiki/False_discovery_rate

    Examples
    --------
    FDR correction of an array of p-values

    >>> from pingouin import fdr
    >>> pvals = [.50, .003, .32, .054, .0003]
    >>> reject, pvals_corr = fdr(pvals, alpha=.05)
    >>> print(reject, pvals_corr)
    [False  True False False  True] [0.5    0.0075 0.4    0.09   0.0015]