Search in sources :

Example 11 with SingularMatrixException

use of org.apache.commons.math3.linear.SingularMatrixException in project tetrad by cmu-phil.

the class SemBicScoreDeterministic method localScore.

/**
 * Calculates the sample likelihood and BIC score for i given its parents in a simple SEM model
 */
public double localScore(int i, int... parents) {
    for (int p : parents) if (forbidden.contains(p))
        return Double.NaN;
    double small = getDeterminismThreshold();
    double s2 = getCovariances().getValue(i, i);
    int p = parents.length;
    TetradMatrix covxx = getSelection(getCovariances(), parents, parents);
    TetradVector covxy = getSelection(getCovariances(), parents, new int[] { i }).getColumn(0);
    try {
        s2 -= covxx.inverse().times(covxy).dotProduct(covxy);
    } catch (SingularMatrixException e) {
        s2 = 0;
    }
    // System.out.println(s2);
    int n = getSampleSize();
    int k = 2 * p + 1;
    if (s2 < small) {
        s2 = 0;
    }
    if (s2 == 0) {
        printDeterminism(i, parents);
        return Double.NaN;
    }
    return -(n) * log(s2) - getPenaltyDiscount() * k * log(n);
}
Also used : TetradVector(edu.cmu.tetrad.util.TetradVector) SingularMatrixException(org.apache.commons.math3.linear.SingularMatrixException) TetradMatrix(edu.cmu.tetrad.util.TetradMatrix)

Example 12 with SingularMatrixException

use of org.apache.commons.math3.linear.SingularMatrixException in project tetrad by cmu-phil.

the class DeltaSextadTest method calcChiSquare.

/**
 * Takes a list of tetrads for the given data set and returns the chi square value for the test. We assume that the
 * tetrads are non-redundant; if not, a matrix exception will be thrown.
 * <p>
 * Calculates the T statistic (Bollen and Ting, p. 161). This is significant if tests as significant using the Chi
 * Square distribution with degrees of freedom equal to the number of nonredundant tetrads tested.
 */
public double calcChiSquare(IntSextad[] sextads) {
    Set<Sigma> boldSigmaSet = new HashSet<>();
    for (IntSextad sextad : sextads) {
        List<Integer> _nodes = sextad.getNodes();
        // 
        for (int k1 = 0; k1 < 3; k1++) {
            for (int k2 = 0; k2 < 3; k2++) {
                boldSigmaSet.add(new Sigma(_nodes.get(k1), _nodes.get(3 + k2)));
            }
        }
    }
    List<Sigma> boldSigma = new ArrayList<>(boldSigmaSet);
    // Need a matrix of variances and covariances of sample covariances.
    TetradMatrix sigma_ss = new TetradMatrix(boldSigma.size(), boldSigma.size());
    for (int i = 0; i < boldSigma.size(); i++) {
        for (int j = i; j < boldSigma.size(); j++) {
            Sigma sigmaef = boldSigma.get(i);
            Sigma sigmagh = boldSigma.get(j);
            int e = sigmaef.getA();
            int f = sigmaef.getB();
            int g = sigmagh.getA();
            int h = sigmagh.getB();
            if (cov != null && cov instanceof CorrelationMatrix) {
                // Assumes multinormality. Using formula 23. (Not implementing formula 22 because that case
                // does not come up.)
                double rr = 0.5 * (r(e, f) * r(g, h)) * (r(e, g) * r(e, g) + r(e, h) * r(e, h) + r(f, g) * r(f, g) + r(f, h) * r(f, h)) + r(e, g) * r(f, h) + r(e, h) * r(f, g) - r(e, f) * (r(f, g) * r(f, h) + r(e, g) * r(e, h)) - r(g, h) * (r(f, g) * r(e, g) + r(f, h) * r(e, h));
                // General.
                double rr2 = r(e, f, g, h) + 0.25 * r(e, f) * r(g, h) * (r(e, e, g, g) * r(f, f, g, g) + r(e, e, h, h) + r(f, f, h, h)) - 0.5 * r(e, f) * (r(e, e, g, h) + r(f, f, g, h)) - 0.5 * r(g, h) * (r(e, f, g, g) + r(e, f, h, h));
                sigma_ss.set(i, j, rr);
                sigma_ss.set(j, i, rr);
            } else if (cov != null && data == null) {
                // Assumes multinormality--see p. 160.
                // double _ss = r(e, g) * r(f, h) + r(e, h) * r(f, g); // + or -? Different advise. + in the code.
                double _ss = r(e, g) * r(f, h) + r(e, h) * r(f, g);
                sigma_ss.set(i, j, _ss);
                sigma_ss.set(j, i, _ss);
            } else {
                double _ss = r(e, f, g, h) - r(e, f) * r(g, h);
                sigma_ss.set(i, j, _ss);
                sigma_ss.set(j, i, _ss);
            }
        }
    }
    // Need a matrix of of population estimates of partial derivatives of tetrads
    // with respect to covariances in boldSigma.
    TetradMatrix del = new TetradMatrix(boldSigma.size(), sextads.length);
    for (int j = 0; j < sextads.length; j++) {
        IntSextad sextad = sextads[j];
        for (int i = 0; i < boldSigma.size(); i++) {
            Sigma sigma = boldSigma.get(i);
            double derivative = getDerivative(sextad, sigma);
            del.set(i, j, derivative);
        }
    }
    // Need a vector of population estimates of the sextads.
    TetradMatrix t = new TetradMatrix(sextads.length, 1);
    for (int i = 0; i < sextads.length; i++) {
        IntSextad sextad = sextads[i];
        List<Integer> nodes = sextad.getNodes();
        TetradMatrix m = new TetradMatrix(3, 3);
        for (int k1 = 0; k1 < 3; k1++) {
            for (int k2 = 0; k2 < 3; k2++) {
                m.set(k1, k2, r(nodes.get(k1), nodes.get(3 + k2)));
            }
        }
        double det = m.det();
        t.set(i, 0, det);
    }
    TetradMatrix sigma_tt = del.transpose().times(sigma_ss).times(del);
    double chisq;
    try {
        chisq = N * t.transpose().times(sigma_tt.inverse()).times(t).get(0, 0);
    } catch (SingularMatrixException e) {
        throw new RuntimeException("Singularity problem.", e);
    }
    return chisq;
}
Also used : TetradMatrix(edu.cmu.tetrad.util.TetradMatrix) SingularMatrixException(org.apache.commons.math3.linear.SingularMatrixException)

