IT++ Logo
random.cpp
Go to the documentation of this file.
00001 
00029 #include <itpp/base/random.h>
00030 #include <itpp/base/itcompat.h>
00031 #include <itpp/base/math/elem_math.h>
00032 
00033 
00034 namespace itpp
00035 {
00036 
00037 // ----------------------------------------------------------------------
00038 // Random_Generator (DSFMT_19937_RNG)
00039 // ----------------------------------------------------------------------
00040 
00041 template <>
00042 bool Random_Generator::initialized = false;
00043 template <>
00044 bool Random_Generator::bigendian = is_bigendian();
00045 template <>
00046 Random_Generator::w128_t Random_Generator::status[N + 1] = { };
00047 template <>
00048 int Random_Generator::idx = 0;
00049 template <>
00050 unsigned int Random_Generator::last_seed = 0U;
00051 #if defined(__SSE2__)
00052 template <>
00053 __m128i Random_Generator::sse2_param_mask = _mm_set_epi32(0, 0, 0, 0);
00054 #endif // __SSE2__
00055 
00056 
00057 // Set the seed of the Global Random Number Generator
00058 void RNG_reset(unsigned int seed)
00059 {
00060   Random_Generator RNG;
00061   RNG.reset(seed);
00062 }
00063 
00064 // Set the seed of the Global Random Number Generator to the same as last time
00065 void RNG_reset()
00066 {
00067   Random_Generator RNG;
00068   RNG.reset();
00069 }
00070 
00071 // Set a random seed for the Global Random Number Generator
00072 void RNG_randomize()
00073 {
00074   Random_Generator RNG;
00075   RNG.randomize();
00076 }
00077 
00078 // Save current full state of generator in memory
00079 void RNG_get_state(ivec &state)
00080 {
00081   Random_Generator RNG;
00082   state = RNG.get_state();
00083 }
00084 
00085 // Resume the state saved in memory
00086 void RNG_set_state(const ivec &state)
00087 {
00088   Random_Generator RNG;
00089   RNG.set_state(state);
00090 }
00091 
00093 // I_Uniform_RNG
00095 
00096 I_Uniform_RNG::I_Uniform_RNG(int min, int max)
00097 {
00098   setup(min, max);
00099 }
00100 
00101 void I_Uniform_RNG::setup(int min, int max)
00102 {
00103   if (min <= max) {
00104     lo = min;
00105     hi = max;
00106   }
00107   else {
00108     lo = max;
00109     hi = min;
00110   }
00111 }
00112 
00113 void I_Uniform_RNG::get_setup(int &min, int &max) const
00114 {
00115   min = lo;
00116   max = hi;
00117 }
00118 
00119 ivec I_Uniform_RNG::operator()(int n)
00120 {
00121   ivec vv(n);
00122 
00123   for (int i=0; i < n; i++)
00124     vv(i) = sample();
00125 
00126   return vv;
00127 }
00128 
00129 imat I_Uniform_RNG::operator()(int h, int w)
00130 {
00131   imat mm(h, w);
00132   int i, j;
00133 
00134   for (i = 0; i < h; i++)
00135     for (j = 0; j < w; j++)
00136       mm(i, j) = sample();
00137 
00138   return mm;
00139 }
00140 
00142 // Uniform_RNG
00144 
00145 Uniform_RNG::Uniform_RNG(double min, double max)
00146 {
00147   setup(min, max);
00148 }
00149 
00150 void Uniform_RNG::setup(double min, double max)
00151 {
00152   if (min <= max) {
00153     lo_bound = min;
00154     hi_bound = max;
00155   }
00156   else {
00157     lo_bound = max;
00158     hi_bound = min;
00159   }
00160 }
00161 
00162 void Uniform_RNG::get_setup(double &min, double &max) const
00163 {
00164   min = lo_bound;
00165   max = hi_bound;
00166 }
00167 
00169 // Exp_RNG
00171 
00172 Exponential_RNG::Exponential_RNG(double lambda)
00173 {
00174   setup(lambda);
00175 }
00176 
00177 vec Exponential_RNG::operator()(int n)
00178 {
00179   vec vv(n);
00180 
00181   for (int i=0; i < n; i++)
00182     vv(i) = sample();
00183 
00184   return vv;
00185 }
00186 
00187 mat Exponential_RNG::operator()(int h, int w)
00188 {
00189   mat mm(h, w);
00190   int i, j;
00191 
00192   for (i = 0; i < h; i++)
00193     for (j = 0; j < w; j++)
00194       mm(i, j) = sample();
00195 
00196   return mm;
00197 }
00198 
00200 // Gamma_RNG
00202 
00203 vec Gamma_RNG::operator()(int n)
00204 {
00205   vec vv(n);
00206   for (int i = 0; i < n; i++)
00207     vv(i) = sample();
00208   return vv;
00209 }
00210 
00211 mat Gamma_RNG::operator()(int r, int c)
00212 {
00213   mat mm(r, c);
00214   for (int i = 0; i < r * c; i++)
00215     mm(i) = sample();
00216   return mm;
00217 }
00218 
00219 double Gamma_RNG::sample()
00220 {
00221   // A copy of rgamma code from the R package, adapted to IT++ by Vasek
00222   // Smidl
00223 
00224   /* Constants : */
00225   const static double sqrt32 = 5.656854;
00226   const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */
00227 
00228   const static double q1 = 0.04166669;
00229   const static double q2 = 0.02083148;
00230   const static double q3 = 0.00801191;
00231   const static double q4 = 0.00144121;
00232   const static double q5 = -7.388e-5;
00233   const static double q6 = 2.4511e-4;
00234   const static double q7 = 2.424e-4;
00235 
00236   const static double a1 = 0.3333333;
00237   const static double a2 = -0.250003;
00238   const static double a3 = 0.2000062;
00239   const static double a4 = -0.1662921;
00240   const static double a5 = 0.1423657;
00241   const static double a6 = -0.1367177;
00242   const static double a7 = 0.1233795;
00243 
00244   /* State variables [FIXME for threading!] :*/
00245   static double aa = 0.;
00246   static double aaa = 0.;
00247   static double s, s2, d;    /* no. 1 (step 1) */
00248   static double q0, b, si, c;/* no. 2 (step 4) */
00249 
00250   double e, p, q, r, t, u, v, w, x, ret_val;
00251   double a = alpha;
00252   double scale = 1.0 / beta;
00253 
00254   it_error_if(!std::isfinite(a) || !std::isfinite(scale) || (a < 0.0)
00255               || (scale <= 0.0), "Gamma_RNG wrong parameters");
00256 
00257   if (a < 1.) { /* GS algorithm for parameters a < 1 */
00258     if (a == 0)
00259       return 0.;
00260     e = 1.0 + exp_m1 * a;
00261     for (;;) { //VS repeat
00262       p = e * RNG.genrand_open_open();
00263       if (p >= 1.0) {
00264         x = -std::log((e - p) / a);
00265         if (-std::log(RNG.genrand_open_close()) >= (1.0 - a) * std::log(x))
00266           break;
00267       }
00268       else {
00269         x = std::exp(std::log(p) / a);
00270         if (-std::log(RNG.genrand_open_close()) >= x)
00271           break;
00272       }
00273     }
00274     return scale * x;
00275   }
00276 
00277   /* --- a >= 1 : GD algorithm --- */
00278 
00279   /* Step 1: Recalculations of s2, s, d if a has changed */
00280   if (a != aa) {
00281     aa = a;
00282     s2 = a - 0.5;
00283     s = std::sqrt(s2);
00284     d = sqrt32 - s * 12.0;
00285   }
00286   /* Step 2: t = standard normal deviate, x = (s,1/2) -normal deviate. */
00287   /* immediate acceptance (i) */
00288   t = NRNG.sample();
00289   x = s + 0.5 * t;
00290   ret_val = x * x;
00291   if (t >= 0.0)
00292     return scale * ret_val;
00293 
00294   /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */
00295   u = RNG.genrand_close_open();
00296   if ((d * u) <= (t * t * t))
00297     return scale * ret_val;
00298 
00299   /* Step 4: recalculations of q0, b, si, c if necessary */
00300   if (a != aaa) {
00301     aaa = a;
00302     r = 1.0 / a;
00303     q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r
00304            + q2) * r + q1) * r;
00305 
00306     /* Approximation depending on size of parameter a */
00307     /* The constants in the expressions for b, si and c */
00308     /* were established by numerical experiments */
00309     if (a <= 3.686) {
00310       b = 0.463 + s + 0.178 * s2;
00311       si = 1.235;
00312       c = 0.195 / s - 0.079 + 0.16 * s;
00313     }
00314     else if (a <= 13.022) {
00315       b = 1.654 + 0.0076 * s2;
00316       si = 1.68 / s + 0.275;
00317       c = 0.062 / s + 0.024;
00318     }
00319     else {
00320       b = 1.77;
00321       si = 0.75;
00322       c = 0.1515 / s;
00323     }
00324   }
00325 
00326   /* Step 5: no quotient test if x not positive */
00327   if (x > 0.0) {
00328     /* Step 6: calculation of v and quotient q */
00329     v = t / (s + s);
00330     if (std::fabs(v) <= 0.25)
00331       q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v
00332                                 + a3) * v + a2) * v + a1) * v;
00333     else
00334       q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
00335 
00336     /* Step 7: quotient acceptance (q) */
00337     if (log(1.0 - u) <= q)
00338       return scale * ret_val;
00339   }
00340 
00341   for (;;) { //VS repeat
00342     /* Step 8: e = standard exponential deviate
00343      *         u =  0,1 -uniform deviate
00344      *         t = (b,si)-double exponential (laplace) sample */
00345     e = -std::log(RNG.genrand_open_close()); //see Exponential_RNG
00346     u = RNG.genrand_open_close();
00347     u = u + u - 1.0;
00348     if (u < 0.0)
00349       t = b - si * e;
00350     else
00351       t = b + si * e;
00352     /* Step 9: rejection if t < tau(1) = -0.71874483771719 */
00353     if (t >= -0.71874483771719) {
00354       /* Step 10:  calculation of v and quotient q */
00355       v = t / (s + s);
00356       if (std::fabs(v) <= 0.25)
00357         q = q0 + 0.5 * t * t *
00358             ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v
00359               + a2) * v + a1) * v;
00360       else
00361         q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
00362       /* Step 11: hat acceptance (h) */
00363       /* (if q not positive go to step 8) */
00364       if (q > 0.0) {
00365         // Try to use w = expm1(q); (Not supported on w32)
00366         w = expm1(q);
00367         /*  ^^^^^ original code had approximation with rel.err < 2e-7 */
00368         /* if t is rejected sample again at step 8 */
00369         if ((c * std::fabs(u)) <= (w * std::exp(e - 0.5 * t * t)))
00370           break;
00371       }
00372     }
00373   } /* repeat .. until `t' is accepted */
00374   x = s + 0.5 * t;
00375   return scale * x * x;
00376 }
00377 
00379 // Normal_RNG
00381 
00382 void Normal_RNG::get_setup(double &meanval, double &variance) const
00383 {
00384   meanval = mean;
00385   variance = sigma * sigma;
00386 }
00387 
00388 // (Ziggurat) tabulated values for the heigt of the Ziggurat levels
00389 const double Normal_RNG::ytab[128] = {
00390   1, 0.963598623011, 0.936280813353, 0.913041104253,
00391   0.892278506696, 0.873239356919, 0.855496407634, 0.838778928349,
00392   0.822902083699, 0.807732738234, 0.793171045519, 0.779139726505,
00393   0.765577436082, 0.752434456248, 0.739669787677, 0.727249120285,
00394   0.715143377413, 0.703327646455, 0.691780377035, 0.68048276891,
00395   0.669418297233, 0.65857233912, 0.647931876189, 0.637485254896,
00396   0.62722199145, 0.617132611532, 0.607208517467, 0.597441877296,
00397   0.587825531465, 0.578352913803, 0.569017984198, 0.559815170911,
00398   0.550739320877, 0.541785656682, 0.532949739145, 0.524227434628,
00399   0.515614886373, 0.507108489253, 0.498704867478, 0.490400854812,
00400   0.482193476986, 0.47407993601, 0.466057596125, 0.458123971214,
00401   0.450276713467, 0.442513603171, 0.434832539473, 0.427231532022,
00402   0.419708693379, 0.41226223212, 0.404890446548, 0.397591718955,
00403   0.390364510382, 0.383207355816, 0.376118859788, 0.369097692334,
00404   0.362142585282, 0.355252328834, 0.348425768415, 0.341661801776,
00405   0.334959376311, 0.328317486588, 0.321735172063, 0.31521151497,
00406   0.308745638367, 0.302336704338, 0.29598391232, 0.289686497571,
00407   0.283443729739, 0.27725491156, 0.271119377649, 0.265036493387,
00408   0.259005653912, 0.253026283183, 0.247097833139, 0.241219782932,
00409   0.235391638239, 0.229612930649, 0.223883217122, 0.218202079518,
00410   0.212569124201, 0.206983981709, 0.201446306496, 0.195955776745,
00411   0.190512094256, 0.185114984406, 0.179764196185, 0.174459502324,
00412   0.169200699492, 0.1639876086, 0.158820075195, 0.153697969964,
00413   0.148621189348, 0.143589656295, 0.138603321143, 0.133662162669,
00414   0.128766189309, 0.123915440582, 0.119109988745, 0.114349940703,
00415   0.10963544023, 0.104966670533, 0.100343857232, 0.0957672718266,
00416   0.0912372357329, 0.0867541250127, 0.082318375932, 0.0779304915295,
00417   0.0735910494266, 0.0693007111742, 0.065060233529, 0.0608704821745,
00418   0.056732448584, 0.05264727098, 0.0486162607163, 0.0446409359769,
00419   0.0407230655415, 0.0368647267386, 0.0330683839378, 0.0293369977411,
00420   0.0256741818288, 0.0220844372634, 0.0185735200577, 0.0151490552854,
00421   0.0118216532614, 0.00860719483079, 0.00553245272614, 0.00265435214565
00422 };
00423 
00424 /*
00425  * (Ziggurat) tabulated values for 2^24 times x[i]/x[i+1], used to accept
00426  * for U*x[i+1]<=x[i] without any floating point operations
00427  */
00428 const unsigned int Normal_RNG::ktab[128] = {
00429   0, 12590644, 14272653, 14988939,
00430   15384584, 15635009, 15807561, 15933577,
00431   16029594, 16105155, 16166147, 16216399,
00432   16258508, 16294295, 16325078, 16351831,
00433   16375291, 16396026, 16414479, 16431002,
00434   16445880, 16459343, 16471578, 16482744,
00435   16492970, 16502368, 16511031, 16519039,
00436   16526459, 16533352, 16539769, 16545755,
00437   16551348, 16556584, 16561493, 16566101,
00438   16570433, 16574511, 16578353, 16581977,
00439   16585398, 16588629, 16591685, 16594575,
00440   16597311, 16599901, 16602354, 16604679,
00441   16606881, 16608968, 16610945, 16612818,
00442   16614592, 16616272, 16617861, 16619363,
00443   16620782, 16622121, 16623383, 16624570,
00444   16625685, 16626730, 16627708, 16628619,
00445   16629465, 16630248, 16630969, 16631628,
00446   16632228, 16632768, 16633248, 16633671,
00447   16634034, 16634340, 16634586, 16634774,
00448   16634903, 16634972, 16634980, 16634926,
00449   16634810, 16634628, 16634381, 16634066,
00450   16633680, 16633222, 16632688, 16632075,
00451   16631380, 16630598, 16629726, 16628757,
00452   16627686, 16626507, 16625212, 16623794,
00453   16622243, 16620548, 16618698, 16616679,
00454   16614476, 16612071, 16609444, 16606571,
00455   16603425, 16599973, 16596178, 16591995,
00456   16587369, 16582237, 16576520, 16570120,
00457   16562917, 16554758, 16545450, 16534739,
00458   16522287, 16507638, 16490152, 16468907,
00459   16442518, 16408804, 16364095, 16301683,
00460   16207738, 16047994, 15704248, 15472926
00461 };
00462 
00463 // (Ziggurat) tabulated values of 2^{-24}*x[i]
00464 const double Normal_RNG::wtab[128] = {
00465   1.62318314817e-08, 2.16291505214e-08, 2.54246305087e-08, 2.84579525938e-08,
00466   3.10340022482e-08, 3.33011726243e-08, 3.53439060345e-08, 3.72152672658e-08,
00467   3.8950989572e-08, 4.05763964764e-08, 4.21101548915e-08, 4.35664624904e-08,
00468   4.49563968336e-08, 4.62887864029e-08, 4.75707945735e-08, 4.88083237257e-08,
00469   5.00063025384e-08, 5.11688950428e-08, 5.22996558616e-08, 5.34016475624e-08,
00470   5.44775307871e-08, 5.55296344581e-08, 5.65600111659e-08, 5.75704813695e-08,
00471   5.85626690412e-08, 5.95380306862e-08, 6.04978791776e-08, 6.14434034901e-08,
00472   6.23756851626e-08, 6.32957121259e-08, 6.42043903937e-08, 6.51025540077e-08,
00473   6.59909735447e-08, 6.68703634341e-08, 6.77413882848e-08, 6.8604668381e-08,
00474   6.94607844804e-08, 7.03102820203e-08, 7.11536748229e-08, 7.1991448372e-08,
00475   7.2824062723e-08, 7.36519550992e-08, 7.44755422158e-08, 7.52952223703e-08,
00476   7.61113773308e-08, 7.69243740467e-08, 7.77345662086e-08, 7.85422956743e-08,
00477   7.93478937793e-08, 8.01516825471e-08, 8.09539758128e-08, 8.17550802699e-08,
00478   8.25552964535e-08, 8.33549196661e-08, 8.41542408569e-08, 8.49535474601e-08,
00479   8.57531242006e-08, 8.65532538723e-08, 8.73542180955e-08, 8.8156298059e-08,
00480   8.89597752521e-08, 8.97649321908e-08, 9.05720531451e-08, 9.138142487e-08,
00481   9.21933373471e-08, 9.30080845407e-08, 9.38259651738e-08, 9.46472835298e-08,
00482   9.54723502847e-08, 9.63014833769e-08, 9.71350089201e-08, 9.79732621669e-08,
00483   9.88165885297e-08, 9.96653446693e-08, 1.00519899658e-07, 1.0138063623e-07,
00484   1.02247952126e-07, 1.03122261554e-07, 1.04003996769e-07, 1.04893609795e-07,
00485   1.05791574313e-07, 1.06698387725e-07, 1.07614573423e-07, 1.08540683296e-07,
00486   1.09477300508e-07, 1.1042504257e-07, 1.11384564771e-07, 1.12356564007e-07,
00487   1.13341783071e-07, 1.14341015475e-07, 1.15355110887e-07, 1.16384981291e-07,
00488   1.17431607977e-07, 1.18496049514e-07, 1.19579450872e-07, 1.20683053909e-07,
00489   1.21808209468e-07, 1.2295639141e-07, 1.24129212952e-07, 1.25328445797e-07,
00490   1.26556042658e-07, 1.27814163916e-07, 1.29105209375e-07, 1.30431856341e-07,
00491   1.31797105598e-07, 1.3320433736e-07, 1.34657379914e-07, 1.36160594606e-07,
00492   1.37718982103e-07, 1.39338316679e-07, 1.41025317971e-07, 1.42787873535e-07,
00493   1.44635331499e-07, 1.4657889173e-07, 1.48632138436e-07, 1.50811780719e-07,
00494   1.53138707402e-07, 1.55639532047e-07, 1.58348931426e-07, 1.61313325908e-07,
00495   1.64596952856e-07, 1.68292495203e-07, 1.72541128694e-07, 1.77574279496e-07,
00496   1.83813550477e-07, 1.92166040885e-07, 2.05295471952e-07, 2.22600839893e-07
00497 };
00498 
00499 // (Ziggurat) position of right-most step
00500 const double Normal_RNG::PARAM_R = 3.44428647676;
00501 
00502 // Get a Normal distributed (0,1) sample
00503 double Normal_RNG::sample()
00504 {
00505   uint32_t u, sign, i, j;
00506   double x, y;
00507 
00508   while (true) {
00509     u = RNG.genrand_uint32();
00510     sign = u & 0x80;            // 1 bit for the sign
00511     i = u & 0x7f;               // 7 bits to choose the step
00512     j = u >> 8;                 // 24 bits for the x-value
00513 
00514     x = j * wtab[i];
00515 
00516     if (j < ktab[i])
00517       break;
00518 
00519     if (i < 127) {
00520       y = ytab[i+1] + (ytab[i] - ytab[i+1]) * RNG.genrand_close_open();
00521     }
00522     else {
00523       x = PARAM_R - std::log(1.0 - RNG.genrand_close_open()) / PARAM_R;
00524       y = std::exp(-PARAM_R * (x - 0.5 * PARAM_R)) * RNG.genrand_close_open();
00525     }
00526 
00527     if (y < std::exp(-0.5 * x * x))
00528       break;
00529   }
00530   return sign ? x : -x;
00531 }
00532 
00533 
00535 // Laplace_RNG
00537 
00538 Laplace_RNG::Laplace_RNG(double meanval, double variance)
00539 {
00540   setup(meanval, variance);
00541 }
00542 
00543 void Laplace_RNG::setup(double meanval, double variance)
00544 {
00545   mean = meanval;
00546   var = variance;
00547   sqrt_12var = std::sqrt(variance / 2.0);
00548 }
00549 
00550 void Laplace_RNG::get_setup(double &meanval, double &variance) const
00551 {
00552   meanval = mean;
00553   variance = var;
00554 }
00555 
00556 
00557 
00558 vec Laplace_RNG::operator()(int n)
00559 {
00560   vec vv(n);
00561 
00562   for (int i=0; i < n; i++)
00563     vv(i) = sample();
00564 
00565   return vv;
00566 }
00567 
00568 mat Laplace_RNG::operator()(int h, int w)
00569 {
00570   mat mm(h, w);
00571   int i, j;
00572 
00573   for (i = 0; i < h; i++)
00574     for (j = 0; j < w; j++)
00575       mm(i, j) = sample();
00576 
00577   return mm;
00578 }
00579 
00581 // AR1_Normal_RNG
00583 
00584 AR1_Normal_RNG::AR1_Normal_RNG(double meanval, double variance, double rho)
00585 {
00586   setup(meanval, variance, rho);
00587 }
00588 
00589 void AR1_Normal_RNG::setup(double meanval, double variance, double rho)
00590 {
00591   mean = meanval;
00592   var = variance;
00593   r = rho;
00594   factr = -2.0 * var * (1.0 - rho * rho);
00595   mem = 0.0;
00596   odd = true;
00597 }
00598 
00599 void AR1_Normal_RNG::get_setup(double &meanval, double &variance,
00600                                double &rho) const
00601 {
00602   meanval = mean;
00603   variance = var;
00604   rho = r;
00605 }
00606 
00607 vec AR1_Normal_RNG::operator()(int n)
00608 {
00609   vec vv(n);
00610 
00611   for (int i=0; i < n; i++)
00612     vv(i) = sample();
00613 
00614   return vv;
00615 }
00616 
00617 mat AR1_Normal_RNG::operator()(int h, int w)
00618 {
00619   mat mm(h, w);
00620   int i, j;
00621 
00622   for (i = 0; i < h; i++)
00623     for (j = 0; j < w; j++)
00624       mm(i, j) = sample();
00625 
00626   return mm;
00627 }
00628 
00629 void AR1_Normal_RNG::reset()
00630 {
00631   mem = 0.0;
00632 }
00633 
00635 // Weibull_RNG
00637 
00638 Weibull_RNG::Weibull_RNG(double lambda, double beta)
00639 {
00640   setup(lambda, beta);
00641 }
00642 
00643 void Weibull_RNG::setup(double lambda, double beta)
00644 {
00645   l = lambda;
00646   b = beta;
00647   mean = tgamma(1.0 + 1.0 / b) / l;
00648   var = tgamma(1.0 + 2.0 / b) / (l * l) - mean;
00649 }
00650 
00651 
00652 vec Weibull_RNG::operator()(int n)
00653 {
00654   vec vv(n);
00655 
00656   for (int i=0; i < n; i++)
00657     vv(i) = sample();
00658 
00659   return vv;
00660 }
00661 
00662 mat Weibull_RNG::operator()(int h, int w)
00663 {
00664   mat mm(h, w);
00665   int i, j;
00666 
00667   for (i = 0; i < h; i++)
00668     for (j = 0; j < w; j++)
00669       mm(i, j) = sample();
00670 
00671   return mm;
00672 }
00673 
00675 // Rayleigh_RNG
00677 
00678 Rayleigh_RNG::Rayleigh_RNG(double sigma)
00679 {
00680   setup(sigma);
00681 }
00682 
00683 vec Rayleigh_RNG::operator()(int n)
00684 {
00685   vec vv(n);
00686 
00687   for (int i=0; i < n; i++)
00688     vv(i) = sample();
00689 
00690   return vv;
00691 }
00692 
00693 mat Rayleigh_RNG::operator()(int h, int w)
00694 {
00695   mat mm(h, w);
00696   int i, j;
00697 
00698   for (i = 0; i < h; i++)
00699     for (j = 0; j < w; j++)
00700       mm(i, j) = sample();
00701 
00702   return mm;
00703 }
00704 
00706 // Rice_RNG
00708 
00709 Rice_RNG::Rice_RNG(double lambda, double beta)
00710 {
00711   setup(lambda, beta);
00712 }
00713 
00714 vec Rice_RNG::operator()(int n)
00715 {
00716   vec vv(n);
00717 
00718   for (int i=0; i < n; i++)
00719     vv(i) = sample();
00720 
00721   return vv;
00722 }
00723 
00724 mat Rice_RNG::operator()(int h, int w)
00725 {
00726   mat mm(h, w);
00727   int i, j;
00728 
00729   for (i = 0; i < h; i++)
00730     for (j = 0; j < w; j++)
00731       mm(i, j) = sample();
00732 
00733   return mm;
00734 }
00735 
00736 } // namespace itpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:30 for IT++ by Doxygen 1.7.4