def rao_blackwell_ledoit_wolf(S, n):
    """Rao-Blackwellized Ledoit-Wolf shrinkaged estimator of the covariance
    matrix.

    Parameters
    ----------
    S : array, shape=(n, n)
        Sample covariance matrix (e.g. estimated with np.cov(X.T))
    n : int
        Number of data points.

    Returns
    -------
    sigma : array, shape=(n, n)
    shrinkage : float

    References
    ----------
    .. [1] Chen, Yilun, Ami Wiesel, and Alfred O. Hero III. "Shrinkage
        estimation of high dimensional covariance matrices" ICASSP (2009)
    """
    p = len(S)
    assert S.shape == (p, p)

    alpha = (n-2)/(n*(n+2))
    beta = ((p+1)*n - 2) / (n*(n+2))

    trace_S2 = np.sum(S*S)  # np.trace(S.dot(S))
    U = ((p * trace_S2 / np.trace(S)**2) - 1)
    rho = min(alpha + beta/U, 1)

    F = (np.trace(S) / p) * np.eye(p)
    return (1-rho)*S + rho*F, rho