Example 13 with SingularMatrixException

use of org.apache.commons.math3.linear.SingularMatrixException in project tetrad by cmu-phil.

the class IndTestFisherZ method partialCorrelation.

private double partialCorrelation(Node x, Node y, List<Node> z) throws SingularMatrixException {
    if (z.isEmpty()) {
        double a = covMatrix.getValue(indexMap.get(x), indexMap.get(y));
        double b = covMatrix.getValue(indexMap.get(x), indexMap.get(x));
        double c = covMatrix.getValue(indexMap.get(y), indexMap.get(y));
        if (b * c == 0)
            throw new SingularMatrixException();
        return -a / Math.sqrt(b * c);
    } else {
        int[] indices = new int[z.size() + 2];
        indices[0] = indexMap.get(x);
        indices[1] = indexMap.get(y);
        for (int i = 0; i < z.size(); i++) indices[i + 2] = indexMap.get(z.get(i));
        TetradMatrix submatrix = covMatrix.getSubmatrix(indices).getMatrix();
        return StatUtils.partialCorrelation(submatrix);
    }
}
Also used : SingularMatrixException(org.apache.commons.math3.linear.SingularMatrixException)

Example 14 with SingularMatrixException

use of org.apache.commons.math3.linear.SingularMatrixException in project tetrad by cmu-phil.

the class IndTestFisherZD method partialCorrelation.

private double partialCorrelation(Node x, Node y, List<Node> z) throws SingularMatrixException {
    if (z.isEmpty()) {
        double a = covMatrix.getValue(indexMap.get(x), indexMap.get(y));
        double b = covMatrix.getValue(indexMap.get(x), indexMap.get(x));
        double c = covMatrix.getValue(indexMap.get(y), indexMap.get(y));
        if (b * c == 0)
            throw new SingularMatrixException();
        return -a / Math.sqrt(b * c);
    } else {
        int[] indices = new int[z.size() + 2];
        indices[0] = indexMap.get(x);
        indices[1] = indexMap.get(y);
        for (int i = 0; i < z.size(); i++) indices[i + 2] = indexMap.get(z.get(i));
        TetradMatrix submatrix = covMatrix.getSubmatrix(indices).getMatrix();
        return StatUtils.partialCorrelation(submatrix);
    }
}
Also used : SingularMatrixException(org.apache.commons.math3.linear.SingularMatrixException) TetradMatrix(edu.cmu.tetrad.util.TetradMatrix)

Aggregations

SingularMatrixException (org.apache.commons.math3.linear.SingularMatrixException)13 RealMatrix (org.apache.commons.math3.linear.RealMatrix)7 ConvergenceException (org.apache.commons.math3.exception.ConvergenceException)4 RealVector (org.apache.commons.math3.linear.RealVector)4 TetradMatrix (edu.cmu.tetrad.util.TetradMatrix)3 Array2DRowRealMatrix (org.apache.commons.math3.linear.Array2DRowRealMatrix)3 ArrayRealVector (org.apache.commons.math3.linear.ArrayRealVector)3 Random (java.util.Random)2 MixtureMultivariateNormalDistribution (org.apache.commons.math3.distribution.MixtureMultivariateNormalDistribution)2 NormalDistribution (org.apache.commons.math3.distribution.NormalDistribution)2 MultivariateNormalMixtureExpectationMaximization (org.apache.commons.math3.distribution.fitting.MultivariateNormalMixtureExpectationMaximization)2 MaxCountExceededException (org.apache.commons.math3.exception.MaxCountExceededException)2 TooManyIterationsException (org.apache.commons.math3.exception.TooManyIterationsException)2 Optimum (org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum)2 LeastSquaresProblem (org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem)2 LevenbergMarquardtOptimizer (org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer)2 DiagonalMatrix (org.apache.commons.math3.linear.DiagonalMatrix)2 RandomGenerator (org.apache.commons.math3.random.RandomGenerator)2 Pair (org.apache.commons.math3.util.Pair)2 Node (edu.cmu.tetrad.graph.Node)1