use of edu.cmu.tetrad.util.TetradVector in project tetrad by cmu-phil.
the class LingamPattern2 method getScore1.
// =============================PRIVATE METHODS=========================//
private Score getScore1(Graph dag, List<TetradMatrix> data, List<Node> variables) {
// System.out.println("Scoring DAG: " + dag);
List<Regression> regressions = new ArrayList<>();
for (TetradMatrix _data : data) {
regressions.add(new RegressionDataset(_data, variables));
}
int totalSampleSize = 0;
for (TetradMatrix _data : data) {
totalSampleSize += _data.rows();
}
int numCols = data.get(0).columns();
List<Node> nodes = dag.getNodes();
double score = 0.0;
double[] pValues = new double[nodes.size()];
TetradMatrix absoluteStandardizedResiduals = new TetradMatrix(totalSampleSize, numCols);
for (int i = 0; i < nodes.size(); i++) {
List<Double> _absoluteStandardizedResiduals = new ArrayList<>();
for (int j = 0; j < data.size(); j++) {
Node _target = nodes.get(i);
List<Node> _regressors = dag.getParents(_target);
Node target = getVariable(variables, _target.getName());
List<Node> regressors = new ArrayList<>();
for (Node _regressor : _regressors) {
Node variable = getVariable(variables, _regressor.getName());
regressors.add(variable);
}
RegressionResult result = regressions.get(j).regress(target, regressors);
TetradVector residualsColumn = result.getResiduals();
DoubleArrayList _absoluteStandardizedResidualsColumn = new DoubleArrayList(residualsColumn.toArray());
double mean = Descriptive.mean(_absoluteStandardizedResidualsColumn);
double std = Descriptive.standardDeviation(Descriptive.variance(_absoluteStandardizedResidualsColumn.size(), Descriptive.sum(_absoluteStandardizedResidualsColumn), Descriptive.sumOfSquares(_absoluteStandardizedResidualsColumn)));
for (int i2 = 0; i2 < _absoluteStandardizedResidualsColumn.size(); i2++) {
_absoluteStandardizedResidualsColumn.set(i2, (_absoluteStandardizedResidualsColumn.get(i2) - mean) / std);
_absoluteStandardizedResidualsColumn.set(i2, Math.abs(_absoluteStandardizedResidualsColumn.get(i2)));
}
for (int k = 0; k < _absoluteStandardizedResidualsColumn.size(); k++) {
_absoluteStandardizedResiduals.add(_absoluteStandardizedResidualsColumn.get(k));
}
}
DoubleArrayList absoluteStandardResidualsList = new DoubleArrayList(absoluteStandardizedResiduals.getColumn(i).toArray());
for (int k = 0; k < _absoluteStandardizedResiduals.size(); k++) {
absoluteStandardizedResiduals.set(k, i, _absoluteStandardizedResiduals.get(k));
}
double _mean = Descriptive.mean(absoluteStandardResidualsList);
double diff = _mean - Math.sqrt(2.0 / Math.PI);
score += diff * diff;
}
for (int j = 0; j < absoluteStandardizedResiduals.columns(); j++) {
double[] x = absoluteStandardizedResiduals.getColumn(j).toArray();
double p = new AndersonDarlingTest(x).getP();
pValues[j] = p;
}
return new Score(score, pValues);
}
use of edu.cmu.tetrad.util.TetradVector in project tetrad by cmu-phil.
the class MVPLikelihood method multipleRegression.
private double multipleRegression(TetradVector Y, TetradMatrix X) {
int n = X.rows();
TetradVector r;
if (X.columns() >= n) {
TetradVector ones = new TetradVector(n);
for (int i = 0; i < n; i++) ones.set(i, 1);
r = ones.scalarMult(ones.dotProduct(Y) / (double) n).minus(Y);
} else {
try {
TetradMatrix Xt = X.transpose();
TetradMatrix XtX = Xt.times(X);
r = X.times(XtX.inverse().times(Xt.times(Y))).minus(Y);
} catch (Exception e) {
TetradVector ones = new TetradVector(n);
for (int i = 0; i < n; i++) ones.set(i, 1);
r = ones.scalarMult(ones.dotProduct(Y) / (double) n).minus(Y);
}
}
double sigma2 = r.dotProduct(r) / n;
double lik;
if (sigma2 < 0) {
TetradVector ones = new TetradVector(n);
for (int i = 0; i < n; i++) ones.set(i, 1);
r = ones.scalarMult(ones.dotProduct(Y) / (double) Math.max(n, 2)).minus(Y);
sigma2 = r.dotProduct(r) / n;
lik = -(n / 2) * (Math.log(2 * Math.PI) + Math.log(sigma2) + 1);
} else if (sigma2 == 0) {
lik = 0;
} else {
lik = -(n / 2) * (Math.log(2 * Math.PI) + Math.log(sigma2) + 1);
}
if (Double.isInfinite(lik) || Double.isNaN(lik)) {
System.out.println(lik);
}
return lik;
}
use of edu.cmu.tetrad.util.TetradVector in project tetrad by cmu-phil.
the class MVPLikelihood method getLik.
public double getLik(int child_index, int[] parents) {
double lik = 0;
Node c = variables.get(child_index);
List<ContinuousVariable> continuous_parents = new ArrayList<>();
List<DiscreteVariable> discrete_parents = new ArrayList<>();
if (c instanceof DiscreteVariable && discretize) {
for (int p : parents) {
Node parent = discreteVariables.get(p);
discrete_parents.add((DiscreteVariable) parent);
}
} else {
for (int p : parents) {
Node parent = variables.get(p);
if (parent instanceof ContinuousVariable) {
continuous_parents.add((ContinuousVariable) parent);
} else {
discrete_parents.add((DiscreteVariable) parent);
}
}
}
int p = continuous_parents.size();
List<List<Integer>> cells = adTree.getCellLeaves(discrete_parents);
int[] continuousCols = new int[p];
for (int j = 0; j < p; j++) continuousCols[j] = nodesHash.get(continuous_parents.get(j));
for (List<Integer> cell : cells) {
// for (int[] cell : cells) {
int r = cell.size();
// int r = cell.length;
if (r > 1) {
double[] mean = new double[p];
double[] var = new double[p];
for (int i = 0; i < p; i++) {
for (int j = 0; j < r; j++) {
mean[i] += continuousData[continuousCols[i]][cell.get(j)];
var[i] += Math.pow(continuousData[continuousCols[i]][cell.get(j)], 2);
}
mean[i] /= r;
var[i] /= r;
var[i] -= Math.pow(mean[i], 2);
var[i] = Math.sqrt(var[i]);
if (Double.isNaN(var[i])) {
System.out.println(var[i]);
}
}
int degree = fDegree;
if (fDegree < 1) {
degree = (int) Math.floor(Math.log(r));
}
TetradMatrix subset = new TetradMatrix(r, p * degree + 1);
for (int i = 0; i < r; i++) {
subset.set(i, p * degree, 1);
for (int j = 0; j < p; j++) {
for (int d = 0; d < degree; d++) {
subset.set(i, p * d + j, Math.pow((continuousData[continuousCols[j]][cell.get(i)] - mean[j]) / var[j], d + 1));
}
}
}
if (c instanceof ContinuousVariable) {
TetradVector target = new TetradVector(r);
for (int i = 0; i < r; i++) {
target.set(i, continuousData[child_index][cell.get(i)]);
// target.set(i, continuousData[child_index][cell[i]]);
}
lik += multipleRegression(target, subset);
} else {
TetradMatrix target = new TetradMatrix(r, ((DiscreteVariable) c).getNumCategories());
for (int i = 0; i < r; i++) {
target.set(i, discreteData[child_index][cell.get(i)], 1);
}
lik += approxMultinomialRegression(target, subset);
}
}
}
return lik;
}
use of edu.cmu.tetrad.util.TetradVector 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;
}
use of edu.cmu.tetrad.util.TetradVector 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;
}
Aggregations