Examples with Min org.apache.commons.math3.stat.descriptive.rank.Min used on opensource projects

Example 6 with Min

use of org.apache.commons.math3.stat.descriptive.rank.Min in project GDSC-SMLM by aherbert.

the class SCMOSLikelihoodWrapperTest method canComputePValue.

private void canComputePValue(BaseNonLinearFunction nlf) {
    System.out.println(nlf.name);
    int n = maxx * maxx;
    double[] a = new double[] { 1 };
    // Simulate sCMOS camera
    nlf.initialise(a);
    RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
    double[] k = Utils.newArray(n, 0, 1.0);
    for (int i = 0; i < n; i++) {
        double u = nlf.eval(i);
        if (u > 0)
            u = rdg.nextPoisson(u);
        k[i] = u * g[i] + rdg.nextGaussian(o[i], sd[i]);
    }
    SCMOSLikelihoodWrapper f = new SCMOSLikelihoodWrapper(nlf, a, k, n, var, g, o);
    double oll = f.computeObservedLikelihood();
    double oll2 = 0;
    double[] op = new double[n];
    for (int j = 0; j < n; j++) {
        op[j] = SCMOSLikelihoodWrapper.likelihood((k[j] - o[j]) / g[j], var[j], g[j], o[j], k[j]);
        oll2 -= Math.log(op[j]);
    }
    System.out.printf("oll=%f, oll2=%f\n", oll, oll2);
    Assert.assertEquals("Observed Log-likelihood", oll2, oll, oll2 * 1e-10);
    double min = Double.POSITIVE_INFINITY;
    double mina = 0;
    for (int i = 5; i <= 15; i++) {
        a[0] = (double) i / 10;
        double ll = f.likelihood(a);
        double llr = f.computeLogLikelihoodRatio(ll);
        BigDecimal product = new BigDecimal(1);
        double ll2 = 0;
        for (int j = 0; j < n; j++) {
            double p1 = SCMOSLikelihoodWrapper.likelihood(nlf.eval(j), var[j], g[j], o[j], k[j]);
            ll2 -= Math.log(p1);
            double ratio = p1 / op[j];
            product = product.multiply(new BigDecimal(ratio));
        }
        double llr2 = -2 * Math.log(product.doubleValue());
        double q = f.computeQValue(ll);
        System.out.printf("a=%f, ll=%f, ll2=%f, llr=%f, llr2=%f, product=%s, p=%f\n", a[0], ll, ll2, llr, llr2, product.round(new MathContext(4)).toString(), q);
        if (min > ll) {
            min = ll;
            mina = a[0];
        }
        // too small to store in a double.
        if (product.doubleValue() > 0)
            Assert.assertEquals("Log-likelihood", llr, llr2, llr * 1e-10);
    }
    Assert.assertEquals("min", 1, mina, 0);
}

Example 7 with Min

use of org.apache.commons.math3.stat.descriptive.rank.Min in project lucene-solr by apache.

the class DescribeEvaluator method evaluate.

public Tuple evaluate(Tuple tuple) throws IOException {
    if (subEvaluators.size() != 1) {
        throw new IOException("describe expects 1 column as a parameters");
    }
    StreamEvaluator colEval = subEvaluators.get(0);
    List<Number> numbers = (List<Number>) colEval.evaluate(tuple);
    DescriptiveStatistics descriptiveStatistics = new DescriptiveStatistics();
    for (Number n : numbers) {
        descriptiveStatistics.addValue(n.doubleValue());
    }
    Map map = new HashMap();
    map.put("max", descriptiveStatistics.getMax());
    map.put("mean", descriptiveStatistics.getMean());
    map.put("min", descriptiveStatistics.getMin());
    map.put("stdev", descriptiveStatistics.getStandardDeviation());
    map.put("sum", descriptiveStatistics.getSum());
    map.put("N", descriptiveStatistics.getN());
    map.put("var", descriptiveStatistics.getVariance());
    map.put("kurtosis", descriptiveStatistics.getKurtosis());
    map.put("skewness", descriptiveStatistics.getSkewness());
    map.put("popVar", descriptiveStatistics.getPopulationVariance());
    map.put("geometricMean", descriptiveStatistics.getGeometricMean());
    map.put("sumsq", descriptiveStatistics.getSumsq());
    return new Tuple(map);
}

Example 8 with Min

