Java RealVector.subtract方法代码示例

import org.apache.commons.math3.linear.RealVector; //导入方法依赖的package包/类
/**
 *
 * @param y     the response vector
 * @param x     the design matrix
 */
private RealMatrix computeBeta(RealVector y, RealMatrix x) {
    if (solver == Solver.QR) {
        return computeBetaQR(y, x);
    } else {
        final int n = x.getRowDimension();
        final int p = x.getColumnDimension();
        final int offset = hasIntercept() ? 1 : 0;
        final RealMatrix xT = x.transpose();
        final RealMatrix xTxInv = new LUDecomposition(xT.multiply(x)).getSolver().getInverse();
        final RealVector betaVector = xTxInv.multiply(xT).operate(y);
        final RealVector residuals = y.subtract(x.operate(betaVector));
        this.rss = residuals.dotProduct(residuals);
        this.errorVariance = rss / (n - p);
        this.stdError = Math.sqrt(errorVariance);
        this.residuals = createResidualsFrame(residuals);
        final RealMatrix covMatrix = xTxInv.scalarMultiply(errorVariance);
        final RealMatrix result = new Array2DRowRealMatrix(p, 2);
        if (hasIntercept()) {
            result.setEntry(0, 0, betaVector.getEntry(0));      //Intercept coefficient
            result.setEntry(0, 1, covMatrix.getEntry(0, 0));    //Intercept variance
        }
        for (int i = 0; i < getRegressors().size(); i++) {
            final int index = i + offset;
            final double variance = covMatrix.getEntry(index, index);
            result.setEntry(index, 1, variance);
            result.setEntry(index, 0, betaVector.getEntry(index));
        }
        return result;
    }
}
 

import org.apache.commons.math3.linear.RealVector; //导入方法依赖的package包/类
/**
 * Computes model parameters and parameter variance using a QR decomposition of the X matrix
 * @param y     the response vector
 * @param x     the design matrix
 */
private RealMatrix computeBetaQR(RealVector y, RealMatrix x) {
    final int n = x.getRowDimension();
    final int p = x.getColumnDimension();
    final int offset = hasIntercept() ? 1 : 0;
    final QRDecomposition decomposition = new QRDecomposition(x, threshold);
    final RealVector betaVector = decomposition.getSolver().solve(y);
    final RealVector residuals = y.subtract(x.operate(betaVector));
    this.rss = residuals.dotProduct(residuals);
    this.errorVariance = rss / (n - p);
    this.stdError = Math.sqrt(errorVariance);
    this.residuals = createResidualsFrame(residuals);
    final RealMatrix rAug = decomposition.getR().getSubMatrix(0, p - 1, 0, p - 1);
    final RealMatrix rInv = new LUDecomposition(rAug).getSolver().getInverse();
    final RealMatrix covMatrix = rInv.multiply(rInv.transpose()).scalarMultiply(errorVariance);
    final RealMatrix result = new Array2DRowRealMatrix(p, 2);
    if (hasIntercept()) {
        result.setEntry(0, 0, betaVector.getEntry(0));      //Intercept coefficient
        result.setEntry(0, 1, covMatrix.getEntry(0, 0));    //Intercept variance
    }
    for (int i = 0; i < getRegressors().size(); i++) {
        final int index = i + offset;
        final double variance = covMatrix.getEntry(index, index);
        result.setEntry(index, 1, variance);
        result.setEntry(index, 0, betaVector.getEntry(index));
    }
    return result;
}
 

import org.apache.commons.math3.linear.RealVector; //导入方法依赖的package包/类
/**
 * Correct the current state estimate with an actual measurement.
 *
 * @param z
 *            the measurement vector
 * @throws NullArgumentException
 *             if the measurement vector is {@code null}
 * @throws DimensionMismatchException
 *             if the dimension of the measurement vector does not fit
 * @throws SingularMatrixException
 *             if the covariance matrix could not be inverted
 */
public void correct(final RealVector z) throws NullArgumentException,
           DimensionMismatchException, SingularMatrixException
{

	// sanity checks
	MathUtils.checkNotNull(z);
	if (z.getDimension() != measurementMatrix.getRowDimension())
	{
		throw new DimensionMismatchException(z.getDimension(),
				measurementMatrix.getRowDimension());
	}

	// S = H * P(k) * H' + R
	RealMatrix s = measurementMatrix.multiply(errorCovariance)
			.multiply(measurementMatrixT)
			.add(measurementModel.getMeasurementNoise());

	// Inn = z(k) - H * xHat(k)-
	RealVector innovation = z.subtract(measurementMatrix
			.operate(stateEstimation));

	// calculate gain matrix
	// K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1
	// K(k) = P(k)- * H' * S^-1

	// instead of calculating the inverse of S we can rearrange the formula,
	// and then solve the linear equation A x X = B with A = S', X = K' and
	// B = (H * P)'

	// K(k) * S = P(k)- * H'
	// S' * K(k)' = H * P(k)-'
	RealMatrix kalmanGain = new CholeskyDecomposition(s).getSolver()
			.solve(measurementMatrix.multiply(errorCovariance.transpose()))
			.transpose();

	// update estimate with measurement z(k)
	// xHat(k) = xHat(k)- + K * Inn
	stateEstimation = stateEstimation.add(kalmanGain.operate(innovation));

	// update covariance of prediction error
	// P(k) = (I - K * H) * P(k)-
	RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain
			.getRowDimension());
	errorCovariance = identity.subtract(
			kalmanGain.multiply(measurementMatrix)).multiply(
			errorCovariance);
}