home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 21 / IOPROG_21.ISO / SOFT / LIBMAT.ZIP / MATRAND.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-16  |  8.8 KB  |  322 lines
  1. /**************************************************/
  2. /*      matrand.c source for matRandom 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 "matrand.hpp"
  17. #include <math.h>
  18.  
  19.  
  20. static const REAL  zero = 0.0, half = 0.5, one = 1.0 , two = 2.0 ,
  21.              pi = 3.141592654 ;
  22.  
  23. /**************************************************************/
  24. /*                       matRandom methods                    */
  25. /**************************************************************/
  26.  
  27. matRandom::matRandom( int newX, int newY, int newZ )
  28. {
  29.     seed( newX, newY, newZ ) ;
  30.     normalOk = FALSE ;
  31.     normal0  = 0.0 ;
  32. } // matRandom
  33.  
  34. matRandom::matRandom( matRandom& ran )
  35. {
  36.     x = ran.x ;
  37.     y = ran.y ;
  38.     z = ran.z ;
  39.     normalOk = ran.normalOk ;
  40.     normal0  = ran.normal0 ;
  41. } // matRandom copy constr
  42.  
  43. void matRandom::operator = ( matRandom& ran )
  44. {
  45.    x = ran.x ;
  46.    y = ran.y ;
  47.    z = ran.z ;
  48.    normalOk = ran.normalOk ;
  49.    normal0  = ran.normal0 ;
  50. } // matRandom =
  51.  
  52. matRandom::~matRandom( void ) {} ;
  53.  
  54. void matRandom::seed( int newX, int newY, int newZ )
  55. {
  56.     static char *mName = "seed" ;
  57.     matFunc func( mName ) ; debugInfo( func ) ;
  58.     if ( newX < 1 || newX > 30000 ||
  59.           newY < 1 || newY > 30000 ||
  60.           newZ < 1 || newZ > 30000 )
  61.       errorExit( mName, NGPAR ) ;
  62.    x = newX ;
  63.    y = newY ;
  64.    z = newZ ;
  65. } // matRandom seed
  66.  
  67. outFile& matRandom::put( outFile& f ) M_CONST
  68. {
  69.    f.putIndex(x).putIndex(y).putIndex(z).newLine() ;
  70.    return f ;
  71. } // matRandom put
  72.  
  73. inFile& matRandom::get( inFile& f )
  74. {
  75.    INDEX x1, y1, z1 ;
  76.    f.getIndex(x1).getIndex(y1).getIndex(z1).nextLine() ;
  77.     x = (INDEX)(x1) ;
  78.    y = (INDEX)(y1) ;
  79.     z = (INDEX)(z1) ;
  80.    return f ;
  81. } // matRandom get
  82.  
  83. outFile& matRandom::info( outFile& f ) M_CONST
  84. {
  85.    if ( matListCtrl() > 3 )
  86.       return f ;
  87.    objectInfo(f) ;
  88.    putField( x, f ) ;
  89.    putField( y, f ) ;
  90.    putField( z, f ) ;
  91.    f.newLine() ;
  92.    return f ;
  93. } // matRandom info
  94.  
  95. REAL matRandom::uniform( void )
  96. /*
  97.    Returns random double between 0 and 1.
  98.    Based on Wichmann & Hill Applied Stats AS183
  99.    1982 Vol 31 188-190.
  100. */
  101. {
  102.    // static char *mName = "uniform" ;
  103.    // matFunc func( mName ) ; debugInfo( func ) ;
  104.    REAL ran;
  105.  
  106.    x = 171 * ( x % 177 ) - 2 * ( x / 177 );
  107.    if ( x < 0 )
  108.       x += 30269;
  109.    y = 172 * ( y % 176 ) - 35 * ( y / 176 );
  110.    if ( y < 0 )
  111.       y += 30307;
  112.    z = 170 * ( z % 178 ) - 63 * ( z / 178 );
  113.    if ( z < 0 )
  114.       z += 30323;
  115.     ran = REAL( x / 30269.0 ) + REAL( y / 30307.0 )
  116.      + REAL( z / 30323.0 ) ;
  117.     ran -= floor( ran );
  118.     return( ran );
  119. } // uniform
  120.  
  121. REAL matRandom::normal( REAL mean, REAL std )
  122. {
  123.     // static char *mName = "normal" ;
  124.     // matFunc func( mName ) ; debugInfo( func ) ;
  125.     REAL normal1, normal2, v1, v2, s ;
  126.     if ( !normalOk ) {
  127.         do {
  128.             v1 = 2.0 * uniform() - 1.0;
  129.             v2 = 2.0 * uniform() - 1.0;
  130.             s  = v1 * v1 + v2 * v2;
  131.       } while ( s >= (REAL) 1.0 );
  132.       normal1 = normal2 = sqrt( (REAL) -2.0 * log(s) / s );
  133.       normal1 = v1 * normal1 ;
  134.       normal2 = v2 * normal2 ;
  135.       normal0 = normal2 ;
  136.       normalOk = TRUE ;
  137.       return mean + std * normal1 ;
  138.    } else {
  139.       normalOk = FALSE ;
  140.       return mean + std * normal0 ;
  141.    } // else
  142. } // normal
  143.  
  144. #include "matspec.hpp"
  145.  
  146. REAL matRandom::expon( REAL mean )
  147. {
  148.    // static char *mName = "expon" ;
  149.    // matFunc func( mName ) ; debugInfo( func ) ;
  150.    return -log( uniform() ) * mean ;
  151. } // expon
  152.  
  153. REAL matRandom::gamma( REAL a )
  154. {
  155.    // See Press et al pp.220-1
  156.    static char *mName = "gamma" ;
  157.    matFunc func( mName ) ; // debugInfo( func ) ;
  158.    REAL gamma, v1, v2, s, t, b, ratio ;
  159.    INDEX i ;
  160.    if ( a < 1 ) {
  161.       errorExit( mName, NGPAR ) ;
  162.    } else if ( a < 6 ) {
  163.       gamma = one ;
  164.       for ( i = 1 ; i <= a ; i++ )
  165.          gamma *= uniform() ;
  166.       gamma = -log( gamma ) ;
  167.    } else {
  168.       do {
  169.          do {
  170.             do {
  171.                v1 = two * uniform() - one ;
  172.                v2 = two * uniform() - one ;
  173.             } while ( v1 * v1 + v2 * v2 > one ) ;
  174.             t = v2 / v1 ;  // ?? small v1 ??
  175.             b = a - 1 ;
  176.             s = sqrt( two * b + one ) ;
  177.             gamma = s * t + b ;
  178.          } while ( gamma <= zero ) ;
  179.          ratio = ( one + t * t ) * exp( b * log( gamma/ b ) - s * t ) ;
  180.       } while ( uniform() > ratio ) ;
  181.    } // else
  182.    return gamma ;
  183. } // matRandom gamma
  184.  
  185.  
  186. REAL matRandom::poisson( REAL mean )
  187. {
  188.    // static char *mName = "poisson" ;
  189.    // matFunc func( mName ) ; debugInfo( func ) ;
  190.    static REAL sq, logm, g, oldm = (-1.0) ;
  191.    // matError error ;
  192.    REAL pois, test, y ;
  193.    if ( mean < 12.0 ) {
  194.       if ( mean != oldm ) {
  195.          oldm = mean ;
  196.          g = exp( -mean ) ;
  197.       } // if
  198.       pois = -1 ;
  199.       test = one ;
  200.       do {
  201.          pois += one ;
  202.          test *= uniform() ;
  203.       } while ( test > g ) ;
  204.    } else {
  205.       if ( mean != oldm ) {
  206.          oldm = mean ;
  207.          sq = sqrt( two * mean ) ;
  208.          logm = log( mean ) ;
  209.          g = mean * logm - logGamma( mean + one ) ;
  210.       } // if
  211.       do {
  212.          do {
  213.             y = tan( pi * uniform() ) ;
  214.             pois = sq * y + mean ;
  215.          } while ( pois < zero ) ;
  216.          pois = floor( pois ) ;
  217.          test = 0.9 * ( one + y * y ) *
  218.                exp( pois * logm - logGamma( pois + one ) - g ) ;
  219.       } while ( uniform() > test ) ;
  220.    } // else
  221.    return pois ;
  222. } // matRandom poisson
  223.  
  224. REAL matRandom::binomial( int n, REAL p )
  225. {
  226.    static char *mName = "binomial" ;
  227.    matFunc func( mName ) ; // debugInfo( func ) ;
  228.    static REAL oldp = (-1.0), q, logp, logpc, olddn, oldg ;
  229.    static int oldn = (-1) ;
  230.    REAL bnl, mean, cmp, g, g1, g2, angle, prob, sq, test, y ;
  231.    // matError error ;
  232.    int j ;
  233.    if ( n <= 0 )
  234.       errorExit( mName, NGPAR ) ;
  235.    prob = ( p <= half ? p : one - p ) ;
  236.    mean = n * prob ;
  237.    if ( n < 25 ) {
  238.       bnl = zero ;
  239.       for ( j = 0; j <= n; j++ )
  240.          if ( uniform() < prob )
  241.             bnl += one ;
  242.    } else if ( mean < one ) {
  243.       g = exp( -mean ) ;
  244.       test = one ;
  245.       for ( j = 0; j <= n && test >= g ; j++ )
  246.          test *= uniform() ;
  247.       bnl = ( j <= n ? j : n ) ;
  248.    } else {
  249.       if ( n != oldn ) {
  250.          oldn = n ;
  251.          olddn = (REAL) oldn ;
  252.          oldg = logGamma( olddn + one ) ;
  253.       } // if
  254.       if ( prob != oldp ) {
  255.          q = one - prob ;
  256.          logp = log( prob ) ;
  257.          logpc = log( q ) ;
  258.          oldp = prob ;
  259.       } // if
  260.       sq = sqrt( two * mean * q ) ;
  261.       do {
  262.          do {
  263.             angle = pi * uniform() ;
  264.             y = tan( angle ) ;
  265.             cmp = sq * y + mean ;
  266.          } while ( cmp < zero || cmp >= ( olddn + one ) ) ;
  267.          cmp = floor( cmp ) ;
  268.          g1 = logGamma( cmp + one ) ;
  269.          g2 = logGamma( olddn - cmp + one ) ;
  270.          test = 1.2 * sq * ( one + y * y ) *
  271.                 exp( oldg - g1 - g2 + cmp * logp
  272.                      + ( olddn - cmp ) * logpc ) ;
  273.       } while ( uniform() > test ) ;
  274.       bnl = cmp ;
  275.    } // else
  276.    if ( prob != p )
  277.       bnl = n - bnl ;
  278.    return bnl ;
  279. } // binomial
  280.  
  281. matrix& matRandom::operator()( matrix& x, distribution dist,
  282.                                REAL p1, REAL p2 )
  283. {
  284.    static char* mName = "op()" ;
  285.    matFunc func( mName ) ; debugInfo( func ) ;
  286.    matError error = NOERROR ;
  287.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  288.    int n ;
  289.    if ( ( dist == GAMMA ) || ( dist == BINOMIAL ) )
  290.       n = (int) p1 ;
  291.    for ( j = 1 ; j <= nc && error == NOERROR ; j++ ) {
  292.       for ( i = 1 ; i <= nr && error == NOERROR ; i++ ) {
  293.          switch( dist ) {
  294.          case UNIFORM :
  295.             x(i,j) = (REAL) uniform() ;
  296.             break ;
  297.          case NORMAL :
  298.             x(i,j) = (REAL) normal( p1, p2 ) ;
  299.             break ;
  300.          case EXPON :
  301.             x(i,j) = (REAL) expon( p1 ) ;
  302.             break ;
  303.          case GAMMA :
  304.             x(i,j) = (REAL) gamma( n ) ;
  305.             break ;
  306.          case POISSON :
  307.             x(i,j) = (REAL) poisson( p1 ) ;
  308.             break ;
  309.          case BINOMIAL :
  310.             x(i,j) = (REAL) binomial( n, p2 ) ;
  311.             break ;
  312.          default :
  313.             matErrNumExit( "rand() ", NRANG ) ;
  314.             break ;
  315.          } // switch
  316.       } // for j
  317.    } // for i
  318.    if ( error )
  319.       errorExit( mName, error ) ;
  320.    return x ;
  321. } // matrandom()
  322.