MillerUpdatingRegression.java
- /*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
- package org.apache.commons.math4.legacy.stat.regression;
- import java.util.Arrays;
- import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
- import org.apache.commons.math4.core.jdkmath.JdkMath;
- import org.apache.commons.numbers.core.Precision;
- /**
- * This class is a concrete implementation of the {@link UpdatingMultipleLinearRegression} interface.
- *
- * <p>The algorithm is described in: <pre>
- * Algorithm AS 274: Least Squares Routines to Supplement Those of Gentleman
- * Author(s): Alan J. Miller
- * Source: Journal of the Royal Statistical Society.
- * Series C (Applied Statistics), Vol. 41, No. 2
- * (1992), pp. 458-478
- * Published by: Blackwell Publishing for the Royal Statistical Society
- * Stable URL: http://www.jstor.org/stable/2347583 </pre>
- *
- * <p>This method for multiple regression forms the solution to the OLS problem
- * by updating the QR decomposition as described by Gentleman.</p>
- *
- * @since 3.0
- */
- public class MillerUpdatingRegression implements UpdatingMultipleLinearRegression {
- /** number of variables in regression. */
- private final int nvars;
- /** diagonals of cross products matrix. */
- private final double[] d;
- /** the elements of the R`Y. */
- private final double[] rhs;
- /** the off diagonal portion of the R matrix. */
- private final double[] r;
- /** the tolerance for each of the variables. */
- private final double[] tol;
- /** residual sum of squares for all nested regressions. */
- private final double[] rss;
- /** order of the regressors. */
- private final int[] vorder;
- /** scratch space for tolerance calc. */
- private final double[] workTolset;
- /** number of observations entered. */
- private long nobs;
- /** sum of squared errors of largest regression. */
- private double sserr;
- /** has rss been called? */
- private boolean rssSet;
- /** has the tolerance setting method been called. */
- private boolean tolSet;
- /** flags for variables with linear dependency problems. */
- private final boolean[] lindep;
- /** singular x values. */
- private final double[] xSing;
- /** workspace for singularity method. */
- private final double[] workSing;
- /** summation of Y variable. */
- private double sumy;
- /** summation of squared Y values. */
- private double sumsqy;
- /** boolean flag whether a regression constant is added. */
- private final boolean hasIntercept;
- /** zero tolerance. */
- private final double epsilon;
- /**
- * Set the default constructor to private access.
- * to prevent inadvertent instantiation
- */
- @SuppressWarnings("unused")
- private MillerUpdatingRegression() {
- this(-1, false, Double.NaN);
- }
- /**
- * This is the augmented constructor for the MillerUpdatingRegression class.
- *
- * @param numberOfVariables number of regressors to expect, not including constant
- * @param includeConstant include a constant automatically
- * @param errorTolerance zero tolerance, how machine zero is determined
- * @throws ModelSpecificationException if {@code numberOfVariables is less than 1}
- */
- public MillerUpdatingRegression(int numberOfVariables, boolean includeConstant, double errorTolerance)
- throws ModelSpecificationException {
- if (numberOfVariables < 1) {
- throw new ModelSpecificationException(LocalizedFormats.NO_REGRESSORS);
- }
- if (includeConstant) {
- this.nvars = numberOfVariables + 1;
- } else {
- this.nvars = numberOfVariables;
- }
- this.hasIntercept = includeConstant;
- this.nobs = 0;
- this.d = new double[this.nvars];
- this.rhs = new double[this.nvars];
- this.r = new double[this.nvars * (this.nvars - 1) / 2];
- this.tol = new double[this.nvars];
- this.rss = new double[this.nvars];
- this.vorder = new int[this.nvars];
- this.xSing = new double[this.nvars];
- this.workSing = new double[this.nvars];
- this.workTolset = new double[this.nvars];
- this.lindep = new boolean[this.nvars];
- for (int i = 0; i < this.nvars; i++) {
- vorder[i] = i;
- }
- if (errorTolerance > 0) {
- this.epsilon = errorTolerance;
- } else {
- this.epsilon = -errorTolerance;
- }
- }
- /**
- * Primary constructor for the MillerUpdatingRegression.
- *
- * @param numberOfVariables maximum number of potential regressors
- * @param includeConstant include a constant automatically
- * @throws ModelSpecificationException if {@code numberOfVariables is less than 1}
- */
- public MillerUpdatingRegression(int numberOfVariables, boolean includeConstant)
- throws ModelSpecificationException {
- this(numberOfVariables, includeConstant, Precision.EPSILON);
- }
- /**
- * A getter method which determines whether a constant is included.
- * @return true regression has an intercept, false no intercept
- */
- @Override
- public boolean hasIntercept() {
- return this.hasIntercept;
- }
- /**
- * Gets the number of observations added to the regression model.
- * @return number of observations
- */
- @Override
- public long getN() {
- return this.nobs;
- }
- /**
- * Adds an observation to the regression model.
- * @param x the array with regressor values
- * @param y the value of dependent variable given these regressors
- * @exception ModelSpecificationException if the length of {@code x} does not equal
- * the number of independent variables in the model
- */
- @Override
- public void addObservation(final double[] x, final double y)
- throws ModelSpecificationException {
- if ((!this.hasIntercept && x.length != nvars) ||
- (this.hasIntercept && x.length + 1 != nvars)) {
- throw new ModelSpecificationException(LocalizedFormats.INVALID_REGRESSION_OBSERVATION,
- x.length, nvars);
- }
- if (!this.hasIntercept) {
- include(Arrays.copyOf(x, x.length), 1.0, y);
- } else {
- final double[] tmp = new double[x.length + 1];
- System.arraycopy(x, 0, tmp, 1, x.length);
- tmp[0] = 1.0;
- include(tmp, 1.0, y);
- }
- ++nobs;
- }
- /**
- * Adds multiple observations to the model.
- * @param x observations on the regressors
- * @param y observations on the regressand
- * @throws ModelSpecificationException if {@code x} is not rectangular, does not match
- * the length of {@code y} or does not contain sufficient data to estimate the model
- */
- @Override
- public void addObservations(double[][] x, double[] y) throws ModelSpecificationException {
- if (x == null || y == null || x.length != y.length) {
- throw new ModelSpecificationException(
- LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
- (x == null) ? 0 : x.length,
- (y == null) ? 0 : y.length);
- }
- if (x.length == 0) { // Must be no y data either
- throw new ModelSpecificationException(
- LocalizedFormats.NO_DATA);
- }
- if (x[0].length + 1 > x.length) {
- throw new ModelSpecificationException(
- LocalizedFormats.NOT_ENOUGH_DATA_FOR_NUMBER_OF_PREDICTORS,
- x.length, x[0].length);
- }
- for (int i = 0; i < x.length; i++) {
- addObservation(x[i], y[i]);
- }
- }
- /**
- * The include method is where the QR decomposition occurs. This statement forms all
- * intermediate data which will be used for all derivative measures.
- * According to the miller paper, note that in the original implementation the x vector
- * is overwritten. In this implementation, the include method is passed a copy of the
- * original data vector so that there is no contamination of the data. Additionally,
- * this method differs slightly from Gentleman's method, in that the assumption is
- * of dense design matrices, there is some advantage in using the original gentleman algorithm
- * on sparse matrices.
- *
- * @param x observations on the regressors
- * @param wi weight of the this observation (-1,1)
- * @param yi observation on the regressand
- */
- private void include(final double[] x, final double wi, final double yi) {
- int nextr = 0;
- double w = wi;
- double y = yi;
- double xi;
- double di;
- double wxi;
- double dpi;
- double xk;
- double wPrev;
- this.rssSet = false;
- sumy = smartAdd(yi, sumy);
- sumsqy = smartAdd(sumsqy, yi * yi);
- for (int i = 0; i < x.length; i++) {
- if (w == 0.0) {
- return;
- }
- xi = x[i];
- if (xi == 0.0) {
- nextr += nvars - i - 1;
- continue;
- }
- di = d[i];
- wxi = w * xi;
- wPrev = w;
- if (di != 0.0) {
- dpi = smartAdd(di, wxi * xi);
- final double tmp = wxi * xi / di;
- if (JdkMath.abs(tmp) > Precision.EPSILON) {
- w = (di * w) / dpi;
- }
- } else {
- dpi = wxi * xi;
- w = 0.0;
- }
- d[i] = dpi;
- for (int k = i + 1; k < nvars; k++) {
- xk = x[k];
- x[k] = smartAdd(xk, -xi * r[nextr]);
- if (di != 0.0) {
- r[nextr] = smartAdd(di * r[nextr], (wPrev * xi) * xk) / dpi;
- } else {
- r[nextr] = xk / xi;
- }
- ++nextr;
- }
- xk = y;
- y = smartAdd(xk, -xi * rhs[i]);
- if (di != 0.0) {
- rhs[i] = smartAdd(di * rhs[i], wxi * xk) / dpi;
- } else {
- rhs[i] = xk / xi;
- }
- }
- sserr = smartAdd(sserr, w * y * y);
- }
- /**
- * Adds to number a and b such that the contamination due to
- * numerical smallness of one addend does not corrupt the sum.
- * @param a - an addend
- * @param b - an addend
- * @return the sum of the a and b
- */
- private double smartAdd(double a, double b) {
- final double aa = JdkMath.abs(a);
- final double ba = JdkMath.abs(b);
- if (aa > ba) {
- final double eps = aa * Precision.EPSILON;
- if (ba > eps) {
- return a + b;
- }
- return a;
- } else {
- final double eps = ba * Precision.EPSILON;
- if (aa > eps) {
- return a + b;
- }
- return b;
- }
- }
- /**
- * As the name suggests, clear wipes the internals and reorders everything in the
- * canonical order.
- */
- @Override
- public void clear() {
- Arrays.fill(this.d, 0.0);
- Arrays.fill(this.rhs, 0.0);
- Arrays.fill(this.r, 0.0);
- Arrays.fill(this.tol, 0.0);
- Arrays.fill(this.rss, 0.0);
- Arrays.fill(this.workTolset, 0.0);
- Arrays.fill(this.workSing, 0.0);
- Arrays.fill(this.xSing, 0.0);
- Arrays.fill(this.lindep, false);
- for (int i = 0; i < nvars; i++) {
- this.vorder[i] = i;
- }
- this.nobs = 0;
- this.sserr = 0.0;
- this.sumy = 0.0;
- this.sumsqy = 0.0;
- this.rssSet = false;
- this.tolSet = false;
- }
- /**
- * This sets up tolerances for singularity testing.
- */
- private void tolset() {
- int pos;
- double total;
- final double eps = this.epsilon;
- for (int i = 0; i < nvars; i++) {
- this.workTolset[i] = JdkMath.sqrt(d[i]);
- }
- tol[0] = eps * this.workTolset[0];
- for (int col = 1; col < nvars; col++) {
- pos = col - 1;
- total = workTolset[col];
- for (int row = 0; row < col; row++) {
- total += JdkMath.abs(r[pos]) * workTolset[row];
- pos += nvars - row - 2;
- }
- tol[col] = eps * total;
- }
- tolSet = true;
- }
- /**
- * The regcf method conducts the linear regression and extracts the
- * parameter vector. Notice that the algorithm can do subset regression
- * with no alteration.
- *
- * @param nreq how many of the regressors to include (either in canonical
- * order, or in the current reordered state)
- * @return an array with the estimated slope coefficients
- * @throws ModelSpecificationException if {@code nreq} is less than 1
- * or greater than the number of independent variables
- */
- private double[] regcf(int nreq) throws ModelSpecificationException {
- int nextr;
- if (nreq < 1) {
- throw new ModelSpecificationException(LocalizedFormats.NO_REGRESSORS);
- }
- if (nreq > this.nvars) {
- throw new ModelSpecificationException(
- LocalizedFormats.TOO_MANY_REGRESSORS, nreq, this.nvars);
- }
- if (!this.tolSet) {
- tolset();
- }
- final double[] ret = new double[nreq];
- boolean rankProblem = false;
- for (int i = nreq - 1; i > -1; i--) {
- if (JdkMath.sqrt(d[i]) < tol[i]) {
- ret[i] = 0.0;
- d[i] = 0.0;
- rankProblem = true;
- } else {
- ret[i] = rhs[i];
- nextr = i * (nvars + nvars - i - 1) / 2;
- for (int j = i + 1; j < nreq; j++) {
- ret[i] = smartAdd(ret[i], -r[nextr] * ret[j]);
- ++nextr;
- }
- }
- }
- if (rankProblem) {
- for (int i = 0; i < nreq; i++) {
- if (this.lindep[i]) {
- ret[i] = Double.NaN;
- }
- }
- }
- return ret;
- }
- /**
- * The method which checks for singularities and then eliminates the offending
- * columns.
- */
- private void singcheck() {
- int pos;
- for (int i = 0; i < nvars; i++) {
- workSing[i] = JdkMath.sqrt(d[i]);
- }
- for (int col = 0; col < nvars; col++) {
- // Set elements within R to zero if they are less than tol(col) in
- // absolute value after being scaled by the square root of their row
- // multiplier
- final double temp = tol[col];
- pos = col - 1;
- for (int row = 0; row < col - 1; row++) {
- if (JdkMath.abs(r[pos]) * workSing[row] < temp) {
- r[pos] = 0.0;
- }
- pos += nvars - row - 2;
- }
- // If diagonal element is near zero, set it to zero, set appropriate
- // element of LINDEP, and use INCLUD to augment the projections in
- // the lower rows of the orthogonalization.
- lindep[col] = false;
- if (workSing[col] < temp) {
- lindep[col] = true;
- if (col < nvars - 1) {
- Arrays.fill(xSing, 0.0);
- int pi = col * (nvars + nvars - col - 1) / 2;
- for (int xi = col + 1; xi < nvars; xi++, pi++) {
- xSing[xi] = r[pi];
- r[pi] = 0.0;
- }
- final double y = rhs[col];
- final double weight = d[col];
- d[col] = 0.0;
- rhs[col] = 0.0;
- this.include(xSing, weight, y);
- } else {
- sserr += d[col] * rhs[col] * rhs[col];
- }
- }
- }
- }
- /**
- * Calculates the sum of squared errors for the full regression.
- * and all subsets in the following manner: <pre>
- * rss[] ={
- * ResidualSumOfSquares_allNvars,
- * ResidualSumOfSquares_FirstNvars-1,
- * ResidualSumOfSquares_FirstNvars-2,
- * ..., ResidualSumOfSquares_FirstVariable} </pre>
- */
- private void ss() {
- double total = sserr;
- rss[nvars - 1] = sserr;
- for (int i = nvars - 1; i > 0; i--) {
- total += d[i] * rhs[i] * rhs[i];
- rss[i - 1] = total;
- }
- rssSet = true;
- }
- /**
- * Calculates the cov matrix assuming only the first nreq variables are
- * included in the calculation. The returned array contains a symmetric
- * matrix stored in lower triangular form. The matrix will have
- * ( nreq + 1 ) * nreq / 2 elements. For illustration <pre>
- * cov =
- * {
- * cov_00,
- * cov_10, cov_11,
- * cov_20, cov_21, cov22,
- * ...
- * } </pre>
- *
- * @param nreq how many of the regressors to include (either in canonical
- * order, or in the current reordered state)
- * @return an array with the variance covariance of the included
- * regressors in lower triangular form
- */
- private double[] cov(int nreq) {
- if (this.nobs <= nreq) {
- return null;
- }
- double rnk = 0.0;
- for (int i = 0; i < nreq; i++) {
- if (!this.lindep[i]) {
- rnk += 1.0;
- }
- }
- final double var = rss[nreq - 1] / (nobs - rnk);
- final double[] rinv = new double[nreq * (nreq - 1) / 2];
- inverse(rinv, nreq);
- final double[] covmat = new double[nreq * (nreq + 1) / 2];
- Arrays.fill(covmat, Double.NaN);
- int pos2;
- int pos1;
- int start = 0;
- double total = 0;
- for (int row = 0; row < nreq; row++) {
- pos2 = start;
- if (!this.lindep[row]) {
- for (int col = row; col < nreq; col++) {
- if (!this.lindep[col]) {
- pos1 = start + col - row;
- if (row == col) {
- total = 1.0 / d[col];
- } else {
- total = rinv[pos1 - 1] / d[col];
- }
- for (int k = col + 1; k < nreq; k++) {
- if (!this.lindep[k]) {
- total += rinv[pos1] * rinv[pos2] / d[k];
- }
- ++pos1;
- ++pos2;
- }
- covmat[ (col + 1) * col / 2 + row] = total * var;
- } else {
- pos2 += nreq - col - 1;
- }
- }
- }
- start += nreq - row - 1;
- }
- return covmat;
- }
- /**
- * This internal method calculates the inverse of the upper-triangular portion
- * of the R matrix.
- * @param rinv the storage for the inverse of r
- * @param nreq how many of the regressors to include (either in canonical
- * order, or in the current reordered state)
- */
- private void inverse(double[] rinv, int nreq) {
- int pos = nreq * (nreq - 1) / 2 - 1;
- int pos1 = -1;
- int pos2 = -1;
- double total = 0.0;
- Arrays.fill(rinv, Double.NaN);
- for (int row = nreq - 1; row > 0; --row) {
- if (!this.lindep[row]) {
- final int start = (row - 1) * (nvars + nvars - row) / 2;
- for (int col = nreq; col > row; --col) {
- pos1 = start;
- pos2 = pos;
- total = 0.0;
- for (int k = row; k < col - 1; k++) {
- pos2 += nreq - k - 1;
- if (!this.lindep[k]) {
- total += -r[pos1] * rinv[pos2];
- }
- ++pos1;
- }
- rinv[pos] = total - r[pos1];
- --pos;
- }
- } else {
- pos -= nreq - row;
- }
- }
- }
- /**
- * In the original algorithm only the partial correlations of the regressors
- * is returned to the user. In this implementation, we have <pre>
- * corr =
- * {
- * corrxx - lower triangular
- * corrxy - bottom row of the matrix
- * }
- * Replaces subroutines PCORR and COR of:
- * ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2 </pre>
- *
- * <p>Calculate partial correlations after the variables in rows
- * 1, 2, ..., IN have been forced into the regression.
- * If IN = 1, and the first row of R represents a constant in the
- * model, then the usual simple correlations are returned.</p>
- *
- * <p>If IN = 0, the value returned in array CORMAT for the correlation
- * of variables Xi & Xj is: <pre>
- * sum ( Xi.Xj ) / Sqrt ( sum (Xi^2) . sum (Xj^2) )</pre>
- *
- * <p>On return, array CORMAT contains the upper triangle of the matrix of
- * partial correlations stored by rows, excluding the 1's on the diagonal.
- * e.g. if IN = 2, the consecutive elements returned are:
- * (3,4) (3,5) ... (3,ncol), (4,5) (4,6) ... (4,ncol), etc.
- * Array YCORR stores the partial correlations with the Y-variable
- * starting with YCORR(IN+1) = partial correlation with the variable in
- * position (IN+1). </p>
- *
- * @param in how many of the regressors to include (either in canonical
- * order, or in the current reordered state)
- * @return an array with the partial correlations of the remainder of
- * regressors with each other and the regressand, in lower triangular form
- */
- public double[] getPartialCorrelations(int in) {
- final double[] output = new double[(nvars - in + 1) * (nvars - in) / 2];
- int pos;
- int pos1;
- int pos2;
- final int rmsOff = -in;
- final int wrkOff = -(in + 1);
- final double[] rms = new double[nvars - in];
- final double[] work = new double[nvars - in - 1];
- double sumxx;
- double sumxy;
- double sumyy;
- final int offXX = (nvars - in) * (nvars - in - 1) / 2;
- if (in < -1 || in >= nvars) {
- return null;
- }
- final int nvm = nvars - 1;
- final int basePos = r.length - (nvm - in) * (nvm - in + 1) / 2;
- if (d[in] > 0.0) {
- rms[in + rmsOff] = 1.0 / JdkMath.sqrt(d[in]);
- }
- for (int col = in + 1; col < nvars; col++) {
- pos = basePos + col - 1 - in;
- sumxx = d[col];
- for (int row = in; row < col; row++) {
- sumxx += d[row] * r[pos] * r[pos];
- pos += nvars - row - 2;
- }
- if (sumxx > 0.0) {
- rms[col + rmsOff] = 1.0 / JdkMath.sqrt(sumxx);
- } else {
- rms[col + rmsOff] = 0.0;
- }
- }
- sumyy = sserr;
- for (int row = in; row < nvars; row++) {
- sumyy += d[row] * rhs[row] * rhs[row];
- }
- if (sumyy > 0.0) {
- sumyy = 1.0 / JdkMath.sqrt(sumyy);
- }
- pos = 0;
- for (int col1 = in; col1 < nvars; col1++) {
- sumxy = 0.0;
- Arrays.fill(work, 0.0);
- pos1 = basePos + col1 - in - 1;
- for (int row = in; row < col1; row++) {
- pos2 = pos1 + 1;
- for (int col2 = col1 + 1; col2 < nvars; col2++) {
- work[col2 + wrkOff] += d[row] * r[pos1] * r[pos2];
- pos2++;
- }
- sumxy += d[row] * r[pos1] * rhs[row];
- pos1 += nvars - row - 2;
- }
- pos2 = pos1 + 1;
- for (int col2 = col1 + 1; col2 < nvars; col2++) {
- work[col2 + wrkOff] += d[col1] * r[pos2];
- ++pos2;
- output[ (col2 - 1 - in) * (col2 - in) / 2 + col1 - in] =
- work[col2 + wrkOff] * rms[col1 + rmsOff] * rms[col2 + rmsOff];
- ++pos;
- }
- sumxy += d[col1] * rhs[col1];
- output[col1 + rmsOff + offXX] = sumxy * rms[col1 + rmsOff] * sumyy;
- }
- return output;
- }
- /**
- * ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2.
- * Move variable from position FROM to position TO in an
- * orthogonal reduction produced by AS75.1.
- *
- * @param from initial position
- * @param to destination
- */
- private void vmove(int from, int to) {
- double d1;
- double d2;
- double x;
- double d1new;
- double d2new;
- double cbar;
- double sbar;
- double y;
- int first;
- int inc;
- int m1;
- int m2;
- int mp1;
- int pos;
- boolean bSkipTo40 = false;
- if (from == to) {
- return;
- }
- if (!this.rssSet) {
- ss();
- }
- int count = 0;
- if (from < to) {
- first = from;
- inc = 1;
- count = to - from;
- } else {
- first = from - 1;
- inc = -1;
- count = from - to;
- }
- int m = first;
- int idx = 0;
- while (idx < count) {
- m1 = m * (nvars + nvars - m - 1) / 2;
- m2 = m1 + nvars - m - 1;
- mp1 = m + 1;
- d1 = d[m];
- d2 = d[mp1];
- // Special cases.
- if (d1 > this.epsilon || d2 > this.epsilon) {
- x = r[m1];
- if (JdkMath.abs(x) * JdkMath.sqrt(d1) < tol[mp1]) {
- x = 0.0;
- }
- if (d1 < this.epsilon || JdkMath.abs(x) < this.epsilon) {
- d[m] = d2;
- d[mp1] = d1;
- r[m1] = 0.0;
- for (int col = m + 2; col < nvars; col++) {
- ++m1;
- x = r[m1];
- r[m1] = r[m2];
- r[m2] = x;
- ++m2;
- }
- x = rhs[m];
- rhs[m] = rhs[mp1];
- rhs[mp1] = x;
- bSkipTo40 = true;
- //break;
- } else if (d2 < this.epsilon) {
- d[m] = d1 * x * x;
- r[m1] = 1.0 / x;
- for (int i = m1 + 1; i < m1 + nvars - m - 1; i++) {
- r[i] /= x;
- }
- rhs[m] /= x;
- bSkipTo40 = true;
- //break;
- }
- if (!bSkipTo40) {
- d1new = d2 + d1 * x * x;
- cbar = d2 / d1new;
- sbar = x * d1 / d1new;
- d2new = d1 * cbar;
- d[m] = d1new;
- d[mp1] = d2new;
- r[m1] = sbar;
- for (int col = m + 2; col < nvars; col++) {
- ++m1;
- y = r[m1];
- r[m1] = cbar * r[m2] + sbar * y;
- r[m2] = y - x * r[m2];
- ++m2;
- }
- y = rhs[m];
- rhs[m] = cbar * rhs[mp1] + sbar * y;
- rhs[mp1] = y - x * rhs[mp1];
- }
- }
- if (m > 0) {
- pos = m;
- for (int row = 0; row < m; row++) {
- x = r[pos];
- r[pos] = r[pos - 1];
- r[pos - 1] = x;
- pos += nvars - row - 2;
- }
- }
- // Adjust variable order (VORDER), the tolerances (TOL) and
- // the vector of residual sums of squares (RSS).
- m1 = vorder[m];
- vorder[m] = vorder[mp1];
- vorder[mp1] = m1;
- x = tol[m];
- tol[m] = tol[mp1];
- tol[mp1] = x;
- rss[m] = rss[mp1] + d[mp1] * rhs[mp1] * rhs[mp1];
- m += inc;
- ++idx;
- }
- }
- /**
- * ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2
- *
- * <p> Re-order the variables in an orthogonal reduction produced by
- * AS75.1 so that the N variables in LIST start at position POS1,
- * though will not necessarily be in the same order as in LIST.
- * Any variables in VORDER before position POS1 are not moved.
- * Auxiliary routine called: VMOVE. </p>
- *
- * <p>This internal method reorders the regressors.</p>
- *
- * @param list the regressors to move
- * @param pos1 where the list will be placed
- * @return -1 error, 0 everything ok
- */
- private int reorderRegressors(int[] list, int pos1) {
- int next;
- int i;
- int l;
- if (list.length < 1 || list.length > nvars + 1 - pos1) {
- return -1;
- }
- next = pos1;
- i = pos1;
- while (i < nvars) {
- l = vorder[i];
- for (int j = 0; j < list.length; j++) {
- if (l == list[j] && i > next) {
- this.vmove(i, next);
- ++next;
- if (next >= list.length + pos1) {
- return 0;
- } else {
- break;
- }
- }
- }
- ++i;
- }
- return 0;
- }
- /**
- * Gets the diagonal of the Hat matrix also known as the leverage matrix.
- *
- * @param rowData returns the diagonal of the hat matrix for this observation
- * @return the diagonal element of the hatmatrix
- */
- public double getDiagonalOfHatMatrix(double[] rowData) {
- double[] wk = new double[this.nvars];
- int pos;
- double total;
- if (rowData.length > nvars) {
- return Double.NaN;
- }
- double[] xrow;
- if (this.hasIntercept) {
- xrow = new double[rowData.length + 1];
- xrow[0] = 1.0;
- System.arraycopy(rowData, 0, xrow, 1, rowData.length);
- } else {
- xrow = rowData;
- }
- double hii = 0.0;
- for (int col = 0; col < xrow.length; col++) {
- if (JdkMath.sqrt(d[col]) < tol[col]) {
- wk[col] = 0.0;
- } else {
- pos = col - 1;
- total = xrow[col];
- for (int row = 0; row < col; row++) {
- total = smartAdd(total, -wk[row] * r[pos]);
- pos += nvars - row - 2;
- }
- wk[col] = total;
- hii = smartAdd(hii, (total * total) / d[col]);
- }
- }
- return hii;
- }
- /**
- * Gets the order of the regressors, useful if some type of reordering
- * has been called. Calling regress with int[]{} args will trigger
- * a reordering.
- *
- * @return int[] with the current order of the regressors
- */
- public int[] getOrderOfRegressors(){
- return Arrays.copyOf(vorder, vorder.length);
- }
- /**
- * Conducts a regression on the data in the model, using all regressors.
- *
- * @return RegressionResults the structure holding all regression results
- * @exception ModelSpecificationException - thrown if number of observations is
- * less than the number of variables
- */
- @Override
- public RegressionResults regress() throws ModelSpecificationException {
- return regress(this.nvars);
- }
- /**
- * Conducts a regression on the data in the model, using a subset of regressors.
- *
- * @param numberOfRegressors many of the regressors to include (either in canonical
- * order, or in the current reordered state)
- * @return RegressionResults the structure holding all regression results
- * @exception ModelSpecificationException - thrown if number of observations is
- * less than the number of variables or number of regressors requested
- * is greater than the regressors in the model
- */
- public RegressionResults regress(int numberOfRegressors) throws ModelSpecificationException {
- if (this.nobs <= numberOfRegressors) {
- throw new ModelSpecificationException(
- LocalizedFormats.NOT_ENOUGH_DATA_FOR_NUMBER_OF_PREDICTORS,
- this.nobs, numberOfRegressors);
- }
- if( numberOfRegressors > this.nvars ){
- throw new ModelSpecificationException(
- LocalizedFormats.TOO_MANY_REGRESSORS, numberOfRegressors, this.nvars);
- }
- tolset();
- singcheck();
- double[] beta = this.regcf(numberOfRegressors);
- ss();
- double[] cov = this.cov(numberOfRegressors);
- int rnk = 0;
- for (int i = 0; i < this.lindep.length; i++) {
- if (!this.lindep[i]) {
- ++rnk;
- }
- }
- boolean needsReorder = false;
- for (int i = 0; i < numberOfRegressors; i++) {
- if (this.vorder[i] != i) {
- needsReorder = true;
- break;
- }
- }
- if (!needsReorder) {
- return new RegressionResults(
- beta, new double[][]{cov}, true, this.nobs, rnk,
- this.sumy, this.sumsqy, this.sserr, this.hasIntercept, false);
- } else {
- double[] betaNew = new double[beta.length];
- double[] covNew = new double[cov.length];
- int[] newIndices = new int[beta.length];
- for (int i = 0; i < nvars; i++) {
- for (int j = 0; j < numberOfRegressors; j++) {
- if (this.vorder[j] == i) {
- betaNew[i] = beta[ j];
- newIndices[i] = j;
- }
- }
- }
- int idx1 = 0;
- int idx2;
- int ii;
- int jj;
- for (int i = 0; i < beta.length; i++) {
- ii = newIndices[i];
- for (int j = 0; j <= i; j++, idx1++) {
- jj = newIndices[j];
- if (ii > jj) {
- idx2 = ii * (ii + 1) / 2 + jj;
- } else {
- idx2 = jj * (jj + 1) / 2 + ii;
- }
- covNew[idx1] = cov[idx2];
- }
- }
- return new RegressionResults(
- betaNew, new double[][]{covNew}, true, this.nobs, rnk,
- this.sumy, this.sumsqy, this.sserr, this.hasIntercept, false);
- }
- }
- /**
- * Conducts a regression on the data in the model, using regressors in array
- * Calling this method will change the internal order of the regressors
- * and care is required in interpreting the hatmatrix.
- *
- * @param variablesToInclude array of variables to include in regression
- * @return RegressionResults the structure holding all regression results
- * @exception ModelSpecificationException - thrown if number of observations is
- * less than the number of variables, the number of regressors requested
- * is greater than the regressors in the model or a regressor index in
- * regressor array does not exist
- */
- @Override
- public RegressionResults regress(int[] variablesToInclude) throws ModelSpecificationException {
- if (variablesToInclude.length > this.nvars) {
- throw new ModelSpecificationException(
- LocalizedFormats.TOO_MANY_REGRESSORS, variablesToInclude.length, this.nvars);
- }
- if (this.nobs <= this.nvars) {
- throw new ModelSpecificationException(
- LocalizedFormats.NOT_ENOUGH_DATA_FOR_NUMBER_OF_PREDICTORS,
- this.nobs, this.nvars);
- }
- Arrays.sort(variablesToInclude);
- int iExclude = 0;
- for (int i = 0; i < variablesToInclude.length; i++) {
- if (i >= this.nvars) {
- throw new ModelSpecificationException(
- LocalizedFormats.INDEX_LARGER_THAN_MAX, i, this.nvars);
- }
- if (i > 0 && variablesToInclude[i] == variablesToInclude[i - 1]) {
- variablesToInclude[i] = -1;
- ++iExclude;
- }
- }
- int[] series;
- if (iExclude > 0) {
- int j = 0;
- series = new int[variablesToInclude.length - iExclude];
- for (int i = 0; i < variablesToInclude.length; i++) {
- if (variablesToInclude[i] > -1) {
- series[j] = variablesToInclude[i];
- ++j;
- }
- }
- } else {
- series = variablesToInclude;
- }
- reorderRegressors(series, 0);
- tolset();
- singcheck();
- double[] beta = this.regcf(series.length);
- ss();
- double[] cov = this.cov(series.length);
- int rnk = 0;
- for (int i = 0; i < this.lindep.length; i++) {
- if (!this.lindep[i]) {
- ++rnk;
- }
- }
- boolean needsReorder = false;
- for (int i = 0; i < series.length; i++) {
- if (this.vorder[i] != series[i]) {
- needsReorder = true;
- break;
- }
- }
- if (!needsReorder) {
- return new RegressionResults(
- beta, new double[][]{cov}, true, this.nobs, rnk,
- this.sumy, this.sumsqy, this.sserr, this.hasIntercept, false);
- } else {
- double[] betaNew = new double[beta.length];
- int[] newIndices = new int[beta.length];
- for (int i = 0; i < series.length; i++) {
- for (int j = 0; j < this.vorder.length; j++) {
- if (this.vorder[j] == series[i]) {
- betaNew[i] = beta[ j];
- newIndices[i] = j;
- }
- }
- }
- double[] covNew = new double[cov.length];
- int idx1 = 0;
- int idx2;
- int ii;
- int jj;
- for (int i = 0; i < beta.length; i++) {
- ii = newIndices[i];
- for (int j = 0; j <= i; j++, idx1++) {
- jj = newIndices[j];
- if (ii > jj) {
- idx2 = ii * (ii + 1) / 2 + jj;
- } else {
- idx2 = jj * (jj + 1) / 2 + ii;
- }
- covNew[idx1] = cov[idx2];
- }
- }
- return new RegressionResults(
- betaNew, new double[][]{covNew}, true, this.nobs, rnk,
- this.sumy, this.sumsqy, this.sserr, this.hasIntercept, false);
- }
- }
- }