Calculates Standard Deviation for Xbar

@param sampleDataCollection
@param nh
@param populationNh
@return
170282
public static double xbarStd(TransposeDataCollection sampleDataCollection, AssociativeArray nh, AssociativeArray populationNh) {
        return Math.sqrt(xbarVariance(sampleDataCollection, nh, populationNh));
    }
1163728
public static int sampleFromDistribution(double[] d, Random random) {
    // sample from the uniform [0,1]
    double r = random.nextDouble();
    // now compare its value to cumulative values to find what interval it falls in
    double total = 0;
    for (int i = 0; i < d.length - 1; i++) {
      if (Double.isNaN(d[i])) {
        throw new RuntimeException("Can't sample from NaN");
      }
      total += d[i];
      if (r < total) {
        return i;
      }
    }
    return d.length - 1; // in case the "double-math" didn't total to exactly 1.0
  }
110418
def sample_id(self, lon):
        ''' 

        '''
        if self.grid == 'WAC':
            sample = np.rint(float(self.SAMPLE_PROJECTION_OFFSET) + 1.0 +
                             (lon * np.pi / 180.0 - float(self.CENTER_LONGITUDE)) *
                             self.A_AXIS_RADIUS *
                             np.cos(self.CENTER_LATITUDE * np.pi / 180.0)
                             / (self.MAP_SCALE * 1e-3))
        else:
            sample = np.rint(float(self.SAMPLE_PROJECTION_OFFSET) + float(self.MAP_RESOLUTION)
                             * (lon - float(self.CENTER_LONGITUDE))) + 1
        return self._control_sample(sample)
542662
public static Statistics statisticsOf(CollectionNumber data) {
        IteratorNumber iterator = data.iterator();
        if (!iterator.hasNext()) {
            return null;
        }
        int count = 0;
        double min = iterator.nextDouble();
        while (Double.isNaN(min)) {
            if (!iterator.hasNext()) {
                return null;
            } else {
                min = iterator.nextDouble();
            }
        }
        double max = min;
        double total = min;
        double totalSquare = min*min;
        count++;
        
        while (iterator.hasNext()) {
            double value = iterator.nextDouble();
            if (!Double.isNaN(value)) {
                if (value > max)
                    max = value;
                if (value < min)
                    min = value;
                total += value;
                totalSquare += value*value;
                count++;
            }
        }
        
        double average = total/count;
        double stdDev = Math.sqrt(totalSquare / count - average * average);
        
        return new StatisticsImpl(count, min, max, average, stdDev);
    }
908449
@Override
	protected double[] probabilities(
		final Seq<Phenotype<G, C>> population,
		final int count
	) {
		assert population != null : "Population must not be null. ";
		assert !population.isEmpty() : "Population is empty.";
		assert count > 0 : "Population to select must be greater than zero. ";

		final double N = population.size();
		final double[] probabilities = new double[population.size()];

		if (N == 1) {
			probabilities[0] = 1;
		} else {
			for (int i = probabilities.length; --i >= 0; ) {
				probabilities[probabilities.length - i - 1] =
					(_nminus + (_nplus - _nminus)*i/(N - 1))/N;
			}
		}

		return probabilities;
	}
423224
public static double pbarStd(double pbar, int sampleN, int populationN) {
        return Math.sqrt(pbarVariance(pbar, sampleN, populationN));
    }
379164
public double randNormal(double mean, double standardDeviation) {
    double x1, x2, w, y1;

    do {
      x1 = 2.0 * randomGenerator.nextDouble() - 1.0;
      x2 = 2.0 * randomGenerator.nextDouble() - 1.0;
      w = x1 * x1 + x2 * x2;
    } while (w >= 1.0);

    w = Math.sqrt((-2.0 * Math.log(w)) / w);
    y1 = x1 * w;
    y1 = y1 * standardDeviation + mean;
    return y1;
  }
