numpy.median

Here are the examples of the python api numpy.median taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.

200 Examples 7

Example 1

Project: bcbio-nextgen
Source File: shared.py
View license
def insert_size_stats(dists):
    """Calcualtes mean/median and MAD from distances, avoiding outliers.

    MAD is the Median Absolute Deviation: http://en.wikipedia.org/wiki/Median_absolute_deviation
    """
    med = numpy.median(dists)
    filter_dists = filter(lambda x: x < med + 10 * med, dists)
    median = numpy.median(filter_dists)
    return {"mean": float(numpy.mean(filter_dists)), "std": float(numpy.std(filter_dists)),
            "median": float(median),
            "mad": float(numpy.median([abs(x - median) for x in filter_dists]))}

Example 2

Project: pypot
Source File: sonar.py
View license
    def _filter(self, data):
        """ Apply a filter to reduce noisy data.

           Return the median value of a heap of data.

        """
        filtered_data = []
        for queue, data in zip(self._raw_data_queues, data):
            queue.append(data)
            filtered_data.append(numpy.median(queue))

        return filtered_data

Example 3

Project: python-oceans
Source File: RPSstuff.py
View license
def gmedian(x, **kw):
    """
    Just like median, except that it skips over bad points.

    """
    xnew = ma.masked_invalid(x)
    return np.median(xnew, **kw)

Example 4

Project: wikidump
Source File: analysis.py
View license
@shelved("length_stats")
def length_stats(prefix):
  doc_lengths = raw_doc_lengths(prefix).values()
  median = n.median(doc_lengths)
  mean   = n.mean(doc_lengths)
  std    = n.std(doc_lengths)
  return median, mean, std

Example 5

Project: cmonkey2
Source File: util.py
View license
def trim_mean(values, trim):
    """returns the trim mean"""
    if not values or len(values) == 0:
        return 0

    values = sorted(values, reverse=True)
    if trim == 0.5:
        return np.median(values)

    k_floor = int(math.floor(trim * len(values)))
    trim_values_floor = values[k_floor:len(values) - k_floor]
    return np.mean(np.ma.masked_array(trim_values_floor,
                                      np.isnan(trim_values_floor)))

Example 6

Project: scipy
Source File: test_binned_statistic.py
View license
    def test_1d_median(self):
        x = self.x
        v = self.v

        stat1, edges1, bc = binned_statistic(x, v, 'median', bins=10)
        stat2, edges2, bc = binned_statistic(x, v, np.median, bins=10)

        assert_allclose(stat1, stat2)
        assert_allclose(edges1, edges2)

Example 7

Project: scipy
Source File: test_binned_statistic.py
View license
    def test_2d_median(self):
        x = self.x
        y = self.y
        v = self.v

        stat1, binx1, biny1, bc = binned_statistic_2d(
            x, y, v, 'median', bins=5)
        stat2, binx2, biny2, bc = binned_statistic_2d(
            x, y, v, np.median, bins=5)

        assert_allclose(stat1, stat2)
        assert_allclose(binx1, binx2)
        assert_allclose(biny1, biny2)

Example 8

Project: scipy
Source File: test_binned_statistic.py
View license
    def test_dd_median(self):
        X = self.X
        v = self.v

        stat1, edges1, bc = binned_statistic_dd(X, v, 'median', bins=3)
        stat2, edges2, bc = binned_statistic_dd(X, v, np.median, bins=3)

        assert_allclose(stat1, stat2)
        assert_allclose(edges1, edges2)

Example 9

Project: theano_pyglm
Source File: fit_latent_network.py
View license
def compute_diff_of_dists(L_true, L_smpls):
    N_smpls = len(L_smpls)

    assert L_true.ndim == 2

    true_dist = compute_distances(L_true)[0]
    inf_dists = compute_distances(L_smpls)

    error_dists = [(inf_dist - true_dist)**2 for inf_dist in inf_dists]
    mse = [error_dist.mean() for error_dist in error_dists]
    medse = [np.median(error_dist.ravel()) for error_dist in error_dists]
    stdse = [error_dist.std() for error_dist in error_dists]

    fig = plt.figure()
    # plt.errorbar(np.arange(N_smpls), mse, yerr=stdse)
    plt.errorbar(np.arange(N_smpls), mse, color='b')
    plt.errorbar(np.arange(N_smpls), medse, color='r')
    plt.show()
    return mse, stdse

