Python numpy.msort方法代码示例

# 需要导入模块: import numpy [as 别名]
# 或者: from numpy import msort [as 别名]
def msort(a):
    """Returns a copy of an array sorted along the first axis.

    Args:
        a (cupy.ndarray): Array to be sorted.

    Returns:
        cupy.ndarray: Array of the same type and shape as ``a``.

    .. note:
        ``cupy.msort(a)``, the CuPy counterpart of ``numpy.msort(a)``, is
        equivalent to ``cupy.sort(a, axis=0)``.

    .. seealso:: :func:`numpy.msort`

    """

    return sort(a, axis=0) 

# 需要导入模块: import numpy [as 别名]
# 或者: from numpy import msort [as 别名]
def msort(a):
    """
    Return a copy of an array sorted along the first axis.

    Parameters
    ----------
    a : array_like
        Array to be sorted.

    Returns
    -------
    sorted_array : ndarray
        Array of the same type and shape as `a`.

    See Also
    --------
    sort

    Notes
    -----
    ``np.msort(a)`` is equivalent to  ``np.sort(a, axis=0)``.

    """
    b = array(a, subok=True, copy=True)
    b.sort(0)
    return b 

# 需要导入模块: import numpy [as 别名]
# 或者: from numpy import msort [as 别名]
def bootstrap_ci(estimator, data, alpha, n_sets=None, sort=numpy.msort, eargs=(), ekwargs={}):
    '''Perform a Monte Carlo bootstrap of a (1-alpha) confidence interval for the given ``estimator``.
    Returns (fhat, ci_lower, ci_upper), where fhat is the result of ``estimator(data, *eargs, **ekwargs)``,
    and ``ci_lower`` and ``ci_upper`` are the lower and upper bounds of the surrounding confidence
    interval, calculated by calling ``estimator(syndata, *eargs, **ekwargs)`` on each synthetic data
    set ``syndata``.  If ``n_sets`` is provided, that is the number of synthetic data sets generated,
    otherwise an appropriate size is selected automatically (see ``calc_mcbs_nsets()``).
    
    ``sort``, if given, is applied to sort the results of calling ``estimator`` on each 
    synthetic data set prior to obtaining the confidence interval. This function must sort
    on the last index.
    
    Individual entries in synthetic data sets are selected by the first index of ``data``, allowing this
    function to be used on arrays of multidimensional data.
    
    Returns (fhat, lb, ub, ub-lb, abs((ub-lb)/fhat), and max(ub-fhat,fhat-lb)) (that is, the estimated value, the
    lower and upper bounds of the confidence interval, the width of the confidence interval, the relative
    width of the confidence interval, and the symmetrized error bar of the confidence interval).'''

    data = numpy.asanyarray(data)
    fhat = numpy.squeeze(estimator(data, *eargs, **ekwargs))
    n_sets = n_sets or calc_mcbs_nsets(alpha)
    fsynth = numpy.empty((n_sets,), dtype=fhat.dtype)
    try:
        return bootstrap_ci_ll(estimator, data, alpha, n_sets or calc_mcbs_nsets(alpha), fsynth, sort, eargs, ekwargs, fhat)
    finally:
        del fsynth 

# 需要导入模块: import numpy [as 别名]
# 或者: from numpy import msort [as 别名]
def test_msort(self):
        self.check(np.msort) 

