Examples with Parameter edu.cmu.tetrad.sem.Parameter used on opensource projects

Example 6 with Parameter

use of edu.cmu.tetrad.sem.Parameter in project tetrad by cmu-phil.

the class PurifyScoreBased method gaussianExpectation.

/**
 * The expectation step for the structural EM algorithm. This is heavily based on "EM Algorithms for ML Factor
 * Analysis", by Rubin and Thayer (Psychometrika, 1982)
 */
private void gaussianExpectation(SemIm semIm) {
    // Get the parameters
    double[][] beta = // latent-to-latent coefficients
    new double[numLatent][numLatent];
    double[][] fi = // latent error terms covariance
    new double[numLatent][numLatent];
    double[][] lambdaI = // observed-to-indicatorcoefficients
    new double[numObserved][numObserved];
    double[][] lambdaL = // latent-to-indicatorcoefficients
    new double[numObserved][numLatent];
    double[][] tau = // measurement error variance
    new double[numObserved][numObserved];
    // structural EM algorithm such as in MimBuildScoreSearch.
    for (int i = 0; i < numLatent; i++) {
        for (int j = 0; j < numLatent; j++) {
            beta[i][j] = 0.;
            fi[i][j] = 0.;
        }
    }
    for (int i = 0; i < numObserved; i++) {
        for (int j = 0; j < numLatent; j++) {
            lambdaL[i][j] = 0.;
        }
    }
    for (int i = 0; i < numObserved; i++) {
        for (int j = 0; j < numObserved; j++) {
            tau[i][j] = 0.;
            lambdaI[i][j] = 0.;
        }
    }
    List parameters = semIm.getFreeParameters();
    double[] paramValues = semIm.getFreeParamValues();
    for (int i = 0; i < parameters.size(); i++) {
        Parameter parameter = (Parameter) parameters.get(i);
        if (parameter.getType() == ParamType.COEF) {
            Node from = parameter.getNodeA();
            Node to = parameter.getNodeB();
            if (to.getNodeType() == NodeType.MEASURED && from.getNodeType() == NodeType.LATENT) {
                // latent-to-indicator edge
                int position1 = (Integer) latentNames.get(from.getName());
                int position2 = (Integer) observableNames.get(to.getName());
                lambdaL[position2][position1] = paramValues[i];
            } else if (to.getNodeType() == NodeType.MEASURED && from.getNodeType() == NodeType.MEASURED) {
                // indicator-to-indicator edge
                int position1 = (Integer) observableNames.get(from.getName());
                int position2 = (Integer) observableNames.get(to.getName());
                lambdaI[position2][position1] = paramValues[i];
            } else if (to.getNodeType() == NodeType.LATENT) {
                // latent-to-latent edge
                int position1 = (Integer) latentNames.get(from.getName());
                int position2 = (Integer) latentNames.get(to.getName());
                beta[position2][position1] = paramValues[i];
            }
        } else if (parameter.getType() == ParamType.VAR) {
            Node exo = parameter.getNodeA();
            if (exo.getNodeType() == NodeType.ERROR) {
                Iterator ci = semIm.getSemPm().getGraph().getChildren(exo).iterator();
                exo = // Assuming error nodes have only one children in SemGraphs...
                (Node) ci.next();
            }
            if (exo.getNodeType() == NodeType.LATENT) {
                fi[((Integer) latentNames.get(exo.getName()))][((Integer) latentNames.get(exo.getName()))] = paramValues[i];
            } else {
                tau[((Integer) observableNames.get(exo.getName()))][((Integer) observableNames.get(exo.getName()))] = paramValues[i];
            }
        } else if (parameter.getType() == ParamType.COVAR) {
            Node exo1 = parameter.getNodeA();
            Node exo2 = parameter.getNodeB();
            // exo1.getNodeType and exo1.getNodeType *should* be error terms of measured variables
            // We will change the pointers to point to their respective indicators
            // Iterator ci = semIm.getEstIm().getGraph().getChildren(exo1)
            // .iterator();
            // exo1 =
            // (Node) ci.next(); //Assuming error nodes have only one children in SemGraphs...
            // ci = semIm.getEstIm().getGraph().getChildren(exo2).iterator();
            // exo2 =
            // (Node) ci.next(); //Assuming error nodes have only one children in SemGraphs...
            exo1 = semIm.getSemPm().getGraph().getVarNode(exo1);
            exo2 = semIm.getSemPm().getGraph().getVarNode(exo2);
            tau[((Integer) observableNames.get(exo1.getName()))][((Integer) observableNames.get(exo2.getName()))] = tau[((Integer) observableNames.get(exo2.getName()))][((Integer) observableNames.get(exo1.getName()))] = paramValues[i];
        }
    }
    // Fill expected sufficiente statistics accordingly to the order of
    // the variables table
    double[][] identity = new double[numLatent][numLatent];
    for (int i = 0; i < numLatent; i++) {
        for (int j = 0; j < numLatent; j++) {
            if (i == j) {
                identity[i][j] = 1.;
            } else {
                identity[i][j] = 0.;
            }
        }
    }
    double[][] identityI = new double[numObserved][numObserved];
    for (int i = 0; i < numObserved; i++) {
        for (int j = 0; j < numObserved; j++) {
            if (i == j) {
                identityI[i][j] = 1.;
            } else {
                identityI[i][j] = 0.;
            }
        }
    }
    double[][] iMinusB = MatrixUtils.inverse(MatrixUtils.subtract(identity, beta));
    double[][] latentImpliedCovar = MatrixUtils.product(iMinusB, MatrixUtils.product(fi, MatrixUtils.transpose(iMinusB)));
    double[][] iMinusI = MatrixUtils.inverse(MatrixUtils.subtract(identityI, lambdaI));
    double[][] indImpliedCovar = MatrixUtils.product(MatrixUtils.product(iMinusI, MatrixUtils.sum(MatrixUtils.product(MatrixUtils.product(lambdaL, latentImpliedCovar), MatrixUtils.transpose(lambdaL)), tau)), MatrixUtils.transpose(iMinusI));
    double[][] loadingLatentCovar = MatrixUtils.product(iMinusI, MatrixUtils.product(lambdaL, latentImpliedCovar));
    double[][] smallDelta = MatrixUtils.product(MatrixUtils.inverse(indImpliedCovar), loadingLatentCovar);
    double[][] bigDelta = MatrixUtils.subtract(latentImpliedCovar, MatrixUtils.product(MatrixUtils.transpose(loadingLatentCovar), smallDelta));
    this.Cyz = MatrixUtils.product(this.Cyy, smallDelta);
    this.Czz = MatrixUtils.sum(MatrixUtils.product(MatrixUtils.transpose(smallDelta), Cyz), bigDelta);
}

