Python scipy.stats.wilcoxon() Examples

def compare_det_test(results_df, detector_name1, detector_name2, experiment=None):
    se1 = get_sensivities(results_df, detector_name1, experiment)
    if len(se1) < 2:
        return 0
    se2 = get_sensivities(results_df, detector_name2, experiment)
    if len(se2) < 2:
        return 0
    l = min(len(se1),len(se2))
    if (l < 20):
        return None
    #print("1:",se1[:l])
    #print("2:",se2[:l])
    try:
        t,p = stats.wilcoxon(se1[:l],se2[:l])
        return p
    except:
        return None 

def main(_):
  stats = FLAGS.stats.split(',')
  root = FLAGS.root
  model_names = FLAGS.model_names.split(',')

  data = {}
  for m in model_names:
    d = load_data(root, m)
    merge_into(data, d)

  print '\nLOADED %d models, %d examples' % (len(data['models']),
                                             len(data['examples']))

  if 'mean' in stats:
    compute_mean(data)
  if 'rank' in stats:
    compute_rank(data)
  if 'diff' in stats:
    compute_diff(data)
  if 'wilcoxon' in stats:
    compute_wilcoxon(data)
  print 

def test_wilcoxon_tie():
    # Regression test for gh-2391.
    # Corresponding R code is:
    #   > result = wilcox.test(rep(0.1, 10), exact=FALSE, correct=FALSE)
    #   > result$p.value
    #   [1] 0.001565402
    #   > result = wilcox.test(rep(0.1, 10), exact=FALSE, correct=TRUE)
    #   > result$p.value
    #   [1] 0.001904195
    stat, p = stats.wilcoxon([0.1] * 10)
    expected_p = 0.001565402
    assert_equal(stat, 0)
    assert_allclose(p, expected_p, rtol=1e-6)

    stat, p = stats.wilcoxon([0.1] * 10, correction=True)
    expected_p = 0.001904195
    assert_equal(stat, 0)
    assert_allclose(p, expected_p, rtol=1e-6) 

def compare_experiments(results_df1,det1):
    se1 = get_sensivities(results_df1, det1, 'sitting')
    if len(se1) < 2:
        return 0
    se2 = get_sensivities(results_df1, det1, 'jogging')
    if len(se2) < 2:
        return 0
    l = min(len(se1),len(se2))
    if (l < 20):
        return None
    try:
        t,p = stats.wilcoxon(se1[:l],se2[:l])
        return p
    except:
        return None 

def wilcoxon_test(v1,v2):# original metric: the smaller the more similar 
    result = stats.wilcoxon(v1, v2).statistic
    if result != result:
        result = 0
    return 1/(math.sqrt(result)+1) 

def compare_cables(results_df1,results_df2,det1,experiment):
    se1 = get_sensivities(results_df1, det1, experiment)
    if len(se1) < 2:
        return 0
    se2 = get_sensivities(results_df2, det1, experiment)
    if len(se2) < 2:
        return 0
    l = min(len(se1),len(se2))
    if (l < 20):
        return None
    try:
        t,p = stats.wilcoxon(se1[:l],se2[:l])
        return p
    except:
        return None 

