List of usage examples for org.apache.commons.math3.analysis.solvers BisectionSolver BisectionSolver
public BisectionSolver()
From source file:jat.core.cm.CRTBP.java
public void findLibrationPoints() { BisectionSolver bs = new BisectionSolver(); UnivariateFunction Lfunction = new L123Func(); double L1 = bs.solve(100, Lfunction, 0, 1.); LibPoints[0] = new Vector3D(L1, 0, 0); C1 = JacobiIntegral(L1, 0, 0, 0, 0, 0); double L2 = bs.solve(100, Lfunction, 1, 2.); LibPoints[1] = new Vector3D(L2, 0, 0); C2 = JacobiIntegral(L2, 0, 0, 0, 0, 0); double L3 = bs.solve(100, Lfunction, -2.0, 0); LibPoints[2] = new Vector3D(L3, 0, 0); C3 = JacobiIntegral(L3, 0, 0, 0, 0, 0); double y45 = Math.sqrt(3.) / 2; LibPoints[3] = new Vector3D(.5 - mu, y45, 0); LibPoints[4] = new Vector3D(.5 - mu, -y45, 0); System.out.println("L1: " + L1 + " C1 " + C1); System.out.println("L2: " + L2 + " C2 " + C2); System.out.println("L3: " + L3 + " C3 " + C3); System.out.println("L4= (" + LibPoints[3].getX() + "," + LibPoints[3].getY() + ")"); System.out.println("L5= (" + LibPoints[4].getX() + "," + LibPoints[4].getY() + ")"); }
From source file:jat.core.cm.CRTBP.java
public void findZeroVelocity() { JacobiFixedx JF = new JacobiFixedx(); BisectionSolver bis = new BisectionSolver(); BrentSolver brs = new BrentSolver(); JF.setC(C);/*from w ww.j a va 2 s.com*/ double x = -1.; double lasty = 0; while ((x += .01) < 1.3) { JF.setx(x); double y = 0; boolean success = false; try { success = true; y = brs.solve(100, JF, lasty - .4, lasty + .4); // y = brs.solve(100, JF, -.1, .8); // double y = bis.solve(100, JF, 0, 2); } catch (NoBracketingException e) { success = false; // System.out.println("exception at x y " + x + " " + y); } finally { if (success) { // System.out.println("x y " + x + " " + y); xzv.add(x); yzv.add(y); lasty = y; //System.out.println("last y " + lasty); } } } }
From source file:it.unibo.alchemist.modelchecker.AlchemistASMC.java
private static int computeMinimum(final double interval, final double confidence) { final UnivariateFunction f = new UnivariateFunction() { @Override//from w w w . j a va2 s .c o m public double value(final double n) { double t; if (Math.ceil(n) == FastMath.floor(n)) { t = new TDistribution((int) n).inverseCumulativeProbability(1 - confidence / 2); } else { double t1 = new TDistribution((int) FastMath.ceil(n)) .inverseCumulativeProbability((1 - confidence / 2)) * (n - Math.floor(n)); double t2 = new TDistribution((int) FastMath.floor(n)) .inverseCumulativeProbability((1 - confidence / 2)) * (Math.ceil(n) - n); t = t1 + t2; } double value = 2 * t / n; return value - interval; } }; final BisectionSolver bs = new BisectionSolver(); return (int) Math.ceil(bs.solve(Integer.MAX_VALUE, f, 1, Integer.MAX_VALUE)); }
From source file:edu.cudenver.bios.distribution.NonCentralFDistribution.java
/** * For this non-central F distribution, F, this function returns the critical value, f, * such that P(F < f).//from w ww .j av a 2 s .c om * * @param probability desired value of P(F < f) * @return critical f such that P(F < f) */ public double inverseCDF(double probability) { if (probability <= 0) return Double.NaN; if (probability >= 1) return Double.POSITIVE_INFINITY; if (method == FMethod.CDF && nonCentrality == 0) { // inverseCdf throws an illegal argument exception when the // non-centrality parameter is non-zero. So we just bisection solve // unless we're really dealing with a central F. Ah, the hazards of // pulling code off of the interwebs. return nonCentralF.inverseCdf(probability); } else { BisectionSolver solver = new BisectionSolver(); NonCentralFQuantileFunction quantFunc = new NonCentralFQuantileFunction(probability); int upperBound = getCDFUpperBound(probability); try { return solver.solve(MAX_ITERATIONS, quantFunc, 0, upperBound); } catch (Exception e) { throw new IllegalArgumentException("Failed to determine F quantile: " + e.getMessage()); } } }
From source file:edu.cudenver.bios.distribution.NonCentralFDistribution.java
/** * For a noncentral F distribution with the specified numerator and denominator * degrees of freedom, returns the noncentrality parameter such that for the given * cumulative probability and value (f), Pr(F < f) = probability * * * @param f value of the random variable * @param probability target cumulative probability * @param ndf numerator degrees of freedom * @param ddf denominator degrees of freedom * @return noncentrality such that Pr(F < f) *//*from w w w . ja v a 2 s . c o m*/ public static double noncentrality(double f, double probability, double ndf, double ddf, double startValue) throws IllegalArgumentException, MathIllegalArgumentException, TooManyEvaluationsException { BisectionSolver solver = new BisectionSolver(); double searchUpperBound = getNoncentralityUpperBound(f, probability, ndf, ddf, startValue); double noncentrality = solver.solve(MAX_ITERATIONS, new NoncentralityFinderCDF(f, ndf, ddf, probability), 0, searchUpperBound); return noncentrality; }
From source file:edu.cudenver.bios.power.OneSampleStudentsTPowerCalculator.java
/** * Calculate sample size for the one sample t test. This function uses a bisection search algorithm * to determine sample size. The actual power is also calculated. * * @see OneSampleStudentsTPower//from w w w .j a v a 2s. c o m * @param alpha type I error * @param mu0 mean under the null hypothesis * @param muA mean under the alternative hypothesis * @param sigma standard deviation * @param sampleSize total sample size * @param isTwoTail if true, power will be calculated for a a two tailed test * @return power object containing sample size results * @throws InvalidAlgorithmParameterException */ protected OneSampleStudentsTPower calculateSampleSize(double alpha, double mu0, double muA, double sigma, double power, boolean isTwoTailed) { /* validate the inputs */ if (sigma < 0) throw new IllegalArgumentException("Invalid standard error: " + sigma); if (power < 0 || power > 1) throw new IllegalArgumentException("Invalid power: " + power); if (alpha <= 0 || alpha >= 1) throw new IllegalArgumentException("Invalid alpha level: " + alpha); BisectionSolver solver = new BisectionSolver(); SampleSizeFunction sampleSizeFunc = new SampleSizeFunction(mu0, muA, sigma, alpha, power, isTwoTailed); try { int sampleSize = (int) Math .ceil(solver.solve(MAX_ITERATIONS, sampleSizeFunc, MIN_SAMPLE_SIZE, MAX_SAMPLE_SIZE)); OneSampleStudentsTPower actualPower = calculatePower(alpha, mu0, muA, sigma, sampleSize, isTwoTailed); actualPower.setNominalPower(power); return actualPower; } catch (Exception e) { throw new IllegalArgumentException("Failed to calculate sample size: " + e.getMessage()); } }
From source file:edu.cudenver.bios.power.glmm.NonCentralityDistribution.java
/** * For this non-centrality distribution, W, this function returns the critical value, w, * such that P(W < w)./*w w w.j a v a2 s.c o m*/ * * @param probability desired value of P(W < w) * @return critical w such that P(W < w) */ public double inverseCDF(double probability) { if (H1 <= 0) return 0; BisectionSolver solver = new BisectionSolver(); NonCentralityQuantileFunction quantFunc = new NonCentralityQuantileFunction(probability); try { return solver.solve(MAX_ITERATIONS, quantFunc, H0, H1); } catch (Exception e) { throw new IllegalArgumentException("Failed to determine non-centrality quantile: " + e.getMessage()); } }
From source file:edu.cudenver.bios.power.GLMMPowerCalculator.java
/** * Find the sample size which achieves the desired power(s) * specified in the input parameters. Uses a bisection search. * * @param params GLMM input parameters//from ww w .j a v a 2 s . c o m * @return sample size */ private GLMMPower getSampleSizeValue(GLMMPowerParameters params, Test test, PowerMethod method, double alpha, double sigmaScale, double betaScale, double targetPower, double quantile) { if (method == PowerMethod.UNCONDITIONAL_POWER) { GLMMPower power = new GLMMPower(test, alpha, targetPower, -1, -1, betaScale, sigmaScale, method, quantile, null); power.setErrorMessage( "We have temporarily disabled sample size calculations using the unconditional power method " + "while we work on computational efficiency. " + "There are two alternatives." + "<ol>" + "<li>" + "You may perform a power calculation for a given sample size, " + "and iterate until you find a sample size and unconditional power that work for your design." + "</li>" + "<li>" + "You may calculate sample size using quantile power " + "calculated at the median power (0.50th quantile) instead of unconditional power. " + "As noted in Glueck and Muller (2003) " + "(see <a href=\"http://samplesizeshop.org/education/related-publications/\">Related Publications</a>), " + "quantile power is a very good approximation for unconditional power." + "</li>" + "</ol>" + "Thank you for your patience while we repair this functionality."); power.setErrorCode(PowerErrorEnum.POWER_METHOD_UNKNOWN); return power; } try { // rescale the beta and sigma matrices and create a test RealMatrix scaledBeta = params.getBeta().scalarMultiply(betaScale, true); RealMatrix scaledSigmaError = params.getSigmaError().scalarMultiply(sigmaScale); GLMMTest glmmTest = GLMMTestFactory.createGLMMTestForPower(test, params.getFApproximationMethod(test), params.getUnivariateCdfMethod(test), params.getUnivariateEpsilonMethod(test), params.getDesignEssence(), params.getXtXInverse(), STARTING_SAMPLE_SIZE, params.getDesignRank(), params.getBetweenSubjectContrast(), params.getWithinSubjectContrast(), params.getTheta(), scaledBeta, scaledSigmaError, (params.getConfidenceIntervalType() != ConfidenceIntervalType.NONE ? params.getSampleSizeForEstimates() - params.getDesignMatrixRankForEstimates() : 0)); // calculate the noncentrality distribution NonCentralityDistribution nonCentralityDist = null; if (method != PowerMethod.CONDITIONAL_POWER) { nonCentralityDist = new NonCentralityDistribution(test, params.getDesignEssence(), params.getXtXInverse(), STARTING_SAMPLE_SIZE, params.getBetweenSubjectContrast(), params.getWithinSubjectContrast(), params.getTheta(), scaledBeta, scaledSigmaError, params.getSigmaGaussianRandom(), params.isNonCentralityCDFExact()); } // Calculate the maximum valid per group N. This avoids multiplication which exceeeds // Integer.MAX_VALUE. Moreover, by limiting to MAX_SAMPLE_SIZE, we avoid calculations // that take many seconds to complete. int designEssenceRows = params.getDesignEssence().getRowDimension(); int maxPerGroupN = Math.min(Integer.MAX_VALUE / designEssenceRows, MAX_SAMPLE_SIZE); /* * find the upper bound on sample size. That is, * find a sample size which produces power greater than or * equal to the desired power. * If an error occurs, we set an error code in the power result */ SampleSizeBound upperBound = getSampleSizeUpperBound(glmmTest, nonCentralityDist, method, targetPower, alpha, quantile, maxPerGroupN); if (upperBound.getError() != null) { GLMMPower power = new GLMMPower(test, alpha, targetPower, upperBound.getActualPower(), upperBound.getSampleSize(), betaScale, sigmaScale, method, quantile, null); switch (upperBound.getError()) { case MAX_SAMPLE_SIZE_EXCEEDED: power.setErrorMessage( // TODO: is the following expression correct? // "The total sample size for this case would exceed " + (maxPerGroupN * designEssenceRows) + ". " "The total sample size for this case would be unreasonably large. " + "For performance reasons, we are not computing its exact value."); power.setErrorCode(PowerErrorEnum.MAX_SAMPLE_SIZE_EXCEEDED); break; case SAMPLE_SIZE_UNDEFINED: power.setErrorMessage(NO_MEAN_DIFFERENCE_MESSAGE); power.setErrorCode(PowerErrorEnum.SAMPLE_SIZE_UNDEFINED); break; case SAMPLE_SIZE_UNDEFINED_DUE_TO_EXCEPTION: power.setErrorMessage("Sample size not well defined."); power.setErrorCode(PowerErrorEnum.SAMPLE_SIZE_UNDEFINED); break; } return power; } /* * find the lower bound on sample size * We search for a lower bound since the F approximations * may be unstable for very small samples */ // TODO: isn't the lower bound just half of the upper bound????? SampleSizeBound lowerBound = getSampleSizeLowerBound(glmmTest, nonCentralityDist, method, upperBound, alpha, quantile); assert lowerBound.getError() == null; /* * At this point we have valid boundaries for searching. * There are two possible scenarios * 1. The upper bound == lower bound. * 2. The upper bound != lower bound and lower bound exceeds required power. * In this case we just take the value at the lower bound. * 3. The upper bound != lower bound and lower bound is less than the required power. * In this case we bisection search */ double calculatedPower; int perGroupSampleSize; if (upperBound.getSampleSize() == lowerBound.getSampleSize()) { // case 1 calculatedPower = upperBound.getActualPower(); perGroupSampleSize = upperBound.getSampleSize(); } else if (lowerBound.getActualPower() >= targetPower) { // case 2 calculatedPower = lowerBound.getActualPower(); perGroupSampleSize = lowerBound.getSampleSize(); } else { // case 3, bisection search time! // create a bisection search function to find the best per group sample size BisectionSolver solver = new BisectionSolver(); SampleSizeFunction sampleSizeFunc = new SampleSizeFunction(glmmTest, nonCentralityDist, method, targetPower, alpha, quantile); double solution = solver.solve(MAX_ITERATIONS, sampleSizeFunc, lowerBound.getSampleSize(), upperBound.getSampleSize()); perGroupSampleSize = (int) Math.rint(solution); // see https://samplesizeshop.atlassian.net/browse/SSS-120 glmmTest.setPerGroupSampleSize(perGroupSampleSize); if (nonCentralityDist != null) nonCentralityDist.setPerGroupSampleSize(perGroupSampleSize); calculatedPower = getPowerByType(glmmTest, nonCentralityDist, method, alpha, quantile); } // build a confidence interval if requested GLMMPowerConfidenceInterval ci; if (params.getConfidenceIntervalType() != ConfidenceIntervalType.NONE) { ci = new GLMMPowerConfidenceInterval(params.getConfidenceIntervalType(), params.getAlphaLowerConfidenceLimit(), params.getAlphaUpperConfidenceLimit(), params.getSampleSizeForEstimates(), params.getDesignMatrixRankForEstimates(), alpha, glmmTest); } else { ci = null; } return new GLMMPower(test, alpha, targetPower, calculatedPower, MatrixUtils.getTotalSampleSize(params.getDesignEssence(), perGroupSampleSize), betaScale, sigmaScale, method, quantile, ci); } catch (PowerException pe) { GLMMPower powerValue = new GLMMPower(test, alpha, targetPower, -1, -1, betaScale, sigmaScale, method, quantile, null); powerValue.setErrorCode(pe.getErrorCode()); powerValue.setErrorMessage(pe.getMessage()); return powerValue; } }
From source file:edu.cudenver.bios.power.GLMMPowerCalculator.java
/** * Perform a bisection search to determine effect size * @param params GLMM input parameters/* w w w . jav a 2 s. c o m*/ * @return detectable difference object * @throws IllegalArgumentException * @throws ArithmeticException */ private GLMMPower getDetectableDifferenceValue(GLMMPowerParameters params, Test test, PowerMethod method, double alpha, double sigmaScale, double targetPower, int sampleSize, double quantile) throws PowerException { RealMatrix scaledBeta = params.getBeta().scalarMultiply(STARTING_BETA_SCALE, true); RealMatrix scaledSigmaError = params.getSigmaError().scalarMultiply(sigmaScale); GLMMTest glmmTest = GLMMTestFactory.createGLMMTestForPower(test, params.getFApproximationMethod(test), params.getUnivariateCdfMethod(test), params.getUnivariateEpsilonMethod(test), params.getDesignEssence(), params.getXtXInverse(), sampleSize, params.getDesignRank(), params.getBetweenSubjectContrast(), params.getWithinSubjectContrast(), params.getTheta(), scaledBeta, scaledSigmaError, (params.getConfidenceIntervalType() != ConfidenceIntervalType.NONE ? params.getSampleSizeForEstimates() - params.getDesignMatrixRankForEstimates() : 0)); NonCentralityDistribution nonCentralityDist = null; if (method != PowerMethod.CONDITIONAL_POWER) { nonCentralityDist = new NonCentralityDistribution(test, params.getDesignEssence(), params.getXtXInverse(), sampleSize, params.getBetweenSubjectContrast(), params.getWithinSubjectContrast(), params.getTheta(), scaledBeta, scaledSigmaError, params.getSigmaGaussianRandom(), params.isNonCentralityCDFExact()); } BisectionSolver solver = new BisectionSolver(); DetectableDifferenceFunction ddFunc = new DetectableDifferenceFunction(glmmTest, nonCentralityDist, method, targetPower, alpha, quantile, params.getBeta()); double lowerBound = 1.0E-10; // need non-zero lower bound or noncentrality dist malfunctions // check if the lower bound beta scale is already greater than the target power scaledBeta = params.getBeta().scalarMultiply(lowerBound, true); glmmTest.setBeta(scaledBeta); if (nonCentralityDist != null) nonCentralityDist.setBeta(scaledBeta); double calculatedPower = getPowerByType(glmmTest, nonCentralityDist, method, alpha, quantile); // find the detectable difference (i.e. beta scale) if the lower bound does not exceed the target power double betaScale = lowerBound; if (calculatedPower < targetPower) { DetectableDifferenceBound upperBound = getDetectableDifferenceUpperBound(glmmTest, nonCentralityDist, params.getBeta(), method, targetPower, alpha, quantile); if (upperBound.getError() != null) { GLMMPower power = new GLMMPower(test, alpha, targetPower, upperBound.getActualPower(), MatrixUtils.getTotalSampleSize(params.getDesignEssence(), sampleSize), upperBound.getBetaScale(), sigmaScale, method, quantile, null); switch (upperBound.getError()) { case MAX_BETA_SCALE_EXCEEDED: power.setErrorMessage("Failed to find valid upper bound on beta scale."); power.setErrorCode(PowerErrorEnum.MAX_BETA_SCALE_EXCEEDED); break; case BETA_SCALE_UNDEFINED: power.setErrorMessage("Beta scale not well defined for no difference."); power.setErrorCode(PowerErrorEnum.BETA_SCALE_UNDEFINED); break; } return power; } betaScale = solver.solve(MAX_ITERATIONS, ddFunc, lowerBound, upperBound.getBetaScale()); } // calculate actual power associated with this beta scale scaledBeta = params.getBeta().scalarMultiply(betaScale, true); glmmTest.setBeta(scaledBeta); if (nonCentralityDist != null) nonCentralityDist.setBeta(scaledBeta); calculatedPower = getPowerByType(glmmTest, nonCentralityDist, method, alpha, quantile); // build a confidence interval if requested GLMMPowerConfidenceInterval ci; if (params.getConfidenceIntervalType() != ConfidenceIntervalType.NONE) { ci = new GLMMPowerConfidenceInterval(params.getConfidenceIntervalType(), params.getAlphaLowerConfidenceLimit(), params.getAlphaUpperConfidenceLimit(), params.getSampleSizeForEstimates(), params.getDesignMatrixRankForEstimates(), alpha, glmmTest); } else { ci = null; } return new GLMMPower(test, alpha, targetPower, calculatedPower, MatrixUtils.getTotalSampleSize(params.getDesignEssence(), sampleSize), betaScale, sigmaScale, method, quantile, ci); }
From source file:objenome.goal.numeric.FindZeros.java
@Override public void solve(Objenome o, List<SetNumericValue> variables) { if (variables.size() == 1) { BisectionSolver solver = new BisectionSolver(); //bind variables values to objenome SetNumericValue var = variables.get(0); double best = solver.solve(1000, d -> { var.setValue(d); return eval(o); }, getMin(var, var), getMax(var, var)); //var.getMin().doubleValue(), var.getMax().doubleValue()); var.setValue(best); } else {/*from w w w . ja va2 s .co m*/ throw new RuntimeException("Unknown how to solve objenome " + o); } }