pacematrix.java

为了下东西 随便发了个 datamining 的源代码

      double dq[] = h1( pj, ks );
      int pk;
      for( int k = ks+1; k < kp; k++ ){
	pk = p[k];
	h2( pj, ks, dq[1], this, pk);
      }
      h2( pj, ks, dq[1], b, 0 ); // for matrix. ???
      A[ks][pj] = dq[0];
      for( int k = ks+1; k < m; k++ )
	A[k][pj] = 0;
    }
    else {          // removing 
      int pj = p[j];
      for( int i = j; i < ks-1; i++ ) 
	p[i] = p[i+1];
      p[ks-1] = pj;
      double [] cs;
      for( int i = j; i < ks-1; i++ ){
	cs = g1( A[i][p[i]], A[i+1][p[i]] );
	for( int l = i; l < kp; l++ ) 
	  g2( cs, i, i+1, p[l] );
	for( int l = 0; l < b.n; l++ )
	  b.g2( cs, i, i+1, l );
      }
    }
  }
    
  /** Solves upper-triangular equation
   *   	R x = b
   *  @param b the response
   *  @param pvt the pivoting vector
   *  @param kp the number of the first columns involved 
   *  @return  void. On output, the solution is stored in b 
   */
  public void  rsolve( PaceMatrix b, IntVector pvt, int kp) {
    if(kp == 0) b.m = 0;
    int i, j, k;
    int [] p = pvt.getArray();
    double s;
    double [][] ba = b.getArray();
    for( k = 0; k < b.n; k++ ) {
      ba[kp-1][k] /= A[kp-1][p[kp-1]];
      for( i = kp - 2; i >= 0; i-- ){
	s = 0;
	for( j = i + 1; j < kp; j++ )
	  s += A[i][p[j]] * ba[j][k];
	ba[i][k] -= s;
	ba[i][k] /= A[i][p[i]];
      }
    } 
    b.m = kp;
  }
    
  /** Returns a new matrix which binds two matrices together with rows. 
   *  @param b  the second matrix
   *  @return the combined matrix
   */
  public PaceMatrix  rbind( PaceMatrix b ){
    if( n != b.n ) 
      throw new IllegalArgumentException("unequal numbers of rows.");
    PaceMatrix c = new PaceMatrix( m + b.m, n );
    c.setMatrix( 0, m - 1, 0, n - 1, this );
    c.setMatrix( m, m + b.m - 1, 0, n - 1, b );
    return c;
  }

  /** Returns a new matrix which binds two matrices with columns.
   *  @param b the second matrix 
   *  @return the combined matrix
   */
  public PaceMatrix  cbind( PaceMatrix b ) {
    if( m != b.m ) 
      throw new IllegalArgumentException("unequal numbers of rows: " + 
					 m + " and " + b.m);
    PaceMatrix c = new PaceMatrix(m, n + b.n);
    c.setMatrix( 0, m - 1, 0, n - 1, this );
    c.setMatrix( 0, m - 1, n, n + b.n - 1, b );
    return c;
  }

  /** Solves the nonnegative linear squares problem. That is, <p>
   *   <center>   min || A x - b||, subject to x >= 0.  </center> <p>
   * 
   *  For algorithm, refer to P161, Chapter 23 of C. L. Lawson and
   *  R. J. Hanson (1974).  "Solving Least Squares
   *  Problems". Prentice-Hall.
   * 	@param b the response
   *  @param pvt vector storing pivoting column indices
   *	@return solution */
  public DoubleVector nnls( PaceMatrix b, IntVector pvt ) {
    int j, t, counter = 0, jm = -1, n = pvt.size();
    double ma, max, alpha, wj;
    int [] p = pvt.getArray();
    DoubleVector x = new DoubleVector( n );
    double [] xA = x.getArray();
    PaceMatrix z = new PaceMatrix(n, 1);
    PaceMatrix bt;
	
    // step 1 
    int kp = 0; // #variables in the positive set P
    while ( true ) {         // step 2 
      if( ++counter > 3*n )  // should never happen
	throw new RuntimeException("Does not converge");
      t = -1;
      max = 0.0;
      bt = new PaceMatrix( b.transpose() );
      for( j = kp; j <= n-1; j++ ) {   // W = A' (b - A x) 
	wj = bt.times( 0, kp, m-1, this, p[j] );
	if( wj > max ) {        // step 4
	  max = wj;
	  t = j;
	}
      }
	    
      // step 3 
      if ( t == -1) break; // optimum achieved 
	    
      // step 5 
      pvt.swap( kp, t );       // move variable from set Z to set P
      kp++;
      xA[kp-1] = 0;
      steplsqr( b, pvt, kp-1, kp-1, true );
      // step 6
      ma = 0;
      while ( ma < 1.5 ) {
	for( j = 0; j <= kp-1; j++ ) z.A[j][0] = b.A[j][0];
	rsolve(z, pvt, kp); 
	ma = 2; jm = -1;
	for( j = 0; j <= kp-1; j++ ) {  // step 7, 8 and 9
	  if( z.A[j][0] <= 0.0 ) { // alpha always between 0 and 1
	    alpha = xA[j] / ( xA[j] - z.A[j][0] ); 
	    if( alpha < ma ) {
	      ma = alpha; jm = j;
	    }
	  }
	}
	if( ma > 1.5 ) 
	  for( j = 0; j <= kp-1; j++ ) xA[j] = z.A[j][0];  // step 7 
	else { 
	  for( j = kp-1; j >= 0; j-- ) { // step 10
	    // Modified to avoid round-off error (which seemingly 
	    // can cause infinite loop).
	    if( j == jm ) { // step 11 
	      xA[j] = 0.0;
	      steplsqr( b, pvt, kp, j, false );
	      kp--;  // move variable from set P to set Z
	    }
	    else xA[j] += ma * ( z.A[j][0] - xA[j] );
	  }
	}
      }
    }
    x.setSize(kp);
    pvt.setSize(kp);
    return x;
  }

  /** Solves the nonnegative least squares problem with equality
   *	constraint. That is, <p>
   *  <center> min ||A x - b||, subject to x >= 0 and c x = d. </center> <p>
   *
   *	@param b the response
   *  @param c coeficients of equality constraints
   *  @param d constants of equality constraints
   *	@param pvt vector storing pivoting column indices */
  public DoubleVector nnlse( PaceMatrix b, PaceMatrix c, PaceMatrix d, 
			     IntVector pvt ) {
    double eps = 1e-10 * Math.max( c.maxAbs(), d.maxAbs() ) /
    Math.max( maxAbs(), b.maxAbs() );
	
    PaceMatrix e = c.rbind( new PaceMatrix( times(eps) ) );
    PaceMatrix f = d.rbind( new PaceMatrix( b.times(eps) ) );

    return e.nnls( f, pvt );
  }

  /** Solves the nonnegative least squares problem with equality
   *	constraint. That is, <p>
   *  <center> min ||A x - b||,  subject to x >= 0 and || x || = 1. </center>
   *  <p>
   *  @param b the response 
   *	@param pvt vector storing pivoting column indices */
  public DoubleVector nnlse1( PaceMatrix b, IntVector pvt ) {
    PaceMatrix c = new PaceMatrix( 1, n, 1 );
    PaceMatrix d = new PaceMatrix( 1, b.n, 1 );
	
    return nnlse(b, c, d, pvt);
  }

  /** Generate matrix with standard-normally distributed random elements
      @param m    Number of rows.
      @param n    Number of colums.
      @return An m-by-n matrix with random elements.  */
  public static Matrix randomNormal( int m, int n ) {
    Random random = new Random();
     
    Matrix A = new Matrix(m,n);
    double[][] X = A.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
	X[i][j] = random.nextGaussian();
      }
    }
    return A;
  }

  public static void  main( String args[] ) {
    System.out.println("================================================" + 
		       "===========");
    System.out.println("To test the pace estimators of linear model\n" + 
		       "coefficients.\n");

    double sd = 2;     // standard deviation of the random error term
    int n = 200;       // total number of observations
    double beta0 = 100;   // intercept
    int k1 = 20;       // number of coefficients of the first cluster
    double beta1 = 0;  // coefficient value of the first cluster
    int k2 = 20;      // number of coefficients of the second cluster
    double beta2 = 5; // coefficient value of the second cluster 
    int k = 1 + k1 + k2;

    DoubleVector beta = new DoubleVector( 1 + k1 + k2 );
    beta.set( 0, beta0 );
    beta.set( 1, k1, beta1 );
    beta.set( k1+1, k1+k2, beta2 );

    System.out.println("The data set contains " + n + 
		       " observations plus " + (k1 + k2) + 
		       " variables.\n\nThe coefficients of the true model"
		       + " are:\n\n" + beta );
	
    System.out.println("\nThe standard deviation of the error term is " + 
		       sd );
	
    System.out.println("===============================================" 
		       + "============");
		
    PaceMatrix X = new PaceMatrix( n, k1+k2+1 );
    X.setMatrix( 0, n-1, 0, 0, 1 );
    X.setMatrix( 0, n-1, 1, k1+k2, random(n, k1+k2) );
	
    PaceMatrix Y = new 
      PaceMatrix( X.times( new PaceMatrix(beta) ).
		  plusEquals( randomNormal(n,1).times(sd) ) );

    IntVector pvt = (IntVector) IntVector.seq(0, k1+k2);

    /*System.out.println( "The OLS estimate (by jama.Matrix.solve()) is:\n\n" + 
      (new PaceMatrix(X.solve(Y))).getColumn(0) );*/
	
    X.lsqrSelection( Y, pvt, 1 );
    X.positiveDiagonal( Y, pvt );

    PaceMatrix sol = (PaceMatrix) Y.clone();
    X.rsolve( sol, pvt, pvt.size() );
    DoubleVector betaHat = sol.getColumn(0).unpivoting( pvt, k );
    System.out.println( "\nThe OLS estimate (through lsqr()) is: \n\n" + 
			betaHat );

    System.out.println( "\nQuadratic loss of the OLS estimate (||X b - X bHat||^2) = " + 
			( new PaceMatrix( X.times( new 
			  PaceMatrix(beta.minus(betaHat)) )))
			.getColumn(0).sum2() );

    System.out.println("=============================================" + 
		       "==============");
    System.out.println("             *** Pace estimation *** \n");
    DoubleVector r = Y.getColumn( pvt.size(), n-1, 0);
    double sde = Math.sqrt(r.sum2() / r.size());
	
    System.out.println( "Estimated standard deviation = " + sde );

    DoubleVector aHat = Y.getColumn( 0, pvt.size()-1, 0).times( 1./sde );
    System.out.println("\naHat = \n" + aHat );
	
    System.out.println("\n========= Based on chi-square mixture ============");

    ChisqMixture d2 = new ChisqMixture();
    int method = MixtureDistribution.NNMMethod;
    DoubleVector AHat = aHat.square();
    d2.fit( AHat, method ); 
    System.out.println( "\nEstimated mixing distribution is:\n" + d2 );
	
    DoubleVector ATilde = d2.pace2( AHat );
    DoubleVector aTilde = ATilde.sqrt().times(aHat.sign());
    PaceMatrix YTilde = new 
      PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
    X.rsolve( YTilde, pvt, pvt.size() );
    DoubleVector betaTilde = 
    YTilde.getColumn(0).unpivoting( pvt, k );
    System.out.println( "\nThe pace2 estimate of coefficients = \n" + 
			betaTilde );
    System.out.println( "Quadratic loss = " + 
			( new PaceMatrix( X.times( new 
			  PaceMatrix(beta.minus(betaTilde)) )))
			.getColumn(0).sum2() );
	
    ATilde = d2.pace4( AHat );
    aTilde = ATilde.sqrt().times(aHat.sign());
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
    X.rsolve( YTilde, pvt, pvt.size() );
    betaTilde = YTilde.getColumn(0).unpivoting( pvt, k );
    System.out.println( "\nThe pace4 estimate of coefficients = \n" + 
			betaTilde );
    System.out.println( "Quadratic loss = " + 
			( new PaceMatrix( X.times( new 
			  PaceMatrix(beta.minus(betaTilde)) )))
			.getColumn(0).sum2() );
	
    ATilde = d2.pace6( AHat );
    aTilde = ATilde.sqrt().times(aHat.sign());
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
    X.rsolve( YTilde, pvt, pvt.size() );
    betaTilde = YTilde.getColumn(0).unpivoting( pvt, k );
    System.out.println( "\nThe pace6 estimate of coefficients = \n" + 
			betaTilde );
    System.out.println( "Quadratic loss = " + 
			( new PaceMatrix( X.times( new 
			  PaceMatrix(beta.minus(betaTilde)) )))
			.getColumn(0).sum2() );
	
    System.out.println("\n========= Based on normal mixture ============");
	
    NormalMixture d = new NormalMixture();
    d.fit( aHat, method ); 
    System.out.println( "\nEstimated mixing distribution is:\n" + d );
	
    aTilde = d.nestedEstimate( aHat );
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
    X.rsolve( YTilde, pvt, pvt.size() );
    betaTilde = YTilde.getColumn(0).unpivoting( pvt, k );
    System.out.println( "The nested estimate of coefficients = \n" + 
			betaTilde );
    System.out.println( "Quadratic loss = " + 
			( new PaceMatrix( X.times( new 
			  PaceMatrix(beta.minus(betaTilde)) )))
			.getColumn(0).sum2() );
	
	
    aTilde = d.subsetEstimate( aHat );
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
    X.rsolve( YTilde, pvt, pvt.size() );
    betaTilde = 
    YTilde.getColumn(0).unpivoting( pvt, k );
    System.out.println( "\nThe subset estimate of coefficients = \n" + 
			betaTilde );
    System.out.println( "Quadratic loss = " + 
			( new PaceMatrix( X.times( new 
			  PaceMatrix(beta.minus(betaTilde)) )))
			.getColumn(0).sum2() );
	
    aTilde = d.empiricalBayesEstimate( aHat );
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
    X.rsolve( YTilde, pvt, pvt.size() );
    betaTilde = YTilde.getColumn(0).unpivoting( pvt, k );
    System.out.println( "\nThe empirical Bayes estimate of coefficients = \n"+
			betaTilde );
	
    System.out.println( "Quadratic loss = " + 
			( new PaceMatrix( X.times( new 
			  PaceMatrix(beta.minus(betaTilde)) )))
			.getColumn(0).sum2() );
	
  }
}