Example 10

Project: MIDAS
Source File: merge_species.py
View license
def compute_stats(args, data):
	species_ids = data.keys()
	stats = dict([(s,{}) for s in species_ids])
	for species_id in species_ids:
		# abundance
		x = data[species_id]['relative_abundance']
		stats[species_id]['median_abundance'] = np.median(x)
		stats[species_id]['mean_abundance'] = np.mean(x)
		# coverage
		x = data[species_id]['coverage']
		stats[species_id]['median_coverage'] = np.median(x)
		stats[species_id]['mean_coverage'] = np.mean(x)
		# prevalence
		stats[species_id]['prevalence'] = prevalence(x, y=args['min_cov'])
	return stats

Example 11

Project: MIDAS
Source File: genes.py
View license
def normalize(args, species, genes):
	""" Count number of bp mapped to each marker gene """
	# compute marker depth
	for gene in genes.values():
		if gene.marker_id is not None:
			species[gene.species_id].markers[gene.marker_id] += gene.depth
	# compute median marker depth
	for sp in species.values():
		sp.marker_coverage = np.median(sp.markers.values())
	# normalize genes by median marker depth
	for gene in genes.values():
		sp = species[gene.species_id]
		if sp.marker_coverage > 0:
			gene.copies = gene.depth/sp.marker_coverage

Example 12

Project: GromacsWrapper
Source File: timeseries.py
View license
def median_histogrammed_function(t, y, **kwargs):
    """Compute median of data *y* in bins along *t*.

    Returns the median-regularised function *F* and the centers of the bins.

    :func:`regularized_function` with *func* = :func:`numpy.median`
    """
    return apply_histogrammed_function(numpy.median, t, y, **kwargs)

Example 13

Project: livetv_mining
Source File: models.py
View license
    @property
    def median(self):
        last30days = datetime.utcnow() - timedelta(days=30)
        online_list = []
        for roomdata in self.hisdata.filter(LiveTVRoomData.create_date > last30days).all():
            online_list.append(roomdata.online)
        return {
            'online': np.median(online_list)
        }

Example 14

Project: mHTM
Source File: plot.py
View license
def compute_err(x, x_med=None, axis=1):
	"""
	Compute the relative errors based off IQR.
	
	@param x: A 2D numpy array containing the data across multiple iterations.
	Each row represents a unique instance of the data.
	
	@param x_med: A 1D numpy array representing the median values.
	
	@param axis: The axis around which to compute the data.
	
	@return: A 2D numpy array containing the -y and +y errors.
	"""
	
	if x_med is None:
		x_med = np.median(x, axis)
	
	return np.array([x_med - np.percentile(x, 25, axis),
		np.percentile(x, 75, axis) - x_med])

Example 15

Project: openstack-neat
Source File: statistics.py
View license
@contract
def mad(data):
    """ Calculate the Median Absolute Deviation from the data.

    :param data: The data to analyze.
     :type data: list(number)

    :return: The calculated MAD.
     :rtype: float
    """
    data_median = median(data)
    return float(median([abs(data_median - x) for x in data]))

Example 16

Project: openstack-neat
Source File: statistics.py
View license
@contract
def tricube_bisquare_weights(data):
    """ Generates a weights according to the tricube bisquare function.

    :param data: The input data.
     :type data: list(float)

    :return: A list of generated weights.
     :rtype: list(float)
    """
    n = len(data)
    s6 = 6 * median(map(abs, data))
    weights = tricube_weights(n)
    weights2 = []
    for i in range(2, n):
            weights2.append(weights[i] * (1 - (data[i] / s6) ** 2) ** 2)
    return [weights2[0], weights2[0]] + weights2

Example 17

Project: tvb-library
Source File: arrays.py
View license
    def _find_summary_info(self):
        "Summarize array contents."
        summary = {"Array type": self.__class__.__name__,
                   "Shape": self.shape,
                   "Maximum": self.value.max(),
                   "Minimum": self.value.min(),
                   "Mean": self.value.mean(),
                   "Median": numpy.median(self.value)}
        return summary

Example 18

