Java Code Examples for org.apache.commons.math.util.FastMath#random()

The following examples show how to use org.apache.commons.math.util.FastMath#random() . You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
Source File: EmpiricalDistributionImpl.java    From astor with GNU General Public License v2.0 6 votes vote down vote up
/**
 * Generates a random value from this distribution.
 *
 * @return the random value.
 * @throws IllegalStateException if the distribution has not been loaded
 */
public double getNextValue() throws IllegalStateException {

    if (!loaded) {
        throw MathRuntimeException.createIllegalStateException(LocalizedFormats.DISTRIBUTION_NOT_LOADED);
    }

    // Start with a uniformly distributed random number in (0,1)
    double x = FastMath.random();

    // Use this to select the bin and generate a Gaussian within the bin
    for (int i = 0; i < binCount; i++) {
       if (x <= upperBounds[i]) {
           SummaryStatistics stats = binStats.get(i);
           if (stats.getN() > 0) {
               if (stats.getStandardDeviation() > 0) {  // more than one obs
                    return randomData.nextGaussian
                        (stats.getMean(),stats.getStandardDeviation());
               } else {
                   return stats.getMean(); // only one obs in bin
               }
           }
       }
    }
    throw new MathRuntimeException(LocalizedFormats.NO_BIN_SELECTED);
}
 
Example 2
Source File: EmpiricalDistributionImpl.java    From astor with GNU General Public License v2.0 6 votes vote down vote up
/**
 * Generates a random value from this distribution.
 *
 * @return the random value.
 * @throws IllegalStateException if the distribution has not been loaded
 */
public double getNextValue() throws IllegalStateException {

    if (!loaded) {
        throw MathRuntimeException.createIllegalStateException(LocalizedFormats.DISTRIBUTION_NOT_LOADED);
    }

    // Start with a uniformly distributed random number in (0,1)
    double x = FastMath.random();

    // Use this to select the bin and generate a Gaussian within the bin
    for (int i = 0; i < binCount; i++) {
       if (x <= upperBounds[i]) {
           SummaryStatistics stats = binStats.get(i);
           if (stats.getN() > 0) {
               if (stats.getStandardDeviation() > 0) {  // more than one obs
                    return randomData.nextGaussian
                        (stats.getMean(),stats.getStandardDeviation());
               } else {
                   return stats.getMean(); // only one obs in bin
               }
           }
       }
    }
    throw new MathRuntimeException(LocalizedFormats.NO_BIN_SELECTED);
}
 
Example 3
Source File: LoessInterpolatorTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
    double dx = 2 * FastMath.PI / xval.length;
    double x = 0;
    for(int i = 0; i < xval.length; ++i) {
        xval[i] = x;
        yval[i] = FastMath.sin(x) + (2 * FastMath.random() - 1) * ynoise;
        x += dx * (1 + (2 * FastMath.random() - 1) * xnoise);
    }
}
 
Example 4
Source File: LoessInterpolatorTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
    double dx = 2 * FastMath.PI / xval.length;
    double x = 0;
    for(int i = 0; i < xval.length; ++i) {
        xval[i] = x;
        yval[i] = FastMath.sin(x) + (2 * FastMath.random() - 1) * ynoise;
        x += dx * (1 + (2 * FastMath.random() - 1) * xnoise);
    }
}
 
Example 5
Source File: LoessInterpolatorTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
    double dx = 2 * FastMath.PI / xval.length;
    double x = 0;
    for(int i = 0; i < xval.length; ++i) {
        xval[i] = x;
        yval[i] = FastMath.sin(x) + (2 * FastMath.random() - 1) * ynoise;
        x += dx * (1 + (2 * FastMath.random() - 1) * xnoise);
    }
}
 
Example 6
Source File: NormalDistributionTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
public void testSetStandardDeviation() throws Exception {
    double sigma = 0.1d + FastMath.random();
    NormalDistribution distribution = (NormalDistribution) getDistribution();
    distribution.setStandardDeviation(sigma);
    assertEquals(sigma, distribution.getStandardDeviation(), 0);
    verifyQuantiles();
    try {
        distribution.setStandardDeviation(0);
        fail("Expecting IllegalArgumentException for sd = 0");
    } catch (IllegalArgumentException ex) {
        // Expected
    }
}
 
