List of usage examples for org.apache.commons.math3.optim MaxEval MaxEval
public MaxEval(int max)
From source file:imagingbook.pub.fd.FourierDescriptor.java
/** * Calculates the 'canonical' start point. This version uses * (a) a coarse search for a global maximum of fp() and subsequently * (b) a numerical optimization using Brent's method * (implemented with Apache Commons Math). */// w ww .j av a2s . com public double getStartPointPhase(int Mp) { Mp = Math.min(Mp, (G.length - 1) / 2); UnivariateFunction fp = new TargetFunction(Mp); // search for the global maximum in coarse steps double cmax = Double.NEGATIVE_INFINITY; int kmax = -1; int K = 25; // number of steps over 180 degrees for (int k = 0; k < K; k++) { final double phi = Math.PI * k / K; // phase to evaluate final double c = fp.value(phi); if (c > cmax) { cmax = c; kmax = k; } } // optimize using previous and next point as the bracket. double minPhi = Math.PI * (kmax - 1) / K; double maxPhi = Math.PI * (kmax + 1) / K; UnivariateOptimizer optimizer = new BrentOptimizer(1E-4, 1E-6); int maxIter = 20; UnivariatePointValuePair result = optimizer.optimize(new MaxEval(maxIter), new UnivariateObjectiveFunction(fp), GoalType.MAXIMIZE, new SearchInterval(minPhi, maxPhi)); double phi0 = result.getPoint(); return phi0; // the canonical start point phase }
From source file:edu.cmu.tetrad.search.Mimbuild2.java
private void optimizeNonMeasureVariancesConditionally(Node[][] indicators, TetradMatrix measurescov, TetradMatrix latentscov, double[][] loadings, int[][] indicatorIndices, double[] delta) { int count = 0; for (int i = 0; i < indicators.length; i++) { for (int j = i; j < indicators.length; j++) { count++;/*from w w w . j a v a 2s.co m*/ } } for (int i = 0; i < indicators.length; i++) { for (int j = 0; j < indicators[i].length; j++) { count++; } } double[] values3 = new double[count]; count = 0; for (int i = 0; i < indicators.length; i++) { for (int j = i; j < indicators.length; j++) { values3[count] = latentscov.get(i, j); count++; } } for (int i = 0; i < indicators.length; i++) { for (int j = 0; j < indicators[i].length; j++) { values3[count] = loadings[i][j]; count++; } } Function2 function2 = new Function2(indicatorIndices, measurescov, loadings, latentscov, delta, count); MultivariateOptimizer search = new PowellOptimizer(1e-7, 1e-7); PointValuePair pair = search.optimize(new InitialGuess(values3), new ObjectiveFunction(function2), GoalType.MINIMIZE, new MaxEval(100000)); minimum = pair.getValue(); }
From source file:com.wwidesigner.optimization.ObjectiveFunctionOptimizer.java
protected static UnivariatePointValuePair runBrent(BrentOptimizer optimizer, BaseObjectiveFunction objective, double[] startPoint) throws TooManyEvaluationsException { UnivariatePointValuePair outcome;/*from w ww.j a v a 2 s. c om*/ outcome = optimizer.optimize(GoalType.MINIMIZE, new UnivariateObjectiveFunction(objective), new MaxEval(objective.getMaxEvaluations()), MaxIter.unlimited(), new SearchInterval(objective.getLowerBounds()[0], objective.getUpperBounds()[0], startPoint[0])); return outcome; }
From source file:edu.cmu.tetrad.search.Mimbuild2.java
private void optimizeMeasureVariancesConditionally(TetradMatrix measurescov, TetradMatrix latentscov, double[][] loadings, int[][] indicatorIndices, double[] delta) { double[] values2 = new double[delta.length]; int count = 0; for (int i = 0; i < delta.length; i++) { values2[count++] = delta[i];//from w w w.j ava 2s . c om } Function2 function2 = new Function2(indicatorIndices, measurescov, loadings, latentscov, delta, count); MultivariateOptimizer search = new PowellOptimizer(1e-7, 1e-7); PointValuePair pair = search.optimize(new InitialGuess(values2), new ObjectiveFunction(function2), GoalType.MINIMIZE, new MaxEval(100000)); minimum = pair.getValue(); }
From source file:gdsc.smlm.ij.plugins.BlinkEstimator.java
/** * Fit the dark time to counts of molecules curve. Only use the first n fitted points. * <p>/* www . j av a 2 s .c o m*/ * Calculates:<br/> * N = The number of photoblinking molecules in the sample<br/> * nBlink = The average number of blinks per flourophore<br/> * tOff = The off-time * * @param td * The dark time * @param ntd * The counts of molecules * @param nFittedPoints * @param log * Write the fitting results to the ImageJ log window * @return The fitted parameters [N, nBlink, tOff], or null if no fit was possible */ public double[] fit(double[] td, double[] ntd, int nFittedPoints, boolean log) { blinkingModel = new BlinkingFunction(); blinkingModel.setLogging(true); for (int i = 0; i < nFittedPoints; i++) blinkingModel.addPoint(td[i], ntd[i]); // Different convergence thresholds seem to have no effect on the resulting fit, only the number of // iterations for convergence double initialStepBoundFactor = 100; double costRelativeTolerance = 1e-6; double parRelativeTolerance = 1e-6; double orthoTolerance = 1e-6; double threshold = Precision.SAFE_MIN; LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(initialStepBoundFactor, costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold); try { double[] obs = blinkingModel.getY(); PointVectorValuePair optimum = optimizer.optimize(new MaxIter(1000), new MaxEval(Integer.MAX_VALUE), new ModelFunctionJacobian(new MultivariateMatrixFunction() { public double[][] value(double[] point) throws IllegalArgumentException { return blinkingModel.jacobian(point); } }), new ModelFunction(blinkingModel), new Target(obs), new Weight(blinkingModel.getWeights()), new InitialGuess(new double[] { ntd[0], 0.1, td[1] })); blinkingModel.setLogging(false); double[] parameters = optimum.getPoint(); double[] exp = optimum.getValue(); double mean = 0; for (double d : obs) mean += d; mean /= obs.length; double ssResiduals = 0, ssTotal = 0; for (int i = 0; i < obs.length; i++) { ssResiduals += (obs[i] - exp[i]) * (obs[i] - exp[i]); ssTotal += (obs[i] - mean) * (obs[i] - mean); } r2 = 1 - ssResiduals / ssTotal; adjustedR2 = getAdjustedCoefficientOfDetermination(ssResiduals, ssTotal, obs.length, parameters.length); if (log) { Utils.log(" Fit %d points. R^2 = %s. Adjusted R^2 = %s", obs.length, Utils.rounded(r2, 4), Utils.rounded(adjustedR2, 4)); Utils.log(" N=%s, nBlink=%s, tOff=%s (%s frames)", Utils.rounded(parameters[0], 4), Utils.rounded(parameters[1], 4), Utils.rounded(parameters[2], 4), Utils.rounded(parameters[2] / msPerFrame, 4)); } return parameters; } catch (TooManyIterationsException e) { if (log) Utils.log(" Failed to fit %d points: Too many iterations (%d)", blinkingModel.size(), optimizer.getIterations()); return null; } catch (ConvergenceException e) { if (log) Utils.log(" Failed to fit %d points", blinkingModel.size()); return null; } }
From source file:com.wwidesigner.optimization.ObjectiveFunctionOptimizer.java
protected static PointValuePair runBobyqa(BOBYQAOptimizer optimizer, BaseObjectiveFunction objective, double[] startPoint, int maxEvaluations) throws TooManyEvaluationsException { PointValuePair outcome;/* w w w .ja va 2s . com*/ EvaluatorInterface originalEvaluator = objective.getEvaluator(); if (objective.isRunTwoStageOptimization()) { objective.setEvaluator(objective.getFirstStageEvaluator()); outcome = optimizer.optimize(GoalType.MINIMIZE, new ObjectiveFunction(objective), new MaxEval(maxEvaluations), MaxIter.unlimited(), new InitialGuess(startPoint), new SimpleBounds(objective.getLowerBounds(), objective.getUpperBounds())); objective.setGeometryPoint(outcome.getPoint()); startPoint = objective.getInitialPoint(); } objective.setEvaluator(originalEvaluator); outcome = optimizer.optimize(GoalType.MINIMIZE, new ObjectiveFunction(objective), new MaxEval(maxEvaluations - optimizer.getEvaluations()), MaxIter.unlimited(), new InitialGuess(startPoint), new SimpleBounds(objective.getLowerBounds(), objective.getUpperBounds())); return outcome; }
From source file:com.itemanalysis.psychometrics.irt.equating.StockingLordMethodTest.java
@Test public void stockingLordTest4() { System.out.println("StockingLordMethod Test 4: Actual Distribution -backwards"); int n = aX.length; LinkedHashMap<String, ItemResponseModel> irmX = new LinkedHashMap<String, ItemResponseModel>(); LinkedHashMap<String, ItemResponseModel> irmY = new LinkedHashMap<String, ItemResponseModel>(); ItemResponseModel irm;/*ww w .j av a2s . c o m*/ for (int i = 0; i < n; i++) { String name = "V" + i; irm = new Irm3PL(aX[i], bX[i], cX[i], 1.0); irmX.put(name, irm); irm = new Irm3PL(aY[i], bY[i], cY[i], 1.0); irmY.put(name, irm); } UserSuppliedDistributionApproximation distX = new UserSuppliedDistributionApproximation(points, xDensity); UserSuppliedDistributionApproximation distY = new UserSuppliedDistributionApproximation(points, yDensity); StockingLordMethod sl = new StockingLordMethod(irmX, irmY, distX, distY, EquatingCriterionType.Q1); sl.setPrecision(4); double[] startValues = { 0, 1 }; int numIterpolationPoints = 2 * 2;//two dimensions A and B BOBYQAOptimizer underlying = new BOBYQAOptimizer(numIterpolationPoints); RandomGenerator g = new JDKRandomGenerator(); RandomVectorGenerator generator = new UncorrelatedRandomVectorGenerator(2, new GaussianRandomGenerator(g)); MultiStartMultivariateOptimizer optimizer = new MultiStartMultivariateOptimizer(underlying, 10, generator); PointValuePair optimum = optimizer.optimize(new MaxEval(1000), new ObjectiveFunction(sl), GoalType.MINIMIZE, SimpleBounds.unbounded(2), new InitialGuess(startValues)); double[] slCoefficients = optimum.getPoint(); sl.setIntercept(slCoefficients[0]); sl.setScale(slCoefficients[1]); System.out.println(" Iterations: " + optimizer.getEvaluations()); System.out.println(" fmin: " + optimum.getValue()); System.out.println(" B = " + slCoefficients[0] + " A = " + slCoefficients[1]); assertEquals(" Intercept test", -0.4834, sl.getIntercept(), 1e-4); assertEquals(" Scale test", 1.0967, sl.getScale(), 1e-4); System.out.println(); }
From source file:edu.cmu.tetrad.search.Mimbuild2.java
private void optimizeAllParamsSimultaneously(Node[][] indicators, TetradMatrix measurescov, TetradMatrix latentscov, double[][] loadings, int[][] indicatorIndices, double[] delta) { double[] values = getAllParams(indicators, latentscov, loadings, delta); Function4 function = new Function4(indicatorIndices, measurescov, loadings, latentscov, delta); MultivariateOptimizer search = new PowellOptimizer(1e-7, 1e-7); PointValuePair pair = search.optimize(new InitialGuess(values), new ObjectiveFunction(function), GoalType.MINIMIZE, new MaxEval(100000)); minimum = pair.getValue();/*from w w w .j a va 2 s.c o m*/ }
From source file:gdsc.smlm.fitting.nonlinear.MaximumLikelihoodFitter.java
public FitStatus fit(int n, double[] y, double[] y_fit, double[] a, double[] a_dev, double[] error, double noise) { numberOfFittedPoints = n;/*from w ww . j ava 2 s . co m*/ LikelihoodWrapper maximumLikelihoodFunction; // We can use different likelihood wrapper functions: switch (likelihoodFunction) { case POISSON_GAMMA_GAUSSIAN: // Poisson-Gamma-Gaussian - EM-CCD data if (alpha > 0 && sigma > 0) { maximumLikelihoodFunction = new PoissonGammaGaussianLikelihoodWrapper(f, a, y, n, alpha, sigma); break; } case POISSON_GAUSSIAN: // Poisson-Gaussian - CCD data if (sigma > 0) { maximumLikelihoodFunction = new PoissonGaussianLikelihoodWrapper(f, a, y, n, sigma); break; } case POISSON: default: // Poisson - most counting data maximumLikelihoodFunction = new PoissonLikelihoodWrapper(f, a, y, n); } // Check if the method requires the gradient but it cannot be computed if (searchMethod.usesGradient && !maximumLikelihoodFunction.canComputeGradient()) { maximumLikelihoodFunction = new PoissonLikelihoodWrapper(f, a, y, n); } try { double[] startPoint = getInitialSolution(a); PointValuePair optimum = null; if (searchMethod == SearchMethod.POWELL || searchMethod == SearchMethod.POWELL_BOUNDED) { // Non-differentiable version using Powell Optimiser // This is as per the method in Numerical Recipes 10.5 (Direction Set (Powell's) method) // I could extend the optimiser and implement bounds on the directions moved. However the mapping // adapter seems to work OK. final boolean basisConvergence = false; // Perhaps these thresholds should be tighter? // The default is to use the sqrt() of the overall tolerance //final double lineRel = FastMath.sqrt(relativeThreshold); //final double lineAbs = FastMath.sqrt(absoluteThreshold); //final double lineRel = relativeThreshold * 1e2; //final double lineAbs = absoluteThreshold * 1e2; // Since we are fitting only a small number of parameters then just use the same tolerance // for each search direction final double lineRel = relativeThreshold; final double lineAbs = absoluteThreshold; CustomPowellOptimizer o = new CustomPowellOptimizer(relativeThreshold, absoluteThreshold, lineRel, lineAbs, null, basisConvergence); OptimizationData maxIterationData = null; if (getMaxIterations() > 0) maxIterationData = new MaxIter(getMaxIterations()); if (searchMethod == SearchMethod.POWELL) { if (powellFunction == null) { // We must map all the parameters into the same range. This is done in the Mortensen MLE // Python code by using the sqrt of the number of photons and background. if (mapGaussian) { Gaussian2DFunction gf = (Gaussian2DFunction) f; // Re-map signal and background using the sqrt int[] indices = gf.gradientIndices(); int[] map = new int[indices.length]; int count = 0; // Background is always first if (indices[0] == Gaussian2DFunction.BACKGROUND) { map[count++] = 0; } // Look for the Signal in multiple peak 2D Gaussians for (int i = 1; i < indices.length; i++) if (indices[i] % 6 == Gaussian2DFunction.SIGNAL) { map[count++] = i; } if (count > 0) { powellFunction = new MappedMultivariateLikelihood(maximumLikelihoodFunction, Arrays.copyOf(map, count)); } } if (powellFunction == null) { powellFunction = new MultivariateLikelihood(maximumLikelihoodFunction); } } // Update the maximum likelihood function in the Powell function wrapper powellFunction.fun = maximumLikelihoodFunction; OptimizationData positionChecker = null; // new org.apache.commons.math3.optim.PositionChecker(relativeThreshold, absoluteThreshold); if (powellFunction.isMapped()) { MappedMultivariateLikelihood adapter = (MappedMultivariateLikelihood) powellFunction; optimum = o.optimize(maxIterationData, new MaxEval(getMaxEvaluations()), new ObjectiveFunction(powellFunction), GoalType.MINIMIZE, new InitialGuess(adapter.map(startPoint)), positionChecker); double[] solution = adapter.unmap(optimum.getPointRef()); optimum = new PointValuePair(solution, optimum.getValue()); } else { optimum = o.optimize(maxIterationData, new MaxEval(getMaxEvaluations()), new ObjectiveFunction(powellFunction), GoalType.MINIMIZE, new InitialGuess(startPoint), positionChecker); } } else { // Try using the mapping adapter for a bounded Powell search MultivariateFunctionMappingAdapter adapter = new MultivariateFunctionMappingAdapter( new MultivariateLikelihood(maximumLikelihoodFunction), lower, upper); optimum = o.optimize(maxIterationData, new MaxEval(getMaxEvaluations()), new ObjectiveFunction(adapter), GoalType.MINIMIZE, new InitialGuess(adapter.boundedToUnbounded(startPoint))); double[] solution = adapter.unboundedToBounded(optimum.getPointRef()); optimum = new PointValuePair(solution, optimum.getValue()); } iterations = o.getIterations(); evaluations = o.getEvaluations(); } else if (searchMethod == SearchMethod.BOBYQA) { // Differentiable approximation using Powell's BOBYQA algorithm. // This is slower than the Powell optimiser and requires a high number of evaluations. int numberOfInterpolationPoints = this.getNumberOfFittedParameters() + 2; BOBYQAOptimizer o = new BOBYQAOptimizer(numberOfInterpolationPoints); optimum = o.optimize(new MaxEval(getMaxEvaluations()), new ObjectiveFunction(new MultivariateLikelihood(maximumLikelihoodFunction)), GoalType.MINIMIZE, new InitialGuess(startPoint), new SimpleBounds(lower, upper)); iterations = o.getIterations(); evaluations = o.getEvaluations(); } else if (searchMethod == SearchMethod.CMAES) { // TODO - Understand why the CMAES optimiser does not fit very well on test data. It appears // to converge too early and the likelihood scores are not as low as the other optimisers. // CMAESOptimiser based on Matlab code: // https://www.lri.fr/~hansen/cmaes.m // Take the defaults from the Matlab documentation double stopFitness = 0; //Double.NEGATIVE_INFINITY; boolean isActiveCMA = true; int diagonalOnly = 0; int checkFeasableCount = 1; RandomGenerator random = new Well19937c(); boolean generateStatistics = false; // The sigma determines the search range for the variables. It should be 1/3 of the initial search region. double[] sigma = new double[lower.length]; for (int i = 0; i < sigma.length; i++) sigma[i] = (upper[i] - lower[i]) / 3; int popSize = (int) (4 + Math.floor(3 * Math.log(sigma.length))); // The CMAES optimiser is random and restarting can overcome problems with quick convergence. // The Apache commons documentations states that convergence should occur between 30N and 300N^2 // function evaluations final int n30 = FastMath.min(sigma.length * sigma.length * 30, getMaxEvaluations() / 2); evaluations = 0; OptimizationData[] data = new OptimizationData[] { new InitialGuess(startPoint), new CMAESOptimizer.PopulationSize(popSize), new MaxEval(getMaxEvaluations()), new CMAESOptimizer.Sigma(sigma), new ObjectiveFunction(new MultivariateLikelihood(maximumLikelihoodFunction)), GoalType.MINIMIZE, new SimpleBounds(lower, upper) }; // Iterate to prevent early convergence int repeat = 0; while (evaluations < n30) { if (repeat++ > 1) { // Update the start point and population size data[0] = new InitialGuess(optimum.getPointRef()); popSize *= 2; data[1] = new CMAESOptimizer.PopulationSize(popSize); } CMAESOptimizer o = new CMAESOptimizer(getMaxIterations(), stopFitness, isActiveCMA, diagonalOnly, checkFeasableCount, random, generateStatistics, new SimpleValueChecker(relativeThreshold, absoluteThreshold)); PointValuePair result = o.optimize(data); iterations += o.getIterations(); evaluations += o.getEvaluations(); //System.out.printf("CMAES [%d] i=%d [%d], e=%d [%d]\n", repeat, o.getIterations(), iterations, // o.getEvaluations(), totalEvaluations); if (optimum == null || result.getValue() < optimum.getValue()) { optimum = result; } } } else if (searchMethod == SearchMethod.BFGS) { // BFGS can use an approximate line search minimisation where as Powell and conjugate gradient // methods require a more accurate line minimisation. The BFGS search does not do a full // minimisation but takes appropriate steps in the direction of the current gradient. // Do not use the convergence checker on the value of the function. Use the convergence on the // point coordinate and gradient //BFGSOptimizer o = new BFGSOptimizer(new SimpleValueChecker(rel, abs)); BFGSOptimizer o = new BFGSOptimizer(); // Configure maximum step length for each dimension using the bounds double[] stepLength = new double[lower.length]; for (int i = 0; i < stepLength.length; i++) { stepLength[i] = (upper[i] - lower[i]) * 0.3333333; if (stepLength[i] <= 0) stepLength[i] = Double.POSITIVE_INFINITY; } // The GoalType is always minimise so no need to pass this in OptimizationData positionChecker = null; //new org.apache.commons.math3.optim.PositionChecker(relativeThreshold, absoluteThreshold); optimum = o.optimize(new MaxEval(getMaxEvaluations()), new ObjectiveFunctionGradient(new MultivariateVectorLikelihood(maximumLikelihoodFunction)), new ObjectiveFunction(new MultivariateLikelihood(maximumLikelihoodFunction)), new InitialGuess(startPoint), new SimpleBounds(lowerConstraint, upperConstraint), new BFGSOptimizer.GradientTolerance(relativeThreshold), positionChecker, new BFGSOptimizer.StepLength(stepLength)); iterations = o.getIterations(); evaluations = o.getEvaluations(); } else { // The line search algorithm often fails. This is due to searching into a region where the // function evaluates to a negative so has been clipped. This means the upper bound of the line // cannot be found. // Note that running it on an easy problem (200 photons with fixed fitting (no background)) the algorithm // does sometimes produces results better than the Powell algorithm but it is slower. BoundedNonLinearConjugateGradientOptimizer o = new BoundedNonLinearConjugateGradientOptimizer( (searchMethod == SearchMethod.CONJUGATE_GRADIENT_FR) ? Formula.FLETCHER_REEVES : Formula.POLAK_RIBIERE, new SimpleValueChecker(relativeThreshold, absoluteThreshold)); // Note: The gradients may become unstable at the edge of the bounds. Or they will not change // direction if the true solution is on the bounds since the gradient will always continue // towards the bounds. This is key to the conjugate gradient method. It searches along a vector // until the direction of the gradient is in the opposite direction (using dot products, i.e. // cosine of angle between them) // NR 10.7 states there is no advantage of the variable metric DFP or BFGS methods over // conjugate gradient methods. So I will try these first. // Try this: // Adapt the conjugate gradient optimiser to use the gradient to pick the search direction // and then for the line minimisation. However if the function is out of bounds then clip the // variables at the bounds and continue. // If the current point is at the bounds and the gradient is to continue out of bounds then // clip the gradient too. // Or: just use the gradient for the search direction then use the line minimisation/rest // as per the Powell optimiser. The bounds should limit the search. // I tried a Bounded conjugate gradient optimiser with clipped variables: // This sometimes works. However when the variables go a long way out of the expected range the gradients // can have vastly different magnitudes. This results in the algorithm stalling since the gradients // can be close to zero and the some of the parameters are no longer adjusted. // Perhaps this can be looked for and the algorithm then gives up and resorts to a Powell optimiser from // the current point. // Changed the bracketing step to very small (default is 1, changed to 0.001). This improves the // performance. The gradient direction is very sensitive to small changes in the coordinates so a // tighter bracketing of the line search helps. // Tried using a non-gradient method for the line search copied from the Powell optimiser: // This also works when the bracketing step is small but the number of iterations is higher. // 24.10.2014: I have tried to get conjugate gradient to work but the gradient function // must not behave suitably for the optimiser. In the current state both methods of using a // Bounded Conjugate Gradient Optimiser perform poorly relative to other optimisers: // Simulated : n=1000, signal=200, x=0.53, y=0.47 // LVM : n=1000, signal=171, x=0.537, y=0.471 (1.003s) // Powell : n=1000, signal=187, x=0.537, y=0.48 (1.238s) // Gradient based PR (constrained): n=858, signal=161, x=0.533, y=0.474 (2.54s) // Gradient based PR (bounded): n=948, signal=161, x=0.533, y=0.473 (2.67s) // Non-gradient based : n=1000, signal=151.47, x=0.535, y=0.474 (1.626s) // The conjugate optimisers are slower, under predict the signal by the most and in the case of // the gradient based optimiser, fail to converge on some problems. This is worse when constrained // fitting is used and not tightly bounded fitting. // I will leave the code in as an option but would not recommend using it. I may remove it in the // future. // Note: It is strange that the non-gradient based line minimisation is more successful. // It may be that the gradient function is not accurate (due to round off error) or that it is // simply wrong when far from the optimum. My JUnit tests only evaluate the function within the // expected range of the answer. // Note the default step size on the Powell optimiser is 1 but the initial directions are unit vectors. // So our bracketing step should be a minimum of 1 / average length of the first gradient vector to prevent // the first step being too large when bracketing. final double gradient[] = new double[startPoint.length]; maximumLikelihoodFunction.likelihood(startPoint, gradient); double l = 0; for (double d : gradient) l += d * d; final double bracketingStep = FastMath.min(0.001, ((l > 1) ? 1.0 / l : 1)); //System.out.printf("Bracketing step = %f (length=%f)\n", bracketingStep, l); o.setUseGradientLineSearch(gradientLineMinimisation); optimum = o.optimize(new MaxEval(getMaxEvaluations()), new ObjectiveFunctionGradient(new MultivariateVectorLikelihood(maximumLikelihoodFunction)), new ObjectiveFunction(new MultivariateLikelihood(maximumLikelihoodFunction)), GoalType.MINIMIZE, new InitialGuess(startPoint), new SimpleBounds(lowerConstraint, upperConstraint), new BoundedNonLinearConjugateGradientOptimizer.BracketingStep(bracketingStep)); iterations = o.getIterations(); evaluations = o.getEvaluations(); //maximumLikelihoodFunction.value(solution, gradient); //System.out.printf("Iter = %d, %g @ %s : %s\n", iterations, ll, Arrays.toString(solution), // Arrays.toString(gradient)); } final double[] solution = optimum.getPointRef(); setSolution(a, solution); //System.out.printf("Iter = %d, Eval = %d, %g @ %s\n", iterations, evaluations, optimum.getValue(), // java.util.Arrays.toString(solution)); // Compute residuals for the FunctionSolver interface if (y_fit == null || y_fit.length < n) y_fit = new double[n]; f.initialise(a); residualSumOfSquares = 0; for (int i = 0; i < n; i++) { y_fit[i] = f.eval(i); final double residual = y[i] - y_fit[i]; residualSumOfSquares += residual * residual; } if (a_dev != null) { // Assume the Maximum Likelihood estimator returns the optimum fit (achieves the Cramer Roa // lower bounds) and so the covariance can be obtained from the Fisher Information Matrix. final int[] gradientIndices = f.gradientIndices(); final int nparams = gradientIndices.length; GradientCalculator calculator = GradientCalculatorFactory.newCalculator(nparams); final double[] I = calculator.fisherInformationDiagonal(n, a, f); for (int i = 0; i < gradientIndices.length; i++) a_dev[gradientIndices[i]] = 1.0 / Math.sqrt(I[i]); } error[0] = NonLinearFit.getError(residualSumOfSquares, noise, n, f.gradientIndices().length); totalSumOfSquares = getSumOfSquares(n, y); } catch (TooManyIterationsException e) { //System.out.printf("Too many iterations = %d\n", e.getMax()); //e.printStackTrace(); return FitStatus.FAILED_TO_CONVERGE; } catch (TooManyEvaluationsException e) { //System.out.printf("Too many evaluations = %d\n", e.getMax()); //e.printStackTrace(); return FitStatus.FAILED_TO_CONVERGE; } catch (ConvergenceException e) { // Occurs when QR decomposition fails - mark as a singular non-linear model (no solution) //System.out.printf("Singular non linear model = %s\n", e.getMessage()); return FitStatus.SINGULAR_NON_LINEAR_MODEL; } catch (BFGSOptimizer.LineSearchRoundoffException e) { //System.out.println("BFGS error: " + e.getMessage()); //e.printStackTrace(); return FitStatus.FAILED_TO_CONVERGE; } catch (Exception e) { //System.out.printf("Unknown error = %s\n", e.getMessage()); e.printStackTrace(); return FitStatus.UNKNOWN; } return FitStatus.OK; }
From source file:com.itemanalysis.psychometrics.irt.equating.HaebaraMethodTest.java
@Test public void haebaraMethodTestMixedFormat() { System.out.println("Haebara Test 5: Mixed format test, symmetric"); LinkedHashMap<String, ItemResponseModel> irmX = new LinkedHashMap<String, ItemResponseModel>(); LinkedHashMap<String, ItemResponseModel> irmY = new LinkedHashMap<String, ItemResponseModel>(); //Form X/*from w w w. j a va2 s . c o m*/ irmX.put("v1", new Irm3PL(0.751335, -0.897391, 0.244001, 1.7)); irmX.put("v2", new Irm3PL(0.955947, -0.811477, 0.242883, 1.7)); irmX.put("v3", new Irm3PL(0.497206, -0.858681, 0.260893, 1.7)); irmX.put("v4", new Irm3PL(0.724000, -0.123911, 0.243497, 1.7)); irmX.put("v5", new Irm3PL(0.865200, 0.205889, 0.319135, 1.7)); irmX.put("v6", new Irm3PL(0.658129, 0.555228, 0.277826, 1.7)); irmX.put("v7", new Irm3PL(1.082118, 0.950549, 0.157979, 1.7)); irmX.put("v8", new Irm3PL(0.988294, 1.377501, 0.084828, 1.7)); irmX.put("v9", new Irm3PL(1.248923, 1.614355, 0.181874, 1.7)); irmX.put("v10", new Irm3PL(1.116682, 2.353932, 0.246856, 1.7)); irmX.put("v11", new Irm3PL(0.438171, 3.217965, 0.309243, 1.7)); irmX.put("v12", new Irm3PL(1.082206, 4.441864, 0.192339, 1.7)); double[] step1 = { 0.0, 1.101266, -1.09327 }; irmX.put("v13", new IrmGPCM(0.269994, step1, 1.7)); double[] step2 = { 0.0, 1.739176, 1.526148 }; irmX.put("v14", new IrmGPCM(0.972506, step2, 1.7)); double[] step3 = { 0.0, 5.566958, 1.362356 }; irmX.put("v15", new IrmGPCM(0.378812, step3, 1.7)); double[] step4 = { 0.0, 0.533540, 2.091335, 0.405283 }; irmX.put("v16", new IrmGPCM(0.537706, step4, 1.7)); double[] step5 = { 0.0, 3.440463, 2.235171, 1.62318 }; irmX.put("v17", new IrmGPCM(0.554506, step5, 1.7)); //Form Y irmY.put("v1", new Irm3PL(0.887276, -1.334798, 0.134406, 1.7)); irmY.put("v2", new Irm3PL(1.184412, -1.129004, 0.237765, 1.7)); irmY.put("v3", new Irm3PL(0.609412, -1.464546, 0.151393, 1.7)); irmY.put("v4", new Irm3PL(0.923812, -0.576435, 0.240097, 1.7)); irmY.put("v5", new Irm3PL(0.822776, -0.476357, 0.192369, 1.7)); irmY.put("v6", new Irm3PL(0.707818, -0.235189, 0.189557, 1.7)); irmY.put("v7", new Irm3PL(1.306976, 0.242986, 0.165553, 1.7)); irmY.put("v8", new Irm3PL(1.295471, 0.598029, 0.090557, 1.7)); irmY.put("v9", new Irm3PL(1.366841, 0.923206, 0.172993, 1.7)); irmY.put("v10", new Irm3PL(1.389624, 1.380666, 0.238008, 1.7)); irmY.put("v11", new Irm3PL(0.293806, 2.028070, 0.203448, 1.7)); irmY.put("v12", new Irm3PL(0.885347, 3.152928, 0.195473, 1.7)); double[] step1Y = { 0.0, 0.399117, -1.38735 }; irmY.put("v13", new IrmGPCM(0.346324, step1Y, 1.7)); double[] step2Y = { 0.0, 0.956014, 0.756514 }; irmY.put("v14", new IrmGPCM(1.252012, step2Y, 1.7)); double[] step3Y = { 0.0, 4.676299, 0.975303 }; irmY.put("v15", new IrmGPCM(0.392282, step3Y, 1.7)); double[] step4Y = { 0.0, 0.042549, 1.104823, -0.118440 }; irmY.put("v16", new IrmGPCM(0.660841, step4Y, 1.7)); double[] step5Y = { 0.0, 2.645241, 1.536046, 0.748514 }; irmY.put("v17", new IrmGPCM(0.669612, step5Y, 1.7)); UniformDistributionApproximation uniform = new UniformDistributionApproximation(-3.0, 3.0, 25); HaebaraMethod hb = new HaebaraMethod(irmX, irmY, uniform, uniform, EquatingCriterionType.Q1Q2); hb.setPrecision(4); double[] startValues = { 0, 1 }; int numIterpolationPoints = 2 * 2;//two dimensions A and B BOBYQAOptimizer underlying = new BOBYQAOptimizer(numIterpolationPoints); RandomGenerator g = new JDKRandomGenerator(); RandomVectorGenerator generator = new UncorrelatedRandomVectorGenerator(2, new GaussianRandomGenerator(g)); MultiStartMultivariateOptimizer optimizer = new MultiStartMultivariateOptimizer(underlying, 10, generator); org.apache.commons.math3.optim.PointValuePair optimum = optimizer.optimize(new MaxEval(1000), new ObjectiveFunction(hb), org.apache.commons.math3.optim.nonlinear.scalar.GoalType.MINIMIZE, SimpleBounds.unbounded(2), new InitialGuess(startValues)); double[] hbCoefficients = optimum.getPoint(); hb.setIntercept(hbCoefficients[0]); hb.setScale(hbCoefficients[1]); System.out.println(" Iterations: " + optimizer.getEvaluations()); System.out.println(" fmin: " + optimum.getValue()); System.out.println(" B = " + hbCoefficients[0] + " A = " + hbCoefficients[1]); assertEquals(" Intercept test", -0.437427, hb.getIntercept(), 1e-4); assertEquals(" Scale test", 0.810364, hb.getScale(), 1e-4); System.out.println(); }