List of usage examples for org.apache.commons.math3.linear RealVector outerProduct
public RealMatrix outerProduct(RealVector v)
From source file:org.rhwlab.BHC.LogNode.java
public double alternative(List<RealVector> data) { int n = data.size(); double nP = n + beta; if (n > maxN) { return 0.0; }//from w ww .j a v a 2s.c o m int D = data.get(0).getDimension(); RealVector xSum = new ArrayRealVector(D); // a vector of zeros RealMatrix XXT = new Array2DRowRealMatrix(D, D); for (RealVector v : data) { xSum = xSum.add(v); RealMatrix v2 = v.outerProduct(v); XXT = XXT.add(v2); } RealMatrix Sp = S.add(XXT); Sp = Sp.add(m.outerProduct(m).scalarMultiply(beta * n / (nP))); Sp = Sp.add(xSum.outerProduct(xSum).scalarMultiply(1.0 / (nP))); Sp = Sp.subtract(m.outerProduct(xSum).add(xSum.outerProduct(m)).scalarMultiply(beta / (nP))); LUDecomposition ed = new LUDecomposition(Sp); if (!ed.getSolver().isNonSingular()) { System.exit(10); } return Utils.eln(ed.getDeterminant()); }
From source file:org.rhwlab.BHC.LogNode.java
public double logMarginalLikelihood(List<RealVector> data) { int n = data.size(); if (n > maxN) { return 0.0; }/*from w w w . ja v a 2s. c o m*/ int D = data.get(0).getDimension(); double rP = beta + n; Double lnrP = Utils.eln(rP); double nuP = nu + n; RealMatrix C = new Array2DRowRealMatrix(D, D); RealVector X = new ArrayRealVector(D); // a vector of zeros for (RealVector v : data) { X = X.add(v); RealMatrix v2 = v.outerProduct(v); C = C.add(v2); } RealVector mP = (m.mapMultiply(beta).add(X)).mapDivide(rP); RealMatrix Sp = C.add(S); RealMatrix rmmP = mP.outerProduct(mP).scalarMultiply(rP); Sp = Sp.add(rmm).subtract(rmmP); LUDecomposition ed = new LUDecomposition(Sp); if (!ed.getSolver().isNonSingular()) { System.exit(10); } double logDetSp = Utils.eln(ed.getDeterminant()); // double logDetSp = alternative(data); Double ret = Utils.elnPow(logPi, -n * D / 2.0); ret = Utils.elnMult(ret, Utils.elnPow(lnBeta, D / 2.0)); ret = Utils.elnMult(ret, -Utils.elnPow(lnrP, D / 2.0)); ret = Utils.elnMult(ret, logdetSnu); ret = Utils.elnMult(ret, -Utils.elnPow(logDetSp, nuP / 2.0)); for (int i = 1; i <= D; ++i) { ret = Utils.elnMult(ret, Gamma.logGamma((nuP + 1 - i) / 2.0)); ret = Utils.elnMult(ret, -Gamma.logGamma((nu + 1 - i) / 2.0)); } if (ThreadedAlgorithm.time < 175) { ret = ret + data.size() * logRSDLikelihood(); } return ret; }
From source file:org.rhwlab.BHCnotused.GaussianGIWPrior.java
@Override public RealMatrix getPrecision() { RealMatrix ret = new Array2DRowRealMatrix(m.getDimension(), m.getDimension()); int n = 0;//from w w w.j a v a 2 s. c om for (LabeledFieldVector v : data) { int label = v.getLabel(); RealVector[] vs = source.getClusterVectors(label); for (int i = 0; i < vs.length; ++i) { RealVector del = vs[i].subtract(getMean()); ret = ret.add(del.outerProduct(del)); ++n; } } ret = ret.scalarMultiply(1.0 / n); LUDecomposition lud = new LUDecomposition(ret); RealMatrix prec = lud.getSolver().getInverse(); return prec; }
From source file:org.rhwlab.BHCnotused.StdGaussianWishartGaussian.java
public void init() { int n = data.size(); int d = m.getDimension(); double rP = r + n; double nuP = nu + n; RealMatrix C = new Array2DRowRealMatrix(d, d); RealVector X = new ArrayRealVector(d); // a vector of zeros for (RealVector v : data) { X = X.add(v);//from w ww . j a va2s .co m RealMatrix v2 = v.outerProduct(v); C = C.add(v2); } RealVector mP = (m.mapMultiply(r).add(X)).mapDivide(rP); RealMatrix Sp = C.add(S); RealMatrix rmmP = mP.outerProduct(mP).scalarMultiply(rP); Sp = Sp.add(rmm).subtract(rmmP); LUDecomposition ed = new LUDecomposition(Sp); double detSp = Math.pow(ed.getDeterminant(), nuP / 2.0); double gamma = 1.0; double gammaP = 1.0; for (int i = 1; i <= d; ++i) { gamma = gamma * Gamma.gamma((nu + 1 - i) / 2.0); gammaP = gammaP * Gamma.gamma((nuP + 1 - i) / 2.0); } double t1 = Math.pow(Math.PI, -n * d / 2.0); double t2 = Math.pow(r / rP, d / 2.0); double t3 = detS / detSp; double t4 = gammaP / gamma; likelihood = t1 * t2 * t3 * t4; }
From source file:org.rhwlab.dispim.datasource.GaussianComponent.java
public RealMatrix precision(RealVector mu) { RealMatrix ret = new Array2DRowRealMatrix(mu.getDimension(), mu.getDimension()); ret.scalarMultiply(0.0); // make sure it is zero int n = 0;//from w w w.j a v a 2 s . co m for (Integer index : indexes) { RealVector x = source.get(index).coords; RealVector del = x.subtract(mu); ret = ret.add(del.outerProduct(del)); ++n; } ret = ret.scalarMultiply(1.0 / n); LUDecomposition lud = new LUDecomposition(ret); RealMatrix prec = lud.getSolver().getInverse(); return prec; }
From source file:org.rhwlab.dispim.datasource.MicroCluster.java
License:asdf
public static RealMatrix precision(List<MicroCluster> data, RealVector mu) { RealMatrix ret = new Array2DRowRealMatrix(mu.getDimension(), mu.getDimension()); RealVector v = new ArrayRealVector(mu.getDimension()); long n = 0;/* www . jav a 2 s . c o m*/ for (MicroCluster micro : data) { for (int p = 0; p < micro.points.length; ++p) { for (int d = 0; d < mu.getDimension(); ++d) { v.setEntry(d, micro.points[p][d]); } RealVector del = v.subtract(mu); ret = ret.add(del.outerProduct(del)); ++n; } } ret = ret.scalarMultiply(1.0 / n); LUDecomposition lud = new LUDecomposition(ret); RealMatrix prec = null; if (lud.getSolver().isNonSingular()) { prec = lud.getSolver().getInverse(); } return prec; }
From source file:org.rhwlab.dispim.datasource.MicroClusterDataSource.java
public RealMatrix getDataVariance() { RealMatrix ret = new Array2DRowRealMatrix(getD(), getD()); for (int k = 0; k < micros.length; ++k) { RealVector v = micros[k].asRealVector(); ret = ret.add(v.outerProduct(v)); }/*w w w .ja v a2 s.c om*/ return ret.scalarMultiply(1.0 / micros.length); }
From source file:org.rhwlab.variationalbayesian.GaussianMixture.java
public void statistics() { // calculate N calculateN();//from ww w. j a v a 2s . c o m // calculate xBar for (int k = 0; k < K; ++k) { if (inModel[k]) { xBar[k].set(0.0); for (int n = 0; n < X.getN(); ++n) { Voxel vox = X.get(n); double intensity = vox.getAdjusted(); RealVector x = vox.coords; // double rnk = (intensity-background)*r.getEntry(n, k); double rnk = intensity * r.getEntry(n, k); xBar[k] = xBar[k].add(x.mapMultiply(rnk)); } xBar[k].mapDivideToSelf(N[k]); } } // calculate the S for (int k = 0; k < K; ++k) { if (inModel[k]) { if (k == 31) { int iusahdfis = 90; } for (int row = 0; row < X.getD(); ++row) { for (int col = 0; col < X.getD(); ++col) { S[k].setEntry(row, col, 0.0); } } for (int n = 0; n < X.getN(); ++n) { Voxel vox = X.get(n); RealVector x = vox.coords; RealVector del = x.subtract(xBar[k]); RealMatrix mat = del.outerProduct(del); double rnk = r.getEntry(n, k); if (rnk > 0) { int iuhsafuid = 0; } // S[k] = S[k].add(mat.scalarMultiply((intensity-background)*rnk)); S[k] = S[k].add(mat.scalarMultiply(vox.getAdjusted() * rnk)); } S[k] = S[k].scalarMultiply(1.0 / N[k]); } } int uisahdfius = 0; }
From source file:org.rhwlab.variationalbayesian.GaussianMixture.java
public void Mstep() { alphaBetaNu();/*from w w w. j a va 2 s . c o m*/ // compute the means of the components for (int k = 0; k < K; ++k) { if (inModel[k]) { RealVector v1 = xBar[k].mapMultiply(N[k]); RealVector v2 = mu0.mapMultiply(beta0); RealVector v3 = v2.add(v1); m[k] = v3.mapDivide(beta[k]); } } // compute the precision matrices for (int k = 0; k < K; ++k) { if (inModel[k]) { RealVector del = xBar[k].subtract(mu0); RealMatrix del2 = del.outerProduct(del); double f = beta0 * N[k] / (beta0 + N[k]); RealMatrix mat = del2.scalarMultiply(f); RealMatrix NS = S[k].scalarMultiply(N[k]); Winverse[k] = W0inverse.add(NS).add(mat); LUDecomposition cd = new LUDecomposition(Winverse[k]); W[k] = cd.getSolver().getInverse(); detWinv[k] = cd.getDeterminant(); detW[k] = 1.0 / cd.getDeterminant(); } } lambdaTilde(); }
From source file:org.rhwlab.variationalbayesian.SuperVoxelGaussianMixture.java
public void statistics() { // calculate N for (int k = 0; k < K; ++k) { N[k] = 0.0;/*w ww. j av a 2 s. c o m*/ for (int n = 0; n < X.getN(); ++n) { SuperVoxel sv = X.get(n); N[k] = N[k] + sv.getVoxels().length * r.getEntry(n, k); } N[k] = N[k] + .000001; } // calculate xBar for (int k = 0; k < K; ++k) { xBar[k].set(0.0); for (int n = 0; n < X.getN(); ++n) { SuperVoxel sv = X.get(n); double rnk = r.getEntry(n, k); RealVector x = sv.getCenter().mapMultiply(rnk * sv.getVoxels().length); xBar[k] = xBar[k].add(x); } xBar[k].mapDivideToSelf(N[k]); } // calculate the S for (int k = 0; k < K; ++k) { for (int row = 0; row < X.getD(); ++row) { for (int col = 0; col < X.getD(); ++col) { S[k].setEntry(row, col, 0.0); } } for (int n = 0; n < X.getN(); ++n) { SuperVoxel sv = X.get(n); for (RealVector vox : sv.getVoxels()) { RealVector del = vox.subtract(xBar[k]); S[k] = S[k].add(del.outerProduct(del).scalarMultiply(r.getEntry(n, k))); } } S[k] = S[k].scalarMultiply(1.0 / N[k]); } }