use of org.apache.commons.math3.stat.descriptive.rank.Min in project GDSC-SMLM by aherbert.

the class PCPALMClusters method fitBinomial.

/**
	 * Fit a zero-truncated Binomial to the cumulative histogram
	 * 
	 * @param histogramData
	 * @return
	 */
private double[] fitBinomial(HistogramData histogramData) {
    // Get the mean and sum of the input histogram
    double mean;
    double sum = 0;
    count = 0;
    for (int i = 0; i < histogramData.histogram[1].length; i++) {
        count += histogramData.histogram[1][i];
        sum += histogramData.histogram[1][i] * i;
    }
    mean = sum / count;
    String name = "Zero-truncated Binomial distribution";
    Utils.log("Mean cluster size = %s", Utils.rounded(mean));
    Utils.log("Fitting cumulative " + name);
    // Convert to a normalised double array for the binomial fitter
    double[] histogram = new double[histogramData.histogram[1].length];
    for (int i = 0; i < histogramData.histogram[1].length; i++) histogram[i] = histogramData.histogram[1][i] / count;
    // Plot the cumulative histogram
    String title = TITLE + " Cumulative Distribution";
    Plot2 plot = null;
    if (showCumulativeHistogram) {
        // Create a cumulative histogram for fitting
        double[] cumulativeHistogram = new double[histogram.length];
        sum = 0;
        for (int i = 0; i < histogram.length; i++) {
            sum += histogram[i];
            cumulativeHistogram[i] = sum;
        }
        double[] values = Utils.newArray(histogram.length, 0.0, 1.0);
        plot = new Plot2(title, "N", "Cumulative Probability", values, cumulativeHistogram);
        plot.setLimits(0, histogram.length - 1, 0, 1.05);
        plot.addPoints(values, cumulativeHistogram, Plot2.CIRCLE);
        Utils.display(title, plot);
    }
    // Do fitting for different N
    double bestSS = Double.POSITIVE_INFINITY;
    double[] parameters = null;
    int worse = 0;
    int N = histogram.length - 1;
    int min = minN;
    final boolean customRange = (minN > 1) || (maxN > 0);
    if (min > N)
        min = N;
    if (maxN > 0 && N > maxN)
        N = maxN;
    Utils.log("Fitting N from %d to %d%s", min, N, (customRange) ? " (custom-range)" : "");
    // Since varying the N should be done in integer steps do this
    // for n=1,2,3,... until the SS peaks then falls off (is worse then the best 
    // score several times in succession)
    BinomialFitter bf = new BinomialFitter(new IJLogger());
    bf.setMaximumLikelihood(maximumLikelihood);
    for (int n = min; n <= N; n++) {
        PointValuePair solution = bf.fitBinomial(histogram, mean, n, true);
        if (solution == null)
            continue;
        double p = solution.getPointRef()[0];
        Utils.log("Fitted %s : N=%d, p=%s. SS=%g", name, n, Utils.rounded(p), solution.getValue());
        if (bestSS > solution.getValue()) {
            bestSS = solution.getValue();
            parameters = new double[] { n, p };
            worse = 0;
        } else if (bestSS < Double.POSITIVE_INFINITY) {
            if (++worse >= 3)
                break;
        }
        if (showCumulativeHistogram)
            addToPlot(n, p, title, plot, new Color((float) n / N, 0, 1f - (float) n / N));
    }
    // Add best it in magenta
    if (showCumulativeHistogram && parameters != null)
        addToPlot((int) parameters[0], parameters[1], title, plot, Color.magenta);
    return parameters;
}

Example 9 with Min

use of org.apache.commons.math3.stat.descriptive.rank.Min in project uPortal by Jasig.

the class JpaStatisticalSummaryTest method testSummaryStatisticsJson.