Example 7 with Parameter

use of edu.cmu.tetrad.sem.Parameter in project tetrad by cmu-phil.

the class PurifyScoreBased method gaussianMaximization.

/*private SemIm getDummyExample()
    {
        this.numObserved = 13;
        this.numLatent = 4;

        ProtoSemGraph newGraph = new ProtoSemGraph();
        Node v[] = new Node[17];
        for (int i = 0; i < 17; i++) {
            if (i < 4)
                v[i] = new GraphNode("L" + (i + 1));
            else
                v[i] = new GraphNode("v" + (i - 3));
            newGraph.addNode(v[i]);
        }
        for (int l = 0; l < numLatent; l++) {
            for (int l2 = l + 1; l2 < numLatent; l2++)
                newGraph.addDirectedEdge(v[l], v[l2]);
            for (int i = 0; i < 3; i++)
                newGraph.addDirectedEdge(v[l], v[l * 3 + i + 4]);
        }
        newGraph.addDirectedEdge(v[3], v[16]);
        newGraph.addDirectedEdge(v[4], v[6]);
        newGraph.addDirectedEdge(v[14], v[15]);
        newGraph.addBidirectedEdge(v[9], v[10]);
        newGraph.addBidirectedEdge(v[7], v[12]);
        newGraph.addBidirectedEdge(v[12], v[16]);

        SemIm realIm = SemIm.newInstance(new SemPm(newGraph));

        DataSet data = new DataSet(realIm.simulateData(1000));
        System.out.println(data.toString());
        System.out.println(new CovarianceMatrix(data));
        System.out.println();

        this.latentNames = new Hashtable();
        this.latentNodes = new ArrayList();
        for (int i = 0; i < numLatent; i++) {
            latentNames.put(v[i].getNode(), new Integer(i));
            latentNodes.add(v[i]);
        }
        this.observableNames = new Hashtable();
        this.measuredNodes = new ArrayList();
        for (int i = numLatent; i < numLatent + numObserved; i++) {
            observableNames.put(v[i].getNode(), new Integer(i - numLatent));
            measuredNodes.add(v[i]);
        }

        double temp[][] = (new CovarianceMatrix(data)).getMatrix();
        Cyy = new double[numObserved][numObserved];
        Czz = new double[numLatent][numLatent];
        Cyz = new double[numObserved][numLatent];
        for (int i = 0; i < numLatent; i++)
            for (int j = 0; j < numLatent; j++)
                Czz[i][j] = temp[i][j];
        for (int i = 0; i < numObserved; i++) {
            for (int j = 0; j < numLatent; j++)
                 Cyz[i][j] = temp[i + numLatent][j];
             for (int j = 0; j < numObserved; j++)
                 Cyy[i][j] = temp[i + numLatent][j + numLatent];
        }


        this.correlatedErrors = new boolean[numObserved][numObserved];
        this.latentParent = new boolean[numObserved][numLatent];
        this.observedParent = new boolean[numObserved][numObserved];

        CovarianceMatrix = new CovarianceMatrix(data);

        //Information for MAG maximization
        this.parents = new int[this.numObserved][];
        this.spouses = new int[this.numObserved][];
        this.nSpouses = new int[this.numObserved];
        this.parentsLat = new int[this.numLatent][];
        this.parentsL = new boolean[this.numObserved][];
        this.parentsCov = new double[this.numObserved][][];
        this.parentsChildCov = new double[this.numObserved][];
        this.parentsLatCov = new double[this.numLatent][][];
        this.parentsChildLatCov = new double[this.numLatent][];
        this.pseudoParentsCov = new double[this.numObserved][][];
        this.pseudoParentsChildCov = new double[this.numObserved][];
        this.covErrors = new double[this.numObserved][this.numObserved];
        this.oldCovErrors = new double[this.numObserved][this.numObserved];
        this.sampleCovErrors = new double[this.numObserved][this.numObserved];
        this.varErrorLatent = new double[this.numLatent];
        this.omega = new double[this.numLatent + this.numObserved - 1][this.numLatent + this.numObserved - 1];
        this.omegaI = new double[this.numLatent + this.numObserved - 1];
        this.selectedInverseOmega = new double[this.numObserved][][];
        this.auxInverseOmega = new double[this.numObserved][][];
        this.parentsResidualsCovar = new double[this.numObserved][][];
        this.iResidualsCovar = new double[this.numObserved + this.numLatent - 1];
        this.betas = new double[this.numObserved][this.numObserved + this.numLatent];
        this.oldBetas = new double[this.numObserved][this.numObserved + this.numLatent];
        this.betasLat = new double[this.numLatent][this.numLatent];

        for (int i = 0; i < numLatent; i++)
             v[i].setNodeType(NodeType.LATENT);
        SemGraph2 semGraph = new SemGraph2(newGraph);
        initializeGaussianEM(semGraph);

        return realIm;
    }*/