652376
public static double harmonicMean(FlatDataCollection flatDataCollection) {
	    int n = 0;
        double harmonicMean = 0.0;
        
        Iterator<Double> it = flatDataCollection.iteratorDouble();
        while(it.hasNext()) {
            Double v = it.next();
            if(v!=null) {
                ++n;
                harmonicMean+=1.0/v;
            }
        }
        
        harmonicMean=n/harmonicMean;
        
        return harmonicMean;
    }
201270
def standard_deviation(numbers):
    """"""
    numbers = list(numbers)
    if not numbers:
        return 0
    mean = sum(numbers) / len(numbers)
    return (sum((n - mean) ** 2 for n in numbers) /
            len(numbers)) ** .5
1194230
protected double computeNNForRealData(final KNNQuery<NumberVector> knnQuery, Relation<NumberVector> relation, final int dim) {
    double w = 0.;
    ModifiableDBIDs dataSampleIds = DBIDUtil.randomSample(relation.getDBIDs(), sampleSize, random);
    for(DBIDIter iter = dataSampleIds.iter(); iter.valid(); iter.advance()) {
      final double kdist = knnQuery.getKNNForDBID(iter, k + 1).getKNNDistance();
      w += MathUtil.powi(kdist, dim);
    }
    return w;
  }
323245
public static FlatDataCollection randomSampling(FlatDataList idList, int n, boolean withReplacement) {
        FlatDataList sampledIds = new FlatDataList();
        
        int populationN = idList.size();
        
        for(int i=0;i<n;) {
            if(withReplacement==false && populationN<=n) {
                /* if replacement is not allowed and we already sampled everything that it can stop */
                break;
            }
            
            int randomPosition = PHPMethods.mt_rand(0, populationN-1);
            
            Object pointID = idList.get(randomPosition);
            
            if(withReplacement==false && sampledIds.contains(pointID)) {
                continue;
            }

            sampledIds.add(pointID);
            ++i;
            
            
        }
    
        return sampledIds.toFlatDataCollection();
    }
657416
def stdev(data, xbar=None):
    
    var = variance(data, xbar)
    try:
        return var.sqrt()
    except AttributeError:
        return math.sqrt(var)
359565
def _calc_q_h0(n, x, h, nt, n_jobs=1, verbose=0, random_state=None):
    
    rng = check_random_state(random_state)
    par, func = parallel_loop(_calc_q_statistic, n_jobs, verbose)
    q = par(func(rng.permutation(x.T).T, h, nt) for _ in range(n))
    
    return np.array(q)
801358
@Override
    public double calculateAnomalyScore(double value) {
        double zScore = (value - mean) / Math.sqrt(variance);
        //Taking absolute value for a more human-readable anomaly score
        return Math.abs(zScore);
    }
