Java Code Examples for org.apache.commons.math.linear.RealMatrix#setEntry()

/**
 * Get the covariance matrix.
 * @return covariance matrix
 */
public RealMatrix getResult() {

    int dimension = sums.length;
    RealMatrix result = MatrixUtils.createRealMatrix(dimension, dimension);

    if (n > 1) {
        double c = 1.0 / (n * (isBiasCorrected ? (n - 1) : n));
        int k = 0;
        for (int i = 0; i < dimension; ++i) {
            for (int j = 0; j <= i; ++j) {
                double e = c * (n * productsSums[k++] - sums[i] * sums[j]);
                result.setEntry(i, j, e);
                result.setEntry(j, i, e);
            }
        }
    }

    return result;

}
 

/**
 * Return sign of values in matrix:
 * val > 0 : 1
 * val = 0 : 0
 * val < 0 : -1
 */
public static RealMatrix sign(RealMatrix matrix) {
    int r = matrix.getRowDimension();
    int c = matrix.getColumnDimension();
    RealMatrix signMatrix = cleanArray2DMatrix(r, c);
    for (int i = 0; i < r; i++) {
        for (int j = 0; j < c; j++) {
            double val = matrix.getEntry(i, j);
            if (val > 0) {
                signMatrix.setEntry(i, j, 1);
            } else if (val < 0) {
                signMatrix.setEntry(i, j, -1);
            }
        }
    }
    return signMatrix;
}
 

/**
 * Calculate cumulative sums across upper right matrix
 * iterative result is entry to left + entry below + orig value - diagonal down (else we double count)
 *
 * @param matrix
 * @param superDiagonal
 * @return
 */
public static RealMatrix sum(RealMatrix matrix, int superDiagonal) {

    int n = Math.min(matrix.getRowDimension(), matrix.getColumnDimension());

    RealMatrix diagMatrix = MatrixTools.cleanArray2DMatrix(matrix.getRowDimension(), matrix.getColumnDimension());
    if (superDiagonal <= 0) {
        diagMatrix = MatrixTools.extractDiagonal(matrix);
        superDiagonal = 1;
    }

    // d = distance from diagonal
    for (int d = superDiagonal; d < n; d++) {
        // i = row, column is i +d;

        for (int i = 0; i < n - d; i++) {
            diagMatrix.setEntry(i, i + d,
                    diagMatrix.getEntry(i, i + d - 1) + diagMatrix.getEntry(i + 1, i + d) +
                            matrix.getEntry(i, i + d) - diagMatrix.getEntry(i + 1, i + d - 1));
        }
    }
    return diagMatrix;
}
 

/**
 * Derives a correlation matrix from a covariance matrix.
 *
 * <p>Uses the formula <br/>
 * <code>r(X,Y) = cov(X,Y)/s(X)s(Y)</code> where
 * <code>r(&middot,&middot;)</code> is the correlation coefficient and
 * <code>s(&middot;)</code> means standard deviation.</p>
 *
 * @param covarianceMatrix the covariance matrix
 * @return correlation matrix
 */
public RealMatrix covarianceToCorrelation(RealMatrix covarianceMatrix) {
    int nVars = covarianceMatrix.getColumnDimension();
    RealMatrix outMatrix = new BlockRealMatrix(nVars, nVars);
    for (int i = 0; i < nVars; i++) {
        double sigma = Math.sqrt(covarianceMatrix.getEntry(i, i));
        outMatrix.setEntry(i, i, 1d);
        for (int j = 0; j < i; j++) {
            double entry = covarianceMatrix.getEntry(i, j) /
                   (sigma * Math.sqrt(covarianceMatrix.getEntry(j, j)));
            outMatrix.setEntry(i, j, entry);
            outMatrix.setEntry(j, i, entry);
        }
    }
    return outMatrix;
}
 

protected RealMatrix createLowerTriangularRealMatrix(double[] data, int dimension) {
    int ptr = 0;
    RealMatrix result = new BlockRealMatrix(dimension, dimension);
    for (int i = 1; i < dimension; i++) {
        for (int j = 0; j < i; j++) {
            result.setEntry(i, j, data[ptr]);
            ptr++;
        }
    }
    return result;
}
 