private double gaussianMaximization(SemIm semIm) {
    semIm.getSemPm().getGraph().setShowErrorTerms(true);
    // Fill matrices with semIm parameters
    for (int i = 0; i < this.numObserved; i++) {
        for (int j = 0; j < this.numObserved + this.numLatent; j++) {
            this.betas[i][j] = 0.;
        }
    }
    for (int i = 0; i < this.numLatent; i++) {
        for (int j = 0; j < this.numLatent; j++) {
            this.betasLat[i][j] = 0.;
        }
    }
    for (int i = 0; i < this.numObserved; i++) {
        for (int j = 0; j < this.numObserved; j++) {
            this.covErrors[i][j] = 0.;
        }
    }
    for (Iterator it = semIm.getFreeParameters().iterator(); it.hasNext(); ) {
        Parameter nextP = (Parameter) it.next();
        if (nextP.getType() == ParamType.COEF) {
            Node node1 = nextP.getNodeA();
            Node node2 = nextP.getNodeB();
            if (node1.getNodeType() == NodeType.LATENT && node2.getNodeType() == NodeType.LATENT) {
                continue;
            }
            Node latent = null, observed = null;
            if (node1.getNodeType() == NodeType.LATENT) {
                latent = node1;
                observed = node2;
            } else if (node2.getNodeType() == NodeType.LATENT) {
                latent = node2;
                observed = node1;
            }
            if (latent != null) {
                int index1 = (Integer) this.latentNames.get(latent.getName());
                int index2 = (Integer) this.observableNames.get(observed.getName());
                this.betas[index2][index1] = semIm.getParamValue(nextP);
            } else {
                int index1 = (Integer) this.observableNames.get(node1.getName());
                int index2 = (Integer) this.observableNames.get(node2.getName());
                if (semIm.getSemPm().getGraph().isParentOf(node1, node2)) {
                    this.betas[index2][this.numLatent + index1] = semIm.getParamValue(nextP);
                } else {
                    this.betas[index1][this.numLatent + index2] = semIm.getParamValue(nextP);
                }
            }
        } else if (nextP.getType() == ParamType.COVAR) {
            Node exo1 = nextP.getNodeA();
            Node exo2 = nextP.getNodeB();
            // exo1.getNodeType and exo1.getNodeType *should* be error terms of measured variables
            // We will change the pointers to point to their respective indicators
            // Iterator ci = semIm.getEstIm().getGraph().getChildren(exo1)
            // .iterator();
            // exo1 =
            // (Node) ci.next(); //Assuming error nodes have only one children in SemGraphs...
            // ci = semIm.getEstIm().getGraph().getChildren(exo2).iterator();
            // exo2 =
            // (Node) ci.next(); //Assuming error nodes have only one children in SemGraphs...
            exo1 = semIm.getSemPm().getGraph().getVarNode(exo1);
            exo2 = semIm.getSemPm().getGraph().getVarNode(exo2);
            int index1 = (Integer) this.observableNames.get(exo1.getName());
            int index2 = (Integer) this.observableNames.get(exo2.getName());
            this.covErrors[index1][index2] = this.covErrors[index2][index1] = semIm.getParamValue(nextP);
        } else if (nextP.getType() == ParamType.VAR) {
            Node exo = nextP.getNodeA();
            if (exo.getNodeType() == NodeType.LATENT) {
                continue;
            }
            // Iterator ci = semIm.getEstIm().getGraph().getChildren(exo)
            // .iterator();
            // exo =
            // (Node) ci.next(); //Assuming error nodes have only one children in SemGraphs...
            exo = semIm.getSemPm().getGraph().getVarNode(exo);
            if (exo.getNodeType() == NodeType.MEASURED) {
                int index = (Integer) this.observableNames.get(exo.getName());
                this.covErrors[index][index] = semIm.getParamValue(nextP);
            }
        }
    }
    // Find estimates for the latent->latent edges and latent variances
    // Assuming latents[0] is always the exogenous node in the latent layer
    this.varErrorLatent[0] = this.Czz[0][0];
    for (int i = 1; i < this.numLatent; i++) {
        for (int j = 0; j < this.parentsLat[i].length; j++) {
            this.parentsChildLatCov[i][j] = this.Czz[i][this.parentsLat[i][j]];
            for (int k = j; k < this.parentsLat[i].length; k++) {
                this.parentsLatCov[i][j][k] = this.Czz[this.parentsLat[i][j]][this.parentsLat[i][k]];
                this.parentsLatCov[i][k][j] = this.parentsLatCov[i][j][k];
            }
        }
        double[] betaL = MatrixUtils.product(MatrixUtils.inverse(this.parentsLatCov[i]), this.parentsChildLatCov[i]);
        this.varErrorLatent[i] = this.Czz[i][i] - MatrixUtils.innerProduct(this.parentsChildLatCov[i], betaL);
        for (int j = 0; j < this.parentsLat[i].length; j++) {
            this.betasLat[i][this.parentsLat[i][j]] = betaL[j];
        }
    }
    // Initialize the covariance matrix for the parents of every observed node
    for (int i = 0; i < this.numObserved; i++) {
        for (int j = 0; j < this.parents[i].length; j++) {
            if (this.parentsL[i][j]) {
                this.parentsChildCov[i][j] = this.Cyz[i][this.parents[i][j]];
            } else {
                this.parentsChildCov[i][j] = this.Cyy[i][this.parents[i][j]];
            }
            for (int k = j; k < this.parents[i].length; k++) {
                if (this.parentsL[i][j] && this.parentsL[i][k]) {
                    this.parentsCov[i][j][k] = this.Czz[this.parents[i][j]][this.parents[i][k]];
                } else if (!this.parentsL[i][j] && this.parentsL[i][k]) {
                    this.parentsCov[i][j][k] = this.Cyz[this.parents[i][j]][this.parents[i][k]];
                } else if (this.parentsL[i][j] && !this.parentsL[i][k]) {
                    this.parentsCov[i][j][k] = this.Cyz[this.parents[i][k]][this.parents[i][j]];
                } else {
                    this.parentsCov[i][j][k] = this.Cyy[this.parents[i][j]][this.parents[i][k]];
                }
                this.parentsCov[i][k][j] = this.parentsCov[i][j][k];
            }
        }
    }
    // ICF algorithm of Drton and Richardson to find estimates for the other edges and variances/covariances
    double change;
    int iter = 0;
    do {
        for (int i = 0; i < this.covErrors.length; i++) {
            for (int j = 0; j < this.covErrors.length; j++) {
                this.oldCovErrors[i][j] = this.covErrors[i][j];
            }
        }
        for (int i = 0; i < this.numObserved; i++) {
            for (int j = 0; j < this.betas[i].length; j++) {
                this.oldBetas[i][j] = this.betas[i][j];
            }
        }
        for (int i = 0; i < this.numObserved; i++) {
            // Build matrix Omega_{-i,-i} as defined in Drton and Richardson (2003)
            for (int ii = 0; ii < this.omega.length; ii++) {
                for (int j = 0; j < this.omega.length; j++) {
                    this.omega[ii][j] = 0.;
                }
            }
            for (int ii = 0; ii < this.numLatent; ii++) {
                this.omegaI[ii] = 0.;
                this.omega[ii][ii] = this.varErrorLatent[ii];
            }
            for (int ii = 0; ii < this.numObserved; ii++) {
                if (ii > i) {
                    this.omegaI[this.numLatent + ii - 1] = this.covErrors[i][ii];
                    this.omega[this.numLatent + ii - 1][this.numLatent + ii - 1] = this.covErrors[ii][ii];
                } else if (ii < i) {
                    this.omegaI[this.numLatent + ii] = this.covErrors[i][ii];
                    this.omega[this.numLatent + ii][this.numLatent + ii] = this.covErrors[ii][ii];
                }
            }
            for (int ii = 0; ii < this.numObserved; ii++) {
                int index_ii;
                if (ii > i) {
                    index_ii = this.numLatent + ii - 1;
                } else if (ii < i) {
                    index_ii = this.numLatent + ii;
                } else {
                    continue;
                }
                for (int j = 0; j < this.nSpouses[ii]; j++) {
                    if (this.spouses[ii][j] > i) {
                        this.omega[index_ii][this.numLatent + this.spouses[ii][j] - 1] = this.covErrors[ii][this.spouses[ii][j]];
                    } else if (this.spouses[ii][j] < i) {
                        this.omega[index_ii][this.numLatent + this.spouses[ii][j]] = this.covErrors[ii][this.spouses[ii][j]];
                    }
                }
            }
            // Find new residuals covariance matrix for every ii != i
            for (int ii = 0; ii < this.numObserved; ii++) {
                if (ii == i) {
                    continue;
                }
                for (int j = ii; j < this.numObserved; j++) {
                    if (j == i) {
                        continue;
                    }
                    this.sampleCovErrors[ii][j] = this.Cyy[ii][j];
                    for (int p = 0; p < this.parents[ii].length; p++) {
                        if (this.parentsL[ii][p]) {
                            this.sampleCovErrors[ii][j] -= this.betas[ii][this.parents[ii][p]] * this.Cyz[j][this.parents[ii][p]];
                        } else {
                            this.sampleCovErrors[ii][j] -= this.betas[ii][this.numLatent + this.parents[ii][p]] * this.Cyy[j][this.parents[ii][p]];
                        }
                    }
                    for (int p = 0; p < this.parents[j].length; p++) {
                        if (this.parentsL[j][p]) {
                            this.sampleCovErrors[ii][j] -= this.betas[j][this.parents[j][p]] * this.Cyz[ii][this.parents[j][p]];
                        } else {
                            this.sampleCovErrors[ii][j] -= this.betas[j][this.numLatent + this.parents[j][p]] * this.Cyy[ii][this.parents[j][p]];
                        }
                    }
                    for (int p1 = 0; p1 < this.parents[ii].length; p1++) {
                        for (int p2 = 0; p2 < this.parents[j].length; p2++) {
                            if (this.parentsL[ii][p1] && this.parentsL[j][p2]) {
                                this.sampleCovErrors[ii][j] += this.betas[ii][this.parents[ii][p1]] * this.betas[j][this.parents[j][p2]] * this.Czz[this.parents[ii][p1]][this.parents[j][p2]];
                            } else if (this.parentsL[ii][p1] && !this.parentsL[j][p2]) {
                                this.sampleCovErrors[ii][j] += this.betas[ii][this.parents[ii][p1]] * this.betas[j][this.numLatent + this.parents[j][p2]] * this.Cyz[this.parents[j][p2]][this.parents[ii][p1]];
                            } else if (!this.parentsL[ii][p1] && this.parentsL[j][p2]) {
                                this.sampleCovErrors[ii][j] += this.betas[ii][this.numLatent + this.parents[ii][p1]] * this.betas[j][this.parents[j][p2]] * this.Cyz[this.parents[ii][p1]][this.parents[j][p2]];
                            } else {
                                this.sampleCovErrors[ii][j] += this.betas[ii][this.numLatent + this.parents[ii][p1]] * this.betas[j][this.numLatent + this.parents[j][p2]] * this.Cyy[this.parents[ii][p1]][this.parents[j][p2]];
                            }
                        }
                    }
                    this.sampleCovErrors[j][ii] = this.sampleCovErrors[ii][j];
                }
            }
            // First, find the covariance of the parents of i and the residuals \epsilon_{-i}
            for (int ii = 0; ii < this.parents[i].length; ii++) {
                // covariance of the parent wrt every residual of latents
                if (this.parentsL[i][ii]) {
                    this.parentsResidualsCovar[i][ii][0] = this.Czz[this.parents[i][ii]][0];
                } else {
                    this.parentsResidualsCovar[i][ii][0] = this.Cyz[this.parents[i][ii]][0];
                }
                for (int j = 1; j < this.numLatent; j++) {
                    if (this.parentsL[i][ii]) {
                        this.parentsResidualsCovar[i][ii][j] = this.Czz[this.parents[i][ii]][j];
                        for (int p = 0; p < this.parentsLat[j].length; p++) {
                            this.parentsResidualsCovar[i][ii][j] -= this.betasLat[j][this.parentsLat[j][p]] * this.Czz[this.parents[i][ii]][this.parentsLat[j][p]];
                        }
                    } else {
                        this.parentsResidualsCovar[i][ii][j] = this.Cyz[this.parents[i][ii]][j];
                        for (int p = 0; p < this.parentsLat[j].length; p++) {
                            this.parentsResidualsCovar[i][ii][j] -= this.betasLat[j][this.parentsLat[j][p]] * this.Cyz[this.parents[i][ii]][this.parentsLat[j][p]];
                        }
                    }
                }
                // covariance of the parent wrt every residual of observables (except for i)
                for (int j = 0; j < this.numObserved; j++) {
                    int index_j;
                    if (j < i) {
                        index_j = this.numLatent + j;
                    } else if (j > i) {
                        index_j = this.numLatent + j - 1;
                    } else {
                        continue;
                    }
                    if (this.parentsL[i][ii]) {
                        this.parentsResidualsCovar[i][ii][index_j] = this.Cyz[j][this.parents[i][ii]];
                        for (int p = 0; p < this.parents[j].length; p++) {
                            if (this.parentsL[j][p]) {
                                this.parentsResidualsCovar[i][ii][index_j] -= this.betas[j][this.parents[j][p]] * this.Czz[this.parents[i][ii]][this.parents[j][p]];
                            } else {
                                this.parentsResidualsCovar[i][ii][index_j] -= this.betas[j][this.numLatent + this.parents[j][p]] * this.Cyz[this.parents[j][p]][this.parents[i][ii]];
                            }
                        }
                    } else {
                        this.parentsResidualsCovar[i][ii][index_j] = this.Cyy[j][this.parents[i][ii]];
                        for (int p = 0; p < this.parents[j].length; p++) {
                            if (this.parentsL[j][p]) {
                                this.parentsResidualsCovar[i][ii][index_j] -= this.betas[j][this.parents[j][p]] * this.Cyz[this.parents[i][ii]][this.parents[j][p]];
                            } else {
                                this.parentsResidualsCovar[i][ii][index_j] -= this.betas[j][this.numLatent + this.parents[j][p]] * this.Cyy[this.parents[j][p]][this.parents[i][ii]];
                            }
                        }
                    }
                }
            }
            // Now, find the covariance of Y_i with respect to everybody else's residuals
            this.iResidualsCovar[0] = // the first latent is exogenous
            this.Cyz[i][0];
            for (int j = 1; j < this.numLatent; j++) {
                this.iResidualsCovar[j] = this.Cyz[i][j];
                for (int p = 0; p < this.parentsLat[j].length; p++) {
                    this.iResidualsCovar[j] -= this.betasLat[j][this.parentsLat[j][p]] * this.Cyz[i][this.parentsLat[j][p]];
                }
            }
            for (int j = 0; j < this.numObserved; j++) {
                int index_j;
                if (j < i) {
                    index_j = this.numLatent + j;
                } else if (j > i) {
                    index_j = this.numLatent + j - 1;
                } else {
                    continue;
                }
                this.iResidualsCovar[index_j] = this.Cyy[i][j];
                for (int p = 0; p < this.parents[j].length; p++) {
                    if (this.parentsL[j][p]) {
                        this.iResidualsCovar[index_j] -= this.betas[j][this.parents[j][p]] * this.Cyz[i][this.parents[j][p]];
                    } else {
                        this.iResidualsCovar[index_j] -= this.betas[j][this.numLatent + this.parents[j][p]] * this.Cyy[i][this.parents[j][p]];
                    }
                }
            }
            // Transform it to get the covariance of parents of i and pseudo-variables Z_sp(i)
            double[][] inverseOmega = MatrixUtils.inverse(this.omega);
            for (int ii = 0; ii < this.nSpouses[i]; ii++) {
                int sp_index;
                if (this.spouses[i][ii] > i) {
                    sp_index = this.numLatent + this.spouses[i][ii] - 1;
                } else {
                    sp_index = this.numLatent + this.spouses[i][ii];
                }
                for (int j = 0; j < this.numLatent + this.numObserved - 1; j++) {
                    this.selectedInverseOmega[i][ii][j] = inverseOmega[sp_index][j];
                }
            }
            for (int ii = 0; ii < this.nSpouses[i]; ii++) {
                for (int j = 0; j < this.numLatent; j++) {
                    this.auxInverseOmega[i][ii][j] = this.selectedInverseOmega[i][ii][j] * this.varErrorLatent[j];
                }
                for (int j = 0; j < this.numObserved; j++) {
                    int index_j;
                    if (j > i) {
                        index_j = this.numLatent + j - 1;
                    } else if (j < i) {
                        index_j = this.numLatent + j;
                    } else {
                        continue;
                    }
                    this.auxInverseOmega[i][ii][index_j] = 0;
                    for (int k = 0; k < this.numObserved; k++) {
                        int index_k;
                        if (k > i) {
                            index_k = this.numLatent + k - 1;
                        } else if (k < i) {
                            index_k = this.numLatent + k;
                        } else {
                            continue;
                        }
                        this.auxInverseOmega[i][ii][index_j] += this.selectedInverseOmega[i][ii][index_k] * this.sampleCovErrors[k][j];
                    }
                }
            }
            for (int ii = 0; ii < this.parents[i].length; ii++) {
                for (int j = ii; j < this.parents[i].length; j++) {
                    this.pseudoParentsCov[i][ii][j] = this.pseudoParentsCov[i][j][ii] = this.parentsCov[i][ii][j];
                }
            }
            for (int ii = 0; ii < this.parents[i].length; ii++) {
                for (int j = 0; j < this.nSpouses[i]; j++) {
                    this.pseudoParentsCov[i][ii][this.parents[i].length + j] = 0.;
                    for (int k = 0; k < this.numLatent + this.numObserved - 1; k++) {
                        this.pseudoParentsCov[i][ii][this.parents[i].length + j] += this.parentsResidualsCovar[i][ii][k] * this.selectedInverseOmega[i][j][k];
                    }
                    this.pseudoParentsCov[i][this.parents[i].length + j][ii] = this.pseudoParentsCov[i][ii][this.parents[i].length + j];
                }
            }
            for (int ii = 0; ii < this.nSpouses[i]; ii++) {
                for (int j = ii; j < this.nSpouses[i]; j++) {
                    this.pseudoParentsCov[i][this.parents[i].length + ii][this.parents[i].length + j] = 0;
                    for (int k = 0; k < this.numLatent + this.numObserved - 1; k++) {
                        this.pseudoParentsCov[i][this.parents[i].length + ii][this.parents[i].length + j] += this.auxInverseOmega[i][ii][k] * this.selectedInverseOmega[i][j][k];
                    }
                    this.pseudoParentsCov[i][this.parents[i].length + j][this.parents[i].length + ii] = this.pseudoParentsCov[i][this.parents[i].length + ii][this.parents[i].length + j];
                    if (this.pseudoParentsCov[i][this.parents[i].length + j][this.parents[i].length + ii] == 0.) {
                        System.out.println("Zero here... Iter = " + iter);
                        /*for (int k = 0; k < this.numLatent + this.numObserved - 1; k++)
                                System.out.println(this.auxInverseOmega[i][ii][k] + " " +  this.selectedInverseOmega[i][j][k]);
                            System.out.println();
                            for (int v = 0; v < this.numLatent + this.numObserved - 1; v++) {
                               for (int k = 0; k < this.numLatent + this.numObserved - 1; k++)
                                   System.out.print(this.omega[v][k] + " ");
                                System.out.println();
                            }
                            System.out.println(semIm.getEstIm().getDag());
                            System.out.println("PARENTS");
                            for (int n = 0; n < this.numObserved; n++) {
                                System.out.print(this.measuredNodes.get(n).toString() + " - ");
                                for (int p = 0; p < this.parents[n].length; p++) {
                                    if (parentsL[n][p])
                                       System.out.print(this.latentNodes.get(parents[n][p]).toString() + " ");
                                    else
                                       System.out.print(this.measuredNodes.get(parents[n][p]).toString() + " ");
                                }
                                System.out.println();
                            }

                            System.out.println("SPOUSES");
                            for (int n = 0; n < this.numObserved; n++) {
                                System.out.print(this.measuredNodes.get(n).toString() + " - ");
                                for (int p = 0; p < this.nSpouses[n]; p++) {
                                    System.out.print(this.measuredNodes.get(spouses[n][p]).toString() + " ");
                                }
                                System.out.println();
                            }
                            System.exit(0);*/
                        iter = 1000;
                        break;
                    }
                }
            }
            // Get the covariance of parents of i and pseudo-variables Z_sp(i) with respect to i
            for (int ii = 0; ii < this.parents[i].length; ii++) {
                this.pseudoParentsChildCov[i][ii] = this.parentsChildCov[i][ii];
            }
            for (int j = 0; j < this.nSpouses[i]; j++) {
                this.pseudoParentsChildCov[i][this.parents[i].length + j] = 0;
                for (int k = 0; k < this.numLatent + this.numObserved - 1; k++) {
                    this.pseudoParentsChildCov[i][this.parents[i].length + j] += selectedInverseOmega[i][j][k] * this.iResidualsCovar[k];
                }
            }
            // Finally, regress Y_i on {parents} union {Z_i}
            // thisI = i;
            double[] params = MatrixUtils.product(MatrixUtils.inverse(this.pseudoParentsCov[i]), this.pseudoParentsChildCov[i]);
            // Update betas and omegas (entries in covErrors)
            for (int j = 0; j < this.parents[i].length; j++) {
                if (this.parentsL[i][j]) {
                    this.betas[i][this.parents[i][j]] = params[j];
                } else {
                    this.betas[i][this.numLatent + this.parents[i][j]] = params[j];
                }
            }
            for (int j = 0; j < this.nSpouses[i]; j++) {
                this.covErrors[i][this.spouses[i][j]] = this.covErrors[this.spouses[i][j]][i] = params[this.parents[i].length + j];
                if (this.spouses[i][j] > i) {
                    this.omegaI[this.numLatent + this.spouses[i][j] - 1] = params[this.parents[i].length + j];
                } else {
                    this.omegaI[this.numLatent + this.spouses[i][j]] = params[this.parents[i].length + j];
                }
            }
            double conditionalVar = this.Cyy[i][i] - MatrixUtils.innerProduct(this.pseudoParentsChildCov[i], params);
            this.covErrors[i][i] = conditionalVar + MatrixUtils.innerProduct(MatrixUtils.product(this.omegaI, inverseOmega), this.omegaI);
        }
        change = 0.;
        for (int i = 0; i < this.covErrors.length; i++) {
            for (int j = i; j < this.covErrors.length; j++) {
                change += Math.abs(this.oldCovErrors[i][j] - this.covErrors[i][j]);
            }
        }
        for (int i = 0; i < this.numObserved; i++) {
            for (int j = 0; j < this.betas[i].length; j++) {
                change += Math.abs(oldBetas[i][j] - this.betas[i][j]);
            }
        }
        iter++;
    // System.out.println("Iteration = " + iter + ", change = " + change);
    } while (iter < 200 && change > 0.01);
    // Now, copy updated parameters back to semIm
    try {
        for (int i = 0; i < this.numObserved; i++) {
            Node node = semIm.getSemPm().getGraph().getNode(this.measuredNodes.get(i).toString());
            // Node nodeErrorTerm = null;
            // for (Node parent : semIm.getEstIm().getGraph().getParents(node)) {
            // if (parent.getNodeType() == NodeType.ERROR) {
            // semIm.setParamValue(parent, parent,
            // this.covErrors[i][i]);
            // nodeErrorTerm = parent;
            // break;
            // }
            // }
            Node nodeErrorTerm = semIm.getSemPm().getGraph().getExogenous(node);
            for (int j = 0; j < this.parents[i].length; j++) {
                Node parent;
                if (this.parentsL[i][j]) {
                    parent = semIm.getSemPm().getGraph().getNode(this.latentNodes.get(this.parents[i][j]).toString());
                } else {
                    parent = semIm.getSemPm().getGraph().getNode(this.measuredNodes.get(this.parents[i][j]).toString());
                }
                if (this.parentsL[i][j]) {
                    semIm.setParamValue(parent, node, this.betas[i][this.parents[i][j]]);
                } else {
                    semIm.setParamValue(parent, node, this.betas[i][this.numLatent + this.parents[i][j]]);
                }
            }
            for (int j = 0; j < this.nSpouses[i]; j++) {
                if (this.spouses[i][j] > i) {
                    Node spouse = semIm.getSemPm().getGraph().getNode(this.measuredNodes.get(this.spouses[i][j]).toString());
                    // Node spouseErrorTerm = null;
                    // for (Iterator it = semIm.getEstIm().getGraph()
                    // .getParents(spouse).iterator(); it.hasNext();) {
                    // Node nextParent = (Node) it.next();
                    // if (nextParent.getNodeType() == NodeType.ERROR) {
                    // spouseErrorTerm = nextParent;
                    // break;
                    // }
                    // }
                    Node spouseErrorTerm = semIm.getSemPm().getGraph().getExogenous(spouse);
                    semIm.setParamValue(nodeErrorTerm, spouseErrorTerm, this.covErrors[i][this.spouses[i][j]]);
                }
            }
        }
        for (int i = 0; i < this.numLatent; i++) {
            Node node = semIm.getSemPm().getGraph().getNode(this.latentNodes.get(i).toString());
            if (semIm.getSemPm().getGraph().getParents(node).size() == 0) {
                semIm.setParamValue(node, node, this.varErrorLatent[i]);
            } else {
                for (Iterator it = semIm.getSemPm().getGraph().getParents(node).iterator(); it.hasNext(); ) {
                    Node nextParent = (Node) it.next();
                    if (nextParent.getNodeType() == NodeType.ERROR) {
                        semIm.setParamValue(nextParent, nextParent, this.varErrorLatent[i]);
                        break;
                    }
                }
                for (int j = 0; j < this.parentsLat[i].length; j++) {
                    Node parent = semIm.getSemPm().getGraph().getNode(this.latentNodes.get(this.parentsLat[i][j]).toString());
                    semIm.setParamValue(parent, node, this.betasLat[i][parentsLat[i][j]]);
                }
            }
        }
        return -semIm.getTruncLL() - 0.5 * semIm.getNumFreeParams() * Math.log(this.covarianceMatrix.getSampleSize());
    } catch (java.lang.IllegalArgumentException e) {
        System.out.println("** Warning: " + e.toString());
        return -Double.MAX_VALUE;
    }
}

