List of usage examples for org.apache.commons.math3.optim.univariate SearchInterval SearchInterval
public SearchInterval(double lo, double hi)
From source file:com.google.cloud.genomics.dataflow.functions.verifybamid.Solver.java
/** * Maximizes a univariate function using a grid search followed by Brent's algorithm. * * @param fn the likelihood function to minimize * @param gridStart the lower bound for the grid search * @param gridEnd the upper bound for the grid search * @param gridStep step size for the grid search * @param relErr relative error tolerance for Brent's algorithm * @param absErr absolute error tolerance for Brent's algorithm * @param maxIter maximum # of iterations to perform in Brent's algorithm * @param maxEval maximum # of Likelihood function evaluations in Brent's algorithm * * @return the value of the parameter that maximizes the function *///w w w .j a v a2 s .com public static double maximize(UnivariateFunction fn, double gridStart, double gridEnd, double gridStep, double relErr, double absErr, int maxIter, int maxEval) { Interval interval = gridSearch(fn, gridStart, gridEnd, gridStep); BrentOptimizer bo = new BrentOptimizer(relErr, absErr); UnivariatePointValuePair max = bo.optimize(new MaxIter(maxIter), new MaxEval(maxEval), new SearchInterval(interval.getInf(), interval.getSup()), new UnivariateObjectiveFunction(fn), GoalType.MAXIMIZE); return max.getPoint(); }
From source file:com.itemanalysis.psychometrics.irt.equating.IrtScaleLinking.java
public void computeCoefficients() { raschFamily = checkRaschModel();//w w w . j av a 2 s.c o m ms.setPrecision(precision); mm.setPrecision(precision); hb.setPrecision(precision); sl.setPrecision(precision); if (raschFamily) { double[] sv = { mm.getIntercept() }; UnivariateOptimizer underlying = new BrentOptimizer(1e-10, 1e-14); JDKRandomGenerator g = new JDKRandomGenerator(); //Haebara method MultiStartUnivariateOptimizer optimizer = new MultiStartUnivariateOptimizer(underlying, 5, g);//Five random starts to Brent optimizer. UnivariatePointValuePair hbPair = optimizer.optimize(new MaxEval(500), new UnivariateObjectiveFunction(hb), GoalType.MINIMIZE, new SearchInterval(-4, 4), new InitialGuess(sv)); hb.setIntercept(hbPair.getPoint()); hb.setScale(1.0); fHB = hbPair.getValue(); //Stocking-Lord method UnivariatePointValuePair slPair = optimizer.optimize(new MaxEval(500), new UnivariateObjectiveFunction(sl), GoalType.MINIMIZE, new SearchInterval(-4, 4), new InitialGuess(sv)); sl.setIntercept(slPair.getPoint()); sl.setScale(1.0); fSL = slPair.getValue(); } else { double[] hbStartValues = { mm.getIntercept(), mm.getScale() }; double[] slStartValues = { mm.getIntercept(), mm.getScale() }; if (useUncmin) { DefaultUncminOptimizer optimizer = new DefaultUncminOptimizer(); try { optimizer.minimize(hb, hbStartValues); double[] param = optimizer.getParameters(); fHB = optimizer.getFunctionValue(); hb.setIntercept(param[0]); if (param.length > 1) { hb.setScale(param[1]); } else { hb.setScale(1.0);//Rasch family of models } optimizer.minimize(sl, slStartValues); param = optimizer.getParameters(); fSL = optimizer.getFunctionValue(); sl.setIntercept(param[0]); if (param.length > 1) { sl.setScale(param[1]); } else { sl.setScale(1.0);//Rasch family of models } } catch (UncminException ex) { ex.printStackTrace(); } } else { 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 hbOptimum = optimizer.optimize(new MaxEval(1000), new ObjectiveFunction(hb), GoalType.MINIMIZE, SimpleBounds.unbounded(2), new InitialGuess(hbStartValues)); double[] hbCoefficients = hbOptimum.getPoint(); hb.setIntercept(hbCoefficients[0]); hb.setScale(hbCoefficients[1]); fHB = hbOptimum.getValue(); PointValuePair slOptimum = optimizer.optimize(new MaxEval(1000), new ObjectiveFunction(sl), GoalType.MINIMIZE, SimpleBounds.unbounded(2), new InitialGuess(slStartValues)); double[] slCoefficients = slOptimum.getPoint(); sl.setIntercept(slCoefficients[0]); sl.setScale(slCoefficients[1]); fSL = slOptimum.getValue(); } } }
From source file:edu.uchc.octane.GaussianFitAstigmatism.java
/** * calculate the Z coordinate based on astigmatism */// w w w .ja v a2 s .co m void calculateZ() { if (calibration_ == null) { z_ = 0; return; } UnivariateFunction func = new UnivariateFunction() { @Override public double value(double z) { double sigmax = FastMath.sqrt(pvp_.getPoint()[3] / 2); double sigmay = FastMath.sqrt(pvp_.getPoint()[4] / 2); double vx = calibration_[0] + (z - calibration_[1]) * (z - calibration_[1]) * calibration_[2] - sigmax; double vy = calibration_[3] + (z - calibration_[4]) * (z - calibration_[4]) * calibration_[5] - sigmay; return vx * vx + vy * vy; } }; BrentOptimizer optimizer = new BrentOptimizer(1e-4, 1e-4); UnivariatePointValuePair upvp = optimizer.optimize(new UnivariateObjectiveFunction(func), GoalType.MINIMIZE, MaxEval.unlimited(), new SearchInterval(2 * z0min_ - z0max_, 2 * z0max_ - z0min_)); if (upvp.getValue() > errTol_) { throw (new ConvergenceException()); } z_ = upvp.getPoint(); }
From source file:com.itemanalysis.psychometrics.irt.estimation.IrtExaminee.java
/** * Maximum likelihood estimate (MLE) of examinee ability. * * @param thetaMin smallest possible ability estimate (lower bound on BrentOptimizer) * @param thetaMax largest possible ability estimate (upper bound on BrentOptimizer) * @return MLE of examinee ability//from w w w . j ava 2s. co m */ public double maximumLikelihoodEstimate(double thetaMin, double thetaMax) { method = EstimationMethod.ML; UnivariateOptimizer optimizer = new BrentOptimizer(1e-10, 1e-14); UnivariatePointValuePair pair = optimizer.optimize(new MaxEval(100), new UnivariateObjectiveFunction(this), GoalType.MAXIMIZE, new SearchInterval(thetaMin, thetaMax)); estimatedTheta = pair.getPoint(); return estimatedTheta; }
From source file:com.itemanalysis.psychometrics.irt.estimation.IrtExaminee.java
/** * Maximum a Posteriori (MAP) estimate of examinee ability using a normal prior * distribution.//from w w w . j ava 2s .c om * * @param mean mean of normal prior distribution * @param sd standard deviation of prior distribution * @param thetaMin smallest possible ability estimate (lower bound on BrentOptimizer) * @param thetaMax largest possible ability estimate (upper bound on BrentOptimizer) * @return MAP estimate of examinee ability */ public double mapEstimate(double mean, double sd, double thetaMin, double thetaMax) { mapPrior = new NormalDistribution(mean, sd); method = EstimationMethod.MAP; UnivariateOptimizer optimizer = new BrentOptimizer(1e-10, 1e-14); UnivariatePointValuePair pair = optimizer.optimize(new MaxEval(100), new UnivariateObjectiveFunction(this), GoalType.MAXIMIZE, new SearchInterval(thetaMin, thetaMax)); estimatedTheta = pair.getPoint(); return estimatedTheta; }
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). *///from ww w . java 2 s . c o m 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:com.itemanalysis.psychometrics.irt.equating.HaebaraMethodTest.java
/** * This test uses Rasch model items. True results from the sirt package in R. * * *//* w w w. ja v a 2s. co m*/ @Test public void raschModelTest() { System.out.println("Haebara test: Rasch model"); LinkedHashMap<String, ItemResponseModel> itemFormX = new LinkedHashMap<String, ItemResponseModel>(); LinkedHashMap<String, ItemResponseModel> itemFormY = new LinkedHashMap<String, ItemResponseModel>(); itemFormX.put("Item1", new Irm3PL(-3.188047976, 1.0)); itemFormX.put("Item2", new Irm3PL(1.031760328, 1.0)); itemFormX.put("Item3", new Irm3PL(0.819040914, 1.0)); itemFormX.put("Item4", new Irm3PL(-2.706947360, 1.0)); itemFormX.put("Item5", new Irm3PL(-0.094527077, 1.0)); itemFormX.put("Item6", new Irm3PL(0.689697135, 1.0)); itemFormX.put("Item7", new Irm3PL(-0.551837153, 1.0)); itemFormX.put("Item8", new Irm3PL(-0.359559276, 1.0)); itemFormY.put("Item1", new Irm3PL(-3.074599226, 1.0)); itemFormY.put("Item2", new Irm3PL(1.012824350, 1.0)); itemFormY.put("Item3", new Irm3PL(0.868538408, 1.0)); itemFormY.put("Item4", new Irm3PL(-2.404483603, 1.0)); itemFormY.put("Item5", new Irm3PL(0.037402866, 1.0)); itemFormY.put("Item6", new Irm3PL(0.700747420, 1.0)); itemFormY.put("Item7", new Irm3PL(-0.602555046, 1.0)); itemFormY.put("Item8", new Irm3PL(-0.350426446, 1.0)); double f = 0; double[] param1 = new double[2]; double[] param2 = new double[2]; UniformDistributionApproximation distX = new UniformDistributionApproximation(-4.0, 4.0, 161);//plink default UniformDistributionApproximation distY = new UniformDistributionApproximation(-4.0, 4.0, 161);//plink default HaebaraMethod hb = new HaebaraMethod(itemFormX, itemFormY, distX, distY, EquatingCriterionType.Q1Q2); hb.setPrecision(6); double[] startValues = { 0 };//Using a start value for the intercept only for the Rasch family of models. DefaultUncminOptimizer optimizer = new DefaultUncminOptimizer(); try { optimizer.minimize(hb, startValues); double[] r = optimizer.getParameters(); param1[0] = r[0]; param1[1] = 1; f = optimizer.getFunctionValue(); hb.setIntercept(param1[0]); hb.setScale(param1[1]); } catch (UncminException ex) { ex.printStackTrace(); } //Check Brent Optimizer values against results from plink. System.out.println(" UNCMIN Optimization"); System.out.println(" Iterations: "); System.out.println(" fmin: " + f); System.out.println(" B = " + hb.getIntercept() + " A = " + hb.getScale()); assertEquals(" Intercept test", 0.065032, hb.getIntercept(), 1e-3);//true result from plink package in R assertEquals(" Scale test", 1.0, hb.getScale(), 1e-4); BrentOptimizer brentOptimizer = new BrentOptimizer(1e-8, 1e-8); UnivariatePointValuePair pair = brentOptimizer.optimize(new MaxEval(200), new UnivariateObjectiveFunction(hb), org.apache.commons.math3.optim.nonlinear.scalar.GoalType.MINIMIZE, new SearchInterval(-3, 3)); param2 = new double[2]; hb.setIntercept(pair.getPoint()); hb.setScale(1.0); f = pair.getValue(); param2[0] = pair.getPoint(); param2[1] = 1; hb.setIntercept(param2[0]); hb.setScale(param2[1]); f = pair.getValue(); //Check Brent results against plink results. System.out.println(); System.out.println(" Brent Optimization"); System.out.println(" Iterations: "); System.out.println(" fmin: " + f); System.out.println(" B = " + hb.getIntercept() + " A = " + hb.getScale()); assertEquals(" Haebara intercept test", 0.06468396, hb.getIntercept(), 1e-4);//True results from sirt package in R assertEquals(" Haebara scale test", 1.0, hb.getScale(), 1e-4); //Compare UMNCMIN and Brent optimizer results. assertEquals(" UNCMIN/Brent intercept test", param1[0], param2[0], 1e-4); assertEquals(" UNCMIN/Brent slope test", param1[1], param2[1], 1e-4); }
From source file:com.itemanalysis.psychometrics.irt.equating.StockingLordMethodTest.java
/** * This test uses Rasch model items./*from w ww.j a va2 s . c om*/ * * */ @Test public void raschModelTest() { System.out.println("Stocking-Lord test: Rasch model"); LinkedHashMap<String, ItemResponseModel> itemFormX = new LinkedHashMap<String, ItemResponseModel>(); LinkedHashMap<String, ItemResponseModel> itemFormY = new LinkedHashMap<String, ItemResponseModel>(); itemFormX.put("Item1", new Irm3PL(-3.188047976, 1.0)); itemFormX.put("Item2", new Irm3PL(1.031760328, 1.0)); itemFormX.put("Item3", new Irm3PL(0.819040914, 1.0)); itemFormX.put("Item4", new Irm3PL(-2.706947360, 1.0)); itemFormX.put("Item5", new Irm3PL(-0.094527077, 1.0)); itemFormX.put("Item6", new Irm3PL(0.689697135, 1.0)); itemFormX.put("Item7", new Irm3PL(-0.551837153, 1.0)); itemFormX.put("Item8", new Irm3PL(-0.359559276, 1.0)); itemFormY.put("Item1", new Irm3PL(-3.074599226, 1.0)); itemFormY.put("Item2", new Irm3PL(1.012824350, 1.0)); itemFormY.put("Item3", new Irm3PL(0.868538408, 1.0)); itemFormY.put("Item4", new Irm3PL(-2.404483603, 1.0)); itemFormY.put("Item5", new Irm3PL(0.037402866, 1.0)); itemFormY.put("Item6", new Irm3PL(0.700747420, 1.0)); itemFormY.put("Item7", new Irm3PL(-0.602555046, 1.0)); itemFormY.put("Item8", new Irm3PL(-0.350426446, 1.0)); double f = 0; double[] param1 = new double[2]; double[] param2 = new double[2]; UniformDistributionApproximation distX = new UniformDistributionApproximation(-4.0, 4.0, 161);//plink default UniformDistributionApproximation distY = new UniformDistributionApproximation(-4.0, 4.0, 161);//plink default StockingLordMethod sl = new StockingLordMethod(itemFormX, itemFormY, distX, distY, EquatingCriterionType.Q1Q2); sl.setPrecision(6); double[] startValues = { 0 };//Using a start value for the intercept only for the Rasch family of models. DefaultUncminOptimizer optimizer = new DefaultUncminOptimizer(); try { optimizer.minimize(sl, startValues); double[] r = optimizer.getParameters(); param1[0] = r[0]; param1[1] = 1; f = optimizer.getFunctionValue(); sl.setIntercept(param1[0]); sl.setScale(param1[1]); } catch (UncminException ex) { ex.printStackTrace(); } //Check Brent Optimizer values against results from plink. System.out.println(" UNCMIN Optimization"); System.out.println(" Iterations: "); System.out.println(" fmin: " + f); System.out.println(" B = " + sl.getIntercept() + " A = " + sl.getScale()); assertEquals(" Intercept test", 0.057566, sl.getIntercept(), 1e-3);//true result from plink package in R assertEquals(" Scale test", 1.0, sl.getScale(), 1e-4); BrentOptimizer brentOptimizer = new BrentOptimizer(1e-8, 1e-8); UnivariatePointValuePair pair = brentOptimizer.optimize(new MaxEval(200), new UnivariateObjectiveFunction(sl), GoalType.MINIMIZE, new SearchInterval(-3, 3)); param2 = new double[2]; sl.setIntercept(pair.getPoint()); sl.setScale(1.0); f = pair.getValue(); param2[0] = pair.getPoint(); param2[1] = 1; sl.setIntercept(param2[0]); sl.setScale(param2[1]); f = pair.getValue(); //Check Brent results against plink results. System.out.println(); System.out.println(" Brent Optimization"); System.out.println(" Iterations: "); System.out.println(" fmin: " + f); System.out.println(" B = " + sl.getIntercept() + " A = " + sl.getScale()); assertEquals(" Intercept test", 0.057566, sl.getIntercept(), 1e-3);//true result from plink package in R assertEquals(" Scale test", 1.0, sl.getScale(), 1e-4); //Compare UMNCMIN and Brent optimizer results. assertEquals(" UNCMIN/Brent intercept test", param1[0], param2[0], 1e-6); assertEquals(" UNCMIN/Brent slope test", param1[1], param2[1], 1e-6); }
From source file:com.itemanalysis.psychometrics.irt.equating.HaebaraMethodTest.java
/** * This test uses a combination of 3PL and PCM items. Results compared to plink * * plink code:/*from w ww. java2 s.c om*/ * * library(plink) * * fX<-matrix(c( * 1, -3.188047976, 0,NA,NA, * 1, 1.031760328, 0,NA,NA, * 1, 0.819040914, 0,NA,NA, * 1, -2.706947360, 0,NA,NA, * 1, -0.094527077, 0,NA,NA, * 1, 0.689697135, 0,NA,NA, * 1, -0.551837153, 0,NA,NA, * 1, -0.359559276, 0,NA,NA, * 1, -1.451470831, -0.146619694, -0.636399040, 0.783018734), * nrow=9, byrow=TRUE) * fX<-as.data.frame(fX) * names(fX)<-c("aparam", "bparam","cparam","s1","s2") * * fY<-matrix(c( * 1,-3.074599226,0,NA,NA, * 1,1.01282435,0,NA,NA, * 1,0.868538408,0,NA,NA, * 1,-2.404483603,0,NA,NA, * 1,0.037402866,0,NA,NA, * 1,0.70074742,0,NA,NA, * 1,-0.602555046,0,NA,NA, * 1,-0.350426446,0,NA,NA, * 1,-1.267744832,-0.185885988,-0.61535623,0.801242218), * nrow=9, byrow=TRUE) * fY<-as.data.frame(fY) * names(fY)<-c("aparam", "bparam","cparam","s1","s2") * * common<-cbind(1:9, 1:9) * cat<-c(rep(2,8),4) * * pmX <- as.poly.mod(9,c("drm","gpcm"),list(1:8,9)) * pmY <- as.poly.mod(9,c("drm","gpcm"),list(1:8,9)) * * pars <- as.irt.pars(list(fx=fX,fy=fY), common, cat=list(fx=cat,fy=cat), * poly.mod=list(pmX,pmY), location=c(TRUE,TRUE)) * * plink(pars, startvals=c(1,0), rescale="SL", base.grp=2, D=1.0, symmetric=TRUE) * * */ @Test public void mixedFormtRaschPCMTest() { System.out.println("Mixed format Haebara test: Rasch and PCM"); LinkedHashMap<String, ItemResponseModel> itemFormX = new LinkedHashMap<String, ItemResponseModel>(); LinkedHashMap<String, ItemResponseModel> itemFormY = new LinkedHashMap<String, ItemResponseModel>(); itemFormX.put("Item1", new Irm3PL(-3.188047976, 1.0)); itemFormX.put("Item2", new Irm3PL(1.031760328, 1.0)); itemFormX.put("Item3", new Irm3PL(0.819040914, 1.0)); itemFormX.put("Item4", new Irm3PL(-2.706947360, 1.0)); itemFormX.put("Item5", new Irm3PL(-0.094527077, 1.0)); itemFormX.put("Item6", new Irm3PL(0.689697135, 1.0)); itemFormX.put("Item7", new Irm3PL(-0.551837153, 1.0)); itemFormX.put("Item8", new Irm3PL(-0.359559276, 1.0)); double[] step1x = { -0.146619694, -0.636399040, 0.783018734 }; itemFormX.put("Item9", new IrmPCM(-1.451470831, step1x, 1.0)); itemFormY.put("Item1", new Irm3PL(-3.074599226, 1.0)); itemFormY.put("Item2", new Irm3PL(1.012824350, 1.0)); itemFormY.put("Item3", new Irm3PL(0.868538408, 1.0)); itemFormY.put("Item4", new Irm3PL(-2.404483603, 1.0)); itemFormY.put("Item5", new Irm3PL(0.037402866, 1.0)); itemFormY.put("Item6", new Irm3PL(0.700747420, 1.0)); itemFormY.put("Item7", new Irm3PL(-0.602555046, 1.0)); itemFormY.put("Item8", new Irm3PL(-0.350426446, 1.0)); double[] step1y = { -0.185885988, -0.61535623, 0.801242218 }; itemFormY.put("Item9", new IrmPCM(-1.267744832, step1y, 1.0)); double f = 0; double[] param1 = new double[2]; double[] param2 = new double[2]; UniformDistributionApproximation distX = new UniformDistributionApproximation(-4.0, 4.0, 161);//plink default UniformDistributionApproximation distY = new UniformDistributionApproximation(-4.0, 4.0, 161);//plink default HaebaraMethod hb = new HaebaraMethod(itemFormX, itemFormY, distX, distY, EquatingCriterionType.Q1Q2); hb.setPrecision(6); double[] startValues = { 0.0 };//Using a start value for the intercept only for the Rasch family of models. DefaultUncminOptimizer optimizer = new DefaultUncminOptimizer(); try { optimizer.minimize(hb, startValues); double[] r = optimizer.getParameters(); param1[0] = r[0]; param1[1] = 1.0; f = optimizer.getFunctionValue(); hb.setIntercept(param1[0]); hb.setScale(param1[1]); } catch (UncminException ex) { ex.printStackTrace(); } //Check UNCMIN values against results from plink. System.out.println(" UNCMIN Optimization"); System.out.println(" Iterations: "); System.out.println(" fmin: " + f); System.out.println(" B = " + hb.getIntercept() + " A = " + hb.getScale()); //TODO there's something wrong with plink. I'm getting different results, but sirt agree with mine (at least for the Rasch model). // assertEquals(" Intercept test", 0.105526, hb.getIntercept(), 1e-4);//True results from plink package in R // assertEquals(" Scale test", 1.0, hb.getScale(), 1e-4); //Check Brent optimizer values against plink BrentOptimizer brentOptimizer = new BrentOptimizer(1e-8, 1e-8); UnivariatePointValuePair pair = brentOptimizer.optimize(new MaxEval(200), new UnivariateObjectiveFunction(hb), org.apache.commons.math3.optim.nonlinear.scalar.GoalType.MINIMIZE, new SearchInterval(-3, 3)); param2[0] = pair.getPoint(); param2[1] = 1.0; hb.setIntercept(param2[0]); hb.setScale(param2[1]); f = pair.getValue(); System.out.println(); System.out.println(" Brent Optimization"); System.out.println(" Iterations: "); System.out.println(" fmin: " + f); System.out.println(" B = " + hb.getIntercept() + " A = " + hb.getScale()); //Check UNCMIN values against results from plink. // assertEquals(" Intercept test", 0.105526, hb.getIntercept(), 1e-4); // assertEquals(" Scale test", 1.0, hb.getScale(), 1e-4); // // // //Compare UMNCMIN and Brent optimizer results. // assertEquals(" UNCMIN/Brent intercept test", param1[0], param2[0], 1e-6); // assertEquals(" UNCMIN/Brent slope test", param1[1], param2[1], 1e-6); }
From source file:com.itemanalysis.psychometrics.irt.equating.StockingLordMethodTest.java
/** * This test uses a combination of 3PL and PCM items. Results compared to plink * * plink code:/*from w w w.j ava 2s .c o m*/ * * library(plink) * * fX<-matrix(c( * 1, -3.188047976, 0,NA,NA, * 1, 1.031760328, 0,NA,NA, * 1, 0.819040914, 0,NA,NA, * 1, -2.706947360, 0,NA,NA, * 1, -0.094527077, 0,NA,NA, * 1, 0.689697135, 0,NA,NA, * 1, -0.551837153, 0,NA,NA, * 1, -0.359559276, 0,NA,NA, * 1, -1.451470831, -0.146619694, -0.636399040, 0.783018734), * nrow=9, byrow=TRUE) * fX<-as.data.frame(fX) * names(fX)<-c("aparam", "bparam","cparam","s1","s2") * * fY<-matrix(c( * 1,-3.074599226,0,NA,NA, * 1,1.01282435,0,NA,NA, * 1,0.868538408,0,NA,NA, * 1,-2.404483603,0,NA,NA, * 1,0.037402866,0,NA,NA, * 1,0.70074742,0,NA,NA, * 1,-0.602555046,0,NA,NA, * 1,-0.350426446,0,NA,NA, * 1,-1.267744832,-0.185885988,-0.61535623,0.801242218), * nrow=9, byrow=TRUE) * fY<-as.data.frame(fY) * names(fY)<-c("aparam", "bparam","cparam","s1","s2") * * common<-cbind(1:9, 1:9) * cat<-c(rep(2,8),4) * * pmX <- as.poly.mod(9,c("drm","gpcm"),list(1:8,9)) * pmY <- as.poly.mod(9,c("drm","gpcm"),list(1:8,9)) * * pars <- as.irt.pars(list(fx=fX,fy=fY), common, cat=list(fx=cat,fy=cat), * poly.mod=list(pmX,pmY), location=c(TRUE,TRUE)) * * plink(pars, startvals=c(1,0), rescale="SL", base.grp=2, D=1.0, symmetric=TRUE) * * */ @Test public void mixedFormtRaschPCMTest() { System.out.println("Mixed format Stocking-Lord test: Rasch and PCM"); LinkedHashMap<String, ItemResponseModel> itemFormX = new LinkedHashMap<String, ItemResponseModel>(); LinkedHashMap<String, ItemResponseModel> itemFormY = new LinkedHashMap<String, ItemResponseModel>(); itemFormX.put("Item1", new Irm3PL(-3.188047976, 1.0)); itemFormX.put("Item2", new Irm3PL(1.031760328, 1.0)); itemFormX.put("Item3", new Irm3PL(0.819040914, 1.0)); itemFormX.put("Item4", new Irm3PL(-2.706947360, 1.0)); itemFormX.put("Item5", new Irm3PL(-0.094527077, 1.0)); itemFormX.put("Item6", new Irm3PL(0.689697135, 1.0)); itemFormX.put("Item7", new Irm3PL(-0.551837153, 1.0)); itemFormX.put("Item8", new Irm3PL(-0.359559276, 1.0)); double[] step1x = { -0.146619694, -0.636399040, 0.783018734 }; itemFormX.put("Item9", new IrmPCM(-1.451470831, step1x, 1.0)); itemFormY.put("Item1", new Irm3PL(-3.074599226, 1.0)); itemFormY.put("Item2", new Irm3PL(1.012824350, 1.0)); itemFormY.put("Item3", new Irm3PL(0.868538408, 1.0)); itemFormY.put("Item4", new Irm3PL(-2.404483603, 1.0)); itemFormY.put("Item5", new Irm3PL(0.037402866, 1.0)); itemFormY.put("Item6", new Irm3PL(0.700747420, 1.0)); itemFormY.put("Item7", new Irm3PL(-0.602555046, 1.0)); itemFormY.put("Item8", new Irm3PL(-0.350426446, 1.0)); double[] step1y = { -0.185885988, -0.61535623, 0.801242218 }; itemFormY.put("Item9", new IrmPCM(-1.267744832, step1y, 1.0)); double f = 0; double[] param1 = new double[2]; double[] param2 = new double[2]; UniformDistributionApproximation distX = new UniformDistributionApproximation(-4.0, 4.0, 161);//plink default UniformDistributionApproximation distY = new UniformDistributionApproximation(-4.0, 4.0, 161);//plink default StockingLordMethod sl = new StockingLordMethod(itemFormX, itemFormY, distX, distY, EquatingCriterionType.Q1Q2); sl.setPrecision(6); double[] startValues = { 0.0 };//Using a start value for the intercept only for the Rasch family of models. DefaultUncminOptimizer optimizer = new DefaultUncminOptimizer(); try { optimizer.minimize(sl, startValues); double[] r = optimizer.getParameters(); param1[0] = r[0]; param1[1] = 1.0; f = optimizer.getFunctionValue(); sl.setIntercept(param1[0]); sl.setScale(param1[1]); } catch (UncminException ex) { ex.printStackTrace(); } //Check UNCMIN values against results from plink. System.out.println(" UNCMIN Optimization"); System.out.println(" Iterations: "); System.out.println(" fmin: " + f); System.out.println(" B = " + sl.getIntercept() + " A = " + sl.getScale()); assertEquals(" Intercept test", 0.101388, sl.getIntercept(), 1e-4);//True results from plink package in R assertEquals(" Scale test", 1.0, sl.getScale(), 1e-4); //Check Brent optimizer values against plink BrentOptimizer brentOptimizer = new BrentOptimizer(1e-8, 1e-8); UnivariatePointValuePair pair = brentOptimizer.optimize(new MaxEval(200), new UnivariateObjectiveFunction(sl), GoalType.MINIMIZE, new SearchInterval(-3, 3)); param2[0] = pair.getPoint(); param2[1] = 1.0; sl.setIntercept(param2[0]); sl.setScale(param2[1]); f = pair.getValue(); System.out.println(); System.out.println(" Brent Optimization"); System.out.println(" Iterations: "); System.out.println(" fmin: " + f); System.out.println(" B = " + sl.getIntercept() + " A = " + sl.getScale()); //Check UNCMIN values against results from plink. assertEquals(" Intercept test", 0.101388, sl.getIntercept(), 1e-4); assertEquals(" Scale test", 1.0, sl.getScale(), 1e-4); //Compare UMNCMIN and Brent optimizer results. assertEquals(" UNCMIN/Brent intercept test", param1[0], param2[0], 1e-6); assertEquals(" UNCMIN/Brent slope test", param1[1], param2[1], 1e-6); }