Examples with TetradMatrix edu.cmu.tetrad.util.TetradMatrix used on opensource projects

Example 56 with TetradMatrix

use of edu.cmu.tetrad.util.TetradMatrix in project tetrad by cmu-phil.

the class Ling method simulateCyclic.

// used to produce dataset if one is not provided as the input to the constructor
private static TetradMatrix simulateCyclic(GraphWithParameters dwp, TetradVector errorCoefficients, int n, Distribution distribution) {
    TetradMatrix reducedForm = reducedForm(dwp);
    TetradMatrix vectors = new TetradMatrix(dwp.getGraph().getNumNodes(), n);
    for (int j = 0; j < n; j++) {
        TetradVector vector = simulateReducedForm(reducedForm, errorCoefficients, distribution);
        vectors.assignColumn(j, vector);
    }
    return vectors;
}

Example 57 with TetradMatrix

use of edu.cmu.tetrad.util.TetradMatrix in project tetrad by cmu-phil.

the class Ling method reducedForm.

// graph matrix is B
// mixing matrix (reduced form) is A
private static TetradMatrix reducedForm(GraphWithParameters graph) {
    TetradMatrix graphMatrix = new TetradMatrix(graph.getGraphMatrix().getDoubleData().toArray());
    int n = graphMatrix.rows();
    // TetradMatrix identityMinusGraphTetradMatrix = TetradMatrixUtils.linearCombination(TetradMatrixUtils.identityTetradMatrix(n), 1, graphTetradMatrix, -1);
    TetradMatrix identityMinusGraphTetradMatrix = TetradMatrix.identity(n).minus(graphMatrix);
    return identityMinusGraphTetradMatrix.inverse();
}

Example 58 with TetradMatrix

use of edu.cmu.tetrad.util.TetradMatrix in project tetrad by cmu-phil.

the class Ling method search.

/**
 * The search method is used to process LiNG. Call search when you want to run the algorithm.
 */
public StoredGraphs search() {
    TetradMatrix A, W;
    StoredGraphs graphs = new StoredGraphs();
    try {
        long sTime = (new Date()).getTime();
        boolean fastIca = true;
        if (fastIca) {
            W = getWFastIca();
            System.out.println("W = " + W);
            // this is the heart of our method:
            graphs = findCandidateModels(dataSet.getVariables(), W, true);
        } else {
            double zeta = 1;
            final List<Mapping> allMappings = createMappings(null, null, dataSet.getNumColumns());
            double[][] _w = estimateW(new TetradMatrix(dataSet.getDoubleData().toArray()), dataSet.getNumColumns(), -zeta, zeta, allMappings);
            W = new TetradMatrix(_w);
            System.out.println("W = " + W);
            // this is the heart of our method:
            graphs = findCandidateModel(dataSet.getVariables(), W, true);
        }
        elapsedTime = (new Date()).getTime() - sTime;
    } catch (Exception e) {
        e.printStackTrace();
    }
    return graphs;
}

Example 59 with TetradMatrix

use of edu.cmu.tetrad.util.TetradMatrix in project tetrad by cmu-phil.

the class Ling method pruneEdgesByResampling.

// ==============================PRIVATE METHODS====================//
/**
 * This is the method used in Patrik's code.
 */
public TetradMatrix pruneEdgesByResampling(TetradMatrix data) {
    TetradMatrix X = new TetradMatrix(data.transpose().toArray());
    int npieces = 10;
    int cols = X.columns();
    int rows = X.rows();
    int piecesize = (int) Math.floor(cols / npieces);
    List<TetradMatrix> bpieces = new ArrayList<>();
    List<TetradVector> diststdpieces = new ArrayList<>();
    List<TetradVector> cpieces = new ArrayList<>();
    for (int p = 0; p < npieces; p++) {
        // % Select subset of data, and permute the variables to the causal order
        // Xp = X(k,((p-1)*piecesize+1):(p*piecesize));
        int p0 = (p) * piecesize;
        int p1 = (p + 1) * piecesize - 1;
        int[] range = range(p0, p1);
        TetradMatrix Xp = X;
        // % Remember to subract out the mean
        // Xpm = mean(Xp,2);
        // Xp = Xp - Xpm*ones(1,size(Xp,2));
        // 
        // % Calculate covariance matrix
        // cov = (Xp*Xp')/size(Xp,2);
        TetradVector Xpm = new TetradVector(rows);
        for (int i = 0; i < rows; i++) {
            double sum = 0.0;
            for (int j = 0; j < Xp.columns(); j++) {
                sum += Xp.get(i, j);
            }
            Xpm.set(i, sum / Xp.columns());
        }
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < Xp.columns(); j++) {
                Xp.set(i, j, Xp.get(i, j) - Xpm.get(i));
            }
        }
        TetradMatrix Xpt = Xp.transpose();
        TetradMatrix cov = Xp.times(Xpt);
        for (int i = 0; i < cov.rows(); i++) {
            for (int j = 0; j < cov.columns(); j++) {
                cov.set(i, j, cov.get(i, j) / Xp.columns());
            }
        }
        // % Do QL decomposition on the inverse square root of cov
        // [Q,L] = tridecomp(cov^(-0.5),'ql');
        boolean posDef = LingUtils.isPositiveDefinite(cov);
        if (!posDef) {
            System.out.println("Covariance matrix is not positive definite.");
        }
        TetradMatrix sqrt = cov.sqrt();
        ;
        TetradMatrix I = TetradMatrix.identity(rows);
        TetradMatrix AI = I.copy();
        TetradMatrix invSqrt = sqrt.inverse();
        QRDecomposition qr = new QRDecomposition(invSqrt.getRealMatrix());
        RealMatrix r = qr.getR();
        // % The estimated disturbance-stds are one over the abs of the diag of L
        // newestdisturbancestd = 1./diag(abs(L));
        TetradVector newestdisturbancestd = new TetradVector(rows);
        for (int t = 0; t < rows; t++) {
            newestdisturbancestd.set(t, 1.0 / Math.abs(r.getEntry(t, t)));
        }
        // 
        for (int s = 0; s < rows; s++) {
            for (int t = 0; t < min(s, cols); t++) {
                r.setEntry(s, t, r.getEntry(s, t) / r.getEntry(s, s));
            }
        }
        // % Calculate corresponding B
        // bnewest = eye(dims)-L;
        TetradMatrix bnewest = TetradMatrix.identity(rows);
        bnewest = bnewest.minus(new TetradMatrix(r));
        TetradVector cnewest = new TetradMatrix(r).times(Xpm);
        bpieces.add(bnewest);
        diststdpieces.add(newestdisturbancestd);
        cpieces.add(cnewest);
    }
    // 
    // for i=1:dims,
    // for j=1:dims,
    // 
    // themean = mean(Bpieces(i,j,:));
    // thestd = std(Bpieces(i,j,:));
    // if abs(themean)<prunefactor*thestd,
    // Bfinal(i,j) = 0;
    // else
    // Bfinal(i,j) = themean;
    // end
    // 
    // end
    // end
    TetradMatrix means = new TetradMatrix(rows, rows);
    TetradMatrix stds = new TetradMatrix(rows, rows);
    TetradMatrix BFinal = new TetradMatrix(rows, rows);
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < rows; j++) {
            double sum = 0.0;
            for (int y = 0; y < npieces; y++) {
                sum += bpieces.get(y).get(i, j);
            }
            double themean = sum / (npieces);
            double sumVar = 0.0;
            for (int y = 0; y < npieces; y++) {
                sumVar += Math.pow((bpieces.get(y).get(i, j)) - themean, 2);
            }
            double thestd = Math.sqrt(sumVar / (npieces));
            means.set(i, j, themean);
            stds.set(i, j, thestd);
            if (Math.abs(themean) < threshold * thestd) {
                // getPruneFactor() * thestd) {
                BFinal.set(i, j, 0);
            } else {
                BFinal.set(i, j, themean);
            }
        }
    }
    return BFinal;
}