Example 7
Source File: LoessInterpolatorTest.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
    double dx = 2 * FastMath.PI / xval.length;
    double x = 0;
    for(int i = 0; i < xval.length; ++i) {
        xval[i] = x;
        yval[i] = FastMath.sin(x) + (2 * FastMath.random() - 1) * ynoise;
        x += dx * (1 + (2 * FastMath.random() - 1) * xnoise);
    }
}
 
Example 8
Source File: CauchyDistributionTest.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
public void testMedian() {
    CauchyDistribution distribution = (CauchyDistribution) getDistribution();
    double expected = FastMath.random();
    distribution.setMedian(expected);
    assertEquals(expected, distribution.getMedian(), 0.0);
}
 
Example 9
Source File: CauchyDistributionTest.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
public void testScale() {
    CauchyDistribution distribution = (CauchyDistribution) getDistribution();
    double expected = FastMath.random();
    distribution.setScale(expected);
    assertEquals(expected, distribution.getScale(), 0.0);
}
 
Example 10
Source File: NormalDistributionTest.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
public void testSetMean() throws Exception {
    double mu = FastMath.random();
    NormalDistribution distribution = (NormalDistribution) getDistribution();
    distribution.setMean(mu);
    verifyQuantiles();
}
 
Example 11
Source File: MullerSolver.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/**
 * Find a real root in the given interval.
 * <p>
 * solve2() differs from solve() in the way it avoids complex operations.
 * Except for the initial [min, max], solve2() does not require bracketing
 * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
 * number arises in the computation, we simply use its modulus as real
 * approximation.</p>
 * <p>
 * Because the interval may not be bracketing, bisection alternative is
 * not applicable here. However in practice our treatment usually works
 * well, especially near real zeros where the imaginary part of complex
 * approximation is often negligible.</p>
 * <p>
 * The formulas here do not use divided differences directly.</p>
 *
 * @param f the function to solve
 * @param min the lower bound for the interval
 * @param max the upper bound for the interval
 * @return the point at which the function value is zero
 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
 * or the solver detects convergence problems otherwise
 * @throws FunctionEvaluationException if an error occurs evaluating the
 * function
 * @throws IllegalArgumentException if any parameters are invalid
 */
