Example usage for org.apache.commons.math3.linear ArrayRealVector outerProduct

Introduction

In this page you can find the example usage for org.apache.commons.math3.linear ArrayRealVector outerProduct.

Prototype

@Override
public RealMatrix outerProduct(RealVector v) 

Source Link

Usage

From source file:edu.cudenver.bios.power.glmm.GLMMTestUnirepGeisserGreenhouse.java

/**
 * Calculate the Geisser-Greenhouse epsilon to correct for violations
 * of sphericity/*from   w w  w  . jav  a2 s .c o  m*/
 */
@Override
protected void calculateEpsilon() {
    super.calculateEpsilon();

    if (this.epsilonMethod == UnivariateEpsilonApproximation.MULLER_BARTON_APPROX) {
        // calculate the expected value of the epsilon estimate
        //  E[f(lambda)] = f(lambda) + g1 / (N - r)
        //  see Muller, Barton (1989) for details
        double b = rankU;
        double g1 = 0;
        for (int i = 0; i < distinctSigmaStarEigenValues.size(); i++) {
            EigenValueMultiplicityPair evmI = distinctSigmaStarEigenValues.get(i);
            double firstDerivative = ((2 * sumLambda) / (b * sumLambdaSquared)
                    - (2 * evmI.eigenValue * sumLambda * sumLambda)
                            / (b * sumLambdaSquared * sumLambdaSquared));

            double secondDerivative = (2 / (b * sumLambdaSquared)
                    - (8 * evmI.eigenValue * sumLambda) / (b * sumLambdaSquared * sumLambdaSquared)
                    + (8 * evmI.eigenValue * evmI.eigenValue * sumLambda * sumLambda)
                            / (b * sumLambdaSquared * sumLambdaSquared * sumLambdaSquared)
                    - (2 * sumLambda * sumLambda) / (b * sumLambdaSquared * sumLambdaSquared));

            // accumulate the first term of g1 (sum over distinct eigen vals of 1st derivative * eigen val ^2 * multiplicity)
            g1 += secondDerivative * evmI.eigenValue * evmI.eigenValue * evmI.multiplicity;
            // loop over elements not equal to current
            for (int j = 0; j < distinctSigmaStarEigenValues.size(); j++) {
                if (i != j) {
                    EigenValueMultiplicityPair evmJ = distinctSigmaStarEigenValues.get(j);
                    // accumulate second term of g1
                    g1 += ((firstDerivative * evmI.eigenValue * evmI.multiplicity * evmJ.eigenValue
                            * evmJ.multiplicity) / (evmI.eigenValue - evmJ.eigenValue));
                }
            }
        }

        expectedEpsilon = epsilonD + g1 / (totalN - rank);

    } else {
        // calculate the expected value of the epsilon estimate
        // see Muller, Edwards, Taylor (2004) for details

        // build a vector with each eigen value repeated per its multiplicity
        ArrayRealVector eigenColumnVector = new ArrayRealVector();
        for (EigenValueMultiplicityPair evmp : distinctSigmaStarEigenValues) {
            // there is probably a more memory-efficient method to do this.
            eigenColumnVector = eigenColumnVector
                    .append(new ArrayRealVector((int) evmp.multiplicity, evmp.eigenValue));
        }

        RealMatrix outerProduct = eigenColumnVector.outerProduct(eigenColumnVector);
        double sum = outerProduct.walkInOptimizedOrder(new SummationVisitor());
        double nu = (totalN - rank);
        double expT1 = (2 * nu * sumLambdaSquared) + (nu * nu * sumLambda * sumLambda);
        double expT2 = (nu * (nu + 1) * sumLambdaSquared) + (nu * sum * sum);
        expectedEpsilon = (1 / (double) rankU) * (expT1 / expT2);
    }
    // ensure that expected value is within bounds 1/b to 1
    if (!Double.isNaN(expectedEpsilon)) {
        if (expectedEpsilon < 1 / rankU) {
            expectedEpsilon = 1 / rankU;
        } else if (expectedEpsilon > 1) {
            expectedEpsilon = 1;
        }
    }
}

---