/**
 * Compute a covariance matrix from a matrix whose columns represent
 * covariates.
 * @param matrix input matrix (must have at least two columns and two rows)
 * @param biasCorrected determines whether or not covariance estimates are bias-corrected
 * @return covariance matrix
 */
protected RealMatrix computeCovarianceMatrix(RealMatrix matrix, boolean biasCorrected) {
    int dimension = matrix.getColumnDimension();
    Variance variance = new Variance(biasCorrected);
    RealMatrix outMatrix = new BlockRealMatrix(dimension, dimension);
    for (int i = 0; i < dimension; i++) {
        for (int j = 0; j < i; j++) {
          double cov = covariance(matrix.getColumn(i), matrix.getColumn(j), biasCorrected);
          outMatrix.setEntry(i, j, cov);
          outMatrix.setEntry(j, i, cov);
        }
        outMatrix.setEntry(i, i, variance.evaluate(matrix.getColumn(i)));
    }
    return outMatrix;
}
 

/**
 * @return matrix flipped Top-Bottom
 */
public static RealMatrix flipTopBottom(RealMatrix matrix) {
    int r = matrix.getRowDimension(), c = matrix.getColumnDimension();
    RealMatrix topBottomFlippedMatrix = cleanArray2DMatrix(r, c);
    for (int i = 0; i < r; i++) {
        for (int j = 0; j < c; j++) {
            topBottomFlippedMatrix.setEntry(r - 1 - i, j, matrix.getEntry(i, j));
        }
    }
    return topBottomFlippedMatrix;
}
 

/**
 * @return matrix flipped Top-Bottom
 */
public static RealMatrix flipTopBottom(RealMatrix matrix) {
    int r = matrix.getRowDimension(), c = matrix.getColumnDimension();
    RealMatrix topBottomFlippedMatrix = cleanArray2DMatrix(r, c);
    for (int i = 0; i < r; i++) {
        for (int j = 0; j < c; j++) {
            topBottomFlippedMatrix.setEntry(r - 1 - i, j, matrix.getEntry(i, j));
        }
    }
    return topBottomFlippedMatrix;
}
 

/**
 * Computes the correlation matrix for the columns of the
 * input matrix.
 * 
 * @param matrix matrix with columns representing variables to correlate
 * @return correlation matrix
 */
public RealMatrix computeCorrelationMatrix(RealMatrix matrix) {
    int nVars = matrix.getColumnDimension();
    RealMatrix outMatrix = new BlockRealMatrix(nVars, nVars);
    for (int i = 0; i < nVars; i++) {
        for (int j = 0; j < i; j++) {
          double corr = correlation(matrix.getColumn(i), matrix.getColumn(j));
          outMatrix.setEntry(i, j, corr);
          outMatrix.setEntry(j, i, corr);
        }
        outMatrix.setEntry(i, i, 1d);
    }
    return outMatrix;
}
 

protected RealMatrix createLowerTriangularRealMatrix(double[] data, int dimension) {
    int ptr = 0;
    RealMatrix result = new BlockRealMatrix(dimension, dimension);
    for (int i = 1; i < dimension; i++) {
        for (int j = 0; j < i; j++) {
            result.setEntry(i, j, data[ptr]);
            ptr++;
        }
    }
    return result;
}
 

/**
 * Compute a covariance matrix from a matrix whose columns represent
 * covariates.
 * @param matrix input matrix (must have at least two columns and two rows)
 * @param biasCorrected determines whether or not covariance estimates are bias-corrected
 * @return covariance matrix
 */
