List of usage examples for org.apache.commons.math3.linear BlockRealMatrix setRowVector
@Override public void setRowVector(final int row, final RealVector vector) throws OutOfRangeException, MatrixDimensionMismatchException
From source file:gamlss.utilities.MatrixFunctions.java
/** * Append rows of the matrices.// w w w . j av a 2s . c om * @param m1 - first matrix * @param m2 - second matrix * @return m1.append.m2 */ public static BlockRealMatrix appendMatricesRows(final BlockRealMatrix m1, final BlockRealMatrix m2) { BlockRealMatrix out = new BlockRealMatrix(m1.getRowDimension(), m1.getColumnDimension() + m2.getColumnDimension()); for (int i = 0; i < m1.getRowDimension(); i++) { out.setRowVector(i, m1.getRowVector(i).append(m2.getRowVector(i))); } return out; }
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 ww .java 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); } }