home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 21 / IOPROG_21.ISO / SOFT / LIBMAT.ZIP / MATOLS.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-16  |  11.6 KB  |  473 lines
  1. /**************************************************/
  2. /*       matols.c source for matols class         */
  3. /**************************************************/
  4.  
  5.  
  6. /**************************************************/
  7. /*            MatClass Source File                */
  8. /*       Copyright of C. R. Birchenhall           */
  9. /*       University of Manchester, UK.            */
  10. /*   MatClass is freeware. This file should be    */
  11. /* made freely available to users of any software */
  12. /* whose creation is wholly or partly dependent   */
  13. /*                on this file.                   */
  14. /**************************************************/
  15.  
  16. #include "matols.hpp"
  17. #include <math.h>
  18.  
  19. /*************************************************************/
  20. /*                     matOls methods                      */
  21. /*************************************************************/
  22.  
  23. void matOls::setNames( void )
  24. {
  25.    Y.name( "olsY" ) ;
  26.    X.name( "olsX" ) ;
  27.    R.name( "olsR" ) ;
  28.    Rinv.name( "olsR-1" ) ;
  29.    beta.name( "olsBta" ) ;
  30.    resid.name( "olsRes" ) ;
  31.    V.name( "olsV" ) ;
  32.    VSqrt.name( "olsVsq" ) ;
  33.    TSS.name( "olsTSS" );
  34.    RSS.name( "olsRSS" ) ;
  35.    YMean.name( "olsYmn" ) ;
  36.    SE.name( "olsSE" ) ;
  37. } // matOls setNames
  38.  
  39. #ifdef __cplusplus
  40. matOls::matOls( void ) : matObject()
  41. #else
  42. matOls::matOls( void ) : ()
  43. #endif
  44. {
  45.    status = 0 ;
  46.    setNames() ;
  47. } // matOls
  48.  
  49. #ifdef __cplusplus
  50. matOls::matOls( const matrix& y, const matrix& x ) : matObject()
  51. #else
  52. matOls::matOls( const matrix& y, const matrix& x ) : ()
  53. #endif
  54. {
  55.    status = 0 ;
  56.    assign( y, x ) ;
  57.    setNames() ;
  58. } // matOls( matrix& )
  59.  
  60. matOls::~matOls( void ) {} // ~matOls
  61.  
  62. #ifdef __cplusplus
  63. matOls::matOls( matOls& ols ) : matObject( ols )
  64. #else
  65. matOls::matOls( matOls& ols ) : ( ols )
  66. #endif
  67. {
  68.    errorExit( "matOls copy", NIMPL ) ;
  69. } // matOls copy constr
  70.  
  71. void matOls::operator = ( matOls& ols )
  72. {
  73.    nm = ols.nm ;
  74.    errorExit( "matOls =", NIMPL ) ;
  75. } // matOls =
  76.  
  77. void matOls::decompose( void )
  78. {
  79.    errorExit( "matOls::decomp", NIMPL ) ;
  80. } // matOls decompose
  81.  
  82. outFile& matOls::info( outFile& f ) M_CONST
  83. {
  84.    errorExit( "matOls::info", NIMPL ) ;
  85.    return f ;
  86. } // matOls::info
  87.  
  88. REAL matOls::varAdd( const matrix& z, INDEX n )
  89. {
  90.    z.errorExit( "matOls::varAdd", NIMPL ) ;
  91.    return (REAL) n ;
  92. } // matOls::varAdd
  93.  
  94. void matOls::initial( void )
  95. {
  96.    static char *mName = "initial" ;
  97.    matrix XMax, XMin ;
  98.    if ( X.nRows() != Y.nRows() )
  99.       errorExit( mName, NEDIM ) ;
  100.    INDEX j ;
  101.    nObs = X.nRows() ;
  102.    nVars = X.nCols() ;
  103.    nDep = Y.nCols() ;
  104.    dof = nObs - nVars ;
  105.    tol = matrixTol() ;
  106.    // constant true iff X(j) is constant for some j
  107.    XMax = X.colMax() ;
  108.    XMin = X.colMin() ;
  109.    for ( constant = 0, j = 1 ; !constant && ( j <= nVars ) ; j++ )
  110.       constant = ( XMax(j) == XMin(j) ) ;
  111.    YMean.colMeanOf( Y ) ;
  112.    TSS.colSqDevOf( Y, YMean ) ;
  113.    RSS.reset( nDep ) ;
  114.    beta.reset( nVars, nDep ) ;
  115.    resid.reset( nObs, nDep ) ;
  116.    R.reset( nVars, nVars ) ;
  117.    Rinv.reset( nVars, nVars ) ;
  118.    V.reset( nVars, nVars ) ;
  119.    VSqrt.reset( nVars ) ;
  120.    SE.reset( nDep ) ;
  121.    status = ASSIGNED ;
  122.    return ;
  123. } // matOls::initial
  124.  
  125. void matOls::assign( const matrix& y, const matrix& x )
  126. {
  127.    static char *mName = "assign" ;
  128.    matFunc func( mName ) ; // debugInfo( func ) ;
  129.    Y = y ;
  130.    X = x ;
  131.    initial() ;
  132. } // matOls::assign
  133.  
  134. void matOls::capture( matrix& y, matrix& x )
  135. {
  136.    static char *mName = "capture" ;
  137.    matFunc func( mName ) ; // debugInfo( func ) ;
  138.    Y.capture(y) ;
  139.    X.capture(x) ;
  140.    initial() ;
  141. } // matOls::capture
  142.  
  143. INDEX matOls::ok( void )
  144. {
  145.    decompose() ;
  146.    return !(status & SINGULAR) ;
  147. } // matOls::ok
  148.  
  149. INDEX matOls::ok( char *mName )
  150. {
  151.    decompose() ;
  152.    if ( status & SINGULAR )
  153.       errorExit( mName, SINGM ) ;
  154.    return OK ;
  155. } // matOls::ok
  156.  
  157. REAL matOls::setTol( REAL newTol )
  158. {
  159.    REAL oldTol = tol ;
  160.    if ( newTol >= 0.0 )
  161.       tol = newTol ;
  162.    return oldTol ;
  163. } // matOls::setTol
  164.  
  165. void matOls::clear( void )
  166. {
  167.    static char *mName = "clear" ;
  168.    matFunc func( mName ) ; debugInfo( func ) ;
  169.    if ( status != 0 ) {
  170.       X.clear() ;
  171.       Y.clear() ;
  172.       R.clear() ;
  173.       Rinv.clear() ;
  174.       beta.clear() ;
  175.       YMean.clear() ;
  176.       RSS.clear() ;
  177.       TSS.clear() ;
  178.       V.clear() ;
  179.       VSqrt.clear() ;
  180.       SE.clear() ;
  181.       status = 0 ;
  182.       constant = 0 ;
  183.    } // if
  184. } // clear
  185.  
  186. charArray& matOls::name( char *newName )
  187. {
  188.    static char *mName = "name" ;
  189.    matFunc func( mName ) ; debugInfo( func ) ;
  190.    if ( newName != 0 )
  191.       nm = newName ;
  192.    return nm ;
  193. } // matOls::name
  194.  
  195. char* matOls::nameStr( void ) M_CONST
  196. {
  197.     return nm.array() ;
  198. } // matOls::name
  199.  
  200. matrix& matOls::coeff( matrix& b )
  201. {
  202.    static char *mName = "coeff" ;
  203.    matFunc func( mName ) ; debugInfo( func ) ;
  204.    if ( ok( mName ) )
  205.       b = beta ;
  206.    return b ;
  207. } // coeff
  208.  
  209. outFile& matOls::olsInfo( outFile& f, char *cName ) M_CONST
  210. {
  211.    if ( matListCtrl() > 4 )
  212.       return f ;
  213.    objectInfo( f ) ;
  214.    putName( nameStr(), f ) ;
  215.    putName( cName, f ) ;
  216.    putField( status, f ) ;
  217.    if ( status & ASSIGNED ) {
  218.      putField( Y.identity(), f ) ;
  219.      putField( X.identity(), f ) ;
  220.    } else
  221.      putField( "NA", f ) ;
  222.    f.newLine() ;
  223.    return f ;
  224. } // matOls::olsInfo
  225.  
  226. outFile& matOls::put( outFile& f ) M_CONST
  227. {
  228.    static char *mName = "put" ;
  229.    matFunc func( mName ) ; debugInfo( func ) ;
  230.    errorExit( mName, NIMPL ) ;
  231.    return f ;
  232. } // matOls::put
  233.  
  234. inFile& matOls::get( inFile& f )
  235. {
  236.    static char *mName = "get" ;
  237.    matFunc func( mName ) ; debugInfo( func ) ;
  238.    errorExit( mName, NIMPL ) ;
  239.    return f ;
  240. } // matOls::get
  241.  
  242. void matOls::formResid( void )
  243. {
  244.    static char *mName = "formRes" ;
  245.    matFunc func( mName ) ; debugInfo( func ) ;
  246.    ok( mName ) ;
  247.    if ( !( status & RESIDUALS ) ) {
  248.       INDEX i ;
  249.       resid = Y ;
  250.       resid -= X * beta ;
  251.       RSS.colSqrOf( resid ) ;
  252.       SE.reset( nDep ) ;
  253.       REAL r ;
  254.       for ( i = 1 ; i <= nDep ; i++ ) {
  255.          r = RSS(i) ;
  256.          SE(i) = sqrt( r / dof ) ;
  257.       } // for
  258.       status |= RESIDUALS ;
  259.    } // if
  260. } // matOls::formResid
  261.  
  262. matrix& matOls::residuals( matrix& res )
  263. {
  264.    static char *mName = "resid" ;
  265.    matFunc func( mName ) ; debugInfo( func ) ;
  266.    formResid() ;
  267.    res = resid ;
  268.    return res ;
  269. } // matOls::residuals
  270.  
  271. matrix& matOls::fitted( matrix& fit )
  272. {
  273.    static char *mName = "fitted" ;
  274.    matFunc func( mName ) ; debugInfo( func ) ;
  275.    formResid() ;
  276.    fit =  Y ;
  277.    fit -= resid ;
  278.    return fit ;
  279. } // matOls::fitted
  280.  
  281. REAL matOls::se( INDEX i )
  282. {
  283.    static char *mName = "se" ;
  284.    matFunc func( mName ) ; debugInfo( func ) ;
  285.    formResid() ;
  286.    return SE(i) ;
  287. } // matOls::se
  288.  
  289. REAL matOls::rss( INDEX i )
  290. {
  291.    static char *mName = "rss" ;
  292.    matFunc func( mName ) ; debugInfo( func ) ;
  293.    formResid() ;
  294.    return RSS(i) ;
  295. } // matOls::rss
  296.  
  297. REAL matOls::rsq( INDEX i )
  298. {
  299.    static char *mName = "rsq" ;
  300.    matFunc func( mName ) ; debugInfo( func ) ;
  301.    formResid() ;
  302.    REAL r = RSS(i), t = TSS(i), y = YMean(i) ;
  303.    if ( constant )
  304.       return 1.0 - r / t ;
  305.    else
  306.       return 1.0 - r / ( t + y * y ) ;
  307. } // matOls::rsq
  308.  
  309. REAL matOls::rBarSq( INDEX i )
  310. {
  311.    static char *mName = "rBarSq" ;
  312.    matFunc func( mName ) ; debugInfo( func ) ;
  313.    formResid() ;
  314.    REAL r = RSS(i), t = TSS(i), y = YMean(i) ;
  315.    if ( constant )
  316.       return 1.0 - ( r * ( nObs - 1 ) ) / ( t * dof ) ;
  317.    else
  318.       return 1.0 - ( r * ( nObs - 1 ) )
  319.                    / ( ( t + y * y ) * dof ) ;
  320. } // matOls::rBarSq
  321.  
  322. REAL matOls::tss( INDEX i )
  323. {
  324.    static char *mName = "tss" ;
  325.    matFunc func( mName ) ; debugInfo( func ) ;
  326.    ok( mName ) ;
  327.    return TSS(i) ;
  328. } // matOls::rss
  329.  
  330. matrix& matOls::stdErr( matrix& std )
  331. {
  332.    std.reset( nVars, nDep ) ;
  333.    INDEX i, j ;
  334.    REAL s ;
  335.    formV() ;
  336.    for ( j = 1 ; j <= nDep ; j++ ) {
  337.       s = SE(j) ;
  338.       for ( i = 1 ; i <= nVars ; i++ )
  339.          std(i,j) = s * VSqrt(i) ;
  340.    } // for
  341.    return std ;
  342. } // matOls::stdErr
  343.  
  344. matrix& matOls::cov( matrix& cv, INDEX i )
  345. {
  346.    static char *mName = "cov" ;
  347.    matFunc func( mName ) ; debugInfo( func ) ;
  348.    formV() ;
  349.    REAL s = SE(i) ;
  350.    cv = V ;
  351.    cv *= s * s ;
  352.    return cv ;
  353. } // matOls::cov
  354.  
  355. REAL matOls::dw( INDEX j )
  356. /*************************
  357.   Durbin-Watson Test
  358. *************************/
  359. {
  360.    static char *mName = "dw" ;
  361.    matFunc func( mName ) ; debugInfo( func ) ;
  362.    DOUBLE e, num = 0.0, r = resid(1,j), denom ;
  363.    INDEX i ;
  364.    formResid() ;
  365.    denom = r * r ;
  366.    for ( i = 2 ; i <= nObs ; i++ ) {
  367.       e = resid(i,j) ;
  368.       denom += e * e ;
  369.       e -= resid(i-1,j) ;
  370.       num += e * e ;
  371.    } // for i
  372.    return num  / denom ;
  373. } // matOls::dw
  374.  
  375. REAL matOls::tTest( const matrix& w, REAL r, INDEX n )
  376. /************************************************
  377.   Return classical t-Value for hypothesis
  378.                 w'beta = r
  379.   where beta is the set of coeff for nth. eqn.
  380. ************************************************/
  381. {
  382.    static char *mName = "tTest" ;
  383.    matFunc func( mName ) ; debugInfo( func ) ;
  384.    DOUBLE num , denom ;
  385.    matrix bCol, v ;
  386.    if ( w.nRows() != nVars )
  387.       errorExit( "ols tTest", NEDIM ) ;
  388.    formV() ;
  389.    // form w'beta - r in num
  390.    bCol = beta.col( n ) ;
  391.    num = w.inner(bCol) - r ;
  392.    bCol.clear() ;
  393.    // form sqrt( w'Vw ) in denom
  394.    v = V * w ;
  395.    denom = sqrt( w.inner(v) ) ;
  396.    return num / ( SE(n) * denom ) ;
  397. } // matOls::tTest
  398.  
  399. REAL matOls::FTest( const matrix& H, const matrix& a, INDEX n )
  400. /***************************************************
  401.    Return classical F-value for the hypotheses
  402.                  H * beta = a
  403.    where beta is set of coeffs for nth eqn.
  404. ***************************************************/
  405. {
  406.    static char *mName = "FTest" ;
  407.    matFunc func( mName ) ; debugInfo( func ) ;
  408.    formV() ;
  409.    INDEX nr = H.nRows() ;
  410.    matrix w( nr ), v( nr ), Z(nr,nr) ;
  411.    // w = H * beta - a
  412.    w = H * beta.col(n) ;
  413.    w -= a ;
  414.    Z = H * ( V.multT( H ) ) ;
  415.    v = w / Z ;
  416.    REAL s = SE(n) ;
  417.    return ( w.inner( v ) / ( H.nRows() * s * s ) ) ;
  418. } // matOls::FTest
  419.  
  420. void matOls::formV( void )
  421. {
  422.    static char *mName = "formV" ;
  423.    matFunc func( mName ) ; debugInfo( func ) ;
  424.    formResid() ;
  425.    if ( !(status & VMATRIX ) ) {
  426.       INDEX i, j, k, n = Rinv.nCols() ;
  427.       DOUBLE sum, r ;
  428.       V.reset(n,n) ;
  429.       VSqrt.reset(n) ;
  430.       for ( i = 1 ; i <= n ; i++ ) {
  431.          for ( j = i ; j <= n ; j++ ) {
  432.             sum = 0.0 ;
  433.             for ( k = j ; k <= n ; k++ ) {
  434.                r = Rinv(i,k) ;
  435.                sum += r * Rinv(j,k) ;
  436.             } // for k
  437.             V(i,j) = (REAL) sum ;
  438.             V(j,i) = sum ;
  439.          } // for j
  440.          r = V(i,i) ;
  441.          VSqrt(i) = sqrt(r) ;
  442.       } // for i
  443.       status |= VMATRIX ;
  444.    } // if
  445. } // matOls::formV
  446.  
  447. REAL matOls::cond( void )
  448. /***************************************
  449.   Estimated condition of X matrix
  450.   based on Frobenius norm.
  451. ***************************************/
  452. {
  453.    static char *mName = "cond" ;
  454.    matFunc func( mName ) ; debugInfo( func ) ;
  455.    ok( mName ) ;
  456.    if ( !(status & CONDITIONED ) ) {
  457.       DOUBLE norm1 = 0.0, norm2 = 0.0 ;
  458.       INDEX i, j ;
  459.       REAL r , ri ;
  460.       for ( i = 1 ; i <= nVars ; i++ ) {
  461.          for ( j = i ; j <= nVars ; j++ ) {
  462.             r = R(i,j) ;
  463.             ri = Rinv(i,j) ;
  464.             norm1 += r * r ;
  465.             norm2 += ri * ri ;
  466.          } // for j
  467.       } // for i
  468.       condition = (REAL) ( sqrt( norm1 ) * sqrt( norm2 ) ) ;
  469.       status |= CONDITIONED ;
  470.    } // if
  471.    return condition ;
  472. } // matOls::cond
  473.