protected RealMatrix computeCovarianceMatrix(RealMatrix matrix, boolean biasCorrected) {
    int dimension = matrix.getColumnDimension();
    Variance variance = new Variance(biasCorrected);
    RealMatrix outMatrix = new BlockRealMatrix(dimension, dimension);
    for (int i = 0; i < dimension; i++) {
        for (int j = 0; j < i; j++) {
          double cov = covariance(matrix.getColumn(i), matrix.getColumn(j), biasCorrected);
          outMatrix.setEntry(i, j, cov);
          outMatrix.setEntry(j, i, cov);
        }
        outMatrix.setEntry(i, i, variance.evaluate(matrix.getColumn(i)));
    }
    return outMatrix;
}
 

/**
 * Returns the values along the diagonal of the matrix
 *
 * @param matrix
 * @return diagonal
 */
public static RealMatrix makeSymmetricMatrix(RealMatrix matrix) {
    RealMatrix symmetricMatrix = extractDiagonal(matrix);
    int n = symmetricMatrix.getRowDimension();
    for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) {
            double val = matrix.getEntry(i, j);
            symmetricMatrix.setEntry(i, j, val);
            symmetricMatrix.setEntry(j, i, val);
        }
    }

    return symmetricMatrix;
}
 

protected void fillUpper(RealMatrix matrix, double diagonalValue) {
    int dimension = matrix.getColumnDimension();
    for (int i = 0; i < dimension; i++) {
        matrix.setEntry(i, i, diagonalValue);
        for (int j = i+1; j < dimension; j++) {
            matrix.setEntry(i, j, matrix.getEntry(j, i));
        }
    }
}
 

public CorrelatedRandomVectorGeneratorTest() {
    mean = new double[] { 0.0, 1.0, -3.0, 2.3 };

    RealMatrix b = MatrixUtils.createRealMatrix(4, 3);
    int counter = 0;
    for (int i = 0; i < b.getRowDimension(); ++i) {
        for (int j = 0; j < b.getColumnDimension(); ++j) {
            b.setEntry(i, j, 1.0 + 0.1 * ++counter);
        }
    }
    RealMatrix bbt = b.multiply(b.transpose());
    covariance = MatrixUtils.createRealMatrix(mean.length, mean.length);
    for (int i = 0; i < covariance.getRowDimension(); ++i) {
        covariance.setEntry(i, i, bbt.getEntry(i, i));
        for (int j = 0; j < covariance.getColumnDimension(); ++j) {
            double s = bbt.getEntry(i, j);
            covariance.setEntry(i, j, s);
            covariance.setEntry(j, i, s);
        }
    }

    RandomGenerator rg = new JDKRandomGenerator();
    rg.setSeed(17399225432l);
    GaussianRandomGenerator rawGenerator = new GaussianRandomGenerator(rg);
    generator = new CorrelatedRandomVectorGenerator(mean,
                                                    covariance,
                                                    1.0e-12 * covariance.getNorm(),
                                                    rawGenerator);
}
 

/**
 * Compute a covariance matrix from a matrix whose columns represent
 * covariates.
 * @param matrix input matrix (must have at least two columns and two rows)
 * @param biasCorrected determines whether or not covariance estimates are bias-corrected
 * @return covariance matrix
 */
protected RealMatrix computeCovarianceMatrix(RealMatrix matrix, boolean biasCorrected) {
    int dimension = matrix.getColumnDimension();
    Variance variance = new Variance(biasCorrected);
    RealMatrix outMatrix = new BlockRealMatrix(dimension, dimension);
    for (int i = 0; i < dimension; i++) {
        for (int j = 0; j < i; j++) {
          double cov = covariance(matrix.getColumn(i), matrix.getColumn(j), biasCorrected);
          outMatrix.setEntry(i, j, cov);
          outMatrix.setEntry(j, i, cov);
        }
        outMatrix.setEntry(i, i, variance.evaluate(matrix.getColumn(i)));
    }
    return outMatrix;
}
 

public static RealMatrix cleanUpNaNs(RealMatrix matrix) {
    for (int r = 0; r < matrix.getRowDimension(); r++)
        for (int c = 0; c < matrix.getColumnDimension(); c++)
            if (Double.isNaN(matrix.getEntry(r, c))) {
                matrix.setEntry(r, c, 0);
            }
    return matrix;
}
 

