code: def rao_blackwell_ledoit_wolf(S, n):
      
      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
docs: 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)