home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 21 / IOPROG_21.ISO / SOFT / LIBMAT.ZIP / MATDIST.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-16  |  11.9 KB  |  374 lines
  1. /**************************************************/
  2. /*  matdist.c source for distribution classes     */
  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 "matspec.hpp"
  17.  
  18. /****************************************************/
  19. /*        Poisson Distribtion Functions             */
  20. /****************************************************/
  21.  
  22. REAL poissonDist( REAL k, REAL mean )
  23. {
  24.    static char* mName = "poissonDist(r)" ;
  25.    matFunc func( mName ) ;
  26.    return 1.0 - incGamma( k, mean ) ;
  27. } // poissonDist
  28.  
  29.  
  30. #ifdef __cplusplus
  31. matrix poissonDist( const matrix& x, REAL mean )
  32. #else
  33. matrix poissonDistM( const matrix& x, REAL mean )
  34. #endif
  35. {
  36.    static char* mName = "poissonDist(m)" ;
  37.    matFunc func( mName ) ; x.debugInfo( func ) ;
  38.    matrix z ;
  39.    z = incGamma( x, mean ) ;
  40.    z.linear( -1.0, z, 1.0 ) ;
  41.    return z ;
  42. } // PoissonDist( const matrix&, REAL )
  43.  
  44.  
  45. /****************************************************/
  46. /*      Chi-squared Distribtion Functions           */
  47. /****************************************************/
  48.  
  49. REAL chi2Dist( REAL r, INDEX dof )
  50. {
  51.    static char* mName = "chi2Dist(r)" ;
  52.    matFunc func( mName ) ;
  53.    return incGamma( 0.5 * r , 0.5 * REAL(dof)  ) ;
  54. } // chi2Dist
  55.  
  56.  
  57. #ifdef __cplusplus
  58. matrix chi2Dist( const matrix& x, INDEX dof )
  59. #else
  60. matrix chi2DistM( const matrix& x, INDEX dof )
  61. #endif
  62. {
  63.    static char* mName = "chi2Dist(m)" ;
  64.    matFunc func( mName ) ; x.debugInfo( func ) ;
  65.    matrix z ;
  66.    z = x ;
  67.    z *= 0.5 ;
  68.    incGamma.par( 0.5 * REAL(dof) ) ;
  69.    return incGamma.transform( z, z ) ;
  70. } // chi2Dist
  71.  
  72.  
  73. /****************************************************/
  74. /*         Student Distribtion Functions            */
  75. /****************************************************/
  76.  
  77. REAL studentDist( REAL r, INDEX dof )
  78. {
  79.    static char* mName = "student(r)" ;
  80.    matFunc func( mName ) ;
  81.    REAL d = REAL(dof), A ;
  82.    if ( r == 0.0 )
  83.       return 0.5 ;
  84.    A = 0.5 * incBeta( d / ( d + r * r ), 0.5 * d, 0.5 ) ;
  85.    if ( r > 0.0 )
  86.       return 1.0 - A ;
  87.    return A ;
  88. } // studentDist REAL
  89.  
  90. #ifdef __cplusplus
  91. matrix studentDist( const matrix& x, INDEX dof )
  92. #else
  93. matrix studentDistM( const matrix& x, INDEX dof )
  94. #endif
  95. {
  96.    static char* mName = "student(m)" ;
  97.    matFunc func( mName ) ; x.debugInfo( func ) ;
  98.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  99.    matrix z( nr, nc ) ;
  100.    for ( j = 1 ; j <= nc ; j++ )
  101.       for ( i = 1 ; i <= nr ; i++ )
  102.          z(i,j) = studentDist( x(i,j), dof ) ;
  103.    return z ;
  104. } // studentDist matrix&
  105.  
  106.  
  107. /****************************************************/
  108. /*             F-Distribtion Functions              */
  109. /****************************************************/
  110.  
  111. REAL FDist( REAL r, INDEX dof1, INDEX dof2 )
  112. {
  113.    static char* mName = "FDist" ;
  114.    matFunc func( mName ) ;
  115.    REAL d1 = REAL(dof1), d2 = REAL(dof2) ;
  116.    if ( r < 0.0 )
  117.       matErrNumExit( mName, BDARG ) ;
  118.    return 1.0 - incBeta( d2 /( d2 + d1 * r ), 0.5 * d2, 0.5 * d1 ) ;
  119. } // FDist
  120.  
  121. #ifdef __cplusplus
  122. matrix FDist( const matrix& x, INDEX dof1, INDEX dof2 )
  123. #else
  124. matrix FDistM( const matrix& x, INDEX dof1, INDEX dof2 )
  125. #endif
  126. {
  127.    static char* mName = "FDist" ;
  128.    matFunc func( mName ) ; x.debugInfo( func ) ;
  129.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  130.    matrix z( nr, nc ) ;
  131.    REAL r ;
  132.    for ( j = 1 ; j <= nc ; j++ )
  133.       for ( i = 1 ; i <= nr ; i++ ) {
  134.          r = x(i,j) ;
  135.          z(i,j) = FDist( r, dof1, dof2 ) ;
  136.       } // for
  137.    return z ;
  138. } // FDist
  139.  
  140. /************************************************************/
  141. /*                      normal class                    */
  142. /************************************************************/
  143.  
  144.  
  145. #ifdef __cplusplus
  146. normalFunc::normalFunc( void ) : matSpecFunc()
  147. #else
  148. normalFunc::normalFunc( void ) : ()
  149. #endif
  150. {
  151.    // !!! these values need more thought !!!
  152.    lower =  8 ;
  153.    upper = 18.66 ;
  154.    calcUpper = 0 ;
  155. } // normalFunc
  156.  
  157. #ifdef __cplusplus
  158. normalFunc::normalFunc( normalFunc& fn )
  159.               : matSpecFunc(fn)
  160. #else
  161. normalFunc::normalFunc( normalFunc& fn )
  162.               : ( fn )
  163. #endif
  164. {
  165.    lower = fn.lower ;
  166.    upper = fn.upper ;
  167. } // normalFunc
  168.  
  169. void normalFunc::operator = ( normalFunc& fn )
  170. {
  171.    matSpecFunc::operator = ( fn ) ;
  172.    lower = fn.lower ;
  173.    upper = fn.upper ;
  174. } // normalFunc
  175.  
  176. normalFunc::~normalFunc( void ) {}
  177.  
  178. outFile& normalFunc::info( outFile& f ) M_CONST
  179. {
  180.    return specInfo( f, "incGam" ) ;
  181. } // info
  182.  
  183. REAL normalFunc::value( REAL r )
  184. /*****************************************************
  185.    See Hill pp.126-129 in Applied Stats Algorithms
  186. ******************************************************/
  187. {
  188.    static REAL a[7]  = { 0.398942280444,  0.399903438504,
  189.                            5.75885480458,  29.8213557808,
  190.                            2.62433121679,  48.6959930692,
  191.                            5.92885724438 } ;
  192.    static REAL b[12] = {   0.398942280385,  3.8052E-8,
  193.                            1.00000615302,   3.98064794E-4,
  194.                            1.98615381364,   0.151679116635,
  195.                            5.29330324926,   4.8385912808,
  196.                           15.1508972451,    0.742380924027,
  197.                           30.789933034,     3.990194417011 } ;
  198.    static char *mName = "value" ;
  199.    matFunc func( mName ) ; debugInfo( func ) ;
  200.    REAL y, z, term ;
  201.    int isUpper = calcUpper ;
  202.    error = NOERROR ;
  203.    if ( r < zero ) {
  204.       isUpper = !isUpper ;
  205.       z = - r ;
  206.    } else
  207.       z = r ;
  208.    if ( z > lower || ( isUpper && z > upper ) )
  209.       return zero ;
  210.    y = half * z * z ;
  211.    if ( z <= 1.28 ) {
  212.       term = half - z * ( a[0] - a[1] * y / ( y + a[2] - a[3] /
  213.                          ( y + a[4] + a[5] / ( y + a[6] ) ) ) ) ;
  214.    } else {
  215.       term = b[0] * exp( -y ) / ( z - b[1] + b[2] / ( z + b[3] + b[4] /
  216.                     ( z - b[5] + b[6] / ( z + b[7] -b[8] /
  217.                     ( z + b[9] + b[10] / ( z + b[11] ) ) ) ) ) ) ;
  218.    } // else
  219.    if ( isUpper )
  220.       return term ;
  221.    // else
  222.    return one - term ;
  223. } // normalFunc value
  224.  
  225. REAL normalFunc::inv( REAL p )
  226. /************************************************
  227.    Adapted from Wichura's PPND16, Algorithm AS241
  228.    Applied Statistics Vol 37 1988 pp 477 - 484
  229. *************************************************/
  230. {
  231.  
  232.    const REAL SPLIT1 = 0.425,
  233.               SPLIT2 = 5.0,
  234.               CONST1 = 0.180625,
  235.               CONST2 = 1.6;
  236.  
  237.    static const REAL a[8] = {
  238.                            3.3871328727963666080E0,
  239.                            1.3314166789178437745E2,
  240.                            1.9715909503065514427E3,
  241.                            1.3731693765509461125E4,
  242.                            4.5921953931549871457E4,
  243.                            6.7265770927008700853E4,
  244.                            3.3430575583588128105E4,
  245.                            2.5090809287301226727E3
  246.                         } ;
  247.  
  248.    static const REAL b[7] = {
  249.                            4.2313330701600911252E1,
  250.                            6.8718700749205790830E2,
  251.                            5.3941960214247511077E3,
  252.                            2.1213794301586595867E4,
  253.                            3.9307895800092710610E4,
  254.                            2.8729085735721942674E4,
  255.                            5.2264952788528545610E3
  256.                         } ;
  257.  
  258.    static const REAL c[8] = {
  259.                            1.42343711074968357734E0,
  260.                            4.63033784615654529590E0,
  261.                            5.76949722146069140550E0,
  262.                            3.64784832476320460504E0,
  263.                            1.27045825245236838258E0,
  264.                            2.41780725177450611770E-1,
  265.                            2.27238449892691845833E-2,
  266.                            7.74545014278341407640E-4
  267.                         } ;
  268.  
  269.    static const REAL d[7] = {
  270.                            2.05319162663775882187E0,
  271.                            1.67638483018380384940E0,
  272.                            6.89767334985100004550E-1,
  273.                            1.48103976427480074590E-1,
  274.                            1.51986665636164571966E-2,
  275.                            5.47593808499534494600E-4,
  276.                            1.05075007164441684324E-9
  277.                         } ;
  278.  
  279.    static REAL e[8] = {
  280.                            6.65790464350110377720E0,
  281.                            5.46378491116411436990E0,
  282.                            1.78482653991729133580E0,
  283.                            2.96560571828504891230E-1,
  284.                            2.65321895265761230930E-2,
  285.                            1.24266094738807843860E-3,
  286.                            2.71155556874348757815E-5,
  287.                            2.01033439929228813265E-7
  288.                         } ;
  289.  
  290.    static const REAL f[7] = {
  291.                            5.99832206555887937690E-1,
  292.                            1.36929880922735805310E-1,
  293.                            1.48753612908506148525E-2,
  294.                            7.86869131145613259100E-4,
  295.                            1.84631831751005468180E-5,
  296.                            1.42151175831644588870E-7,
  297.                            2.04426310338993978564E-15
  298.                         } ;
  299.  
  300.    REAL q = p - half, r, x ;
  301.    error = NOERROR ;
  302.    if ( fabs( q ) < SPLIT1 ) {
  303.       r = CONST1 - q * q ;
  304.       return q * ((((((( a[7] * r + a[6] ) * r + a[5] ) * r + a[4] ) * r
  305.                         + a[3] ) * r + a[2] ) * r + a[1] ) * r + a[0] ) /
  306.                  ((((((( b[6] * r + b[5] ) * r + b[4] ) * r + b[3] ) * r
  307.                         + b[2] ) * r + b[1] ) * r + b[0] ) * r + one ) ;
  308.    } else {
  309.       if ( q < zero )
  310.          r = p ;
  311.       else
  312.          r = one - p ;
  313.       if ( r <= zero ) {
  314.          error = BDARG ;
  315.          return zero ;
  316.       } // if
  317.       r = sqrt( -log( r ) ) ;
  318.       if ( r <= SPLIT2 ) {
  319.          r -= CONST2 ;
  320.          x = ((((((( c[7] * r + c[6] ) * r + c[5] ) * r + c[4] ) * r
  321.                      + c[3] ) * r + c[2] ) * r + c[1] ) * r + c[0] ) /
  322.                ((((((( d[6] * r + d[5] ) * r + d[4] ) * r + d[3] ) * r
  323.                      + d[2] ) * r + d[1] ) * r + d[0] ) * r + one ) ;
  324.       } else {
  325.          r -= SPLIT2 ;
  326.          x =  ((((((( e[7] * r + e[6] ) * r + e[5] ) * r + e[4] ) * r
  327.                      + e[3] ) * r + e[2] ) * r + e[1] ) * r + e[0] ) /
  328.                ((((((( f[6] * r + f[5] ) * r + f[4] ) * r + f[3] ) * r
  329.                      + f[2] ) * r + f[1] ) * r + f[0] ) * r + one ) ;
  330.       } // else
  331.       if ( q < zero )
  332.          x = -x ;
  333.       return x ;
  334.    } // else
  335. } // normalFunc inv
  336.  
  337. matrix normalFunc::inv( const matrix& x )
  338. {
  339.    static char *mName = "inv" ;
  340.    matFunc func( mName ) ; debugInfo( func ) ;
  341.    return matSpecFunc::inv(x) ;
  342. } // normalFunc inv( matrix )
  343.  
  344. normalFunc& normalFunc::par( REAL newLower, REAL newUpper )
  345. {
  346.    static char *mName = "par" ;
  347.    matFunc func( mName ) ; debugInfo( func ) ;
  348.    if ( newLower <= zero || newUpper <= zero)
  349.       errorExit( "par", NGPAR ) ;
  350.    lower = newLower ;
  351.    upper = newUpper ;
  352.    return *this ;
  353. } // normalFunc par
  354.  
  355. matrix normalFunc::operator () ( const matrix& x )
  356. {
  357.    static char *mName = "normal(m)" ;
  358.    matFunc func( mName ) ; debugInfo( func ) ;
  359.    INDEX nr = x.nRows(), nc = x.nCols() ;
  360.    matrix z(nr,nc) ;
  361.    transform(z,x) ;
  362.    return z ;
  363. } // normal op ()
  364.  
  365. REAL normalFunc::operator () ( REAL r )
  366. {
  367.    static char *mName = "normal(r)" ;
  368.    matFunc func( mName ) ; debugInfo( func ) ;
  369.    return value(r) ;
  370. } // normal op ()
  371.  
  372. normalFunc normal ;
  373.  
  374.