def t_test_raw(
        x0,
        x1,
):
    """
    Perform two-sided t-test allowing for unequal group variances (Welch's t-test) on raw data
    along second axis (for each gene).

    The t-test assumes normally distributed data. This function computes
    the necessary statistics and calls t_test_moments().

    :param x0: np.array (observations x genes)
        Observations in first group by gene
    :param x1:  np.array (observations x genes)
        Observations in second group by gene.
    """
    axis = 1
    if x0.shape[1] != x1.shape[1]:
        raise ValueError('stats.t_test_raw(): x0 and x1 have to contain the same number of genes')
    # Reshape into 2D array if only one test is performed.
    if np.ndim(x0) == 1:
        x0 = x0.reshape([x0.shape[0], 1])
        x1 = x1.reshape([x1.shape[0], 1])
    if np.any(x0.shape[axis] != x1.shape[axis]):
        raise ValueError(
            'stats.wilcoxon(): the first axis (number of tests) is not allowed to differ between x0 and x1!')

    mu0 = np.asarray(np.mean(x0, axis=0)).flatten()
    var0 = np.asarray(np.var(x0, axis=0)).flatten()
    mu1 = np.asarray(np.mean(x1, axis=0)).flatten()
    var1 = np.asarray(np.var(x1, axis=0)).flatten()
    n0 = x0.shape[0]
    n1 = x1.shape[0]

    pval = t_test_moments(mu0=mu0, mu1=mu1, var0=var0, var1=var1, n0=n0, n1=n1)
    return pval 

def compute_wilcoxon(data):
  print '\nWILCOXON SIGNED-RANK TEST\n'
  # We take the first model as the basis for each comparison.
  ssim = np.array(data['ssim'])
  psnr = np.array(data['psnr'])
  for i, m in enumerate(data['models']):
    if i == 0:
      print '    [differences from %s]' % m
      continue
    ssim_v, ssim_p = wilcoxon(ssim[:, i], ssim[:, 0])
    psnr_v, psnr_p = wilcoxon(psnr[:, i], psnr[:, 0])
    print '%30s    ssim %.3f, p %.1e    psnr %.2f, p %.1e' % (m, ssim_v, ssim_p,
                                                              psnr_v, psnr_p) 

def score_event(nx, ny, ncommon, barcode_frequencies, resamples=100):
    """
    perform a resampling test, based on the number of fragments in each region
    and the barcode frequency distribution (the latter because some barcodes
    have more DNA in them than others)
    """

    samples = []
    for _ in range(resamples):
        s1 = numpy.random.choice(len(barcode_frequencies), nx, replace=False, p=barcode_frequencies)
        s2 = numpy.random.choice(len(barcode_frequencies), ny, replace=False, p=barcode_frequencies)
        common = len(set(s1).intersection(set(s2)))
        samples.append(common)

    # make a one-sided one-sample Wilcoxon test
    statistic, pvalue = stats.wilcoxon(numpy.array(samples)-ncommon)
    if numpy.mean(samples) > ncommon:
        pvalue = 1 - pvalue/2.0
    else:
        pvalue = pvalue/2.0

    # result = r["wilcox.test"](samples, mu=ncommon, alternative="less")
    # pvalue2 = result.rx2("p.value")[0]

    # print "::::::", nx, ny, ncommon, (numpy.array(samples)>ncommon).sum(), pvalue, pvalue2
    return pvalue 

def getActiveTF(self,dTD,dTG,dMb):
                # the idea : the targets for given TFs are signiciantly differenent between the fromNode and toNode
                HGL=[item.upper() for item in self.GL]
                pv_cut=1e-1
                ATF=[]

                FE=self.fromNode.getAvgEx()
                TE=self.toNode.getAvgEx()
                for i in dMb:
                        tfi=dMb[i]
                        tfi=[HGL.index(item) for item in tfi]

                        fe=[FE[item] for item in tfi]
                        te=[TE[item] for item in tfi]

                        #fe=[self.fromNode.E[item] for item in tfi]
                        #te=[self.toNode.E[item] for item in tfi]
                        #fe=reduce(lambda x,y:x+y, [[item.E[j] for item in self.fromNode.cells] for j in tfi])
                        #te=reduce(lambda x,y:x+y, [[item.E[j] for item in self.toNode.cells] for j in tfi])
                        #pdb.set_trace()
                        pv=wilcoxon(fe,te)[-1]
                        if pv<pv_cut:
                                ATF.append([pv,i])
                ATF.sort()
                return [item[1] for item in ATF]

        #-------------------------------------------------------------------
        # regresion model for each path 