View license
  def benchmarkCreateLocalCluster(self):
    deltas = []
    iters = 5
    for _ in range(iters):
      start_time = time.time()
      create_local_cluster(num_workers=1, num_ps=10)
      end_time = time.time()
      deltas.append(end_time - start_time)

    median_deltas = np.median(deltas)
    print(
        "\n\nbenchmark_create_local_cluster_1_worker_10_ps.  "
        "iterations: %d, median wall time: %g\n\n" % (iters, median_deltas))
    self.report_benchmark(
        iters=iters,
        wall_time=median_deltas,
        name="benchmark_create_local_cluster_1_worker_10_ps")

Example 19

View license
  def benchmarkCreateLocalCluster(self):
    deltas = []
    iters = 5
    for _ in range(iters):
      start_time = time.time()
      create_local_cluster(num_workers=1, num_ps=10)
      end_time = time.time()
      deltas.append(end_time - start_time)

    median_deltas = np.median(deltas)
    print(
        "\n\nbenchmark_create_local_cluster_1_worker_10_ps.  "
        "iterations: %d, median wall time: %g\n\n" % (iters, median_deltas))
    self.report_benchmark(
        iters=iters,
        wall_time=median_deltas,
        name="benchmark_create_local_cluster_1_worker_10_ps")

Example 20

Project: Sushi
Source File: sushi.py
View license
def running_median(values, window_size):
    if window_size % 2 != 1:
        raise SushiError('Median window size should be odd')
    half_window = window_size // 2
    medians = []
    items_count = len(values)
    for idx in xrange(items_count):
        radius = min(half_window, idx, items_count-idx-1)
        med = np.median(values[idx-radius:idx+radius+1])
        medians.append(med)
    return medians

Example 21

Project: rf
Source File: imaging.py
View license
def _get_geoaxes(crs=None, latlons=None):
    """Return cartopy geoaxis"""
    if crs is None:
        from cartopy.crs import AzimuthalEquidistant
        latlon0 = np.median(latlons, axis=0)
        crs = AzimuthalEquidistant(*latlon0[::-1])
    return plt.axes(projection=crs)

Example 22

