work family of fitting examples

Author: JohnMount

Count/expFit.tsv

@@ -15,6 +15,11 @@ public final class BalanceBasedJacobian implements BalanceJacobianCalc {
 	public BalanceBasedJacobian(final Link link) {
 		this.link = link;
 	}
+	
+	@Override
+	public String toString() {
+		return link.toString();
+	}
 
 	@Override
 	public BalanceJacobianCoef calc(final Obs obs, final double[] beta) {

Count/src/com/mzlabs/fit/BalanceBasedJacobian.java

@@ -18,6 +18,11 @@ public BalanceJacobianCoef calc(Obs obs, double[] beta) {
 		return new BalanceJacobianCoef(balCoef,jacCoef);
 	}
 
+	@Override
+	public String toString() {
+		return "PoissonRegression(log-link)";
+	}
+	
 	@Override
 	public double evalEst(double[] beta, double[] x) {
 		return Math.exp(Obs.dot(beta,x));

Count/src/com/mzlabs/fit/DirectPoissonJacobian.java

@@ -1,7 +1,10 @@
 package com.mzlabs.fit;
 
+import java.util.ArrayList;
 import java.util.Arrays;
+import java.util.Random;
 
+import com.winvector.linalg.LinalgFactory;
 import com.winvector.linalg.colt.ColtMatrix;
 
 public final class ExpectationAndSqLoss implements VectorFnWithJacobian {
@@ -24,7 +27,7 @@ public static double dotNE(double[] beta, double[] x) {
 	}
 
 	@Override
-	public double evalEst(double[] beta, double[] x) {
+	public double evalEst(final double[] beta, final double[] x) {
 		return link.invLink(dotNE(beta,x));
 	}
 
@@ -86,4 +89,104 @@ public void balanceAndJacobian(final Iterable<Obs> obs, final double[] beta,
 	public int dim(final Obs obs) {
 		return obs.x.length + 1;
 	}
+	
+	@Override
+	public String toString() {
+		return "ExpectationAndSquareLoss(" + link + ")";
+	}
+	
+	public static void main(final String[] args) {
+		// build some example data
+		final Random rand = new Random(343406L);
+		final ArrayList<Obs> obs = new ArrayList<Obs>();
+		for(int i=0;i<200;++i) {
+			final double[] x = new double[] {1, rand.nextGaussian(), rand.nextGaussian()};
+			final double y = Math.max(0.1,Math.exp(x[0] + 2.0+x[1] + 3.0+x[2]) + 10*rand.nextDouble());
+			obs.add(new Obs(x,y,1.0));
+		}
+		// provision fitters
+		final ExpectationAndSqLoss fn = new ExpectationAndSqLoss(LinkBasedGradHess.logLink);
+		final int elIndex = 2;
+		final NewtonFitter[] fitters = { 
+				new NewtonFitter(new SquareLossOfExp()),
+				new NewtonFitter(DirectPoissonJacobian.poissonLink),
+				new NewtonFitter(fn)  // elIndex == 2
+		};
+		final int dim = obs.get(0).x.length;
+		final int nTrain = obs.size()/2;
+		final LinearFitter lf = new LinearFitter(dim);
+		// scan data
+		for(int i=0;i<nTrain;++i) {
+			final Obs obsi = obs.get(i);
+			lf.addObservation(obsi.x,Math.log(obsi.y),1.0);
+			for(final NewtonFitter fitter: fitters) {
+				fitter.addObservation(obsi.x,obsi.y,1.0);
+			}
+		}
+		// solve
+		final double[] lfSoln = lf.solve();
+		final double[][] fSoln = new double[fitters.length][];
+		for(int j=0;j<fitters.length;++j) {
+			fSoln[j] = fitters[j].solve();
+		}
+		// check balance condition
+		final int pdim = fSoln[elIndex].length;
+		final double[] balance = new double[pdim];
+		final LinalgFactory<ColtMatrix> factory = ColtMatrix.factory;
+		final ColtMatrix jacobian = factory.newMatrix(pdim,pdim,false);
+		fn.balanceAndJacobian(obs.subList(0,nTrain), fSoln[elIndex], balance, jacobian);
+		double balanceCheck = 0.0;
+		for(int i=0;i<nTrain;++i) {
+			final Obs obsi = obs.get(i);
+			final double llfit = fn.evalEst(fSoln[elIndex],obsi.x);
+			balanceCheck += obsi.wt*(obsi.y-llfit);
+		}
+		for(final double bi: balance) {
+			if(Math.abs(bi)>1.0e-3) {
+				throw new IllegalStateException("didn't balance");
+			}
+		}
+		if(Math.abs(balanceCheck)>1.0e-4) {
+			throw new IllegalStateException("didn't balance");
+		}
+		// print data and estimates
+		for(int j=1;j<dim;++j) {
+			System.out.print("x"+j + "\t");
+		}
+		System.out.print("y" + "\t" + "isTrain" + "\t" + "logYest");
+		for(int j=0;j<fitters.length;++j) {
+			System.out.print("\t" + fitters[j].fn);
+		}
+		System.out.println();
+		for(int i=0;i<obs.size();++i) {
+			final Obs obsi = obs.get(i);
+			for(int j=1;j<dim;++j) {
+				System.out.print(obsi.x[j] + "\t");
+			}
+			System.out.print(obsi.y + "\t" + (i<nTrain?1:0) + "\t" + Math.exp(lf.evalEst(lfSoln,obsi.x)));
+			for(int j=0;j<fitters.length;++j) {
+				System.out.print("\t" + fitters[j].evalEst(fSoln[j],obsi.x));
+			}
+			System.out.println();
+		}
+		/**
+		 R: steps
+		 	
+		 	library(ggplot2)
+ 			library(reshape2)
+ 			d <- read.table('expFit.tsv',sep='\t',stringsAsFactors=FALSE,header=TRUE)
+ 			ests <- c('logYest','SquareLossOfExp','GLM.PoissonRegression.log.link..','ExpectationAndSquareLoss.log.link.')
+ 			dTrain <- subset(d,isTrain==1)
+ 			dTest <- subset(d,isTrain==0)
+ 			for(v in ests) {
+ 			   print(paste(v,sum(dTrain$y-dTrain[,v]),sum((dTrain$y-dTrain[,v])^2)))
+ 			}
+ 			for(v in ests) {
+ 			   print(paste(v,sum(dTest$y-dTest[,v]),sum((dTest$y-dTest[,v])^2)))
+ 			}
+			dplot <- melt(subset(d,isTrain==1),id.vars=c('x1','x2','isTrain','y'),variable.name='estimate')
+			ggplot(data=dplot,aes(x=value,y=y,color=estimate,shape=estimate)) + geom_point()
+			ggplot(data=dplot,aes(x=value,y=y,color=estimate,shape=estimate)) + geom_point() + scale_x_log10() + scale_y_log10()
+		 */
+	}
 }

Count/src/com/mzlabs/fit/ExpectationAndSqLoss.java

@@ -8,8 +8,10 @@ public interface Fitter {
 	 * @param y
 	 * @param wt weight of observation (set to 1.0 in many cases)
 	 */
-	public abstract void addObservation(final double[] x, final double y,
+	void addObservation(final double[] x, final double y,
 			final double wt);
 
-	public abstract double[] solve();
+	double[] solve();
+	
+	double evalEst(double[] beta, double[] x);
 }

Count/src/com/mzlabs/fit/Fitter.java

@@ -11,6 +11,11 @@ public GLMModel(final BalanceJacobianCalc balanceJacobianCalc) {
 		this.balanceJacobianCalc = balanceJacobianCalc;
 	}
 
+	@Override
+	public String toString() {
+		return "GLM(" + balanceJacobianCalc.toString() + ")";
+	}
+	
 	@Override
 	public void balanceAndJacobian(final Iterable<Obs> obs, final double[] beta,
 			final double[] balance, final ColtMatrix jacobian) {

Count/src/com/mzlabs/fit/GLMModel.java

@@ -61,10 +61,8 @@ public double[] solve() {
 		return soln;
 	}
 
-	/* (non-Javadoc)
-	 * @see com.mzlabs.count.util.Fitter#predict(double[], double[])
-	 */
-	public double predict(final double[] soln, final double[] x) {
+	@Override
+	public double evalEst(final double[] soln, final double[] x) {
 		final int n = soln.length;
 		if((n!=x.length)||(n!=soln.length)) {
 			throw new IllegalArgumentException();

Count/src/com/mzlabs/fit/LinearFitter.java

@@ -15,6 +15,11 @@ public LinkBasedGradHess(final ProbYGivenZ prob, final Link link) {
 		this.link = link;
 	}
 
+	@Override
+	public String toString() {
+		return link.toString();
+	}
+	
 	private static double sq(final double z) {
 		return z*z;
 	}
@@ -77,6 +82,10 @@ public double invLink(final double z) {
 		public double heuristicLink(double y) {
 			return Math.log(Math.abs(y)+1.0); // near log(y), but well behaved
 		}
+		@Override
+		public String toString() {
+			return "log-link";
+		}
 	};
 
 

Count/src/com/mzlabs/fit/LinkBasedGradHess.java

@@ -117,4 +117,9 @@ public double[] solve() {
 		}
 		return beta;
 	}
+
+	@Override
+	public double evalEst(final double[] beta, final double[] x) {
+		return fn.evalEst(beta, x);
+	}
 }

Count/src/com/mzlabs/fit/NewtonFitter.java

@@ -19,6 +19,11 @@ public double evalEst(final double[] beta, final double[] x) {
 	public double heuristicLink(final double y) {
 		return Math.log(Math.abs(y)+1.0); // near log(y), but well behaved
 	}
+	
+	@Override
+	public String toString() {
+		return "SquareLossOfExp";
+	}
 
 
 	@Override