def wilcoxon_similar_distribution(self, column,
                                       pvalue_threshold=0.05,
                                       num_rounds=3):
        p_value = permutation_test(
            self.new_data[column],
            self.historical_data[column],
            method="approximate",
            num_rounds=num_rounds,
            func=lambda x, y: stats.wilcoxon(x, y).statistic,
            seed=0)
        if p_value < pvalue_threshold:
            return False
        return True 

def test_wilcoxon_result_attributes():
    x = np.array([120, 114, 181, 188, 180, 146, 121, 191, 132, 113, 127, 112])
    y = np.array([133, 143, 119, 189, 112, 199, 198, 113, 115, 121, 142, 187])
    res = stats.wilcoxon(x, y, correction=False)
    attributes = ('statistic', 'pvalue')
    check_named_results(res, attributes) 

def test_accuracy_wilcoxon():
    freq = [1, 4, 16, 15, 8, 4, 5, 1, 2]
    nums = range(-4, 5)
    x = np.concatenate([[u] * v for u, v in zip(nums, freq)])
    y = np.zeros(x.size)

    T, p = stats.wilcoxon(x, y, "pratt")
    assert_allclose(T, 423)
    assert_allclose(p, 0.00197547303533107)

    T, p = stats.wilcoxon(x, y, "zsplit")
    assert_allclose(T, 441)
    assert_allclose(p, 0.0032145343172473055)

    T, p = stats.wilcoxon(x, y, "wilcox")
    assert_allclose(T, 327)
    assert_allclose(p, 0.00641346115861)

    # Test the 'correction' option, using values computed in R with:
    # > wilcox.test(x, y, paired=TRUE, exact=FALSE, correct={FALSE,TRUE})
    x = np.array([120, 114, 181, 188, 180, 146, 121, 191, 132, 113, 127, 112])
    y = np.array([133, 143, 119, 189, 112, 199, 198, 113, 115, 121, 142, 187])
    T, p = stats.wilcoxon(x, y, correction=False)
    assert_equal(T, 34)
    assert_allclose(p, 0.6948866, rtol=1e-6)
    T, p = stats.wilcoxon(x, y, correction=True)
    assert_equal(T, 34)
    assert_allclose(p, 0.7240817, rtol=1e-6) 

def test_wilcoxon_arg_type():
    # Should be able to accept list as arguments.
    # Address issue 6070.
    arr = [1, 2, 3, 0, -1, 3, 1, 2, 1, 1, 2]

    _ = stats.wilcoxon(arr, zero_method="pratt")
    _ = stats.wilcoxon(arr, zero_method="zsplit")
    _ = stats.wilcoxon(arr, zero_method="wilcox") 

def test_wilcoxon_bad_arg():
    # Raise ValueError when two args of different lengths are given or
    # zero_method is unknown.
    assert_raises(ValueError, stats.wilcoxon, [1], [1, 2])
    assert_raises(ValueError, stats.wilcoxon, [1, 2], [1, 2], "dummy") 

def test_accuracy_wilcoxon():
    freq = [1, 4, 16, 15, 8, 4, 5, 1, 2]
    nums = range(-4, 5)
    x = np.concatenate([[u] * v for u, v in zip(nums, freq)])
    y = np.zeros(x.size)

    T, p = stats.wilcoxon(x, y, "pratt")
    assert_allclose(T, 423)
    assert_allclose(p, 0.00197547303533107)

    T, p = stats.wilcoxon(x, y, "zsplit")
    assert_allclose(T, 441)
    assert_allclose(p, 0.0032145343172473055)

    T, p = stats.wilcoxon(x, y, "wilcox")
    assert_allclose(T, 327)
    assert_allclose(p, 0.00641346115861)

    # Test the 'correction' option, using values computed in R with:
    # > wilcox.test(x, y, paired=TRUE, exact=FALSE, correct={FALSE,TRUE})
    x = np.array([120, 114, 181, 188, 180, 146, 121, 191, 132, 113, 127, 112])
    y = np.array([133, 143, 119, 189, 112, 199, 198, 113, 115, 121, 142, 187])
    T, p = stats.wilcoxon(x, y, correction=False)
    assert_equal(T, 34)
    assert_allclose(p, 0.6948866, rtol=1e-6)
    T, p = stats.wilcoxon(x, y, correction=True)
    assert_equal(T, 34)
    assert_allclose(p, 0.7240817, rtol=1e-6) 

