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

Example 16 with TetradMatrix

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

the class KernelUtils method constructH.

/**
 * Constructs the projection matrix on 1/m
 *
 * @param m the sample size
 */
public static TetradMatrix constructH(int m) {
    TetradMatrix H = new TetradMatrix(m, m);
    double od = -1.0 / (double) m;
    double d = od + 1;
    for (int i = 0; i < m; i++) {
        for (int j = i; j < m; j++) {
            if (i == j) {
                H.set(i, j, d);
            } else {
                H.set(i, j, od);
            }
        }
    }
    return H;
}

Example 17 with TetradMatrix

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

the class KernelUtils method constructGramMatrix.

/**
 * Constructs Gram matrix for a given vector valued sample. The set of kernels corresponds to the variables in the
 * set. The output matrix is the tensor product of Gram matrices for each variable.
 *
 * @param kernels the kernels for each variable
 * @param dataset the dataset containing each variable
 * @param nodes   the variables to construct the Gram matrix for
 */
public static TetradMatrix constructGramMatrix(List<Kernel> kernels, DataSet dataset, List<Node> nodes) {
    int m = dataset.getNumRows();
    TetradMatrix gram = new TetradMatrix(m, m);
    for (int k = 0; k < nodes.size(); k++) {
        Node node = nodes.get(k);
        int col = dataset.getColumn(node);
        Kernel kernel = kernels.get(k);
        for (int i = 0; i < m; i++) {
            for (int j = i; j < m; j++) {
                double keval = kernel.eval(dataset.getDouble(i, col), dataset.getDouble(j, col));
                if (k == 0) {
                    gram.set(i, j, keval);
                } else {
                    keval *= gram.get(i, j);
                    gram.set(i, j, keval);
                }
            }
        }
    }
    return gram;
}

Example 18 with TetradMatrix

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

the class KernelUtils method constructCentralizedGramMatrix.

/**
 * Constructs the centralized Gram matrix for a given vector valued sample. The set of kernels corresponds to the
 * variables in the set. The output matrix is the tensor product of Gram matrices for each variable.
 *
 * @param kernels the kernels for each variable
 * @param dataset the dataset containing each variable
 * @param nodes   the variables to construct the Gram matrix for
 */
public static TetradMatrix constructCentralizedGramMatrix(List<Kernel> kernels, DataSet dataset, List<Node> nodes) {
    int m = dataset.getNumRows();
    TetradMatrix gram = constructGramMatrix(kernels, dataset, nodes);
    TetradMatrix H = constructH(m);
    TetradMatrix KH = gram.times(H);
    TetradMatrix HKH = H.times(KH);
    return HKH;
}

Example 19 with TetradMatrix

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

the class SemEstimatorGibbs method estimate.

// ==============================PUBLIC METHODS=========================//
/**
 * Runs the estimator on the data and SemPm passed in through the
 * constructor.
 */
