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
Generated on Sat Jul 9 2011 15:21:30 for IT++ by Doxygen 1.7.4