873369
def compute_offset_nuth(dh, slope, aspect, min_count=100, remove_outliers=True, plot=True):
    
    import scipy.optimize as optimization

    if dh.count() < min_count:
        sys.exit("Not enough dh samples")
    if slope.count() < min_count:
        sys.exit("Not enough slope/aspect samples")

    #mean_dh = dh.mean()
    #mean_slope = slope.mean()
    #c_seed = (mean_dh/np.tan(np.deg2rad(mean_slope))) 
    med_dh = malib.fast_median(dh)
    med_slope = malib.fast_median(slope)
    c_seed = (med_dh/np.tan(np.deg2rad(med_slope))) 

    x0 = np.array([0.0, 0.0, c_seed])
  
    print("Computing common mask")
    common_mask = ~(malib.common_mask([dh, aspect, slope]))

    #Prepare x and y data
    xdata = aspect[common_mask].data
    ydata = (dh[common_mask]/np.tan(np.deg2rad(slope[common_mask]))).data

    print("Initial sample count:")
    print(ydata.size)

    if remove_outliers:
        print("Removing outliers")
        #print("Absolute dz filter: %0.2f" % max_dz)
        #diff = np.ma.masked_greater(diff, max_dz)
        #print(diff.count())

        #Outlier dz filter
        f = 3
        sigma, u = (ydata.std(), ydata.mean())
        #sigma, u = malib.mad(ydata, return_med=True)
        rmin = u - f*sigma
        rmax = u + f*sigma
        print("3-sigma filter: %0.2f - %0.2f" % (rmin, rmax))
        idx = (ydata >= rmin) & (ydata <= rmax)
        xdata = xdata[idx]
        ydata = ydata[idx]
        print(ydata.size)

    #Generate synthetic data to test curve_fit
    #xdata = np.arange(0,360,0.01)
    #ydata = f(xdata, 20.0, 130.0, -3.0) + 20*np.random.normal(size=len(xdata))
    
    #Limit sample size
    #n = 10000
    #idx = random.sample(range(xdata.size), n)
    #xdata = xdata[idx]
    #ydata = ydata[idx]

    #Compute robust statistics for 1-degree bins
    nbins = 360
    bin_range = (0., 360.)
    bin_width = 1.0
    bin_count, bin_edges, bin_centers = malib.bin_stats(xdata, ydata, stat='count', nbins=nbins, bin_range=bin_range)
    bin_med, bin_edges, bin_centers = malib.bin_stats(xdata, ydata, stat='median', nbins=nbins, bin_range=bin_range)
    #Needed to estimate sigma for weighted lsq
    #bin_mad, bin_edges, bin_centers = malib.bin_stats(xdata, ydata, stat=malib.mad, nbins=nbins, bin_range=bin_range)
    #Started implementing this for more generic binning, needs testing
    #bin_count, x_bin_edges, y_bin_edges = malib.get_2dhist(xdata, ydata, \
    #        xlim=bin_range, nbins=(nbins, nbins), stat='count')

    """
    #Mask bins in grid directions, can potentially contain biased stats
    #Especially true for SGM algorithm
    #badbins = [0, 90, 180, 270, 360]
    badbins = [0, 45, 90, 135, 180, 225, 270, 315, 360]
    bin_stat = np.ma.masked_where(np.around(bin_edges[:-1]) % 45 == 0, bin_stat)
    bin_edges = np.ma.masked_where(np.around(bin_edges[:-1]) % 45 == 0, bin_edges)
    """

    #Remove any bins with only a few points
    min_bin_sample_count = 9
    idx = (bin_count.filled(0) >= min_bin_sample_count) 
    bin_count = bin_count[idx].data
    bin_med = bin_med[idx].data
    #bin_mad = bin_mad[idx].data
    bin_centers = bin_centers[idx]

    fit = None
    fit_fig = None

    #Want a good distribution of bins, at least 1/4 to 1/2 of sinusoid, to ensure good fit
    #Need at least 3 valid bins to fit 3 parameters in nuth_func
    #min_bin_count = 3
    min_bin_count = 90 
    
    #Not going to help if we have a step function between two plateaus, but better than nothing
    #Calculate bin aspect spread
    bin_ptp = np.cos(np.radians(bin_centers)).ptp()
    min_bin_ptp = 1.0 

    #Should iterate here, if not enough bins, increase bin width
    if len(bin_med) >= min_bin_count and bin_ptp >= min_bin_ptp:

        print("Computing fit")
        #Unweighted fit
        fit = optimization.curve_fit(nuth_func, bin_centers, bin_med, x0)[0]

        #Weight by observed spread in each bin 
        #sigma = bin_mad
        #fit = optimization.curve_fit(nuth_func, bin_centers, bin_med, x0, sigma, absolute_sigma=True)[0]

        #Weight by bin count
        #sigma = bin_count.max()/bin_count
        #fit = optimization.curve_fit(nuth_func, bin_centers, bin_med, x0, sigma, absolute_sigma=False)[0]

        print(fit)

        if plot:
            print("Generating Nuth and Kaab plot")
            bin_idx = np.digitize(xdata, bin_edges)
            output = []
            for i in np.arange(1, len(bin_edges)):
                output.append(ydata[bin_idx==i])
            #flierprops={'marker':'.'}
            lw = 0.25
            whiskerprops={'linewidth':lw}
            capprops={'linewidth':lw}
            boxprops={'facecolor':'k', 'linewidth':0}
            medianprops={'marker':'o', 'ms':1, 'color':'r'}
            fit_fig, ax = plt.subplots(figsize=(6,6))
            #widths = (bin_width/2.0)
            widths = 2.5*(bin_count/bin_count.max())
            #widths = bin_count/np.percentile(bin_count, 50)
            #Stride
            s=3
            #This is inefficient, but we have list of arrays with different length, need to filter
            #Reduntant with earlier filter, should refactor
            bp = ax.boxplot(np.array(output)[idx][::s], positions=bin_centers[::s], widths=widths[::s], showfliers=False, \
                    patch_artist=True, boxprops=boxprops, whiskerprops=whiskerprops, capprops=capprops, \
                    medianprops=medianprops)
            bin_ticks = [0, 45, 90, 135, 180, 225, 270, 315, 360]
            ax.set_xticks(bin_ticks)
            ax.set_xticklabels(bin_ticks)
            """
            #Can pull out medians from boxplot
            #We are computing multiple times, inefficient
            bp_bin_med = []
            for medline in bp['medians']:
                bp_bin_med.append(medline.get_ydata()[0])
            """

            #Plot the fit
            f_a = nuth_func(bin_centers, fit[0], fit[1], fit[2])
            nuth_func_str = r'$y=%0.2f*cos(%0.2f-x)+%0.2f$' % tuple(fit)
            ax.plot(bin_centers, f_a, 'b', label=nuth_func_str)

            ax.set_xlabel('Aspect (deg)')
            ax.set_ylabel('dh/tan(slope) (m)')
            ax.axhline(color='gray', linewidth=0.5)

            ax.set_xlim(*bin_range)
            ylim = ax.get_ylim()
            abs_ylim = np.max(np.abs(ylim))
            #abs_ylim = np.max(np.abs([ydata.min(), ydata.max()]))
            #pad = 0.2 * abs_ylim 
            pad = 0
            ylim = (-abs_ylim - pad, abs_ylim + pad)
            minylim = (-10,10)
            if ylim[0] > minylim[0]:
                ylim = minylim
            ax.set_ylim(*ylim)
            ax.legend(prop={'size':8})

    return fit, fit_fig
