Changeset 266

Show
Ignore:
Timestamp:
02/09/09 23:15:51 (17 years ago)
Author:
smidl
Message:

switch to svn version of itpp

Location:
bdm
Files:
2 modified

Legend:

Unmodified
Added
Removed

  • bdm/itpp_ext.cpp

    r263 r266  
    6161//Gamma 
    6262 
    63         Gamma_RNG::Gamma_RNG ( double a, double b ) { 
    64                 setup ( a,b ); 
    65         } 
     63//      Gamma_RNG::Gamma_RNG ( double a, double b ) { 
     64//              setup ( a,b ); 
     65//      } 
    6666 
    6767 
    … …  
    8585                return temp; 
    8686        } 
    87 #define log std::log 
    88 #define exp std::exp 
    89 #define sqrt std::sqrt 
    90 #define R_FINITE std::isfinite 
    91  
    92         double  Gamma_RNG::sample() { 
    93                 //A copy of rgamma code from the R package!! 
    94                 // 
    95  
    96                 /* Constants : */ 
    97                 const static double sqrt32 = 5.656854; 
    98                 const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */ 
    99  
    100                 /* Coefficients q[k] - for q0 = sum(q[k]*a^(-k)) 
    101                  * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k) 
    102                  * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k) 
    103                  */ 
    104                 const static double q1 = 0.04166669; 
    105                 const static double q2 = 0.02083148; 
    106                 const static double q3 = 0.00801191; 
    107                 const static double q4 = 0.00144121; 
    108                 const static double q5 = -7.388e-5; 
    109                 const static double q6 = 2.4511e-4; 
    110                 const static double q7 = 2.424e-4; 
    111  
    112                 const static double a1 = 0.3333333; 
    113                 const static double a2 = -0.250003; 
    114                 const static double a3 = 0.2000062; 
    115                 const static double a4 = -0.1662921; 
    116                 const static double a5 = 0.1423657; 
    117                 const static double a6 = -0.1367177; 
    118                 const static double a7 = 0.1233795; 
    119  
    120                 /* State variables [FIXME for threading!] :*/ 
    121                 static double aa = 0.; 
    122                 static double aaa = 0.; 
    123                 static double s, s2, d;    /* no. 1 (step 1) */ 
    124                 static double q0, b, si, c;/* no. 2 (step 4) */ 
    125  
    126                 double e, p, q, r, t, u, v, w, x, ret_val; 
    127                 double a=alpha; 
    128                 double scale=1.0/beta; 
    129  
    130                 if ( !R_FINITE ( a ) || !R_FINITE ( scale ) || a < 0.0 || scale <= 0.0 ) 
    131                         {it_error ( "Gamma_RNG wrong parameters" );} 
    132  
    133                 if ( a < 1. ) { /* GS algorithm for parameters a < 1 */ 
    134                         if ( a == 0 ) 
    135                                 return 0.; 
    136                         e = 1.0 + exp_m1 * a; 
    137                         for ( ;; ) {  //VS repeat 
    138                                 p = e * unif_rand(); 
    139                                 if ( p >= 1.0 ) { 
    140                                         x = -log ( ( e - p ) / a ); 
    141                                         if ( exp_rand() >= ( 1.0 - a ) * log ( x ) ) 
    142                                                 break; 
    143                                 } 
    144                                 else { 
    145                                         x = exp ( log ( p ) / a ); 
    146                                         if ( exp_rand() >= x ) 
    147                                                 break; 
    148                                 } 
    149                         } 
    150                         return scale * x; 
    151                 } 
    152  
    153                 /* --- a >= 1 : GD algorithm --- */ 
    154  
    155                 /* Step 1: Recalculations of s2, s, d if a has changed */ 
    156                 if ( a != aa ) { 
    157                         aa = a; 
    158                         s2 = a - 0.5; 
    159                         s = sqrt ( s2 ); 
    160                         d = sqrt32 - s * 12.0; 
    161                 } 
    162                 /* Step 2: t = standard normal deviate, 
    163                            x = (s,1/2) -normal deviate. */ 
    164  
    165                 /* immediate acceptance (i) */ 
    166                 t = norm_rand(); 
    167                 x = s + 0.5 * t; 
    168                 ret_val = x * x; 
    169                 if ( t >= 0.0 ) 
    170                         return scale * ret_val; 
    171  
    172                 /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */ 
    173                 u = unif_rand(); 
    174                 if ( ( d * u ) <= ( t * t * t ) ) 
    175                         return scale * ret_val; 
    176  
    177                 /* Step 4: recalculations of q0, b, si, c if necessary */ 
    178  
    179                 if ( a != aaa ) { 
    180                         aaa = a; 
    181                         r = 1.0 / a; 
    182                         q0 = ( ( ( ( ( ( q7 * r + q6 ) * r + q5 ) * r + q4 ) * r + q3 ) * r 
    183                                  + q2 ) * r + q1 ) * r; 
    184  
    185                         /* Approximation depending on size of parameter a */ 
    186                         /* The constants in the expressions for b, si and c */ 
    187                         /* were established by numerical experiments */ 
    188  
    189                         if ( a <= 3.686 ) { 
    190                                 b = 0.463 + s + 0.178 * s2; 
    191                                 si = 1.235; 
    192                                 c = 0.195 / s - 0.079 + 0.16 * s; 
    193                         } 
    194                         else if ( a <= 13.022 ) { 
    195                                 b = 1.654 + 0.0076 * s2; 
    196                                 si = 1.68 / s + 0.275; 
    197                                 c = 0.062 / s + 0.024; 
    198                         } 
    199                         else { 
    200                                 b = 1.77; 
    201                                 si = 0.75; 
    202                                 c = 0.1515 / s; 
    203                         } 
    204                 } 
    205                 /* Step 5: no quotient test if x not positive */ 
    206  
    207                 if ( x > 0.0 ) { 
    208                         /* Step 6: calculation of v and quotient q */ 
    209                         v = t / ( s + s ); 
    210                         if ( fabs ( v ) <= 0.25 ) 
    211                                 q = q0 + 0.5 * t * t * ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v 
    212                                                              + a3 ) * v + a2 ) * v + a1 ) * v; 
    213                         else 
    214                                 q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v ); 
    215  
    216  
    217                         /* Step 7: quotient acceptance (q) */ 
    218                         if ( log ( 1.0 - u ) <= q ) 
    219                                 return scale * ret_val; 
    220                 } 
    221  
    222                 for ( ;; ) { //VS repeat 
    223                         /* Step 8: e = standard exponential deviate 
    224                          *      u =  0,1 -uniform deviate 
    225                          *      t = (b,si)-double exponential (laplace) sample */ 
    226                         e = exp_rand(); 
    227                         u = unif_rand(); 
    228                         u = u + u - 1.0; 
    229                         if ( u < 0.0 ) 
    230                                 t = b - si * e; 
    231                         else 
    232                                 t = b + si * e; 
    233                         /* Step  9:  rejection if t < tau(1) = -0.71874483771719 */ 
    234                         if ( t >= -0.71874483771719 ) { 
    235                                 /* Step 10:      calculation of v and quotient q */ 
    236                                 v = t / ( s + s ); 
    237                                 if ( fabs ( v ) <= 0.25 ) 
    238                                         q = q0 + 0.5 * t * t * 
    239                                             ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v + a3 ) * v 
    240                                                 + a2 ) * v + a1 ) * v; 
    241                                 else 
    242                                         q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v ); 
    243                                 /* Step 11:      hat acceptance (h) */ 
    244                                 /* (if q not positive go to step 8) */ 
    245                                 if ( q > 0.0 ) { 
    246                                         // TODO: w = expm1(q); 
    247                                         w = exp ( q )-1; 
    248                                         /*  ^^^^^ original code had approximation with rel.err < 2e-7 */ 
    249                                         /* if t is rejected sample again at step 8 */ 
    250                                         if ( ( c * fabs ( u ) ) <= ( w * exp ( e - 0.5 * t * t ) ) ) 
    251                                                 break; 
    252                                 } 
    253                         } 
    254                 } /* repeat .. until  `t' is accepted */ 
    255                 x = s + 0.5 * t; 
    256                 return scale * x * x; 
    257         } 
     87// #define log std::log 
     88// #define exp std::exp 
     89// #define sqrt std::sqrt 
     90// #define R_FINITE std::isfinite 
     91//  
     92//      double  Gamma_RNG::sample() { 
     93//              //A copy of rgamma code from the R package!! 
     94//              // 
     95//  
     96//              /* Constants : */ 
     97//              const static double sqrt32 = 5.656854; 
     98//              const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */ 
     99//  
     100//              /* Coefficients q[k] - for q0 = sum(q[k]*a^(-k)) 
     101//               * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k) 
     102//               * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k) 
     103//               */ 
     104//              const static double q1 = 0.04166669; 
     105//              const static double q2 = 0.02083148; 
     106//              const static double q3 = 0.00801191; 
     107//              const static double q4 = 0.00144121; 
     108//              const static double q5 = -7.388e-5; 
     109//              const static double q6 = 2.4511e-4; 
     110//              const static double q7 = 2.424e-4; 
     111//  
     112//              const static double a1 = 0.3333333; 
     113//              const static double a2 = -0.250003; 
     114//              const static double a3 = 0.2000062; 
     115//              const static double a4 = -0.1662921; 
     116//              const static double a5 = 0.1423657; 
     117//              const static double a6 = -0.1367177; 
     118//              const static double a7 = 0.1233795; 
     119//  
     120//              /* State variables [FIXME for threading!] :*/ 
     121//              static double aa = 0.; 
     122//              static double aaa = 0.; 
     123//              static double s, s2, d;    /* no. 1 (step 1) */ 
     124//              static double q0, b, si, c;/* no. 2 (step 4) */ 
     125//  
     126//              double e, p, q, r, t, u, v, w, x, ret_val; 
     127//              double a=alpha; 
     128//              double scale=1.0/beta; 
     129//  
     130//              if ( !R_FINITE ( a ) || !R_FINITE ( scale ) || a < 0.0 || scale <= 0.0 ) 
     131//                      {it_error ( "Gamma_RNG wrong parameters" );} 
     132//  
     133//              if ( a < 1. ) { /* GS algorithm for parameters a < 1 */ 
     134//                      if ( a == 0 ) 
     135//                              return 0.; 
     136//                      e = 1.0 + exp_m1 * a; 
     137//                      for ( ;; ) {  //VS repeat 
     138//                              p = e * unif_rand(); 
     139//                              if ( p >= 1.0 ) { 
     140//                                      x = -log ( ( e - p ) / a ); 
     141//                                      if ( exp_rand() >= ( 1.0 - a ) * log ( x ) ) 
     142//                                              break; 
     143//                              } 
     144//                              else { 
     145//                                      x = exp ( log ( p ) / a ); 
     146//                                      if ( exp_rand() >= x ) 
     147//                                              break; 
     148//                              } 
     149//                      } 
     150//                      return scale * x; 
     151//              } 
     152//  
     153//              /* --- a >= 1 : GD algorithm --- */ 
     154//  
     155//              /* Step 1: Recalculations of s2, s, d if a has changed */ 
     156//              if ( a != aa ) { 
     157//                      aa = a; 
     158//                      s2 = a - 0.5; 
     159//                      s = sqrt ( s2 ); 
     160//                      d = sqrt32 - s * 12.0; 
     161//              } 
     162//              /* Step 2: t = standard normal deviate, 
     163//                         x = (s,1/2) -normal deviate. */ 
     164//  
     165//              /* immediate acceptance (i) */ 
     166//              t = norm_rand(); 
     167//              x = s + 0.5 * t; 
     168//              ret_val = x * x; 
     169//              if ( t >= 0.0 ) 
     170//                      return scale * ret_val; 
     171//  
     172//              /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */ 
     173//              u = unif_rand(); 
     174//              if ( ( d * u ) <= ( t * t * t ) ) 
     175//                      return scale * ret_val; 
     176//  
     177//              /* Step 4: recalculations of q0, b, si, c if necessary */ 
     178//  
     179//              if ( a != aaa ) { 
     180//                      aaa = a; 
     181//                      r = 1.0 / a; 
     182//                      q0 = ( ( ( ( ( ( q7 * r + q6 ) * r + q5 ) * r + q4 ) * r + q3 ) * r 
     183//                               + q2 ) * r + q1 ) * r; 
     184//  
     185//                      /* Approximation depending on size of parameter a */ 
     186//                      /* The constants in the expressions for b, si and c */ 
     187//                      /* were established by numerical experiments */ 
     188//  
     189//                      if ( a <= 3.686 ) { 
     190//                              b = 0.463 + s + 0.178 * s2; 
     191//                              si = 1.235; 
     192//                              c = 0.195 / s - 0.079 + 0.16 * s; 
     193//                      } 
     194//                      else if ( a <= 13.022 ) { 
     195//                              b = 1.654 + 0.0076 * s2; 
     196//                              si = 1.68 / s + 0.275; 
     197//                              c = 0.062 / s + 0.024; 
     198//                      } 
     199//                      else { 
     200//                              b = 1.77; 
     201//                              si = 0.75; 
     202//                              c = 0.1515 / s; 
     203//                      } 
     204//              } 
     205//              /* Step 5: no quotient test if x not positive */ 
     206//  
     207//              if ( x > 0.0 ) { 
     208//                      /* Step 6: calculation of v and quotient q */ 
     209//                      v = t / ( s + s ); 
     210//                      if ( fabs ( v ) <= 0.25 ) 
     211//                              q = q0 + 0.5 * t * t * ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v 
     212//                                                           + a3 ) * v + a2 ) * v + a1 ) * v; 
     213//                      else 
     214//                              q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v ); 
     215//  
     216//  
     217//                      /* Step 7: quotient acceptance (q) */ 
     218//                      if ( log ( 1.0 - u ) <= q ) 
     219//                              return scale * ret_val; 
     220//              } 
     221//  
     222//              for ( ;; ) { //VS repeat 
     223//                      /* Step 8: e = standard exponential deviate 
     224//                       *      u =  0,1 -uniform deviate 
     225//                       *      t = (b,si)-double exponential (laplace) sample */ 
     226//                      e = exp_rand(); 
     227//                      u = unif_rand(); 
     228//                      u = u + u - 1.0; 
     229//                      if ( u < 0.0 ) 
     230//                              t = b - si * e; 
     231//                      else 
     232//                              t = b + si * e; 
     233//                      /* Step  9:  rejection if t < tau(1) = -0.71874483771719 */ 
     234//                      if ( t >= -0.71874483771719 ) { 
     235//                              /* Step 10:      calculation of v and quotient q */ 
     236//                              v = t / ( s + s ); 
     237//                              if ( fabs ( v ) <= 0.25 ) 
     238//                                      q = q0 + 0.5 * t * t * 
     239//                                          ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v + a3 ) * v 
     240//                                              + a2 ) * v + a1 ) * v; 
     241//                              else 
     242//                                      q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v ); 
     243//                              /* Step 11:      hat acceptance (h) */ 
     244//                              /* (if q not positive go to step 8) */ 
     245//                              if ( q > 0.0 ) { 
     246//                                      // TODO: w = expm1(q); 
     247//                                      w = exp ( q )-1; 
     248//                                      /*  ^^^^^ original code had approximation with rel.err < 2e-7 */ 
     249//                                      /* if t is rejected sample again at step 8 */ 
     250//                                      if ( ( c * fabs ( u ) ) <= ( w * exp ( e - 0.5 * t * t ) ) ) 
     251//                                              break; 
     252//                              } 
     253//                      } 
     254//              } /* repeat .. until  `t' is accepted */ 
     255//              x = s + 0.5 * t; 
     256//              return scale * x * x; 
     257//      } 
    258258 
    259259 

  • bdm/itpp_ext.h

    r263 r266  
    3838        } 
    3939 
     40#ifdef TTT 
    4041        /*! 
    4142          \brief Gamma distribution 
    … …  
    7374 
    7475        }; 
    75  
     76#endif 
    7677        bool qr ( const mat &A, mat &R ); 
    7778