public double solve2(final UnivariateRealFunction f,
                     final double min, final double max)
    throws MaxIterationsExceededException, FunctionEvaluationException {

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = f.value(x0);
    double x1 = max;
    double y1 = f.value(x1);
    double x2 = 0.5 * (x0 + x1);
    double y2 = f.value(x2);

    // check for zeros before verifying bracketing
    if (y0 == 0.0) { return min; }
    if (y1 == 0.0) { return max; }
    verifyBracketing(min, max, f);

    double oldx = Double.POSITIVE_INFINITY;
    for (int i = 1; i <= maximalIterationCount; ++i) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = f.value(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance) {
            setResult(x, i);
            return result;
        }
        if (FastMath.abs(y) <= functionValueAccuracy) {
            setResult(x, i);
            return result;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
    throw new MaxIterationsExceededException(maximalIterationCount);
}
 
Example 12
Source File: IntegerDistributionAbstractTest.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/**
 * Verifies that floating point arguments are correctly handled by
 * cumulativeProbablility(-,-)
 * JIRA: MATH-184
 */
public void testFloatingPointArguments() throws Exception {
    for (int i = 0; i < cumulativeTestPoints.length; i++) {
        double arg = cumulativeTestPoints[i];
        assertEquals(
                "Incorrect cumulative probability value returned for " +
                cumulativeTestPoints[i],
                cumulativeTestValues[i],
                distribution.cumulativeProbability(arg), tolerance);
        if (i < cumulativeTestPoints.length - 1) {
            double arg2 = cumulativeTestPoints[i + 1];
            assertEquals("Inconsistent probability for discrete range " +
                    "[ " + arg + "," + arg2 + " ]",
               distribution.cumulativeProbability(
                       cumulativeTestPoints[i],
                       cumulativeTestPoints[i + 1]),
               distribution.cumulativeProbability(arg, arg2), tolerance);
            arg = arg - FastMath.random();
            arg2 = arg2 + FastMath.random();
            assertEquals("Inconsistent probability for discrete range " +
                    "[ " + arg + "," + arg2 + " ]",
               distribution.cumulativeProbability(
                       cumulativeTestPoints[i],
                       cumulativeTestPoints[i + 1]),
               distribution.cumulativeProbability(arg, arg2), tolerance);
        }
    }
    int one = 1;
    int ten = 10;
    int two = 2;
    double oned = one;
    double twod = two;
    double tend = ten;
    assertEquals(distribution.cumulativeProbability(one, two),
            distribution.cumulativeProbability(oned, twod), tolerance);
    assertEquals(distribution.cumulativeProbability(one, two),
            distribution.cumulativeProbability(oned - tolerance,
                    twod + 0.9), tolerance);
    assertEquals(distribution.cumulativeProbability(two, ten),
            distribution.cumulativeProbability(twod, tend), tolerance);
    assertEquals(distribution.cumulativeProbability(two, ten),
            distribution.cumulativeProbability(twod - tolerance,
                    tend + 0.9), tolerance);
}
 
Example 13
Source File: WeibullDistributionTest.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
public void testAlpha() {
    WeibullDistribution distribution = (WeibullDistribution) getDistribution();
    double expected = FastMath.random();
    distribution.setShape(expected);
    assertEquals(expected, distribution.getShape(), 0.0);
}
 
Example 14
Source File: WeibullDistributionTest.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
public void testBeta() {
    WeibullDistribution distribution = (WeibullDistribution) getDistribution();
    double expected = FastMath.random();
    distribution.setScale(expected);
    assertEquals(expected, distribution.getScale(), 0.0);
}
 
Example 15
Source File: MullerSolver2.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/**
 * {@inheritDoc}
 */
@Override
protected double doSolve() {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}
 
Example 16
Source File: IntegerDistributionAbstractTest.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/**
 * Verifies that floating point arguments are correctly handled by
 * cumulativeProbablility(-,-)
 * JIRA: MATH-184
 */
@Test
public void testFloatingPointArguments() throws Exception {
    for (int i = 0; i < cumulativeTestPoints.length; i++) {
        double arg = cumulativeTestPoints[i];
        Assert.assertEquals(
                "Incorrect cumulative probability value returned for " +
                cumulativeTestPoints[i],
                cumulativeTestValues[i],
                distribution.cumulativeProbability(arg), tolerance);
        if (i < cumulativeTestPoints.length - 1) {
            double arg2 = cumulativeTestPoints[i + 1];
            Assert.assertEquals("Inconsistent probability for discrete range " +
                    "[ " + arg + "," + arg2 + " ]",
               distribution.cumulativeProbability(
                       cumulativeTestPoints[i],
                       cumulativeTestPoints[i + 1]),
               distribution.cumulativeProbability(arg, arg2), tolerance);
            arg = arg - FastMath.random();
            arg2 = arg2 + FastMath.random();
            Assert.assertEquals("Inconsistent probability for discrete range " +
                    "[ " + arg + "," + arg2 + " ]",
               distribution.cumulativeProbability(
                       cumulativeTestPoints[i],
                       cumulativeTestPoints[i + 1]),
               distribution.cumulativeProbability(arg, arg2), tolerance);
        }
    }
    int one = 1;
    int ten = 10;
    int two = 2;
    double oned = one;
    double twod = two;
    double tend = ten;
    Assert.assertEquals(distribution.cumulativeProbability(one, two),
            distribution.cumulativeProbability(oned, twod), tolerance);
    Assert.assertEquals(distribution.cumulativeProbability(one, two),
            distribution.cumulativeProbability(oned - tolerance,
                    twod + 0.9), tolerance);
    Assert.assertEquals(distribution.cumulativeProbability(two, ten),
            distribution.cumulativeProbability(twod, tend), tolerance);
    Assert.assertEquals(distribution.cumulativeProbability(two, ten),
            distribution.cumulativeProbability(twod - tolerance,
                    tend + 0.9), tolerance);
}
 
Example 17
Source File: IntegerDistributionAbstractTest.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/**
 * Verifies that floating point arguments are correctly handled by
 * cumulativeProbablility(-,-)
 * JIRA: MATH-184
 */
public void testFloatingPointArguments() throws Exception {
    for (int i = 0; i < cumulativeTestPoints.length; i++) {
        double arg = cumulativeTestPoints[i];
        assertEquals(
                "Incorrect cumulative probability value returned for " +
                cumulativeTestPoints[i],
                cumulativeTestValues[i],
                distribution.cumulativeProbability(arg), tolerance);
        if (i < cumulativeTestPoints.length - 1) {
            double arg2 = cumulativeTestPoints[i + 1];
            assertEquals("Inconsistent probability for discrete range " +
                    "[ " + arg + "," + arg2 + " ]",
               distribution.cumulativeProbability(
                       cumulativeTestPoints[i],
                       cumulativeTestPoints[i + 1]),
               distribution.cumulativeProbability(arg, arg2), tolerance);
            arg = arg - FastMath.random();
            arg2 = arg2 + FastMath.random();
            assertEquals("Inconsistent probability for discrete range " +
                    "[ " + arg + "," + arg2 + " ]",
               distribution.cumulativeProbability(
                       cumulativeTestPoints[i],
                       cumulativeTestPoints[i + 1]),
               distribution.cumulativeProbability(arg, arg2), tolerance);
        }
    }
    int one = 1;
    int ten = 10;
    int two = 2;
    double oned = one;
    double twod = two;
    double tend = ten;
    assertEquals(distribution.cumulativeProbability(one, two),
            distribution.cumulativeProbability(oned, twod), tolerance);
    assertEquals(distribution.cumulativeProbability(one, two),
            distribution.cumulativeProbability(oned - tolerance,
                    twod + 0.9), tolerance);
    assertEquals(distribution.cumulativeProbability(two, ten),
            distribution.cumulativeProbability(twod, tend), tolerance);
    assertEquals(distribution.cumulativeProbability(two, ten),
            distribution.cumulativeProbability(twod - tolerance,
                    tend + 0.9), tolerance);
}
 
Example 18
Source File: MullerSolver2.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/**
 * {@inheritDoc}
 */
@Override
protected double doSolve() {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}
 
Example 19
Source File: IntegerDistributionAbstractTest.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/**
 * Verifies that floating point arguments are correctly handled by
 * cumulativeProbablility(-,-)
 * JIRA: MATH-184
 */
@Test
public void testFloatingPointArguments() throws Exception {
    for (int i = 0; i < cumulativeTestPoints.length; i++) {
        double arg = cumulativeTestPoints[i];
        Assert.assertEquals(
                "Incorrect cumulative probability value returned for " +
                cumulativeTestPoints[i],
                cumulativeTestValues[i],
                distribution.cumulativeProbability(arg), tolerance);
        if (i < cumulativeTestPoints.length - 1) {
            double arg2 = cumulativeTestPoints[i + 1];
            Assert.assertEquals("Inconsistent probability for discrete range " +
                    "[ " + arg + "," + arg2 + " ]",
               distribution.cumulativeProbability(
                       cumulativeTestPoints[i],
                       cumulativeTestPoints[i + 1]),
               distribution.cumulativeProbability(arg, arg2), tolerance);
            arg = arg - FastMath.random();
            arg2 = arg2 + FastMath.random();
            Assert.assertEquals("Inconsistent probability for discrete range " +
                    "[ " + arg + "," + arg2 + " ]",
               distribution.cumulativeProbability(
                       cumulativeTestPoints[i],
                       cumulativeTestPoints[i + 1]),
               distribution.cumulativeProbability(arg, arg2), tolerance);
        }
    }
    int one = 1;
    int ten = 10;
    int two = 2;
    double oned = one;
    double twod = two;
    double tend = ten;
    Assert.assertEquals(distribution.cumulativeProbability(one, two),
            distribution.cumulativeProbability(oned, twod), tolerance);
    Assert.assertEquals(distribution.cumulativeProbability(one, two),
            distribution.cumulativeProbability(oned - tolerance,
                    twod + 0.9), tolerance);
    Assert.assertEquals(distribution.cumulativeProbability(two, ten),
            distribution.cumulativeProbability(twod, tend), tolerance);
    Assert.assertEquals(distribution.cumulativeProbability(two, ten),
            distribution.cumulativeProbability(twod - tolerance,
                    tend + 0.9), tolerance);
}
 
Example 20
Source File: MullerSolver2.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/**
 * {@inheritDoc}
 */
@Override
protected double doSolve() {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}