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

Example 41 with TetradMatrix

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

the class CovMatrixSumWrapper method calcSum.

// public CovMatrixSumWrapper(SemEstimatorWrapper wrapper1, DataWrapper wrapper2) {
// if (wrapper1 == null || wrapper2 == null) {
// throw new NullPointerException("The data must not be null");
// }
// 
// DataModel model2 = wrapper2.getSelectedDataModel();
// 
// if (!(model2 instanceof ICovarianceMatrix)) {
// throw new IllegalArgumentException("Expecting corrariance matrices.");
// }
// 
// TetradMatrix corr1 = wrapper1.getEstimatedSemIm().getImplCovarMeas();
// TetradMatrix corr2 = ((ICovarianceMatrix) model2).getMatrix();
// 
// TetradMatrix corr3 = calcSum(corr1, corr2);
// 
// ICovarianceMatrix corrWrapper = new CovarianceMatrix(model2.getVariable(), corr3,
// ((ICovarianceMatrix) model2).getSampleSize());
// 
// setDataModel(corrWrapper);
// setSourceGraph(wrapper2.getSourceGraph());
// LogDataUtils.logDataModelList("Difference of matrices.", getDataModelList());
// 
// }
// public CovMatrixSumWrapper(SemImWrapper wrapper1, DataWrapper wrapper2) {
// try {
// if (wrapper1 == null || wrapper2 == null) {
// throw new NullPointerException("The data must not be null");
// }
// 
// DataModel model2 = wrapper2.getSelectedDataModel();
// 
// if (!(model2 instanceof ICovarianceMatrix)) {
// throw new IllegalArgumentException("Expecting corrariance matrices.");
// }
// 
// TetradMatrix corr1 = wrapper1.getSemIm().getImplCovarMeas();
// TetradMatrix corr2 = ((ICovarianceMatrix) model2).getMatrix();
// 
// TetradMatrix corr3 = calcSum(corr1, corr2);
// 
// ICovarianceMatrix corrWrapper = new CovarianceMatrix(model2.getVariable(), corr3,
// ((ICovarianceMatrix) model2).getSampleSize());
// 
// setDataModel(corrWrapper);
// setSourceGraph(wrapper2.getSourceGraph());
// LogDataUtils.logDataModelList("Difference of matrices.", getDataModelList());
// } catch (Exception e) {
// e.printStackTrace();
// throw new RuntimeException(e);
// }
// 
// }
private TetradMatrix calcSum(TetradMatrix corr1, TetradMatrix corr2) {
    if (corr1.rows() != corr2.rows()) {
        throw new IllegalArgumentException("Covariance matrices must be the same size.");
    }
    TetradMatrix corr3 = new TetradMatrix(corr2.rows(), corr2.rows());
    for (int i = 0; i < corr3.rows(); i++) {
        for (int j = 0; j < corr3.rows(); j++) {
            double v = corr1.get(i, j) + corr2.get(i, j);
            corr3.set(i, j, v);
            corr3.set(j, i, v);
        }
    }
    return corr3;
}

Example 42 with TetradMatrix

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

the class SubsetSelectedVariablesWrapper method createCovarianceModel.

private static DataModel createCovarianceModel(ICovarianceMatrix data) {
    int numSelected = 0;
    for (Node node : data.getVariables()) {
        if (data.isSelected(node)) {
            numSelected++;
        }
    }
    int[] selectedIndices = new int[numSelected];
    String[] nodeNames = new String[numSelected];
    int index = -1;
    for (int i = 0; i < data.getVariables().size(); i++) {
        Node node = data.getVariables().get(i);
        if (data.isSelected(node)) {
            ++index;
            selectedIndices[index] = i;
            nodeNames[index] = node.getName();
        }
    }
    TetradMatrix matrix = data.getMatrix();
    TetradMatrix newMatrix = matrix.getSelection(selectedIndices, selectedIndices).copy();
    ICovarianceMatrix newCov = new CovarianceMatrix(DataUtils.createContinuousVariables(nodeNames), newMatrix, matrix.rows());
    return newCov;
}

Example 43 with TetradMatrix

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

the class FastIca method findComponents.

/**
 * Runs the Fast ICA algorithm (following the R version) and returns the
 * list of result items that the R version returns.
 *
 * @return this list, as an FastIca.IcaResult object.
 */
