List of usage examples for org.apache.commons.math3.linear BlockRealMatrix copy
@Override
public BlockRealMatrix copy()
From source file:gamlss.algorithm.AdditiveFit.java
/** * Performs a smoothing process - creates an approximating function * that attempts to capture important patterns in the data, while * leaving out noise or other fine-scale structures/rapid phenomena. * @param xMat - design matrix/*from w ww . j av a 2 s .co m*/ * @param y - response variable * @param sMat - smoother matrix * @param distParameter - distribution parameter * @return matrix of smoothed values */ public BlockRealMatrix fitSmoother(final ArrayRealVector y, final ArrayRealVector weights, final BlockRealMatrix xMat, final BlockRealMatrix sMat, final int distParameter) { this.s = sMat.copy(); this.w = weights; this.whichDistParameter = distParameter; //residuals <- as.vector(y - s %*% rep(1, ncol(s))) tempV = new ArrayRealVector(s.getColumnDimension()); tempV.set(1.0); tempV = (ArrayRealVector) s.operate(tempV); residuals = MatrixFunctions.vecSub(y, tempV); tempV = null; //fit <- list(fitted.values = 0) fittedValues = new ArrayRealVector(residuals.getDimension()); //rss <- weighted.mean(residuals^2, w) //rss = calculateRss(residuals, w); //tempArr = null; //rssold <- rss * 10 //rssOld = rss*10; //nit <- 0 nit = 0; //df <- rep(NA, length(who)) //lambda <- rep(NA, length(who)) //double[] df = new double[s.getColumnDimension()]; //var <- s var = s.copy(); //if(trace) cat("\nADDITIVE iter rss/n term\n") if (Controls.BACKFIT_TRACE) { System.out.println("ADDITIVE iter rss/n term"); } //ndig <- -log10(tol) + 1 //double ndig = -FastMath.log10(tol)+1; //RATIO <- tol + 1 ratio = Controls.BACKFIT_TOL + 1; //while(RATIO > tol & nit < maxit) while (ratio > Controls.BACKFIT_TOL & nit < Controls.BACKFIT_CYCLES) { //rssold <- rss //rssOld = rss; //nit <- nit + 1 nit++; //z <- residuals + fit$fitted.values z = residuals.add(fittedValues); //org.apache.commons.lang3.time.StopWatch watch = new org.apache.commons.lang3.time.StopWatch(); //watch.start(); //fit <- lm.wfit(x, z, w, method="qr", ...) wls.setNoIntercept(Controls.NO_INTERCEPT[whichDistParameter - 1]); wls.newSampleData(z, xMat.copy(), w.copy()); wls.setNoIntercept(false); //residuals <- fit$residuals fittedValues = (ArrayRealVector) wls.calculateFittedValues(Controls.IS_SVD); //watch.stop(); //System.out.println(watch.getNanoTime()/(1e+9)); //watch.reset(); //residuals = z.subtract(fittedValues); //residuals = (ArrayRealVector) wls.calculateResiduals(); residuals = z.subtract(fittedValues); //rss <- weighted.mean(residuals^2, w) //rss = calculateRss(residuals, w); //if(trace) cat(nit, " ", //format(round(rss, ndig)), " Parametric -- lm.wfit\n", sep = "") if (Controls.BACKFIT_TRACE) { //System.out.println(" " + nit + " " + //rss + " Parametric -- lm.wfit"); } //deltaf <- 0 deltaf = 0; //for(j in seq(names.calls)) for (colNum = 0; colNum < s.getColumnDimension(); colNum++) { //old <- s[, j] old = (ArrayRealVector) s.getColumnVector(colNum); //z <- residuals + s[, j] z = residuals.add(old); //fit.call <- eval(smooth.calls[[j]]) //residuals <- as.double(fit.call$residuals) residuals = smoother.solve(this); //if(length(residuals) != n) //stop(paste(names.calls[j], "returns a vector //of the wrong length")) if (residuals.getDimension() != y.getDimension()) { System.err.println(colNum + " column of matrix s has a" + " vector of the wrong length"); } //s[, j] <- z - residual s.setColumnVector(colNum, z.subtract(residuals)); //if (length(fit.call$lambda)>1) //{for cases where there are multiple lambdas //ambda[j] <- fit.call$lambda[1] //coefSmo[[j]] <- if(is.null(fit.call$coefSmo)) //0 else fit.call$coefSmo //deltaf <- deltaf + weighted.mean((s[, j] - old)^2, w) tempV = MatrixFunctions.vecSub((ArrayRealVector) s.getColumnVector(colNum), old); deltaf = deltaf + mean.evaluate(tempV.ebeMultiply(tempV).getDataRef(), w.getDataRef()); tempV = null; //rss <- weighted.mean(residuals^2, w) //rss = calculateRss(residuals, w); // if(trace) if (Controls.BACKFIT_TRACE) { //cat(" ", nit, " ", format(round(rss, ndig)), //" Nonparametric -- ", names.calls[j], "\n", sep = "") //System.out.println(" " + nit +" " + rss + " //Nonparametric " + "pb(column "+ i+")"); } //df[j] <- fit.call$nl.df //df[i] = pb.getEdf(); //if(se) if (Controls.IS_SE) { //var[, j] <- fit.call$var var.setColumnVector(colNum, smoother.getVar()); } } //RATIO <- sqrt(deltaf/sum(w * apply(s, 1, sum)^2)) tempD = 0.0; tempM = new BlockRealMatrix(s.getRowDimension(), 1); for (int j = 0; j < s.getRowDimension(); j++) { for (int k = 0; k < s.getColumnDimension(); k++) { tempD = tempD + s.getEntry(j, k); } tempM.setEntry(j, 0, tempD); tempD = 0.0; } size = tempM.getRowDimension(); for (int j = 0; j < size; j++) { tempD = tempD + tempM.getEntry(j, 0) * tempM.getEntry(j, 0) * w.getEntry(j); } ratio = FastMath.sqrt(deltaf / tempD); tempD = null; tempM = null; //if(trace) //cat("Relative change in functions:", //format(round(RATIO, ndig)), "\n") if (Controls.BACKFIT_TRACE) { System.out.println("Relative change in functions: " + ratio); } } //if(nit == maxit && maxit > 1) //warning(paste("additive.fit convergence not //obtained in ", maxit, " iterations")) if (ratio > Controls.BACKFIT_TOL) { System.out.println( "AdditiveFit convergence is not obtained in " + Controls.BACKFIT_CYCLES + " iterations"); } //fit$fitted.values <- y - residuals fittedValues = y.subtract(residuals); //rl <- c(fit, list(smooth = s, nl.df = //sum(df), lambda=lambda, coefSmo=coefSmo)) return s; }