# 需要导入模块: import numpy [as 别名]
# 或者: from numpy import msort [as 别名]
def bootstrap_ci(estimator, data, alpha, n_sets=None, args=(), kwargs={}, sort=numpy.msort, extended_output = False):
    '''Perform a Monte Carlo bootstrap of a (1-alpha) confidence interval for the given ``estimator``.
    Returns (fhat, ci_lower, ci_upper), where fhat is the result of ``estimator(data, *args, **kwargs)``,
    and ``ci_lower`` and ``ci_upper`` are the lower and upper bounds of the surrounding confidence
    interval, calculated by calling ``estimator(syndata, *args, **kwargs)`` on each synthetic data
    set ``syndata``.  If ``n_sets`` is provided, that is the number of synthetic data sets generated,
    otherwise an appropriate size is selected automatically (see ``get_bssize()``).
    
    ``sort``, if given, is applied to sort the results of calling ``estimator`` on each 
    synthetic data set prior to obtaining the confidence interval.
    
    Individual entries in synthetic data sets are selected by the first index of ``data``, allowing this
    function to be used on arrays of multidimensional data.
    
    If ``extended_output`` is True (by default not), instead of returning (fhat, lb, ub), this function returns
    (fhat, lb, ub, ub-lb, abs((ub-lb)/fhat), and max(ub-fhat,fhat-lb)) (that is, the estimated value, the
    lower and upper bounds of the confidence interval, the width of the confidence interval, the relative
    width of the confidence interval, and the symmetrized error bar of the confidence interval).'''
    
    data = numpy.asanyarray(data)
    
    fhat = estimator(data, *args, **kwargs)
    
    try:
        estimator_shape = fhat.shape
    except AttributeError:
        estimator_shape = ()
        
    try:
        estimator_dtype = fhat.dtype
    except AttributeError:
        estimator_dtype = type(fhat) 
        
    dlen = len(data)
    n_sets = n_sets or get_bssize(alpha)
    
    f_synth = numpy.empty((n_sets,) + estimator_shape, dtype=estimator_dtype)
    
    for i in range(0, n_sets):
        indices = numpy.random.randint(dlen, size=(dlen,))
        f_synth[i] = estimator(data[indices], *args, **kwargs)
        
    f_synth_sorted = sort(f_synth)
    lbi = int(math.floor(n_sets*alpha/2))
    ubi = int(math.ceil(n_sets*(1-alpha/2)))
    lb = f_synth_sorted[lbi]
    ub = f_synth_sorted[ubi]
    
    try:
        if extended_output:
            return (fhat, lb, ub, ub-lb, abs((ub-lb)/fhat) if fhat else 0, max(ub-fhat,fhat-lb))
        else:
            return (fhat, lb, ub)
    finally:
        # Do a little explicit memory management
        del f_synth, f_synth_sorted 

# 需要导入模块: import numpy [as 别名]
# 或者: from numpy import msort [as 别名]
def derive_threshold(auc_mtx: pd.DataFrame, regulon_name: str, seed=None, method: str = 'hdt') -> float:
    '''
    Derive threshold on the AUC values of the given regulon to binarize the cells in two clusters: "on" versus "off"
    state of the regulator.

    :param auc_mtx: The dataframe with the AUC values for all cells and regulons (n_cells x n_regulons).
    :param regulon_name: the name of the regulon for which to predict the threshold.
    :param method: The method to use to decide if the distribution of AUC values for the given regulon is not unimodel.
        Can be either Hartigan's Dip Test (HDT) or Bayesian Information Content (BIC). The former method performs better
        but takes considerable more time to execute (40min for 350 regulons). The BIC compares the BIC for two Gaussian
        Mixture Models: single versus two components.
    :return: The threshold on the AUC values.
    '''
    assert auc_mtx is not None and not auc_mtx.empty
    assert regulon_name in auc_mtx.columns
    assert method in {'hdt', 'bic'}

    data = auc_mtx[regulon_name].values

    if seed:
        np.random.seed(seed=seed)

    def isbimodal(data, method):
        if method == 'hdt':
            # Use Hartigan's dip statistic to decide if distribution deviates from unimodality.
            _, pval, _ = diptst(np.msort(data))
            return (pval is not None) and (pval <= 0.05)
        else:
            # Compare Bayesian Information Content of two Gaussian Mixture Models.
            X = data.reshape(-1, 1)
            gmm2 = mixture.GaussianMixture(n_components=2, covariance_type='full', random_state=seed).fit(X)
            gmm1 = mixture.GaussianMixture(n_components=1, covariance_type='full', random_state=seed).fit(X)
            return gmm2.bic(X) <= gmm1.bic(X)

    if not isbimodal(data, method):
        # For a unimodal distribution the threshold is set as mean plus two standard deviations.
        return data.mean() + 2.0*data.std()
    else:
        # Fit a two component Gaussian Mixture model on the AUC distribution using an Expectation-Maximization algorithm
        # to identify the peaks in the distribution.
        gmm2 = mixture.GaussianMixture(n_components=2, covariance_type='full', random_state=seed).fit(data.reshape(-1, 1))
        # For a bimodal distribution the threshold is defined as the "trough" in between the two peaks.
        # This is solved as a minimization problem on the kernel smoothed density.
        return minimize_scalar(fun=stats.gaussian_kde(data),
                               bounds=sorted(gmm2.means_),
                               method='bounded').x[0]