/**
 * Compute a covariance matrix from a matrix whose columns represent
 * covariates.
 * @param matrix input matrix (must have at least two columns and two rows)
 * @param biasCorrected determines whether or not covariance estimates are bias-corrected
 * @return covariance matrix
 */
protected RealMatrix computeCovarianceMatrix(RealMatrix matrix, boolean biasCorrected) {
    int dimension = matrix.getColumnDimension();
    Variance variance = new Variance(biasCorrected);
    RealMatrix outMatrix = new BlockRealMatrix(dimension, dimension);
    for (int i = 0; i < dimension; i++) {
        for (int j = 0; j < i; j++) {
          double cov = covariance(matrix.getColumn(i), matrix.getColumn(j), biasCorrected);
          outMatrix.setEntry(i, j, cov);
          outMatrix.setEntry(j, i, cov);
        }
        outMatrix.setEntry(i, i, variance.evaluate(matrix.getColumn(i)));
    }
    return outMatrix;
}
 

protected void fillUpper(RealMatrix matrix, double diagonalValue) {
    int dimension = matrix.getColumnDimension();
    for (int i = 0; i < dimension; i++) {
        matrix.setEntry(i, i, diagonalValue);
        for (int j = i+1; j < dimension; j++) {
            matrix.setEntry(i, j, matrix.getEntry(j, i));
        }
    }
}
 

/**
 * Generate an error covariance matrix and sample data representing models
 * with this error structure. Then verify that GLS estimated coefficients,
 * on average, perform better than OLS.
 */
@Test
public void testGLSEfficiency() throws Exception {
    RandomGenerator rg = new JDKRandomGenerator();
    rg.setSeed(200);  // Seed has been selected to generate non-trivial covariance
    
    // Assume model has 16 observations (will use Longley data).  Start by generating
    // non-constant variances for the 16 error terms.
    final int nObs = 16;
    double[] sigma = new double[nObs];
    for (int i = 0; i < nObs; i++) {
        sigma[i] = 10 * rg.nextDouble();
    }
    
    // Now generate 1000 error vectors to use to estimate the covariance matrix
    // Columns are draws on N(0, sigma[col])
    final int numSeeds = 1000;
    RealMatrix errorSeeds = MatrixUtils.createRealMatrix(numSeeds, nObs);
    for (int i = 0; i < numSeeds; i++) {
        for (int j = 0; j < nObs; j++) {
            errorSeeds.setEntry(i, j, rg.nextGaussian() * sigma[j]);
        }
    }
    
    // Get covariance matrix for columns
    RealMatrix cov = (new Covariance(errorSeeds)).getCovarianceMatrix();
      
    // Create a CorrelatedRandomVectorGenerator to use to generate correlated errors
    GaussianRandomGenerator rawGenerator = new GaussianRandomGenerator(rg);
    double[] errorMeans = new double[nObs];  // Counting on init to 0 here
    CorrelatedRandomVectorGenerator gen = new CorrelatedRandomVectorGenerator(errorMeans, cov,
     1.0e-12 * cov.getNorm(), rawGenerator);
    
    // Now start generating models.  Use Longley X matrix on LHS
    // and Longley OLS beta vector as "true" beta.  Generate
    // Y values by XB + u where u is a CorrelatedRandomVector generated
    // from cov.
    OLSMultipleLinearRegression ols = new OLSMultipleLinearRegression();
    ols.newSampleData(longley, nObs, 6);
    final RealVector b = ols.calculateBeta().copy();
    final RealMatrix x = ols.X.copy();
    
    // Create a GLS model to reuse
    GLSMultipleLinearRegression gls = new GLSMultipleLinearRegression();
    gls.newSampleData(longley, nObs, 6);
    gls.newCovarianceData(cov.getData());
    
    // Create aggregators for stats measuring model performance
    DescriptiveStatistics olsBetaStats = new DescriptiveStatistics();
    DescriptiveStatistics glsBetaStats = new DescriptiveStatistics();
    
    // Generate Y vectors for 10000 models, estimate GLS and OLS and
    // Verify that OLS estimates are better
    final int nModels = 10000;
    for (int i = 0; i < nModels; i++) {
        
        // Generate y = xb + u with u cov
        RealVector u = MatrixUtils.createRealVector(gen.nextVector());
        double[] y = u.add(x.operate(b)).getData();
        
        // Estimate OLS parameters
        ols.newYSampleData(y);
        RealVector olsBeta = ols.calculateBeta();
        
        // Estimate GLS parameters
        gls.newYSampleData(y);
        RealVector glsBeta = gls.calculateBeta();
        
        // Record deviations from "true" beta
        double dist = olsBeta.getDistance(b);
        olsBetaStats.addValue(dist * dist);
        dist = glsBeta.getDistance(b);
        glsBetaStats.addValue(dist * dist);
        
    }
    
    // Verify that GLS is on average more efficient, lower variance
    assert(olsBetaStats.getMean() > 1.5 * glsBetaStats.getMean());
    assert(olsBetaStats.getStandardDeviation() > glsBetaStats.getStandardDeviation());  
}
 