public void estimate() {
    // dogibbs in pascal
    // In the comments, getgibsprefs, PRIORINIT, GIBBSINIT, FORMAPPROXDIST,
    // DRAWFROMAPPROX refer to procedure in the Pascal version from which
    // this was adapted.  The same is true of the private methods such
    // as brent, neglogpost, etc.
    // Initialize method variables
    List<Parameter> parameters = this.semPm.getParameters();
    int numParameters = parameters.size();
    double[][] parameterCovariances = new double[numParameters][numParameters];
    this.parameterMeans = new double[numParameters];
    this.paramConstraints = new ParamConstraint[numParameters];
    TetradMatrix data = new TetradMatrix(parameters.size(), numIterations / 50);
    // PRIORINIT
    if (flatPrior) {
        // this is used to construct the prior covariance matrix, means
        for (int i = 0; i < numParameters; i++) {
            Parameter param = parameters.get(i);
            this.parameterMeans[i] = (param.isFixed()) ? 0.0 : this.priorVariance;
            // Default parameter constraints.  The user should have the
            // option to change these via the GUI
            paramConstraints[i] = // ParamType.VAR = 'Error Variance'
            (param.getType() == ParamType.VAR) ? new ParamConstraint(this.startIm, param, ParamConstraintType.GT, 0.0) : new ParamConstraint(this.startIm, param, ParamConstraintType.NONE, 0.0);
            for (int j = 0; j < numParameters; j++) {
                parameterCovariances[i][j] = (i == j && !param.isFixed()) ? this.priorVariance : 0.0;
            }
        }
        this.priorCov = new TetradMatrix(parameterCovariances);
    } else {
        System.out.println("Informative Prior. Exiting.");
        return;
    }
    // END PRIORINIT
    // GIBBSINIT
    SemIm posteriorIm = new SemIm(this.startIm);
    List postFreeParams = posteriorIm.getFreeParameters();
    System.out.println("entering main loop");
    for (int iter = 1; iter <= this.numIterations; iter++) {
        System.out.println(iter);
        for (int param = 0; param < postFreeParams.size(); param++) {
            Parameter p = parameters.get(param);
            ParamConstraint constraint = this.paramConstraints[param];
            if (!p.isFixed()) {
                // FORMAPPROXDIST begin
                double number = (constraint.getParam2() == null) ? constraint.getNumber() : this.startIm.getParamValue(constraint.getParam2());
                double ax, bx, cx;
                // Mark - these constraints follow pascal code
                if (constraint.getType() == ParamConstraintType.NONE) {
                    ax = -500.0;
                    bx = 0.0;
                    cx = 500.0;
                } else if (constraint.getType() == ParamConstraintType.GT) {
                    ax = number;
                    cx = number + 500.0;
                    bx = (ax + cx) / 2.0;
                } else if (constraint.getType() == ParamConstraintType.LT) {
                    cx = number;
                    ax = number - 500.0;
                    bx = (ax + cx) / 2.0;
                } else if (constraint.getType() == ParamConstraintType.EQ) {
                    bx = number;
                    ax = number - 500.0;
                    cx = number + 500.0;
                } else {
                    ax = -500.0;
                    bx = 0.0;
                    cx = 500.0;
                }
                double[] mean = new double[1];
                // dmean is the density at the mean
                double dmean = -brent(param, ax, bx, cx, this.tolerance, mean, parameters);
                double gap = 0.005;
                double denom;
                do {
                    gap = 2.0 * gap;
                    int gapThreshold = 1;
                    double minDenom = 0.01;
                    if (gap > gapThreshold) {
                        denom = minDenom;
                        break;
                    }
                    System.out.println(p.getNodeA() + " " + p.getNodeA().getNodeType());
                    System.out.println(p.getNodeB() + " " + p.getNodeB().getNodeType());
                    double dmeanplus = neglogpost(param, mean[0] + gap, parameters);
                    denom = dmean + dmeanplus;
                    if (denom < minDenom)
                        denom = minDenom;
                // System.out.println("gap = "+gap+"; denom = "+denom+"; dmean = "+dmean+"; dmeanplus = "+dmeanplus);
                } while (denom < 0.0);
                double vr = (this.stretch1 * 0.5 * gap * gap) / denom;
                // System.out.println("vr = "+vr+" param = "+param);
                // FORMAPPROXDIST end
                // DRAWFROMAPPROX begin
                boolean realdraw = false;
                double rj = 0.0, accept = 0.0, cand = 0.0;
                while (!realdraw || rj <= accept) {
                    cand = mean[0] + Math.max(RandomUtil.getInstance().nextNormal(0, 1) * Math.sqrt(vr), 0);
                    realdraw = (constraint.wouldBeSatisfied(cand));
                    if (realdraw) {
                        // System.out.println("dcand start");
                        double dcand = -1.0 * neglogpost(param, cand, parameters);
                        // System.out.println("dcand end");
                        double numer = dcand - dmean;
                        double denom1 = (-1.0 * Math.sqrt(cand - mean[0]) / (2.0 * vr)) - Math.log(this.stretch2);
                        rj = numer - denom1;
                        accept = Math.log(RandomUtil.getInstance().nextDouble());
                        int rejectionThreshold = 5;
                        if (rj > rejectionThreshold) {
                            // System.out.println("rj = "+rj);
                            rj = rejectionThreshold;
                        }
                    }
                }
                // DRAWFROMAPPROX end
                // System.out.println("end of iteration");
                // UPDATEPARM
                Parameter ppost = (Parameter) postFreeParams.get(param);
                if (ppost.isFixed())
                    posteriorIm.setFixedParamValue(ppost, cand);
                else
                    posteriorIm.setParamValue(ppost, cand);
            // UPDATEPARM end
            }
        }
        int subsampleStride = 50;
        if (iter % subsampleStride == 0 && iter > 0) {
            for (int i = 0; i < numParameters; i++) {
                Parameter ppost = (posteriorIm.getSemPm()).getParameters().get(i);
                data.set(i, iter / subsampleStride - 1, posteriorIm.getParamValue(ppost));
            }
        }
    }
    dataSet = data;
    this.estimatedSem = posteriorIm;
// setMeans(posteriorIm, data);
}

Example 20 with TetradMatrix

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

the class SemOptimizerEm method updateMatrices.

private void updateMatrices() {
    TetradMatrix impliedCovar = semIm.getImplCovar(true);
    for (int i = 0; i < numObserved; i++) {
        for (int j = i; j < numObserved; j++) {
            yCovModel.set(i, j, impliedCovar.get(idxObserved[i], idxObserved[j]));
            yCovModel.set(j, i, impliedCovar.get(idxObserved[i], idxObserved[j]));
        }
        for (int j = 0; j < numLatent; j++) {
            yzCovModel.set(i, j, impliedCovar.get(idxObserved[i], idxLatent[j]));
        }
    }
    for (int i = 0; i < numLatent; i++) {
        for (int j = i; j < numLatent; j++) {
            zCovModel.set(i, j, impliedCovar.get(idxLatent[i], idxLatent[j]));
            zCovModel.set(j, i, impliedCovar.get(idxLatent[i], idxLatent[j]));
        }
    }
}