public IcaResult findComponents() {
    int n = X.rows();
    int p = X.columns();
    if (numComponents > Math.min(n, p)) {
        TetradLogger.getInstance().log("info", "Requested number of components is too large.");
        TetradLogger.getInstance().log("info", "Reset to " + Math.min(n, p));
        numComponents = Math.min(n, p);
    }
    if (wInit == null) {
        wInit = new TetradMatrix(numComponents, numComponents);
        for (int i = 0; i < wInit.rows(); i++) {
            for (int j = 0; j < wInit.columns(); j++) {
                wInit.set(i, j, RandomUtil.getInstance().nextNormal(0, 1));
            }
        }
    } else if (wInit.rows() != wInit.columns()) {
        throw new IllegalArgumentException("wInit is the wrong size.");
    }
    if (verbose) {
        TetradLogger.getInstance().log("info", "Centering");
    }
    X = center(X);
    if (colNorm) {
        X = scale(X);
    }
    X = X.transpose();
    if (verbose) {
        TetradLogger.getInstance().log("info", "Whitening");
    }
    TetradMatrix V = X.times(X.transpose()).scalarMult(1.0 / n);
    // v.scalarMult(1.0 / n);
    SingularValueDecomposition s = new SingularValueDecomposition(V.getRealMatrix());
    TetradMatrix D = new TetradMatrix(s.getS());
    TetradMatrix U = new TetradMatrix(s.getU());
    for (int i = 0; i < D.rows(); i++) {
        D.set(i, i, 1.0 / Math.sqrt(D.get(i, i)));
    }
    TetradMatrix K = D.times(U.transpose());
    // This SVD gives -U from R's SVD.
    K = K.scalarMult(-1);
    K = K.getPart(0, numComponents - 1, 0, p - 1);
    TetradMatrix X1 = K.times(X);
    TetradMatrix b;
    if (algorithmType == DEFLATION) {
        b = icaDeflation(X1, numComponents, tolerance, function, alpha, maxIterations, verbose, wInit);
    } else if (algorithmType == PARALLEL) {
        b = icaParallel(X1, numComponents, tolerance, function, alpha, maxIterations, verbose, wInit);
    } else {
        throw new IllegalStateException();
    }
    TetradMatrix w = b.times(K);
    TetradMatrix S = w.times(X);
    TetradMatrix A = w.transpose().times(w.times(w.transpose()).inverse());
    return new IcaResult(X.transpose(), K.transpose(), b.transpose(), A.transpose(), S.transpose());
}

Example 44 with TetradMatrix

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

the class FastIca method icaDeflation.