def test_wilcoxon_bad_arg():
    """Raise ValueError when two args of different lengths are given or
       zero_method is unknwon"""
    assert_raises(ValueError, stats.wilcoxon, [1], [1,2])
    assert_raises(ValueError, stats.wilcoxon, [1,2], [1,2], "dummy") 

def symtest(alpha, axis=None):
    """
    Non-parametric test for symmetry around the median. Works by performing a
    Wilcoxon sign rank test on the differences to the median.

    H0: the population is symmetrical around the median
    HA: the population is not symmetrical around the median


    :param alpha: sample of angles in radian
    :param axis:  compute along this dimension, default is None
                  if axis=None, array is raveled
    :return pval: two-tailed p-value
    :return T:    test statistics of underlying wilcoxon test


    References: [Zar2009]_
    """

    m = descriptive.median(alpha, axis=axis)

    d = np.angle(np.exp(1j * m[np.newaxis]) / np.exp(1j * alpha))

    if axis is not None:
        oshape = d.shape[1:]
        d2 = d.reshape((d.shape[0], int(np.prod(d.shape[1:]))))
        T, pval = map(lambda x: np.asarray(x).reshape(
            oshape), zip(*[stats.wilcoxon(dd) for dd in d2.T]))
    else:
        T, pval = stats.wilcoxon(d)

    return pval, T 

def pearsonr_test(v1,v2):#original metric: the larger the more similar
    result = stats.mstats.pearsonr(v1,v2)[0]
    result = stats.wilcoxon(v1, v2).statistic
    if result != result:
        result = -1
    return sigmoid(result) 

def spearmanr_test(v1,v2):#original metric: the larger the more similar 
    result = stats.mstats.spearmanr(v1,v2).correlation
    result = stats.wilcoxon(v1, v2).statistic
    if result != result:
        result = -1
    return sigmoid(result) 

def entroy_test(v1,v2):#original metric: the smaller the more similar
    result = stats.entropy(v1,v2)
    result = stats.wilcoxon(v1, v2).statistic
    if result != result:
        result = 0
    return result 

def get_valid_tests(equal_var, independent, normal, num_samples):
    '''
    Get valid tests given number of samples and statistical characterization of
    samples:

    Equal variance
    Indepenence
    Normality
    '''
    if num_samples == 1:
        valid_tests = {
            'chisquare': stats.chisquare,
            'power_divergence': stats.power_divergence,
            'kstest': stats.kstest
        }
        if normal:
            valid_tests['input']['one_sample_ttest'] = stats.ttest_1samp

    elif num_samples == 2:
        if independent:
            valid_tests = {
                'mannwhitneyu': stats.mannwhitneyu,
                'kruskal': stats.kruskal,
                'ks_2samp': stats.ks_2samp
            }
            if normal:
                valid_tests['two_sample_ttest'] = stats.ttest_ind
                if equal_var:
                    valid_tests['f_oneway'] = stats.f_oneway
        else:
            valid_tests = {
                'two_sample_ks': stats.ks_2samp,
                'wilcoxon': stats.wilcoxon
            }
            if normal:
                valid_tests['two_sample_related_ttest'] = stats.ttest_rel

    elif num_samples >= 3:
        if independent:
            valid_tests = {
                'kruskal': stats.kruskal
            }
            if normal and equal_var:
                valid_tests['f_oneway'] = stats.f_oneway

        else:
            valid_tests['friedmanchisquare'] = stats.friedmanchisquare

    return valid_tests 

