get rid of implicit dc term in fitters

JohnMount · JohnMount · commit 9bd524ba789d · 2014-09-26T16:32:24.000-07:00
diff --git a/Count/src/com/mzlabs/count/ctab/CTab.java b/Count/src/com/mzlabs/count/ctab/CTab.java
@@ -215,7 +215,6 @@ public static void main(final String[] args) {
 		System.out.println("n" + "\t" + "total" + "\t" + "target" + "\t" + "count" + "\t" + "date" + "\t" + "cacheSizes" + "\t" + "tableFinishTimeEst");
 		for(int n=1;n<=9;++n) {
 			final CTab ctab = new CTab(n,true);
-			//final NewtonFitter lf = new NewtonFitter(new SquareLossOfExp());
 			final NewtonFitter lf = new NewtonFitter(DirectPoissonJacobian.poissonLink);
 			final int tLast = (n*n-3*n+2)/2;
 			for(int total=0;total<=tLast;++total) {
@@ -226,14 +225,14 @@ public static void main(final String[] args) {
 				long remainingTimeEstMS = 10000;
 				if(total>2) { 
 					// simplistic model: time ~ exp(a + b*size)
-					final double[] x = { total };
+					final double[] x = { 1, total };
 					final double y = 10000.0+curTime.getTime() - startTime.getTime();
 					lf.addObservation(x,y,1.0);
 					if(total>6) {
 						final double[] beta = lf.solve();
 						double timeEstMS = 0.0;
 						for(int j=total+1;j<=tLast;++j) {
-							final double predict = lf.link.evalEst(beta,new double[] {j});
+							final double predict = lf.link.evalEst(beta,new double[] {1,j});
 							timeEstMS += predict;
 						}
 						remainingTimeEstMS = (long)Math.ceil(timeEstMS);
diff --git a/Count/src/com/mzlabs/fit/BalanceJacobianCalc.java b/Count/src/com/mzlabs/fit/BalanceJacobianCalc.java
@@ -4,15 +4,15 @@ public interface BalanceJacobianCalc {
 	/**
 	 * 
 	 * @param obs
-	 * @param beta beta.length==obs.x.length+1
+	 * @param beta beta.length==obs.x.length
 	 * @return balance eqns and Jacobian coefficients of underlying loss function for parametes beta with respect to datum obs
 	 */
 	BalanceJacobianCoef calc(Obs obs, double[] beta);
 	
 	/**
 	 * 
 	 * @param beta
-	 * @param x x.length==beta.length-1 (last beta is the coefficient matchng an implicit constant term by convention)
+	 * @param x x.length==beta.length (last beta is the coefficient matchng an implicit constant term by convention)
 	 * @return the estimate of y given parameters beta and datum x
 	 */
 	double evalEst(double[] beta, double[] x);
diff --git a/Count/src/com/mzlabs/fit/GLMModel.java b/Count/src/com/mzlabs/fit/GLMModel.java
@@ -24,10 +24,10 @@ public void balanceAndJacobian(final Iterable<Obs> obs, final double[] beta,
 		for(final Obs obsi: obs) {
 			final BalanceJacobianCoef ghc = balanceJacobianCalc.calc(obsi,beta);
 			for(int i=0;i<dim;++i) {
-				final double xi = i<dim-1?obsi.x[i]:1.0;
+				final double xi = obsi.x[i];
 				balance[i] += ghc.balanceCoef*xi;
 				for(int j=0;j<dim;++j) {
-					final double xj = j<dim-1?obsi.x[j]:1.0;
+					final double xj = obsi.x[j];
 					final double hij = jacobian.get(i,j);
 					jacobian.set(i,j,hij+xi*xj*ghc.jacobianCoef);
 				}
diff --git a/Count/src/com/mzlabs/fit/LinearFitter.java b/Count/src/com/mzlabs/fit/LinearFitter.java
@@ -19,24 +19,24 @@ public final class LinearFitter implements Fitter {
 	 */
 	public LinearFitter(final int n) {
 		final LinalgFactory<ColtMatrix> factory = ColtMatrix.factory;
-		xTx = factory.newMatrix(n+1,n+1,false);
-		xTy = new double[n+1];
+		xTx = factory.newMatrix(n,n,false);
+		xTy = new double[n];
 	}
 
 	/* (non-Javadoc)
 	 * @see com.mzlabs.count.util.Fitter#addObservation(double[], double, double)
 	 */
 	@Override
 	public void addObservation(final double[] x, final double y, final double wt) {
-		final int n = xTx.rows()-1;
+		final int n = xTx.rows();
 		if(n!=x.length) {
 			throw new IllegalArgumentException();
 		}
-		for(int i=0;i<=n;++i) {
-			final double xi = i<n?x[i]:1.0;
+		for(int i=0;i<n;++i) {
+			final double xi = x[i];
 			xTy[i] += wt*xi*y;
-			for(int j=0;j<=n;++j) {
-				final double xj = j<n?x[j]:1.0;
+			for(int j=0;j<n;++j) {
+				final double xj = x[j];
 				xTx.set(i,j,xTx.get(i, j)+wt*xi*xj);
 			}
 		}
@@ -47,15 +47,15 @@ public void addObservation(final double[] x, final double y, final double wt) {
 	 */
 	@Override
 	public double[] solve() {
-		final int n = xTx.rows()-1;
+		final int n = xTx.rows();
 		final double epsilon = 1.0e-5;
 		final double[] xTxii = new double[n+1];
-		for(int i=0;i<=n;++i) {
+		for(int i=0;i<n;++i) {
 			xTxii[i] = xTx.get(i,i);
 			xTx.set(i,i,xTxii[i]+epsilon);  // Ridge term
 		}
 		final double[] soln = xTx.solve(xTy);
-		for(int i=0;i<=n;++i) {
+		for(int i=0;i<n;++i) {
 			xTx.set(i,i,xTxii[i]);
 		}
 		return soln;
@@ -65,14 +65,13 @@ public double[] solve() {
 	 * @see com.mzlabs.count.util.Fitter#predict(double[], double[])
 	 */
 	public double predict(final double[] soln, final double[] x) {
-		final int n = soln.length-1;
-		if((n!=x.length)||(n+1!=soln.length)) {
+		final int n = soln.length;
+		if((n!=x.length)||(n!=soln.length)) {
 			throw new IllegalArgumentException();
 		}
 		double sum = 0.0;
-		for(int i=0;i<=n;++i) {
-			final double xi = i<n?x[i]:1.0;
-			sum += xi*soln[i];
+		for(int i=0;i<n;++i) {
+			sum += x[i]*soln[i];
 		}
 		return sum;
 	}
diff --git a/Count/src/com/mzlabs/fit/NewtonFitter.java b/Count/src/com/mzlabs/fit/NewtonFitter.java
@@ -34,9 +34,9 @@ public void addObservation(final double[] x, final double y, final double wt) {
 	@Override
 	public double[] solve() {
 		final LinalgFactory<ColtMatrix> factory = ColtMatrix.factory;
-		final int dim = obs.get(0).x.length+1;
+		final int dim = obs.get(0).x.length;
 		// roughly: often solving y ~ f(b.x), so start at f^-1(y) ~ b.x
-		final Fitter sf = new LinearFitter(dim-1);
+		final Fitter sf = new LinearFitter(dim);
 		for(final Obs obsi: obs) {
 			sf.addObservation(obsi.x, link.heuristicLink(obsi.y), obsi.wt);
 		}
diff --git a/Count/src/com/mzlabs/fit/Obs.java b/Count/src/com/mzlabs/fit/Obs.java
@@ -15,26 +15,25 @@ public Obs(final double[] x, final double y, final double wt) {
 	
 	/**
 	 * 
-	 * @param soln soln.length==x.length+1
+	 * @param soln soln.length==x.length
 	 * @param x
 	 * @return
 	 */
 	public static double dot(final double soln[], final double[] x) {
 		final int n = x.length;
-		if(soln.length!=n+1) {
+		if(soln.length!=n) {
 			throw new IllegalArgumentException();
 		}
 		double sum = 0.0;
-		for(int i=0;i<=n;++i) {
-			final double xi = i<n?x[i]:1.0;
-			sum += xi*soln[i];
+		for(int i=0;i<n;++i) {
+			sum += x[i]*soln[i];
 		}
-		return sum;		
+		return sum;
 	}
 	
 	/**
 	 * 
-	 * @param soln.length==x.length+1
+	 * @param soln.length==x.length
 	 * @return
 	 */
 	public double dot(final double[] soln) {
diff --git a/Count/src/com/mzlabs/fit/SquareLossOfExp.java b/Count/src/com/mzlabs/fit/SquareLossOfExp.java
@@ -37,10 +37,10 @@ public void balanceAndJacobian(final Iterable<Obs> obs, final double[] beta,
 			final double gradCoef = -2*diff*ebx*obsi.wt;
 			final double hessCoef = -2*(obsi.y-2*ebx)*ebx*obsi.wt;
 			for(int i=0;i<dim;++i) {
-				final double xi = i<dim-1?obsi.x[i]:1.0;
+				final double xi = obsi.x[i];
 				balance[i] += gradCoef*xi;
 				for(int j=0;j<dim;++j) {
-					final double xj = j<dim-1?obsi.x[j]:1.0;
+					final double xj = obsi.x[j];
 					final double hij = jacobian.get(i,j);
 					jacobian.set(i,j,hij+xi*xj*hessCoef);
 				}
diff --git a/Count/src/com/mzlabs/fit/VectorFnWithJacobian.java b/Count/src/com/mzlabs/fit/VectorFnWithJacobian.java
@@ -11,7 +11,7 @@ public interface VectorFnWithJacobian {
 	/**
 	 * 
 	 * @param beta
-	 * @param x x.length==beta.length-1 (last beta is the coefficient matchng an implicit constant term by convention)
+	 * @param x x.length==beta.length (last beta is the coefficient matchng an implicit constant term by convention)
 	 * @return the estimate of y given parameters beta and datum x
 	 */
 	double evalEst(double[] beta, double[] x);
diff --git a/Count/tests/com/mzlabs/count/util/TestLogLinFitter.java b/Count/tests/com/mzlabs/count/util/TestLogLinFitter.java
@@ -19,14 +19,14 @@
 public class TestLogLinFitter {
 	@Test
 	public void testLFit() {
-		final LinearFitter lf = new LinearFitter(1);
+		final LinearFitter lf = new LinearFitter(2);
 		final NewtonFitter llf = new NewtonFitter(new SquareLossOfExp());
 		final Random rand = new Random(343406L);
 		final ArrayList<Obs> obs = new ArrayList<Obs>();
 		for(int i=1;i<7;++i) {
 			final double y = Math.exp(2.0*i);
 			for(int j=0;j<10;++j) {
-				final double[] x = new double[] {i};
+				final double[] x = new double[] {1,i};
 				final double yObserved = y*(1+0.3*rand.nextGaussian());
 				llf.addObservation(x,yObserved,1.0);
 				lf.addObservation(x,Math.log(Math.max(1.0,yObserved)),1.0);
@@ -55,7 +55,7 @@ public void testLFit() {
 	@Test
 	public void testPLinks() {
 		final double y = 1.55528;
-		final double[] x = { 5.0 };
+		final double[] x = { 1.0, 5.0 };
 		final Obs obs = new Obs(x,y,1.0);
 		final double[] beta = { 0.2 , -0.1};
 		final BalanceJacobianCoef lpgh = LinkBasedGradHess.poissonGradHess.calc(obs, beta);
@@ -72,7 +72,7 @@ public void testPFit() {
 		final Random rand = new Random(343406L);
 		for(int i=1;i<=5;++i) {
 			final double y = Math.exp(0.4*i) + rand.nextGaussian();
-			final double[] x = new double[] {i};
+			final double[] x = new double[] {1,i};
 			llf.addObservation(x,y,1.0);
 			obs.add(new Obs(x,y,1.0));
 		}
@@ -89,8 +89,7 @@ public void testPFit() {
 			final double[] x = obsi.x;
 			final double llfit = llf.link.evalEst(llsoln,x);
 			for(int i=0;i<dim;++i) {
-				final double xi = i<dim-1?obsi.x[i]:1.0;
-				sums[i] += obsi.wt*xi*(y-llfit);
+				sums[i] += obsi.wt*x[i]*(y-llfit);
 			}
 		}
 		for(final double si: sums) {

Original file line number	Diff line number	Diff line change
`@@ -11,7 +11,7 @@ public interface VectorFnWithJacobian {`
`11`	`11`	`/**`
`12`	`12`	`*`
`13`	`13`	`* @param beta`
`14`		`- * @param x x.length==beta.length-1 (last beta is the coefficient matchng an implicit constant term by convention)`
	`14`	`+ * @param x x.length==beta.length (last beta is the coefficient matchng an implicit constant term by convention)`
`15`	`15`	`* @return the estimate of y given parameters beta and datum x`
`16`	`16`	`*/`
`17`	`17`	`double evalEst(double[] beta, double[] x);`