Example 60 with TetradMatrix

use of edu.cmu.tetrad.util.TetradMatrix in project tetrad by cmu-phil.

the class LingamPattern method search.

public Graph search() {
    long initialTime = System.currentTimeMillis();
    Graph _pattern = GraphUtils.bidirectedToUndirected(getPattern());
    TetradLogger.getInstance().log("info", "Making list of all dags in pattern...");
    List<Graph> dags = SearchGraphUtils.getAllGraphsByDirectingUndirectedEdges(_pattern);
    TetradLogger.getInstance().log("normalityTests", "Anderson Darling P value for Variables\n");
    NumberFormat nf = new DecimalFormat("0.0000");
    if (dags.isEmpty()) {
        return null;
    }
    TetradMatrix data = getDataSet().getDoubleData();
    List<Node> variables = getDataSet().getVariables();
    if (dags.size() == 0) {
        throw new IllegalArgumentException("The data set is empty.");
    }
    // Check that all the daga and the data contain the same variables.
    List<Score> scores = new ArrayList<>();
    for (Graph dag : dags) {
        scores.add(getScore(dag, data, variables));
    }
    double maxScore = 0.0;
    int maxj = -1;
    for (int j = 0; j < dags.size(); j++) {
        double _score = scores.get(j).score;
        if (_score > maxScore) {
            maxScore = _score;
            maxj = j;
        }
    }
    Graph dag = dags.get(maxj);
    this.bestDag = new EdgeListGraph(dags.get(maxj));
    this.pValues = scores.get(maxj).pvals;
    TetradLogger.getInstance().log("graph", "winning dag = " + dag);
    TetradLogger.getInstance().log("normalityTests", "Anderson Darling P value for Residuals\n");
    for (int j = 0; j < getDataSet().getNumColumns(); j++) {
        TetradLogger.getInstance().log("normalityTests", getDataSet().getVariable(j) + ": " + nf.format(scores.get(maxj).pvals[j]));
    }
    // System.out.println();
    Graph ngDagPattern = SearchGraphUtils.patternFromDag(dag);
    List<Node> nodes = ngDagPattern.getNodes();
    for (Edge edge : ngDagPattern.getEdges()) {
        Node node1 = edge.getNode1();
        Node node2 = edge.getNode2();
        double p1 = getPValues()[nodes.indexOf(node1)];
        double p2 = getPValues()[nodes.indexOf(node2)];
        boolean node1Nongaussian = p1 < getAlpha();
        boolean node2Nongaussian = p2 < getAlpha();
        if (node1Nongaussian || node2Nongaussian) {
            if (!Edges.isUndirectedEdge(edge)) {
                continue;
            }
            ngDagPattern.removeEdge(edge);
            ngDagPattern.addEdge(dag.getEdge(node1, node2));
            if (node1Nongaussian) {
                TetradLogger.getInstance().log("edgeOrientations", node1 + " nongaussian ");
            }
            if (node2Nongaussian) {
                TetradLogger.getInstance().log("edgeOrientations", node2 + " nongaussian ");
            }
            TetradLogger.getInstance().log("nongaussianOrientations", "Nongaussian orientation: " + dag.getEdge(node1, node2));
        }
    }
    // System.out.println();
    // 
    // System.out.println("Applying Meek rules.");
    // System.out.println();
    new MeekRules().orientImplied(ngDagPattern);
    this.ngDagPattern = ngDagPattern;
    TetradLogger.getInstance().log("graph", "Returning: " + ngDagPattern);
    return ngDagPattern;
}