Example 8 with Parameter

use of edu.cmu.tetrad.sem.Parameter in project tetrad by cmu-phil.

the class SemPmGraphicalEditor method beginNodeEdit.

private void beginNodeEdit(Node node) {
    Parameter parameter = getNodeParameter(node);
    if (parameter == null) {
        throw new IllegalStateException("There is no variance parameter in " + "model for node " + node + ".");
    }
    ParameterEditor paramEditor = new ParameterEditor(parameter, semPm());
    paramEditor.addPropertyChangeListener(new PropertyChangeListener() {

        public void propertyChange(PropertyChangeEvent evt) {
            System.out.println("** " + evt.getPropertyName());
            firePropertyChange(evt.getPropertyName(), null, null);
        }
    });
    int ret = JOptionPane.showOptionDialog(workbench(), paramEditor, "Parameter Properties", JOptionPane.OK_CANCEL_OPTION, JOptionPane.PLAIN_MESSAGE, null, null, null);
    if (ret == JOptionPane.OK_OPTION) {
        parameter.setName(paramEditor.getParamName());
        parameter.setFixed(paramEditor.isFixed());
        parameter.setInitializedRandomly(paramEditor.isInitializedRandomly());
        parameter.setStartingValue(paramEditor.getStartingValue());
        resetLabels();
    }
}

