2
* This program is free software; you can redistribute it and/or modify
3
* it under the terms of the GNU General Public License as published by
4
* the Free Software Foundation; either version 2 of the License, or (at
5
* your option) any later version.
7
* This program is distributed in the hope that it will be useful, but
8
* WITHOUT ANY WARRANTY; without even the implied warranty of
9
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
10
* General Public License for more details.
12
* You should have received a copy of the GNU General Public License
13
* along with this program; if not, write to the Free Software
14
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
18
* Copyright (C) 2002 University of Waikato, Hamilton, New Zealand
22
package weka.classifiers.functions.pace;
24
import weka.core.matrix.DoubleVector;
25
import weka.core.matrix.FlexibleDecimalFormat;
26
import weka.core.matrix.IntVector;
27
import weka.core.matrix.Matrix;
28
import weka.core.matrix.Maths;
30
import java.util.Random;
31
import java.text.DecimalFormat;
34
* Class for matrix manipulation used for pace regression. <p>
38
* Wang, Y. (2000). "A new approach to fitting linear models in high
39
* dimensional spaces." PhD Thesis. Department of Computer Science,
40
* University of Waikato, New Zealand. <p>
42
* Wang, Y. and Witten, I. H. (2002). "Modeling for optimal probability
43
* prediction." Proceedings of ICML'2002. Sydney. <p>
45
* @author Yong Wang (yongwang@cs.waikato.ac.nz)
46
* @version $Revision: 1.5 $
48
public class PaceMatrix
51
/** for serialization */
52
static final long serialVersionUID = 2699925616857843973L;
54
/* ------------------------
56
* ------------------------ */
58
/** Construct an m-by-n PACE matrix of zeros.
59
@param m Number of rows.
60
@param n Number of colums.
62
public PaceMatrix( int m, int n ) {
66
/** Construct an m-by-n constant PACE matrix.
67
@param m Number of rows.
68
@param n Number of colums.
69
@param s Fill the matrix with this scalar value.
71
public PaceMatrix( int m, int n, double s ) {
75
/** Construct a PACE matrix from a 2-D array.
76
@param A Two-dimensional array of doubles.
77
@throws IllegalArgumentException All rows must have the same length
79
public PaceMatrix( double[][] A ) {
83
/** Construct a PACE matrix quickly without checking arguments.
84
@param A Two-dimensional array of doubles.
85
@param m Number of rows.
86
@param n Number of colums.
88
public PaceMatrix( double[][] A, int m, int n ) {
92
/** Construct a PaceMatrix from a one-dimensional packed array
93
@param vals One-dimensional array of doubles, packed by columns (ala Fortran).
94
@param m Number of rows.
95
@throws IllegalArgumentException Array length must be a multiple of m.
97
public PaceMatrix( double vals[], int m ) {
101
/** Construct a PaceMatrix with a single column from a DoubleVector
102
@param v DoubleVector
104
public PaceMatrix( DoubleVector v ) {
106
setMatrix( 0, v.size()-1, 0, v );
109
/** Construct a PaceMatrix from a Matrix
112
public PaceMatrix( Matrix X ) {
113
super( X.getRowDimension(), X.getColumnDimension() );
117
/* ------------------------
119
* ------------------------ */
121
/** Set the row dimenion of the matrix
122
* @param rowDimension the row dimension
124
public void setRowDimension( int rowDimension )
129
/** Set the column dimenion of the matrix
130
* @param columnDimension the column dimension
132
public void setColumnDimension( int columnDimension )
138
* Clone the PaceMatrix object.
142
public Object clone () {
143
PaceMatrix X = new PaceMatrix(m,n);
144
double[][] C = X.getArray();
145
for (int i = 0; i < m; i++) {
146
for (int j = 0; j < n; j++) {
153
/** Add a value to an element and reset the element
154
* @param i the row number of the element
155
* @param j the column number of the element
156
* @param s the double value to be added with
158
public void setPlus(int i, int j, double s) {
162
/** Multiply a value with an element and reset the element
163
* @param i the row number of the element
164
* @param j the column number of the element
165
* @param s the double value to be multiplied with
167
public void setTimes(int i, int j, double s) {
171
/** Set the submatrix A[i0:i1][j0:j1] with a same value
172
* @param i0 the index of the first element of the column
173
* @param i1 the index of the last element of the column
174
* @param j0 the index of the first column
175
* @param j1 the index of the last column
176
* @param s the value to be set to
178
public void setMatrix( int i0, int i1, int j0, int j1, double s ) {
180
for( int i = i0; i <= i1; i++ ) {
181
for( int j = j0; j <= j1; j++ ) {
185
} catch( ArrayIndexOutOfBoundsException e ) {
186
throw new ArrayIndexOutOfBoundsException( "Index out of bounds" );
190
/** Set the submatrix A[i0:i1][j] with the values stored in a
192
* @param i0 the index of the first element of the column
193
* @param i1 the index of the last element of the column
194
* @param j the index of the column
195
* @param v the vector that stores the values*/
196
public void setMatrix( int i0, int i1, int j, DoubleVector v ) {
197
for( int i = i0; i <= i1; i++ ) {
198
A[i][j] = v.get(i-i0);
202
/** Set the whole matrix from a 1-D array
203
* @param v 1-D array of doubles
204
* @param columnFirst Whether to fill the column first or the row.
205
* @throws ArrayIndexOutOfBoundsException Submatrix indices
207
public void setMatrix ( double[] v, boolean columnFirst ) {
209
if( v.length != m * n )
210
throw new IllegalArgumentException("sizes not match.");
213
for( i = 0; i < m; i++ ) {
214
for( j = 0; j < n; j++ ) {
221
for( j = 0; j < n; j++ ) {
222
for( i = 0; i < m; i++ ){
229
} catch( ArrayIndexOutOfBoundsException e ) {
230
throw new ArrayIndexOutOfBoundsException( "Submatrix indices" );
234
/** Returns the maximum absolute value of all elements
235
@return the maximum value
237
public double maxAbs () {
238
double ma = Math.abs(A[0][0]);
239
for (int j = 0; j < n; j++) {
240
for (int i = 0; i < m; i++) {
241
ma = Math.max(ma, Math.abs(A[i][j]));
247
/** Returns the maximum absolute value of some elements of a column,
248
that is, the elements of A[i0:i1][j].
249
@param i0 the index of the first element of the column
250
@param i1 the index of the last element of the column
251
@param j the index of the column
252
@return the maximum value */
253
public double maxAbs ( int i0, int i1, int j ) {
254
double m = Math.abs(A[i0][j]);
255
for (int i = i0+1; i <= i1; i++) {
256
m = Math.max(m, Math.abs(A[i][j]));
261
/** Returns the minimum absolute value of some elements of a column,
262
that is, the elements of A[i0:i1][j].
263
@param i0 the index of the first element of the column
264
@param i1 the index of the last element of the column
265
@param column the index of the column
266
@return the minimum value
268
public double minAbs ( int i0, int i1, int column ) {
269
double m = Math.abs(A[i0][column]);
270
for (int i = i0+1; i <= i1; i++) {
271
m = Math.min(m, Math.abs(A[i][column]));
276
/** Check if the matrix is empty
277
* @return true if the matrix is empty
279
public boolean isEmpty(){
280
if(m == 0 || n == 0) return true;
281
if(A == null) return true;
285
/** Return a DoubleVector that stores a column of the matrix
286
* @param j the index of the column
289
public DoubleVector getColumn( int j ) {
290
DoubleVector v = new DoubleVector( m );
291
double [] a = v.getArray();
292
for(int i = 0; i < m; i++)
297
/** Return a DoubleVector that stores some elements of a column of the
299
* @param i0 the index of the first element of the column
300
* @param i1 the index of the last element of the column
301
* @param j the index of the column
302
* @return the DoubleVector
304
public DoubleVector getColumn( int i0, int i1, int j ) {
305
DoubleVector v = new DoubleVector( i1-i0+1 );
306
double [] a = v.getArray();
308
for( int i = i0; i <= i1; i++ ) {
316
/** Multiplication between a row (or part of a row) of the first matrix
317
* and a column (or part or a column) of the second matrix.
318
* @param i the index of the row in the first matrix
319
* @param j0 the index of the first column in the first matrix
320
* @param j1 the index of the last column in the first matrix
321
* @param B the second matrix
322
* @param l the index of the column in the second matrix
323
* @return the result of the multiplication
325
public double times( int i, int j0, int j1, PaceMatrix B, int l ) {
327
for(int j = j0; j <= j1; j++ ) {
328
s += A[i][j] * B.A[j][l];
333
/** Decimal format for converting a matrix into a string
334
* @return the default decimal format
336
protected DecimalFormat [] format() {
337
return format(0, m-1, 0, n-1, 7, false );
340
/** Decimal format for converting a matrix into a string
341
* @param digits the number of digits
342
* @return the decimal format
344
protected DecimalFormat [] format( int digits ) {
345
return format(0, m-1, 0, n-1, digits, false);
348
/** Decimal format for converting a matrix into a string
349
* @param digits the number of digits
351
* @return the decimal format
353
protected DecimalFormat [] format( int digits, boolean trailing ) {
354
return format(0, m-1, 0, n-1, digits, trailing);
357
/** Decimal format for converting a matrix into a string
361
* @param digits the number of digits
363
* @return the decimal format
365
protected DecimalFormat format(int i0, int i1, int j, int digits,
367
FlexibleDecimalFormat df = new FlexibleDecimalFormat(digits, trailing);
369
for(int i = i0; i <= i1; i ++ )
370
df.update( A[i][j] );
374
/** Decimal format for converting a matrix into a string
380
* @param digits the number of digits
381
* @return the decimal format
383
protected DecimalFormat [] format(int i0, int i1, int j0, int j1,
384
int digits, boolean trailing) {
385
DecimalFormat [] f = new DecimalFormat[j1-j0+1];
386
for( int j = j0; j <= j1; j++ ) {
387
f[j] = format(i0, i1, j, digits, trailing);
393
* Converts matrix to string
395
* @return the matrix as string
397
public String toString() {
398
return toString( 5, false );
402
* Converts matrix to string
404
* @param digits number of digits after decimal point
405
* @param trailing true if trailing zeros are padded
406
* @return the matrix as string
408
public String toString( int digits, boolean trailing ) {
410
if( isEmpty() ) return "null matrix";
412
StringBuffer text = new StringBuffer();
413
DecimalFormat [] nf = format( digits, trailing );
419
int [] nCols = new int[n];
421
for( int j = 0; j < n; j++ ) {
422
lenNumber = nf[j].format( A[0][j]).length();
423
if( count + 1 + lenNumber > width -1 ) {
424
nCols[nk++] = numCols;
428
count += 1 + lenNumber;
434
for( int k = 0; k < n; ) {
435
for( int i = 0; i < m; i++ ) {
436
for( int j = k; j < k + nCols[nk]; j++)
437
text.append( " " + nf[j].format( A[i][j]) );
445
return text.toString();
448
/** Squared sum of a column or row in a matrix
449
* @param j the index of the column or row
450
* @param i0 the index of the first element
451
* @param i1 the index of the last element
452
* @param col if true, sum over a column; otherwise, over a row
453
* @return the squared sum
455
public double sum2( int j, int i0, int i1, boolean col ) {
457
if( col ) { // column
458
for( int i = i0; i <= i1; i++ )
459
s2 += A[i][j] * A[i][j];
462
for( int i = i0; i <= i1; i++ )
463
s2 += A[j][i] * A[j][i];
468
/** Squared sum of columns or rows of a matrix
469
* @param col if true, sum over columns; otherwise, over rows
470
* @return the squared sum
472
public double[] sum2( boolean col ) {
475
double [] s2 = new double[l];
476
for( int i = 0; i < l; i++ )
477
s2[i] = sum2( i, 0, p-1, col );
481
/** Constructs single Householder transformation for a column
483
@param j the index of the column
484
@param k the index of the row
487
public double [] h1( int j, int k ) {
488
double dq[] = new double[2];
489
double s2 = sum2(j, k, m-1, true);
490
dq[0] = A[k][j] >= 0 ? - Math.sqrt( s2 ) : Math.sqrt( s2 );
492
dq[1] = A[k][j] * dq[0];
496
/** Performs single Householder transformation on one column of a matrix
498
@param j the index of the column
499
@param k the index of the row
500
@param q q = - u'u/2; must be negative
501
@param b the matrix to be transformed
502
@param l the column of the matrix b
504
public void h2( int j, int k, double q, PaceMatrix b, int l ) {
506
for( int i = k; i < m; i++ )
507
s += A[i][j] * b.A[i][l];
509
for( int i = k; i < m; i++ )
510
b.A[i][l] += alpha * A[i][j];
513
/** Constructs the Givens rotation
516
* @return a double array that stores the cosine and sine values
518
public double [] g1( double a, double b ) {
519
double cs[] = new double[2];
520
double r = Maths.hypot(a, b);
532
/** Performs the Givens rotation
533
* @param cs a array storing the cosine and sine values
534
* @param i0 the index of the row of the first element
535
* @param i1 the index of the row of the second element
536
* @param j the index of the column
538
public void g2( double cs[], int i0, int i1, int j ){
539
double w = cs[0] * A[i0][j] + cs[1] * A[i1][j];
540
A[i1][j] = - cs[1] * A[i0][j] + cs[0] * A[i1][j];
544
/** Forward ordering of columns in terms of response explanation. On
545
* input, matrices A and b are already QR-transformed. The indices of
546
* transformed columns are stored in the pivoting vector.
548
*@param b the PaceMatrix b
549
*@param pvt the pivoting vector
550
*@param k0 the first k0 columns (in pvt) of A are not to be changed
552
public void forward( PaceMatrix b, IntVector pvt, int k0 ) {
553
for( int j = k0; j < Math.min(pvt.size(), m); j++ ) {
554
steplsqr( b, pvt, j, mostExplainingColumn(b, pvt, j), true );
558
/** Returns the index of the column that has the largest (squared)
559
* response, when each of columns pvt[ks:] is moved to become the
560
* ks-th column. On input, A and b are both QR-transformed.
563
* @param pvt pivoting index of A
564
* @param ks columns pvt[ks:] of A are to be tested
565
* @return the index of the column
567
public int mostExplainingColumn( PaceMatrix b, IntVector pvt, int ks ) {
569
int [] p = pvt.getArray();
570
double ma = columnResponseExplanation( b, pvt, ks, ks );
572
for( int i = ks+1; i < pvt.size(); i++ ) {
573
val = columnResponseExplanation( b, pvt, i, ks );
582
/** Backward ordering of columns in terms of response explanation. On
583
* input, matrices A and b are already QR-transformed. The indices of
584
* transformed columns are stored in the pivoting vector.
586
* A and b must have the same number of rows, being the (pseudo-)rank.
588
* @param b PaceMatrix b
589
* @param pvt pivoting vector
590
* @param ks number of QR-transformed columns; psuedo-rank of A
591
* @param k0 first k0 columns in pvt[] are not to be ordered.
593
public void backward( PaceMatrix b, IntVector pvt, int ks, int k0 ) {
594
for( int j = ks; j > k0; j-- ) {
595
steplsqr( b, pvt, j, leastExplainingColumn(b, pvt, j, k0), false );
599
/** Returns the index of the column that has the smallest (squared)
600
* response, when the column is moved to become the (ks-1)-th
601
* column. On input, A and b are both QR-transformed.
604
* @param pvt pivoting index of A
605
* @param ks psudo-rank of A
606
* @param k0 A[][pvt[0:(k0-1)]] are excluded from the testing.
607
* @return the index of the column
609
public int leastExplainingColumn( PaceMatrix b, IntVector pvt, int ks,
612
int [] p = pvt.getArray();
613
double mi = columnResponseExplanation( b, pvt, ks-1, ks );
615
for( int i = k0; i < ks - 1; i++ ) {
616
val = columnResponseExplanation( b, pvt, i, ks );
625
/** Returns the squared ks-th response value if the j-th column becomes
626
* the ks-th after orthogonal transformation. A[][pvt[ks:j]] (or
627
* A[][pvt[j:ks]], if ks > j) and b[] are already QR-transformed
628
* on input and will remain unchanged on output.
630
* More generally, it returns the inner product of the corresponding
631
* row vector of the response PaceMatrix. (To be implemented.)
633
*@param b PaceMatrix b
634
*@param pvt pivoting vector
635
*@param j the column A[pvt[j]][] is to be moved
636
*@param ks the target column A[pvt[ks]][]
637
*@return the squared response value */
638
public double columnResponseExplanation( PaceMatrix b, IntVector pvt,
642
* If j == ks - 1, returns the squared ks-th response directly.
644
* If j > ks -1, returns the ks-th response after
645
* Householder-transforming the j-th column and the response.
647
* If j < ks - 1, returns the ks-th response after a sequence of
648
* Givens rotations starting from the j-th row. */
651
double [] xxx = new double[n];
652
int [] p = pvt.getArray();
655
if( j == ks -1 ) val = b.A[j][0];
656
else if( j > ks - 1 ) {
657
int jm = Math.min(n-1, j);
658
DoubleVector u = getColumn(ks,jm,p[j]);
659
DoubleVector v = b.getColumn(ks,jm,0);
660
val = v.innerProduct(u) / u.norm2();
663
for( k = j+1; k < ks; k++ ) // make a copy of A[j][]
667
for( k = j+1; k < ks; k++ ) {
668
cs = g1( xxx[k], A[k][p[k]] );
669
for( l = k+1; l < ks; l++ )
670
xxx[l] = - cs[1] * xxx[l] + cs[0] * A[k][p[l]];
671
val = - cs[1] * val + cs[0] * b.A[k][0];
674
return val * val; // or inner product in later implementation???
678
* QR transformation for a least squares problem<br/>
680
* implicitly both A and b are transformed. pvt.size() is the psuedo-rank of
683
* @param b PaceMatrix b
684
* @param pvt pivoting vector
685
* @param k0 the first k0 columns of A (indexed by pvt) are pre-chosen.
686
* (But subject to rank examination.)
688
* For example, the constant term may be reserved, in which
691
public void lsqr( PaceMatrix b, IntVector pvt, int k0 ) {
692
final double TINY = 1e-15;
693
int [] p = pvt.getArray();
694
int ks = 0; // psuedo-rank
695
for(int j = 0; j < k0; j++ ) // k0 pre-chosen columns
696
if( sum2(p[j],ks,m-1,true) > TINY ){ // large diagonal element
697
steplsqr(b, pvt, ks, j, true);
700
else { // collinear column
702
pvt.setSize(pvt.size()-1);
707
// initial QR transformation
708
for(int j = k0; j < Math.min( pvt.size(), m ); j++ ) {
709
if( sum2(p[j], ks, m-1, true) > TINY ) {
710
steplsqr(b, pvt, ks, j, true);
713
else { // collinear column
715
pvt.setSize(pvt.size()-1);
720
b.m = m = ks; // reset number of rows
724
/** QR transformation for a least squares problem <br/>
726
* implicitly both A and b are transformed. pvt.size() is the psuedo-rank of A.
728
* @param b PaceMatrix b
729
* @param pvt pivoting vector
730
* @param k0 the first k0 columns of A (indexed by pvt) are pre-chosen.
731
* (But subject to rank examination.)
733
* For example, the constant term may be reserved, in which
736
public void lsqrSelection( PaceMatrix b, IntVector pvt, int k0 ) {
737
int numObs = m; // number of instances
738
int numXs = pvt.size();
742
if( numXs > 200 || numXs > numObs ) { // too many columns.
745
backward(b, pvt, pvt.size(), k0);
749
* Sets all diagonal elements to be positive (or nonnegative) without
750
* changing the least squares solution
751
* @param Y the response
752
* @param pvt the pivoted column index
754
public void positiveDiagonal( PaceMatrix Y, IntVector pvt ) {
756
int [] p = pvt.getArray();
757
for( int i = 0; i < pvt.size(); i++ ) {
758
if( A[i][p[i]] < 0.0 ) {
759
for( int j = i; j < pvt.size(); j++ )
760
A[i][p[j]] = - A[i][p[j]];
761
Y.A[i][0] = - Y.A[i][0];
766
/** Stepwise least squares QR-decomposition of the problem
768
@param b PaceMatrix b
769
@param pvt pivoting vector
770
@param ks number of transformed columns
771
@param j pvt[j], the column to adjoin or delete
772
@param adjoin to adjoin if true; otherwise, to delete */
773
public void steplsqr( PaceMatrix b, IntVector pvt, int ks, int j,
775
final int kp = pvt.size(); // number of columns under consideration
776
int [] p = pvt.getArray();
778
if( adjoin ) { // adjoining
781
double dq[] = h1( pj, ks );
783
for( int k = ks+1; k < kp; k++ ){
785
h2( pj, ks, dq[1], this, pk);
787
h2( pj, ks, dq[1], b, 0 ); // for matrix. ???
789
for( int k = ks+1; k < m; k++ )
794
for( int i = j; i < ks-1; i++ )
798
for( int i = j; i < ks-1; i++ ){
799
cs = g1( A[i][p[i]], A[i+1][p[i]] );
800
for( int l = i; l < kp; l++ )
801
g2( cs, i, i+1, p[l] );
802
for( int l = 0; l < b.n; l++ )
803
b.g2( cs, i, i+1, l );
808
/** Solves upper-triangular equation <br/>
810
* On output, the solution is stored in b
811
* @param b the response
812
* @param pvt the pivoting vector
813
* @param kp the number of the first columns involved
815
public void rsolve( PaceMatrix b, IntVector pvt, int kp) {
818
int [] p = pvt.getArray();
820
double [][] ba = b.getArray();
821
for( k = 0; k < b.n; k++ ) {
822
ba[kp-1][k] /= A[kp-1][p[kp-1]];
823
for( i = kp - 2; i >= 0; i-- ){
825
for( j = i + 1; j < kp; j++ )
826
s += A[i][p[j]] * ba[j][k];
828
ba[i][k] /= A[i][p[i]];
834
/** Returns a new matrix which binds two matrices together with rows.
835
* @param b the second matrix
836
* @return the combined matrix
838
public PaceMatrix rbind( PaceMatrix b ){
840
throw new IllegalArgumentException("unequal numbers of rows.");
841
PaceMatrix c = new PaceMatrix( m + b.m, n );
842
c.setMatrix( 0, m - 1, 0, n - 1, this );
843
c.setMatrix( m, m + b.m - 1, 0, n - 1, b );
847
/** Returns a new matrix which binds two matrices with columns.
848
* @param b the second matrix
849
* @return the combined matrix
851
public PaceMatrix cbind( PaceMatrix b ) {
853
throw new IllegalArgumentException("unequal numbers of rows: " +
855
PaceMatrix c = new PaceMatrix(m, n + b.n);
856
c.setMatrix( 0, m - 1, 0, n - 1, this );
857
c.setMatrix( 0, m - 1, n, n + b.n - 1, b );
861
/** Solves the nonnegative linear squares problem. That is, <p>
862
* <center> min || A x - b||, subject to x >= 0. </center> <p>
864
* For algorithm, refer to P161, Chapter 23 of C. L. Lawson and
865
* R. J. Hanson (1974). "Solving Least Squares
866
* Problems". Prentice-Hall.
867
* @param b the response
868
* @param pvt vector storing pivoting column indices
869
* @return solution */
870
public DoubleVector nnls( PaceMatrix b, IntVector pvt ) {
871
int j, t, counter = 0, jm = -1, n = pvt.size();
872
double ma, max, alpha, wj;
873
int [] p = pvt.getArray();
874
DoubleVector x = new DoubleVector( n );
875
double [] xA = x.getArray();
876
PaceMatrix z = new PaceMatrix(n, 1);
880
int kp = 0; // #variables in the positive set P
881
while ( true ) { // step 2
882
if( ++counter > 3*n ) // should never happen
883
throw new RuntimeException("Does not converge");
886
bt = new PaceMatrix( b.transpose() );
887
for( j = kp; j <= n-1; j++ ) { // W = A' (b - A x)
888
wj = bt.times( 0, kp, m-1, this, p[j] );
889
if( wj > max ) { // step 4
896
if ( t == -1) break; // optimum achieved
899
pvt.swap( kp, t ); // move variable from set Z to set P
902
steplsqr( b, pvt, kp-1, kp-1, true );
906
for( j = 0; j <= kp-1; j++ ) z.A[j][0] = b.A[j][0];
909
for( j = 0; j <= kp-1; j++ ) { // step 7, 8 and 9
910
if( z.A[j][0] <= 0.0 ) { // alpha always between 0 and 1
911
alpha = xA[j] / ( xA[j] - z.A[j][0] );
918
for( j = 0; j <= kp-1; j++ ) xA[j] = z.A[j][0]; // step 7
920
for( j = kp-1; j >= 0; j-- ) { // step 10
921
// Modified to avoid round-off error (which seemingly
922
// can cause infinite loop).
923
if( j == jm ) { // step 11
925
steplsqr( b, pvt, kp, j, false );
926
kp--; // move variable from set P to set Z
928
else xA[j] += ma * ( z.A[j][0] - xA[j] );
938
/** Solves the nonnegative least squares problem with equality
939
* constraint. That is, <p>
940
* <center> min ||A x - b||, subject to x >= 0 and c x = d. </center> <p>
942
* @param b the response
943
* @param c coeficients of equality constraints
944
* @param d constants of equality constraints
945
* @param pvt vector storing pivoting column indices
946
* @return the solution
948
public DoubleVector nnlse( PaceMatrix b, PaceMatrix c, PaceMatrix d,
950
double eps = 1e-10 * Math.max( c.maxAbs(), d.maxAbs() ) /
951
Math.max( maxAbs(), b.maxAbs() );
953
PaceMatrix e = c.rbind( new PaceMatrix( times(eps) ) );
954
PaceMatrix f = d.rbind( new PaceMatrix( b.times(eps) ) );
956
return e.nnls( f, pvt );
959
/** Solves the nonnegative least squares problem with equality
960
* constraint. That is, <p>
961
* <center> min ||A x - b||, subject to x >= 0 and || x || = 1. </center>
963
* @param b the response
964
* @param pvt vector storing pivoting column indices
965
* @return the solution
967
public DoubleVector nnlse1( PaceMatrix b, IntVector pvt ) {
968
PaceMatrix c = new PaceMatrix( 1, n, 1 );
969
PaceMatrix d = new PaceMatrix( 1, b.n, 1 );
971
return nnlse(b, c, d, pvt);
974
/** Generate matrix with standard-normally distributed random elements
975
@param m Number of rows.
976
@param n Number of colums.
977
@return An m-by-n matrix with random elements. */
978
public static Matrix randomNormal( int m, int n ) {
979
Random random = new Random();
981
Matrix A = new Matrix(m,n);
982
double[][] X = A.getArray();
983
for (int i = 0; i < m; i++) {
984
for (int j = 0; j < n; j++) {
985
X[i][j] = random.nextGaussian();
994
* @param args the commandline arguments - ignored
996
public static void main( String args[] ) {
997
System.out.println("================================================" +
999
System.out.println("To test the pace estimators of linear model\n" +
1002
double sd = 2; // standard deviation of the random error term
1003
int n = 200; // total number of observations
1004
double beta0 = 100; // intercept
1005
int k1 = 20; // number of coefficients of the first cluster
1006
double beta1 = 0; // coefficient value of the first cluster
1007
int k2 = 20; // number of coefficients of the second cluster
1008
double beta2 = 5; // coefficient value of the second cluster
1009
int k = 1 + k1 + k2;
1011
DoubleVector beta = new DoubleVector( 1 + k1 + k2 );
1012
beta.set( 0, beta0 );
1013
beta.set( 1, k1, beta1 );
1014
beta.set( k1+1, k1+k2, beta2 );
1016
System.out.println("The data set contains " + n +
1017
" observations plus " + (k1 + k2) +
1018
" variables.\n\nThe coefficients of the true model"
1019
+ " are:\n\n" + beta );
1021
System.out.println("\nThe standard deviation of the error term is " +
1024
System.out.println("==============================================="
1027
PaceMatrix X = new PaceMatrix( n, k1+k2+1 );
1028
X.setMatrix( 0, n-1, 0, 0, 1 );
1029
X.setMatrix( 0, n-1, 1, k1+k2, random(n, k1+k2) );
1032
PaceMatrix( X.times( new PaceMatrix(beta) ).
1033
plusEquals( randomNormal(n,1).times(sd) ) );
1035
IntVector pvt = (IntVector) IntVector.seq(0, k1+k2);
1037
/*System.out.println( "The OLS estimate (by jama.Matrix.solve()) is:\n\n" +
1038
(new PaceMatrix(X.solve(Y))).getColumn(0) );*/
1040
X.lsqrSelection( Y, pvt, 1 );
1041
X.positiveDiagonal( Y, pvt );
1043
PaceMatrix sol = (PaceMatrix) Y.clone();
1044
X.rsolve( sol, pvt, pvt.size() );
1045
DoubleVector betaHat = sol.getColumn(0).unpivoting( pvt, k );
1046
System.out.println( "\nThe OLS estimate (through lsqr()) is: \n\n" +
1049
System.out.println( "\nQuadratic loss of the OLS estimate (||X b - X bHat||^2) = " +
1050
( new PaceMatrix( X.times( new
1051
PaceMatrix(beta.minus(betaHat)) )))
1052
.getColumn(0).sum2() );
1054
System.out.println("=============================================" +
1056
System.out.println(" *** Pace estimation *** \n");
1057
DoubleVector r = Y.getColumn( pvt.size(), n-1, 0);
1058
double sde = Math.sqrt(r.sum2() / r.size());
1060
System.out.println( "Estimated standard deviation = " + sde );
1062
DoubleVector aHat = Y.getColumn( 0, pvt.size()-1, 0).times( 1./sde );
1063
System.out.println("\naHat = \n" + aHat );
1065
System.out.println("\n========= Based on chi-square mixture ============");
1067
ChisqMixture d2 = new ChisqMixture();
1068
int method = MixtureDistribution.NNMMethod;
1069
DoubleVector AHat = aHat.square();
1070
d2.fit( AHat, method );
1071
System.out.println( "\nEstimated mixing distribution is:\n" + d2 );
1073
DoubleVector ATilde = d2.pace2( AHat );
1074
DoubleVector aTilde = ATilde.sqrt().times(aHat.sign());
1075
PaceMatrix YTilde = new
1076
PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
1077
X.rsolve( YTilde, pvt, pvt.size() );
1078
DoubleVector betaTilde =
1079
YTilde.getColumn(0).unpivoting( pvt, k );
1080
System.out.println( "\nThe pace2 estimate of coefficients = \n" +
1082
System.out.println( "Quadratic loss = " +
1083
( new PaceMatrix( X.times( new
1084
PaceMatrix(beta.minus(betaTilde)) )))
1085
.getColumn(0).sum2() );
1087
ATilde = d2.pace4( AHat );
1088
aTilde = ATilde.sqrt().times(aHat.sign());
1089
YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
1090
X.rsolve( YTilde, pvt, pvt.size() );
1091
betaTilde = YTilde.getColumn(0).unpivoting( pvt, k );
1092
System.out.println( "\nThe pace4 estimate of coefficients = \n" +
1094
System.out.println( "Quadratic loss = " +
1095
( new PaceMatrix( X.times( new
1096
PaceMatrix(beta.minus(betaTilde)) )))
1097
.getColumn(0).sum2() );
1099
ATilde = d2.pace6( AHat );
1100
aTilde = ATilde.sqrt().times(aHat.sign());
1101
YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
1102
X.rsolve( YTilde, pvt, pvt.size() );
1103
betaTilde = YTilde.getColumn(0).unpivoting( pvt, k );
1104
System.out.println( "\nThe pace6 estimate of coefficients = \n" +
1106
System.out.println( "Quadratic loss = " +
1107
( new PaceMatrix( X.times( new
1108
PaceMatrix(beta.minus(betaTilde)) )))
1109
.getColumn(0).sum2() );
1111
System.out.println("\n========= Based on normal mixture ============");
1113
NormalMixture d = new NormalMixture();
1114
d.fit( aHat, method );
1115
System.out.println( "\nEstimated mixing distribution is:\n" + d );
1117
aTilde = d.nestedEstimate( aHat );
1118
YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
1119
X.rsolve( YTilde, pvt, pvt.size() );
1120
betaTilde = YTilde.getColumn(0).unpivoting( pvt, k );
1121
System.out.println( "The nested estimate of coefficients = \n" +
1123
System.out.println( "Quadratic loss = " +
1124
( new PaceMatrix( X.times( new
1125
PaceMatrix(beta.minus(betaTilde)) )))
1126
.getColumn(0).sum2() );
1129
aTilde = d.subsetEstimate( aHat );
1130
YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
1131
X.rsolve( YTilde, pvt, pvt.size() );
1133
YTilde.getColumn(0).unpivoting( pvt, k );
1134
System.out.println( "\nThe subset estimate of coefficients = \n" +
1136
System.out.println( "Quadratic loss = " +
1137
( new PaceMatrix( X.times( new
1138
PaceMatrix(beta.minus(betaTilde)) )))
1139
.getColumn(0).sum2() );
1141
aTilde = d.empiricalBayesEstimate( aHat );
1142
YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times( sde ));
1143
X.rsolve( YTilde, pvt, pvt.size() );
1144
betaTilde = YTilde.getColumn(0).unpivoting( pvt, k );
1145
System.out.println( "\nThe empirical Bayes estimate of coefficients = \n"+
1148
System.out.println( "Quadratic loss = " +
1149
( new PaceMatrix( X.times( new
1150
PaceMatrix(beta.minus(betaTilde)) )))
1151
.getColumn(0).sum2() );
added
removed
weka/classifiers/functions/pace/PaceMatrix.java