76478
public final double sampleSigma(int n)
  {
    synchronized (_lock) {
      long count = _count.get();
      long lastCount = _lastStdCount;
      _lastStdCount = count;

      double sum = _sum.get();
      double lastSum = _lastStdSum;
      _lastStdSum = sum;

      double sumSquare = _sumSquare;
      _sumSquare = 0;

      if (count == lastCount)
        return 0;

      double avg = (sum - lastSum) / (count - lastCount);
      double part = (count - lastCount) * sumSquare - sum * sum;

      if (part < 0)
        part = 0;

      double std = Math.sqrt(part) / (count - lastCount);

      return _scale * (avg + n * std);
    }
  }
805773
public static double getPvalue(FlatDataCollection flatDataCollection, double median) {
        int n=0;
        AssociativeArray Di = new AssociativeArray();
        Iterator<Double> it = flatDataCollection.iteratorDouble();
        while(it.hasNext()) {
            double delta=it.next()-median;

            if(delta==0.0) {
                continue; //don't count it at all
            }

            String key="+";
            if(delta<0) {
                key="-";
            }
            Di.put(key+String.valueOf(n), Math.abs(delta));
            ++n;
        }
        if(n<=0) {
            throw new IllegalArgumentException("The provided collection can't be empty.");
        }

        //converts the values of the table with its Ranks
        Ranks.getRanksFromValues(Di);
        double W=0.0;
        for(Map.Entry<Object, Object> entry : Di.entrySet()) {
            String key = entry.getKey().toString();
            Double rank = TypeInference.toDouble(entry.getValue());
        
            if(key.charAt(0)=='+') {
                W+=rank;
            }
        }

        double pvalue= scoreToPvalue(W, n);

        return pvalue;
    }