List of usage examples for org.apache.commons.math3.analysis.solvers BrentSolver BrentSolver
public BrentSolver(double relativeAccuracy, double absoluteAccuracy)
From source file:jat.examples.CommonsMath.RootFinder.java
/** * @param args//from ww w . j a v a2 s. c o m */ public static void main(String[] args) { final double relativeAccuracy = 1.0e-12; final double absoluteAccuracy = 1.0e-8; UnivariateFunction function = new MyFunction(); UnivariateSolver nonBracketing = new BrentSolver(relativeAccuracy, absoluteAccuracy); // double baseRoot = bs.solve(100, function, -2.0, 0); System.out.println("Roots of f(x)=x^2-2: "); double baseRoot; baseRoot = nonBracketing.solve(100, function, -2.0, 0); System.out.println("root1: " + baseRoot); baseRoot = nonBracketing.solve(100, function, 0, 2.0); System.out.println("root2: " + baseRoot); // System.out.println("rs: " + Math.sqrt(2.)); }
From source file:gamlss.distributions.GT.java
/** * This is the Skew t type3 distribution with supplied link * function for each of the distribution parameters. * @param muLink - link function for mu distribution parameter * @param sigmaLink - link function for sigma distribution parameter * @param nuLink - link function for nu distribution parameter * @param tauLink - link function for tau distribution parameter *//* w ww. j a v a2 s. co m*/ public GT(final int muLink, final int sigmaLink, final int nuLink, final int tauLink) { distributionParameterLink.put(DistributionSettings.MU, MakeLinkFunction.checkLink(DistributionSettings.GT, muLink)); distributionParameterLink.put(DistributionSettings.SIGMA, MakeLinkFunction.checkLink(DistributionSettings.GT, sigmaLink)); distributionParameterLink.put(DistributionSettings.NU, MakeLinkFunction.checkLink(DistributionSettings.GT, nuLink)); distributionParameterLink.put(DistributionSettings.TAU, MakeLinkFunction.checkLink(DistributionSettings.GT, tauLink)); integrator = new RombergIntegrator(); function = new IntegratingFunction(); interval = new double[2]; uniRootSolver = new BrentSolver(1.0e-12, 1.0e-8); uniRootObj = new UniRootObjFunction(); }
From source file:gamlss.distributions.ST1.java
/** This is the Skew t (Azzalini type 1) distribution with * supplied link function for each of the distribution parameters. * @param muLink - link function for mu distribution parameter * @param sigmaLink - link function for sigma distribution parameter * @param nuLink - link function for nu distribution parameter * @param tauLink - link function for tau distribution parameter *///from w w w . j a v a2s . c o m public ST1(final int muLink, final int sigmaLink, final int nuLink, final int tauLink) { distributionParameterLink.put(DistributionSettings.MU, MakeLinkFunction.checkLink(DistributionSettings.ST1, muLink)); distributionParameterLink.put(DistributionSettings.SIGMA, MakeLinkFunction.checkLink(DistributionSettings.ST1, sigmaLink)); distributionParameterLink.put(DistributionSettings.NU, MakeLinkFunction.checkLink(DistributionSettings.ST1, nuLink)); distributionParameterLink.put(DistributionSettings.TAU, MakeLinkFunction.checkLink(DistributionSettings.ST1, tauLink)); tdDist = new TDistr(); noDist = new NormalDistribution(); integrator = new LegendreGaussIntegrator(2, LegendreGaussIntegrator.DEFAULT_RELATIVE_ACCURACY, LegendreGaussIntegrator.DEFAULT_ABSOLUTE_ACCURACY); function = new IntegratingFunction(); interval = new double[2]; uniRootSolver = new BrentSolver(1.0e-12, 1.0e-8); uniRootObj = new UniRootObjFunction(); }
From source file:gamlss.smoothing.PB.java
/** * The main fitting method, initiate appropriate smoothing * function according to incoming parameters. * @param additiveFit -object of AdditiveFoit class * @return reiduals/*from w w w.j a v a 2s .c o m*/ */ //gamlss.pb <- function(x, y, w, xeval = NULL, ...) public ArrayRealVector solve(final AdditiveFit additiveFit) { Double lambda = Controls.LAMBDAS_USER_DEFINED; colNum = additiveFit.getColNum(); y = additiveFit.getZ(); w = additiveFit.getW(); whichDistParameter = additiveFit.getWhichDistParameter(); basisM = (BlockRealMatrix) getBasisM().get(colNum); penaltyM = (BlockRealMatrix) getPenaltyM().get(colNum); df = (Integer) getDfValues().get(colNum); //n <- nrow(X) # the no of observations n = basisM.getRowDimension(); //lambdaS <- get(startLambdaName, envir=gamlss.env) //geting the starting value lambdaS = getLambdas().get(whichDistParameter); //if (lambdaS>=1e+07) lambda <- 1e+07 if (lambdaS >= 1e+07) { lambda = 1e+07; } //if (lambdaS<=1e-07) lambda <- 1e-07 if (lambdaS <= 1e-07) { lambda = 1e-07; } //if (is.null(df)&&!is.null(lambda)||!is.null(df)&&!is.null(lambda)) if (lambda != null) { //fit <- regpen(y, X, w, lambda, D) beta = regpen(lambda); //fv <- X %*% fit$beta fv = (ArrayRealVector) basisM.operate(beta); //else if (is.null(df)&&is.null(lambda)) } else if (df == null) { //lambda <- lambdaS lambda = lambdaS; //switch(control$c," ML"={ switch (Controls.SMOOTH_METHOD) { case DistributionSettings.GAIC: fv = functionGAIC(lambda); break; case DistributionSettings.ML: fv = functionML(lambda); break; case DistributionSettings.ML1: fv = functionML1(lambda); break; case DistributionSettings.EM: System.err.println("EM has not been implemented yet"); break; case DistributionSettings.GCV: fv = functionGCV(lambda); break; default: System.err.println( " Cannot recognise the " + "smoothing method or it has" + "not been implemented yet"); } } else { //QR <- qr(sqrt(w)*X) //Rinv <- solve(qr.R(QR)) rM = (BlockRealMatrix) new QRDecomposition( MatrixFunctions.multVectorMatrix(MatrixFunctions.sqrtVec(w), basisM)).getR(); rM = rM.getSubMatrix(0, rM.getColumnDimension() - 1, 0, rM.getColumnDimension() - 1); rMinv = (BlockRealMatrix) new QRDecomposition(rM).getSolver().getInverse(); //S <- t(D)%*%D sM = penaltyM.transpose().multiply(penaltyM); //UDU <- eigen(t(Rinv)%*%S%*%Rinv) uduV = new EigenDecomposition(rMinv.transpose().multiply(sM).multiply(rMinv)).getRealEigenvalues(); //lambda <- if (sign(edf1_df(0))==sign(edf1_df(100000))) 100000 //in case they have the some sign edfTemp1 = edf1_df(0, df); edfTemp2 = edf1_df(100000.0, df); if (FastMath.signum(edfTemp1) == FastMath.signum(edfTemp2)) { lambda = 100000.0; } else { //else uniroot(edf1_df, c(0,100000))$root uniRootObj.setDf(df); final double relativeAccuracy = 1.0e-12; final double absoluteAccuracy = 1.0e-8; UnivariateSolver uniRootSolver = new BrentSolver(relativeAccuracy, absoluteAccuracy); lambda = uniRootSolver.solve(1000, uniRootObj, 0.0, 100000.0); } //fit <- regpen(y, X, w, lambda, D) beta = regpen(lambda); fv = (ArrayRealVector) basisM.operate(beta); } if (Controls.IF_NEED_THIS) { //waug <- as.vector(c(w, rep(1,nrow(D)))) waug = w.append(MatrixFunctions.repV(1.0, penaltyM.getRowDimension())); //xaug <- as.matrix(rbind(X,sqrt(lambda)*D)) xaug = MatrixFunctions.appendMatricesRows(basisM, (BlockRealMatrix) penaltyM.scalarMultiply(FastMath.sqrt(lambda))); //lev <- hat(sqrt(waug)*xaug,intercept=FALSE)[1:n] //get the hat matrix lev = (ArrayRealVector) MatrixFunctions.getMainDiagonal(new BlockRealMatrix(MatrixFunctions .calculateHat(MatrixFunctions.multVectorMatrix(MatrixFunctions.sqrtVec(waug), xaug)).getData())) .getSubVector(0, n); //lev <- (lev-.hat.WX(w,x)) //subtract the linear since is already fitted lev = lev.subtract( hatWX((ArrayRealVector) getSmootherMatrices().get(whichDistParameter).getColumnVector(colNum))); // var <- lev/w //the variance of the smootherz var = lev.ebeDivide(w); } //if (is.null(xeval)) # if no prediction if (Controls.XEVAL_USER_DEFINED == null) { // Residuals return y.subtract(fv); //else # for prediction } else { //ll <- dim(as.matrix(attr(x,"X")))[1] //nx <- as.matrix(attr(x,"X"))[seq(length(y)+1,ll),] tempArr = ArithmeticSeries.getSeries(y.getDimension() + 1, basisM.getRowDimension(), 1); BlockRealMatrix nx = new BlockRealMatrix(tempArr.length, basisM.getColumnDimension()); for (int i = 0; i < tempArr.length; i++) { nx.setRowVector(i, basisM.getRowVector((int) tempArr[i])); } //pred <- drop(nx %*% fit$beta) ArrayRealVector pred = (ArrayRealVector) nx.operate(beta); // Residuals return y.subtract(fv); } }