Example 9 with Parameter

use of edu.cmu.tetrad.sem.Parameter in project tetrad by cmu-phil.

the class SemPmGraphicalEditor method beginEdgeEdit.

// ========================PRIVATE PROTECTED METHODS======================//
private void beginEdgeEdit(Edge edge) {
    Parameter parameter = getEdgeParameter(edge);
    ParameterEditor paramEditor = new ParameterEditor(parameter, semPm());
    paramEditor.addPropertyChangeListener(new PropertyChangeListener() {

        public void propertyChange(PropertyChangeEvent evt) {
            firePropertyChange("modelChanged", null, null);
        }
    });
    int ret = JOptionPane.showOptionDialog(workbench(), paramEditor, "Parameter Properties", JOptionPane.OK_CANCEL_OPTION, JOptionPane.PLAIN_MESSAGE, null, null, null);
    if (ret == JOptionPane.OK_OPTION) {
        parameter.setName(paramEditor.getParamName());
        parameter.setFixed(paramEditor.isFixed());
        parameter.setInitializedRandomly(paramEditor.isInitializedRandomly());
        parameter.setStartingValue(paramEditor.getStartingValue());
        resetLabels();
    }
}

Example 10 with Parameter

use of edu.cmu.tetrad.sem.Parameter in project tetrad by cmu-phil.

the class SemPmGraphicalEditor method getEquationOfNode.

private String getEquationOfNode(Node node) {
    String eqn = node.getName() + " = B0_" + node.getName();
    SemGraph semGraph = semPm().getGraph();
    List parentNodes = semGraph.getParents(node);
    for (Object parentNodeObj : parentNodes) {
        Node parentNode = (Node) parentNodeObj;
        // Parameter edgeParam = semPm().getEdgeParameter(
        // semGraph.getEdge(parentNode, node));
        Parameter edgeParam = getEdgeParameter(semGraph.getDirectedEdge(parentNode, node));
        if (edgeParam != null) {
            eqn = eqn + " + " + edgeParam.getName() + "*" + parentNode;
        }
    }
    eqn = eqn + " + " + semPm().getGraph().getExogenous(node);
    return eqn;
}