def pearson_multi(X):
    """
    Compute Pearson correlation across a set of samples.

    Parameters
    ----------
    X: `numpy.ndarray` or `scipy.sparse.csr_matrix`
        Matrix with rows corresponding to all of the samples and columns
        corresponding to features that define the correlation matrices.

    Returns
    -------
    corr: `np.ndarray`
        Dense Pearson coexpression matrix with dimension
        `(X.shape[1], X.shape[1])`.
    """
    from scipy.special import betainc
    from sklearn.covariance import EmpiricalCovariance
    from sklearn.preprocessing import StandardScaler
    X = StandardScaler().fit_transform(X)
    corr = EmpiricalCovariance(store_precision=False, assume_centered=True
        ).fit(X).covariance_
    corr[corr > 1.0] = 1.0
    corr[corr < -1.0] = -1.0
    return corr