// @Test
// public void testJpaStatisticalSummary() {
// final long id = this.executeInTransaction(new Callable<Long>() {
// @Override
// public Long call() throws Exception {
// final JpaStatisticalSummary jpaStatisticalSummary = new
// JpaStatisticalSummary();
// 
// final Random r = new Random(0);
// for (int i = 0; i < 10; i++) {
// final int nextInt = r.nextInt(100000000);
// jpaStatisticalSummary.addValue(nextInt);
// }
// 
// getEntityManager().persist(jpaStatisticalSummary);
// 
// System.out.println(jpaStatisticalSummary);
// 
// return jpaStatisticalSummary.getStatSummaryId();
// }
// });
// 
// System.out.println(id);
// 
// this.executeInTransaction(new CallableWithoutResult() {
// @Override
// protected void callWithoutResult() {
// final JpaStatisticalSummary jpaStatisticalSummary =
// getEntityManager().find(JpaStatisticalSummary.class, id);
// 
// System.out.println(jpaStatisticalSummary);
// }
// });
// }
@Ignore
@Test
public void testSummaryStatisticsJson() throws Exception {
    final SecondMoment secondMoment = new SecondMoment();
    final Sum sum = new Sum();
    final SumOfSquares sumsq = new SumOfSquares();
    final Min min = new Min();
    final Max max = new Max();
    final SumOfLogs sumLog = new SumOfLogs();
    final Random r = new Random(0);
    for (int i = 0; i < 10; i++) {
        final int nextInt = r.nextInt(100000000);
        secondMoment.increment(nextInt);
        sum.increment(nextInt);
        sumsq.increment(nextInt);
        min.increment(nextInt);
        max.increment(nextInt);
        sumLog.increment(nextInt);
    }
    testStorelessUnivariateStatistic(secondMoment, 7.513432791665536E15);
    testStorelessUnivariateStatistic(sum, 6.01312177E8);
    testStorelessUnivariateStatistic(sumsq, 4.3671066212513456E16);
    testStorelessUnivariateStatistic(min, 2116447.0);
    testStorelessUnivariateStatistic(max, 8.5505948E7);
    testStorelessUnivariateStatistic(sumLog, 175.91713800250577);
}

Example 10 with Min

use of org.apache.commons.math3.stat.descriptive.rank.Min in project gatk-protected by broadinstitute.

the class PosteriorSummaryUtils method calculatePosteriorMode.

/**
     * Given a list of posterior samples, returns an estimate of the posterior mode (using
     * mllib kernel density estimation in {@link KernelDensity} and {@link BrentOptimizer}).
     * Note that estimate may be poor if number of samples is small (resulting in poor kernel density estimation),
     * or if posterior is not unimodal (or is sufficiently pathological otherwise). If the samples contain
     * {@link Double#NaN}, {@link Double#NaN} will be returned.
     * @param samples   posterior samples, cannot be {@code null} and number of samples must be greater than 0
     * @param ctx       {@link JavaSparkContext} used by {@link KernelDensity} for mllib kernel density estimation
     */
public static double calculatePosteriorMode(final List<Double> samples, final JavaSparkContext ctx) {
    Utils.nonNull(samples);
    Utils.validateArg(samples.size() > 0, "Number of samples must be greater than zero.");
    //calculate sample min, max, mean, and standard deviation
    final double sampleMin = Collections.min(samples);
    final double sampleMax = Collections.max(samples);
    final double sampleMean = new Mean().evaluate(Doubles.toArray(samples));
    final double sampleStandardDeviation = new StandardDeviation().evaluate(Doubles.toArray(samples));
    //if samples are all the same or contain NaN, can simply return mean
    if (sampleStandardDeviation == 0. || Double.isNaN(sampleMean)) {
        return sampleMean;
    }
    //use Silverman's rule to set bandwidth for kernel density estimation from sample standard deviation
    //see https://en.wikipedia.org/wiki/Kernel_density_estimation#Practical_estimation_of_the_bandwidth
    final double bandwidth = SILVERMANS_RULE_CONSTANT * sampleStandardDeviation * Math.pow(samples.size(), SILVERMANS_RULE_EXPONENT);
    //use kernel density estimation to approximate posterior from samples
    final KernelDensity pdf = new KernelDensity().setSample(ctx.parallelize(samples, 1)).setBandwidth(bandwidth);
    //use Brent optimization to find mode (i.e., maximum) of kernel-density-estimated posterior
    final BrentOptimizer optimizer = new BrentOptimizer(RELATIVE_TOLERANCE, RELATIVE_TOLERANCE * (sampleMax - sampleMin));
    final UnivariateObjectiveFunction objective = new UnivariateObjectiveFunction(f -> pdf.estimate(new double[] { f })[0]);
    //search for mode within sample range, start near sample mean
    final SearchInterval searchInterval = new SearchInterval(sampleMin, sampleMax, sampleMean);
    return optimizer.optimize(objective, GoalType.MAXIMIZE, searchInterval, BRENT_MAX_EVAL).getPoint();
}