def wilcoxon_test(self, metric_name=None):
        """http://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test
        `Wilcoxon signed-rank test
        <https://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test>`_.
        Tests whether two  related paired samples come from the same
        distribution.
        In particular, it tests whether the distribution of the differences
        x-y is symmetric about zero
        """
        self._check_is_evaluated()
        metric_name = self._validate_metric_name(metric_name)
        metrics_per_estimator_dataset = \
            self._get_metrics_per_estimator_dataset(
                metric_name)

        wilcoxon_df = pd.DataFrame()
        values = np.array([])
        prod = itertools.product(metrics_per_estimator_dataset.keys(),
                                 repeat=2)
        for p in prod:
            estim_1 = p[0]
            estim_2 = p[1]
            w, p_val = stats.wilcoxon(metrics_per_estimator_dataset[p[0]],
                                      metrics_per_estimator_dataset[p[1]])

            w_test = {
                "estimator_1": estim_1,
                "estimator_2": estim_2,
                "statistic": w,
                "p_val": p_val
            }

            wilcoxon_df = wilcoxon_df.append(w_test, ignore_index=True)
            values = np.append(values, w)
            values = np.append(values, p_val)

        index = wilcoxon_df["estimator_1"].unique()
        values_names = ["statistic", "p_val"]
        col_idx = pd.MultiIndex.from_product([index, values_names])
        values_reshaped = values.reshape(len(index),
                                         len(values_names) * len(index))

        values_df_multiindex = pd.DataFrame(values_reshaped, index=index,
                                            columns=col_idx)

        return wilcoxon_df, values_df_multiindex 

def responsive_units(spike_times, spike_clusters, event_times,
                     pre_time=[0.5, 0], post_time=[0, 0.5], alpha=0.05):
    """
    Determine responsive neurons by doing a Wilcoxon Signed-Rank test between a baseline period
    before a certain task event (e.g. stimulus onset) and a period after the task event.

    Parameters
    ----------
    spike_times : 1D array
        spike times (in seconds)
    spike_clusters : 1D array
        cluster ids corresponding to each event in `spikes`
    event_times : 1D array
        times (in seconds) of the events from the two groups
    pre_time : two-element array
        time (in seconds) preceding the event to get the baseline (e.g. [0.5, 0.2] would be a
        window starting 0.5 seconds before the event and ending at 0.2 seconds before the event)
    post_time : two-element array
        time (in seconds) to follow the event times
    alpha : float
        alpha to use for statistical significance

    Returns
    -------
    significant_units : ndarray
        an array with the indices of clusters that are significatly modulated
    stats : 1D array
        the statistic of the test that was performed
    p_values : ndarray
        the p-values of all the clusters
    cluster_ids : ndarray
        cluster ids of the p-values
    """

    # Get spike counts for baseline and event timewindow
    baseline_times = np.column_stack(((event_times - pre_time[0]), (event_times - pre_time[1])))
    baseline_counts, cluster_ids = _get_spike_counts_in_bins(spike_times, spike_clusters,
                                                             baseline_times)
    times = np.column_stack(((event_times + post_time[0]), (event_times + post_time[1])))
    spike_counts, cluster_ids = _get_spike_counts_in_bins(spike_times, spike_clusters, times)

    # Do statistics
    p_values = np.empty(spike_counts.shape[0])
    stats = np.empty(spike_counts.shape[0])
    for i in range(spike_counts.shape[0]):
        if np.sum(baseline_counts[i, :] - spike_counts[i, :]) == 0:
            p_values[i] = 1
            stats[i] = 0
        else:
            stats[i], p_values[i] = wilcoxon(baseline_counts[i, :], spike_counts[i, :])

    # Perform FDR correction for multiple testing
    sig_units, p_values, _, _ = multipletests(p_values, alpha)
    significant_units = cluster_ids[sig_units]

    return significant_units, stats, p_values, cluster_ids