/** 
 * Solve an estimation problem using a least squares criterion.
 *
 * <p>This method set the unbound parameters of the given problem
 * starting from their current values through several iterations. At
 * each step, the unbound parameters are changed in order to
 * minimize a weighted least square criterion based on the
 * measurements of the problem.</p>
 *
 * <p>The iterations are stopped either when the criterion goes
 * below a physical threshold under which improvement are considered
 * useless or when the algorithm is unable to improve it (even if it
 * is still high). The first condition that is met stops the
 * iterations. If the convergence it not reached before the maximum
 * number of iterations, an {@link EstimationException} is
 * thrown.</p>
 *
 * @param problem estimation problem to solve
 * @exception EstimationException if the problem cannot be solved
 *
 * @see EstimationProblem
 *
 */
@Override
public void estimate(EstimationProblem problem)
throws EstimationException {

    initializeEstimate(problem);

    // work matrices
    double[] grad             = new double[parameters.length];
    ArrayRealVector bDecrement = new ArrayRealVector(parameters.length);
    double[] bDecrementData   = bDecrement.getDataRef();
    RealMatrix wGradGradT     = MatrixUtils.createRealMatrix(parameters.length, parameters.length);

    // iterate until convergence is reached
    double previous = Double.POSITIVE_INFINITY;
    do {

        // build the linear problem
        incrementJacobianEvaluationsCounter();
        RealVector b = new ArrayRealVector(parameters.length);
        RealMatrix a = MatrixUtils.createRealMatrix(parameters.length, parameters.length);
        for (int i = 0; i < measurements.length; ++i) {
            if (! measurements [i].isIgnored()) {

                double weight   = measurements[i].getWeight();
                double residual = measurements[i].getResidual();

                // compute the normal equation
                for (int j = 0; j < parameters.length; ++j) {
                    grad[j] = measurements[i].getPartial(parameters[j]);
                    bDecrementData[j] = weight * residual * grad[j];
                }

                // build the contribution matrix for measurement i
                for (int k = 0; k < parameters.length; ++k) {
                    double gk = grad[k];
                    for (int l = 0; l < parameters.length; ++l) {
                        wGradGradT.setEntry(k, l, weight * gk * grad[l]);
                    }
                }

                // update the matrices
                a = a.add(wGradGradT);
                b = b.add(bDecrement);

            }
        }

        try {

            // solve the linearized least squares problem
            RealVector dX = new LUDecompositionImpl(a).getSolver().solve(b);

            // update the estimated parameters
            for (int i = 0; i < parameters.length; ++i) {
                parameters[i].setEstimate(parameters[i].getEstimate() + dX.getEntry(i));
            }

        } catch(InvalidMatrixException e) {
            throw new EstimationException("unable to solve: singular problem");
        }


        previous = cost;
        updateResidualsAndCost();

    } while ((getCostEvaluations() < 2) ||
             (Math.abs(previous - cost) > (cost * steadyStateThreshold) &&
              (Math.abs(cost) > convergence)));

}