// ==============================PRIVATE METHODS==========================//
private TetradMatrix icaDeflation(TetradMatrix X, int numComponents, double tolerance, int function, double alpha, int maxIterations, boolean verbose, TetradMatrix wInit) {
    if (verbose && function == LOGCOSH) {
        TetradLogger.getInstance().log("info", "Deflation FastIca using lgcosh approx. to neg-entropy function");
    }
    if (verbose && function == EXP) {
        TetradLogger.getInstance().log("info", "Deflation FastIca using exponential approx. to neg-entropy function");
    }
    int p = X.columns();
    TetradMatrix W = new TetradMatrix(numComponents, numComponents);
    for (int i = 0; i < numComponents; i++) {
        if (verbose) {
            TetradLogger.getInstance().log("fastIcaDetails", "Component " + (i + 1));
        }
        TetradVector w = wInit.getRow(i);
        if (i > 0) {
            TetradVector t = w.like();
            for (int u = 0; u < i; u++) {
                double k = 0.0;
                for (int j = 0; j < numComponents; j++) {
                    k += w.get(j) * W.get(u, j);
                }
                for (int j = 0; j < numComponents; j++) {
                    t.set(j, t.get(j) + k * W.get(u, j));
                }
            }
            for (int j = 0; j < numComponents; j++) {
                w.set(j, w.get(j) - t.get(j));
            }
        }
        double rms = rms(w);
        for (int j = 0; j < numComponents; j++) {
            w.set(j, w.get(j) / rms);
        }
        int it = 0;
        double _tolerance = Double.POSITIVE_INFINITY;
        if (function == LOGCOSH) {
            while (_tolerance > tolerance && ++it <= maxIterations) {
                TetradVector wx = X.transpose().times(w);
                TetradVector gwx0 = new TetradVector(p);
                for (int j = 0; j < p; j++) {
                    gwx0.set(j, Math.tanh(alpha * wx.get(j)));
                }
                TetradMatrix gwx = new TetradMatrix(numComponents, p);
                for (int _i = 0; _i < numComponents; _i++) {
                    for (int j = 0; j < p; j++) {
                        gwx.set(_i, j, gwx0.get(j));
                    }
                }
                TetradMatrix xgwx = new TetradMatrix(numComponents, p);
                for (int _i = 0; _i < numComponents; _i++) {
                    for (int j = 0; j < p; j++) {
                        xgwx.set(_i, j, X.get(_i, j) * gwx0.get(j));
                    }
                }
                TetradVector v1 = new TetradVector(numComponents);
                for (int k = 0; k < numComponents; k++) {
                    v1.set(k, mean(xgwx.getRow(k)));
                }
                TetradVector g_wx = new TetradVector(p);
                for (int k = 0; k < p; k++) {
                    double tmp1 = Math.tanh(alpha * wx.get(k));
                    g_wx.set(k, alpha * (1.0 - tmp1 * tmp1));
                }
                TetradVector v2 = w.copy();
                double meanGwx = mean(g_wx);
                v2 = v2.scalarMult(meanGwx);
                TetradVector w1 = v1.copy();
                // w1.assign(v2, PlusMult.plusMult(-1));
                w1 = w1.minus(v2);
                if (i > 0) {
                    TetradVector t = w1.like();
                    for (int u = 0; u < i; u++) {
                        double k = 0.0;
                        for (int j = 0; j < numComponents; j++) {
                            k += w1.get(j) * W.get(u, j);
                        }
                        for (int j = 0; j < numComponents; j++) {
                            t.set(j, t.get(j) + k * W.get(u, j));
                        }
                    }
                    for (int j = 0; j < numComponents; j++) {
                        w1.set(j, w1.get(j) - t.get(j));
                    }
                }
                double _rms = rms(w1);
                for (int k = 0; k < numComponents; k++) {
                    w1.set(k, w1.get(k) / _rms);
                }
                _tolerance = 0.0;
                for (int k = 0; k < numComponents; k++) {
                    _tolerance += w1.get(k) * w.get(k);
                }
                _tolerance = Math.abs(Math.abs(_tolerance) - 1.0);
                if (verbose) {
                    TetradLogger.getInstance().log("fastIcaDetails", "Iteration " + it + " tol = " + _tolerance);
                }
                w = w1;
            }
        } else if (function == EXP) {
            while (_tolerance > tolerance && ++it <= maxIterations) {
                TetradVector wx = X.transpose().times(w);
                TetradVector gwx0 = new TetradVector(p);
                for (int j = 0; j < p; j++) {
                    gwx0.set(j, wx.get(j) * Math.exp(-(wx.get(j) * wx.get(j)) / 2));
                }
                TetradMatrix gwx = new TetradMatrix(numComponents, p);
                for (int _i = 0; _i < numComponents; _i++) {
                    for (int j = 0; j < p; j++) {
                        gwx.set(_i, j, gwx0.get(j));
                    }
                }
                TetradMatrix xgwx = new TetradMatrix(numComponents, p);
                for (int _i = 0; _i < numComponents; _i++) {
                    for (int j = 0; j < p; j++) {
                        xgwx.set(_i, j, X.get(_i, j) * gwx0.get(j));
                    }
                }
                TetradVector v1 = new TetradVector(numComponents);
                for (int k = 0; k < numComponents; k++) {
                    v1.set(k, mean(xgwx.getRow(k)));
                }
                TetradVector g_wx = new TetradVector(p);
                for (int j = 0; j < p; j++) {
                    g_wx.set(j, (1.0 - wx.get(j) * wx.get(j)) * Math.exp(-(wx.get(j) * wx.get(j)) / 2));
                }
                TetradVector v2 = w.copy();
                double meanGwx = mean(g_wx);
                TetradVector w1 = v2.scalarMult(meanGwx).minus(v2);
                if (i > 0) {
                    TetradVector t = w1.like();
                    for (int u = 0; u < i; u++) {
                        double k = 0.0;
                        for (int j = 0; j < numComponents; j++) {
                            k += w1.get(j) * W.get(u, j);
                        }
                        for (int j = 0; j < numComponents; j++) {
                            t.set(j, t.get(j) + k * W.get(u, j));
                        }
                    }
                    for (int j = 0; j < numComponents; j++) {
                        w1.set(j, w1.get(j) - t.get(j));
                    }
                }
                double _rms = rms(w1);
                for (int k = 0; k < numComponents; k++) {
                    w1.set(k, w1.get(k) / _rms);
                }
                _tolerance = 0.0;
                for (int k = 0; k < numComponents; k++) {
                    _tolerance += w1.get(k) * w.get(k);
                }
                _tolerance = Math.abs(Math.abs(_tolerance) - 1.0);
                if (verbose) {
                    TetradLogger.getInstance().log("fastIcaDetails", "Iteration " + it + " tol = " + _tolerance);
                }
                w = w1;
            }
        }
        W.assignRow(i, w);
    }
    return W;
}

