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);
}
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;
}
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);
}
}
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);
}
}
Aggregations