Project: rf
Source File: profile.py
View license
def _find_box(latlon, boxes, crs=None):
    """Return the box which encloses the coordinates."""
    import cartopy.crs as ccrs
    from shapely.geometry import Point
    if crs is None:
        latlons = [boxes[len(boxes)//2]['latlon']]
        latlon0 = np.median(latlons, axis=0)
        crs = ccrs.AzimuthalEquidistant(*latlon0[::-1])
    pc = ccrs.PlateCarree()
    p = crs.project_geometry(Point(*latlon[::-1]), pc)
    for box in boxes:
        poly = crs.project_geometry(box['poly'], pc)
        if p.within(poly):
            return box

Example 23

Project: msaf
Source File: main.py
View license
def get_beats(y, sr, hop_length):

    odf = librosa.onset.onset_strength(y=y, sr=sr, hop_length=hop_length, aggregate=np.median, n_mels=128)
    bpm, beats = librosa.beat.beat_track(onset_envelope=odf, sr=sr, hop_length=hop_length, start_bpm=120)

    if bpm < MIN_TEMPO:
        bpm, beats = librosa.beat.beat_track(onset_envelope=odf, sr=sr, hop_length=hop_length, bpm=2*bpm)

    return bpm, beats

Example 24

Project: qiime
Source File: make_distance_boxplots.py
View license
def _sort_distributions(plot_data, plot_labels, plot_colors, sort_type):
    """Sort plot data, labels, and colors.

    Can sort in order of increasing median or alphabetically (based on plot
    labels).

    """
    if sort_type == 'median':
        sort_key = lambda item: np.median(item[0])
    elif sort_type == 'alphabetical':
        sort_key = lambda item: item[1]
    else:
        raise ValueError("Invalid sort type '%s'." % sort_type)

    # Taken from http://stackoverflow.com/a/9764364
    return zip(*sorted(zip(plot_data, plot_labels, plot_colors), key=sort_key))

Example 25

Project: orange
Source File: quickmenu.py
View license
    def sizeHint(self):
        """
        Size hint is the maximum width and median height of the widgets
        contained in the stack.

        """
        default_size = QSize(200, 400)
        widget_hints = [default_size]
        for i in range(self.count()):
            hint = self.widget(i).sizeHint()
            widget_hints.append(hint)

        width = max([s.width() for s in widget_hints])
        # Take the median for the height
        height = numpy.median([s.height() for s in widget_hints])

        return QSize(width, int(height))

Example 26

Project: tsfresh
Source File: feature_calculators.py
View license
@set_property("fctype", "aggregate")
def median(x):
    """
    Returns the median of x

    :param x: the time series to calculate the feature of
    :type x: pandas.Series
    :return: the value of this feature
    :return type: float
    """
    return np.median(x)

Example 27

Project: pipeline
Source File: dynamicEnhancer.py
View license
def getMedianSignal(enhancerFile,name,dataFile):

    '''
    returns the median enhancer signal of a file
    '''
    dataDict = pipeline_dfci.loadDataTable(dataFile)
    enhancerTable = utils.parseTable(enhancerFile,'\t')

    backgroundName = dataDict[name]['background']
    if dataDict.has_key(backgroundName):
        enhancerVector = [float(line[6]) - float(line[7]) for line in enhancerTable[6:]]
    else:
        enhancerVector = [float(line[6]) for line in enhancerTable[6:]]
            

    median= numpy.median(enhancerVector)

    return median

Example 28

Project: combined-pvalues
Source File: stepsize.py
View license
def stepsize(bed_files, col):

    D1 = []
    for bed_file in bed_files:
        for _, chromlist in groupby(bediter(bed_file, col), itemgetter('chrom')):
            L = list(chromlist)

            last_start = 0
            for i, ibed in enumerate(L):
                assert ibed['start'] >= last_start
                # look around ibed. nearest could be up or down-stream
                if i + 2 == len(L): break
                D1.append(L[i + 1]['start'] - ibed['start'])
        # round up to the nearest 10
    return int(round(np.median(D1) + 5, -1))

Example 29

Project: ldsc
Source File: munge_sumstats.py
View license
def check_median(x, expected_median, tolerance, name):
    '''Check that median(x) is within tolerance of expected_median.'''
    m = np.median(x)
    if np.abs(m - expected_median) > tolerance:
        msg = 'WARNING: median value of {F} is {V} (should be close to {M}). This column may be mislabeled.'
        raise ValueError(msg.format(F=name, M=expected_median, V=round(m, 2)))
    else:
        msg = 'Median value of {F} was {C}, which seems sensible.'.format(
            C=m, F=name)

    return msg

Example 30

View license
def test_scatter_res_raw():
    """Test feature that measures scatter of Lomb-Scargle residuals."""
    times, values, errors = irregular_random()
    lomb_model = lomb_scargle.lomb_scargle_model(times, values, errors)
    residuals = values - lomb_model['freq_fits'][0]['model']
    resid_mad = np.median(np.abs(residuals - np.median(residuals)))
    value_mad = np.median(np.abs(values - np.median(values)))
    f = generate_features(times, values, errors, ['scatter_res_raw'])
    npt.assert_allclose(f['scatter_res_raw'], resid_mad / value_mad, atol=3e-2)

Example 31

Project: UMI-tools
Source File: dedup.py
View license
    def _get_best_percentile(self, cluster, counts):
        ''' return all UMIs with counts >1% of the
        median counts in the cluster '''

        if len(cluster) == 1:
            return list(cluster)
        else:
            threshold = np.median(counts.values())/100
            return [read for read in cluster if counts[read] > threshold]

Example 32

Project: bcbio-nextgen
Source File: __init__.py
View license
def estimate_fragment_size(bam_file, nreads=5000):
    """
    estimate median fragment size of a SAM/BAM file
    """
    with open_samfile(bam_file) as bam_handle:
        reads = tz.itertoolz.take(nreads, bam_handle)
        # it would be good to skip spliced paired reads.
        lengths = [x.template_length for x in reads if x.template_length > 0]
    if not lengths:
        return 0
    return int(numpy.median(lengths))

Example 33

Project: bcbio-nextgen
Source File: coverage.py
View license
def _average_genome_coverage(data, bam_file):
    total = sum([c.size for c in ref.file_contigs(dd.get_ref_file(data), data["config"])])
    read_counts = sambamba.number_of_mapped_reads(data, bam_file, keep_dups=False)
    with pysam.Samfile(bam_file, "rb") as pysam_bam:
        read_size = np.median(list(itertools.islice((a.query_length for a in pysam_bam.fetch()), 1e5)))
    avg_cov = float(read_counts * read_size) / total
    return avg_cov

Example 34

Project: bcbio-nextgen-vm
Source File: devel.py
View license
def _calculate_common_memory(kvs):
    """Get the median memory specification, in megabytes.
    """
    mems = []
    for key, val in kvs:
        cur_val, _ = _get_cur_mem(key, val)
        mems.append(cur_val)
    return numpy.median(mems)

Example 35

Project: urlparse4
Source File: urls.py
View license
def benchmark(name, func, debug=False):
    times = []
    for n in range(0, REPEATS):
        for i, url in enumerate(URLS):
            u = url.strip()
            if debug:
                print u
            t = clock()
            func(u)
            times.append(clock() - t)

    row = [name, sum(times), mean(times), median(times), percentile(times, 90)]
    print row
    data.append(row)

Example 36

Project: tract_querier
Source File: operations.py
View license
@tract_math_operation('<scalar>: calculates median of a scalar quantity for each tract')
def scalar_tract_median(optional_flags, tractography, scalar):
    try:
        tracts = tractography.original_tracts_data()[scalar]
        result = OrderedDict((
            ('tract file', []),
            ('median %s' % scalar, []),
        ))
        for i, t in enumerate(tracts):
            result['tract file'].append('Tract %04d' % i)
            result['median %s' % scalar].append(float(numpy.median(t)))

        return result
    except KeyError:
        raise ValueError("Tractography does not contain this scalar data")

Example 37

Project: george
Source File: models.py
View license
    def __init__(self, *args, **kwargs):
        kwargs["median"] = False
        super(GPModel, self).__init__(*args, **kwargs)

        # Normalize the fluxes.
        mu = np.median(self.f)
        self.f /= mu
        self.fe /= mu

        # Set up the GP model.
        self.kernel = 1e-6 * kernels.Matern32Kernel(3.0)
        self.gp = george.GP(self.kernel, mean=_mean_function(self.system))
        self.gp.compute(self.t, self.fe)

Example 38

Project: mmfeat
Source File: mmspace.py
View license
    def setMethodType(self, methodType):
        if methodType in ['dfmm', 'dwmm']:
            if not hasattr(self.visSpace, 'dispersions'):
                raise ValueError('Dispersion-dependent method selected but no dispersions loaded.')
            if methodType == 'dfmm':
                self.filterThreshold = np.median(self.visSpace.dispersions.values())
        self.methodType = methodType

Example 39

Project: pysparnn
Source File: pysparnn_utils.py
View license
def knn_benchmark(search_index, docs, answers, n_trials=1000, docs_per_query=50):
    # Bootstrap-ish code to measure the time and accuracy
    times = []
    recalls = []
    for i in range(n_trials):
        query_index = random.sample(range(len(docs)), docs_per_query)
        time, recall = time_it(search_index, docs, query_index, answers)
        time = time / docs_per_query
        times.append(time)
        recalls.append(recall)
    return np.median(times), np.median(recalls)    

Example 40

Project: TADbit
Source File: tad_cmo.py
View license
def matrix2binnary_contacts(tad1, tad2):
    cutoff = median(tad1)
    contacts1 = []
    contacts2 = []
    for i in xrange(len(tad1)):
        for j in xrange(i, len(tad1)):
            if tad1[i][j] > cutoff:
                contacts1.append((i, j))
    cutoff = median(tad2)
    for i in xrange(len(tad2)):
        for j in xrange(i, len(tad2)):
            if tad2[i][j] > cutoff:
                contacts2.append((i, j))
    return contacts1, contacts2

Example 41

Project: oq-engine
Source File: scenario_test.py
View license
    def medians(self, case):
        [gmfa] = self.execute(case.__file__, 'job.ini').values()
        median = {imt: [] for imt in self.calc.oqparam.imtls}
        for imti, imt in enumerate(self.calc.oqparam.imtls):
            gmfa_by_imt = get_array(gmfa, imti=imti)
            for sid in self.calc.sitecol.sids:
                gmvs = get_array(gmfa_by_imt, sid=sid)['gmv']
                median[imt].append(numpy.median(gmvs))
        return median

Example 42

Project: oq-risklib
Source File: scenario_test.py
View license
    def medians(self, case):
        [gmfa] = self.execute(case.__file__, 'job.ini').values()
        median = {imt: [] for imt in self.calc.oqparam.imtls}
        for imti, imt in enumerate(self.calc.oqparam.imtls):
            gmfa_by_imt = get_array(gmfa, imti=imti)
            for sid in self.calc.sitecol.sids:
                gmvs = get_array(gmfa_by_imt, sid=sid)['gmv']
                median[imt].append(numpy.median(gmvs))
        return median

Example 43

Project: mtpy
Source File: pek2dforward.py
View license
    def force_isotropy(self):
        """
        force isotropy at depths shallower than anisotropy_min_depth. Clears
        up some bins for resistivity - these are limited
        
        """

        rnew = 1.*self.resistivity
        rmedian = 1.*np.median(self.resistivity[:,:,:2],axis=2)
        
        for j in range(len(rnew)):
            for k in range(len(rnew[j])):
                for i in range(2):
                    if self.blockcentres_z[j] < self.anisotropy_min_depth:
                        rnew[j,k,i] = rmedian[j,k]

        self.resistivity = rnew

Example 44

Project: PrettyPandas
Source File: styler.py
View license
    def median(self, title="Median", **kwargs):
        """Add a median summary to this table.

        :param title: Title to be displayed.
        :param kwargs: Keyword arguments passed to ``numpy.median``.
        """
        return self.summary(np.median, title, **kwargs)

Example 45

Project: vcsi
Source File: vcsi.py
View license
    def avg9x(self, matrix, percentage=0.05):
        """Computes the median of the top n% highest values.
        By default, takes the top 5%
        """
        xs = matrix.flatten()
        srt = sorted(xs, reverse=True)
        length = int(math.floor(percentage * len(srt)))

        matrix_subset = srt[:length]
        return numpy.median(matrix_subset)

Example 46

View license
def auto_canny(image, sigma=0.25):
    # compute the median of the single channel pixel intensities
    v = np.median(image)

    # apply automatic Canny edge detection using the computed median
    lower = int(max(0, (1.0 - sigma) * v))
    upper = int(min(255, (1.0 + sigma) * v))
    edged = cv2.Canny(image, lower, upper)
    # return the edged image
    return edged

Example 47

Project: lhcb_trigger_ml
Source File: losses.py
View license
    def update_tree_leaf(self, leaf, indices_in_leaf,
                         X, y, y_pred, sample_weight, update_mask, residual):
        if self.use_median:
            residual = residual[indices_in_leaf]
            return numpy.median(residual)
        else:
            return numpy.clip(y_pred[indices_in_leaf[0]], -10, 10)

Example 48

Project: lhcb_trigger_ml
Source File: losses.py
View license
    def update_tree_leaf(self, leaf, indices_in_leaf,
                         X, y, y_pred, sample_weight, update_mask, residual):
        if self.use_median:
            residual = residual[indices_in_leaf]
            return numpy.median(residual)
        else:
            return numpy.clip(y_pred[indices_in_leaf[0]], -10, 10)

Example 49

Project: nmrglue
Source File: proc_bl.py
View license
def _find_noise_sd(sd_set, ratio):
    '''Calculate the median m1 from SDset. exclude the elements greater
    than 2m1from SDset and recalculate the median m2. Repeat until
    m2/m1 converge(sd_set)'''
    m1 = np.median(sd_set)
    S = sd_set <= 2.0 * m1
    tmp = S * sd_set
    sd_set = tmp[tmp != 0]
    m2 = np.median(sd_set)
    while m2/m1 < ratio:
        m1 = np.median(sd_set)
        S = sd_set <= 2.0 * m1
        tmp = S * sd_set
        sd_set = tmp[tmp != 0]
    return m2

Example 50

Project: nmrglue
Source File: proc_bl.py
View license
def _find_noise_sd(sd_set, ratio):
    '''Calculate the median m1 from SDset. exclude the elements greater
    than 2m1from SDset and recalculate the median m2. Repeat until
    m2/m1 converge(sd_set)'''
    m1 = np.median(sd_set)
    S = sd_set <= 2.0 * m1
    tmp = S * sd_set
    sd_set = tmp[tmp != 0]
    m2 = np.median(sd_set)
    while m2/m1 < ratio:
        m1 = np.median(sd_set)
        S = sd_set <= 2.0 * m1
        tmp = S * sd_set
        sd_set = tmp[tmp != 0]
    return m2