Example 45 with TetradMatrix

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

the class DeltaTetradTest2 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(Tetrad... tetrads) {
    this.df = tetrads.length;
    // Need a list of symbolic covariances--i.e. covariances that appear in tetrads.
    Set<Sigma> boldSigmaSet = new LinkedHashSet<>();
    List<Sigma> boldSigma = new ArrayList<>();
    for (Tetrad tetrad : tetrads) {
        boldSigmaSet.add(new Sigma(tetrad.getI(), tetrad.getK()));
        boldSigmaSet.add(new Sigma(tetrad.getI(), tetrad.getL()));
        boldSigmaSet.add(new Sigma(tetrad.getJ(), tetrad.getK()));
        boldSigmaSet.add(new Sigma(tetrad.getJ(), tetrad.getL()));
    }
    for (Sigma sigma : boldSigmaSet) {
        boldSigma.add(sigma);
    }
    // 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 = 0; j < boldSigma.size(); j++) {
            Sigma sigmaef = boldSigma.get(i);
            Sigma sigmagh = boldSigma.get(j);
            Node e = sigmaef.getA();
            Node f = sigmaef.getB();
            Node g = sigmagh.getA();
            Node 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 * (sxy(e, f) * sxy(g, h)) * (sxy(e, g) * sxy(e, g) + sxy(e, h) * sxy(e, h) + sxy(f, g) * sxy(f, g) + sxy(f, h) * sxy(f, h)) + sxy(e, g) * sxy(f, h) + sxy(e, h) * sxy(f, g) - sxy(e, f) * (sxy(f, g) * sxy(f, h) + sxy(e, g) * sxy(e, h)) - sxy(g, h) * (sxy(f, g) * sxy(e, g) + sxy(f, h) * sxy(e, h));
                sigma_ss.set(i, j, rr);
            } else if (cov != null && dataSet == null) {
                // Assumes multinormality--see p. 160.
                // + or -? Different advise. + in the code.
                double _ss = sxy(e, g) * sxy(f, h) - sxy(e, h) * sxy(f, g);
                sigma_ss.set(i, j, _ss);
            } else {
                double _ss = sxyzw(e, f, g, h) - sxy(e, f) * sxy(g, h);
                sigma_ss.set(i, j, _ss);
            }
        }
    }
    // Need a matrix of of population estimates of partial derivatives of tetrads
    // with respect to covariances in boldSigma.w
    TetradMatrix del = new TetradMatrix(boldSigma.size(), tetrads.length);
    for (int i = 0; i < boldSigma.size(); i++) {
        for (int j = 0; j < tetrads.length; j++) {
            Sigma sigma = boldSigma.get(i);
            Tetrad tetrad = tetrads[j];
            Node e = tetrad.getI();
            Node f = tetrad.getJ();
            Node g = tetrad.getK();
            Node h = tetrad.getL();
            double derivative = getDerivative(e, f, g, h, sigma.getA(), sigma.getB());
            del.set(i, j, derivative);
        }
    }
    // Need a vector of population estimates of the tetrads.
    TetradMatrix t = new TetradMatrix(tetrads.length, 1);
    for (int i = 0; i < tetrads.length; i++) {
        Tetrad tetrad = tetrads[i];
        Node e = tetrad.getI();
        Node f = tetrad.getJ();
        Node g = tetrad.getK();
        Node h = tetrad.getL();
        double d1 = sxy(e, f);
        double d2 = sxy(g, h);
        double d3 = sxy(e, g);
        double d4 = sxy(f, h);
        double value = d1 * d2 - d3 * d4;
        t.set(i, 0, value);
    }
    // Now multiply to get Sigma_tt
    TetradMatrix w1 = del.transpose().times(sigma_ss);
    TetradMatrix sigma_tt = w1.times(del);
    // And now invert and multiply to get T.
    TetradMatrix v0 = sigma_tt.inverse();
    TetradMatrix v1 = t.transpose().times(v0);
    TetradMatrix v2 = v1.times(t);
    double chisq = N * v2.get(0, 0);
    this.chisq = chisq;
    return chisq;
}