Mercurial > octave
annotate src/DLD-FUNCTIONS/rand.cc @ 7533:ff52243af934
save state separately for each MT random number generator
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Tue, 26 Feb 2008 05:28:59 -0500 |
parents | 1c7b3e1fda19 |
children | 08125419efcb |
rev | line source |
---|---|
2928 | 1 /* |
2 | |
7017 | 3 Copyright (C) 1996, 1997, 1998, 1999, 2000, 2002, 2003, 2005, 2006, |
4 2007 John W. Eaton | |
2928 | 5 |
6 This file is part of Octave. | |
7 | |
8 Octave is free software; you can redistribute it and/or modify it | |
9 under the terms of the GNU General Public License as published by the | |
7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
11 option) any later version. | |
2928 | 12 |
13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
16 for more details. | |
17 | |
18 You should have received a copy of the GNU General Public License | |
7016 | 19 along with Octave; see the file COPYING. If not, see |
20 <http://www.gnu.org/licenses/>. | |
2928 | 21 |
22 */ | |
23 | |
24 #ifdef HAVE_CONFIG_H | |
25 #include <config.h> | |
26 #endif | |
27 | |
28 #include <ctime> | |
29 | |
30 #include <string> | |
31 | |
32 #include "f77-fcn.h" | |
33 #include "lo-mappers.h" | |
4307 | 34 #include "oct-rand.h" |
4153 | 35 #include "quit.h" |
2928 | 36 |
37 #include "defun-dld.h" | |
38 #include "error.h" | |
39 #include "gripes.h" | |
40 #include "oct-obj.h" | |
41 #include "unwind-prot.h" | |
42 #include "utils.h" | |
43 | |
6437 | 44 /* |
45 %!shared __random_statistical_tests__ | |
46 %! % Flag whether the statistical tests should be run in "make check" or not | |
47 %! __random_statistical_tests__ = 0; | |
48 */ | |
49 | |
4307 | 50 static octave_value |
5730 | 51 do_rand (const octave_value_list& args, int nargin, const char *fcn, |
52 bool additional_arg = false) | |
2928 | 53 { |
4307 | 54 octave_value retval; |
5730 | 55 NDArray a; |
56 int idx = 0; | |
57 dim_vector dims; | |
2928 | 58 |
5730 | 59 if (additional_arg) |
60 { | |
61 if (nargin == 0) | |
62 { | |
63 error ("%s: expecting at least one argument", fcn); | |
64 goto done; | |
65 } | |
66 else if (args(0).is_string()) | |
67 additional_arg = false; | |
68 else | |
69 { | |
70 a = args(0).array_value (); | |
71 if (error_state) | |
72 { | |
73 error ("%s: expecting scalar or matrix arguments", fcn); | |
74 goto done; | |
75 } | |
76 idx++; | |
77 nargin--; | |
78 } | |
79 } | |
2928 | 80 |
4543 | 81 switch (nargin) |
2928 | 82 { |
4543 | 83 case 0: |
84 { | |
5730 | 85 if (additional_arg) |
86 dims = a.dims (); | |
87 else | |
88 { | |
89 dims.resize (2); | |
4543 | 90 |
5730 | 91 dims(0) = 1; |
92 dims(1) = 1; | |
93 } | |
4543 | 94 goto gen_matrix; |
95 } | |
96 break; | |
2928 | 97 |
4543 | 98 case 1: |
99 { | |
5730 | 100 octave_value tmp = args(idx); |
4543 | 101 |
102 if (tmp.is_string ()) | |
103 { | |
104 std::string s_arg = tmp.string_value (); | |
2928 | 105 |
4543 | 106 if (s_arg == "dist") |
107 { | |
108 retval = octave_rand::distribution (); | |
109 } | |
5730 | 110 else if (s_arg == "seed") |
4543 | 111 { |
112 retval = octave_rand::seed (); | |
113 } | |
5730 | 114 else if (s_arg == "state" || s_arg == "twister") |
115 { | |
7533
ff52243af934
save state separately for each MT random number generator
John W. Eaton <jwe@octave.org>
parents:
7421
diff
changeset
|
116 retval = octave_rand::state (fcn); |
5730 | 117 } |
4543 | 118 else if (s_arg == "uniform") |
119 { | |
120 octave_rand::uniform_distribution (); | |
121 } | |
122 else if (s_arg == "normal") | |
123 { | |
124 octave_rand::normal_distribution (); | |
125 } | |
5730 | 126 else if (s_arg == "exponential") |
127 { | |
128 octave_rand::exponential_distribution (); | |
129 } | |
130 else if (s_arg == "poisson") | |
131 { | |
132 octave_rand::poisson_distribution (); | |
133 } | |
134 else if (s_arg == "gamma") | |
135 { | |
136 octave_rand::gamma_distribution (); | |
137 } | |
4543 | 138 else |
4664 | 139 error ("%s: unrecognized string argument", fcn); |
4543 | 140 } |
141 else if (tmp.is_scalar_type ()) | |
142 { | |
143 double dval = tmp.double_value (); | |
2928 | 144 |
4543 | 145 if (xisnan (dval)) |
146 { | |
4664 | 147 error ("%s: NaN is invalid a matrix dimension", fcn); |
4543 | 148 } |
149 else | |
150 { | |
151 dims.resize (2); | |
152 | |
5275 | 153 dims(0) = NINTbig (tmp.double_value ()); |
154 dims(1) = NINTbig (tmp.double_value ()); | |
2928 | 155 |
4543 | 156 if (! error_state) |
157 goto gen_matrix; | |
158 } | |
159 } | |
160 else if (tmp.is_range ()) | |
161 { | |
162 Range r = tmp.range_value (); | |
163 | |
164 if (r.all_elements_are_ints ()) | |
165 { | |
5275 | 166 octave_idx_type n = r.nelem (); |
4543 | 167 |
168 dims.resize (n); | |
169 | |
5275 | 170 octave_idx_type base = NINTbig (r.base ()); |
171 octave_idx_type incr = NINTbig (r.inc ()); | |
172 octave_idx_type lim = NINTbig (r.limit ()); | |
2928 | 173 |
4543 | 174 if (base < 0 || lim < 0) |
4664 | 175 error ("%s: all dimensions must be nonnegative", fcn); |
4543 | 176 else |
177 { | |
5275 | 178 for (octave_idx_type i = 0; i < n; i++) |
4543 | 179 { |
180 dims(i) = base; | |
181 base += incr; | |
182 } | |
2928 | 183 |
4543 | 184 goto gen_matrix; |
185 } | |
186 } | |
187 else | |
4664 | 188 error ("%s: expecting all elements of range to be integers", |
189 fcn); | |
4543 | 190 } |
191 else if (tmp.is_matrix_type ()) | |
192 { | |
193 Array<int> iv = tmp.int_vector_value (true); | |
194 | |
195 if (! error_state) | |
196 { | |
5275 | 197 octave_idx_type len = iv.length (); |
2928 | 198 |
4543 | 199 dims.resize (len); |
200 | |
5275 | 201 for (octave_idx_type i = 0; i < len; i++) |
4543 | 202 { |
5275 | 203 octave_idx_type elt = iv(i); |
4543 | 204 |
205 if (elt < 0) | |
206 { | |
4664 | 207 error ("%s: all dimensions must be nonnegative", fcn); |
4543 | 208 goto done; |
209 } | |
210 | |
211 dims(i) = iv(i); | |
212 } | |
2928 | 213 |
4543 | 214 goto gen_matrix; |
215 } | |
216 else | |
4664 | 217 error ("%s: expecting integer vector", fcn); |
4543 | 218 } |
219 else | |
220 { | |
221 gripe_wrong_type_arg ("rand", tmp); | |
222 return retval; | |
223 } | |
224 } | |
225 break; | |
226 | |
227 default: | |
228 { | |
5730 | 229 octave_value tmp = args(idx); |
4543 | 230 |
231 if (nargin == 2 && tmp.is_string ()) | |
232 { | |
5164 | 233 std::string ts = tmp.string_value (); |
234 | |
5730 | 235 if (ts == "seed") |
4543 | 236 { |
5782 | 237 if (args(idx+1).is_real_scalar ()) |
238 { | |
239 double d = args(idx+1).double_value (); | |
2928 | 240 |
5782 | 241 if (! error_state) |
242 octave_rand::seed (d); | |
243 } | |
244 else | |
245 error ("%s: seed must be a real scalar", fcn); | |
4543 | 246 } |
5730 | 247 else if (ts == "state" || ts == "twister") |
248 { | |
249 ColumnVector s = | |
250 ColumnVector (args(idx+1).vector_value(false, true)); | |
251 | |
252 if (! error_state) | |
7533
ff52243af934
save state separately for each MT random number generator
John W. Eaton <jwe@octave.org>
parents:
7421
diff
changeset
|
253 octave_rand::state (s, fcn); |
5730 | 254 } |
4543 | 255 else |
4664 | 256 error ("%s: unrecognized string argument", fcn); |
4543 | 257 } |
258 else | |
259 { | |
260 dims.resize (nargin); | |
261 | |
262 for (int i = 0; i < nargin; i++) | |
263 { | |
5760 | 264 dims(i) = args(idx+i).int_value (); |
4543 | 265 |
266 if (error_state) | |
267 { | |
4664 | 268 error ("%s: expecting integer arguments", fcn); |
4543 | 269 goto done; |
270 } | |
271 } | |
272 | |
273 goto gen_matrix; | |
274 } | |
275 } | |
276 break; | |
2928 | 277 } |
278 | |
4543 | 279 done: |
2928 | 280 |
281 return retval; | |
282 | |
283 gen_matrix: | |
284 | |
5355 | 285 dims.chop_trailing_singletons (); |
286 | |
5730 | 287 if (additional_arg) |
288 { | |
289 if (a.length() == 1) | |
290 return octave_rand::nd_array (dims, a(0)); | |
291 else | |
292 { | |
293 if (a.dims() != dims) | |
294 { | |
295 error ("%s: mismatch in argument size", fcn); | |
296 return retval; | |
297 } | |
298 octave_idx_type len = a.length (); | |
299 NDArray m (dims); | |
300 double *v = m.fortran_vec (); | |
301 for (octave_idx_type i = 0; i < len; i++) | |
302 v[i] = octave_rand::scalar (a(i)); | |
303 return m; | |
304 } | |
305 } | |
306 else | |
307 return octave_rand::nd_array (dims); | |
2928 | 308 } |
309 | |
4665 | 310 DEFUN_DLD (rand, args, , |
3369 | 311 "-*- texinfo -*-\n\ |
312 @deftypefn {Loadable Function} {} rand (@var{x})\n\ | |
313 @deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\ | |
5730 | 314 @deftypefnx {Loadable Function} {} rand (\"state\", @var{x})\n\ |
315 @deftypefnx {Loadable Function} {} rand (\"seed\", @var{x})\n\ | |
3369 | 316 Return a matrix with random elements uniformly distributed on the\n\ |
317 interval (0, 1). The arguments are handled the same as the arguments\n\ | |
5730 | 318 for @code{eye}.\n\ |
319 \n\ | |
320 You can query the state of the random number generator using the\n\ | |
3369 | 321 form\n\ |
2928 | 322 \n\ |
3369 | 323 @example\n\ |
5730 | 324 v = rand (\"state\")\n\ |
325 @end example\n\ | |
326 \n\ | |
327 This returns a column vector @var{v} of length 625. Later, you can\n\ | |
328 restore the random number generator to the state @var{v}\n\ | |
329 using the form\n\ | |
330 \n\ | |
331 @example\n\ | |
332 rand (\"state\", v)\n\ | |
3369 | 333 @end example\n\ |
334 \n\ | |
335 @noindent\n\ | |
5730 | 336 You may also initialize the state vector from an arbitrary vector of\n\ |
337 length <= 625 for @var{v}. This new state will be a hash based on the\n\ | |
5798 | 338 value of @var{v}, not @var{v} itself.\n\ |
5730 | 339 \n\ |
340 By default, the generator is initialized from @code{/dev/urandom} if it is\n\ | |
341 available, otherwise from cpu time, wall clock time and the current\n\ | |
342 fraction of a second.\n\ | |
343 \n\ | |
7096 | 344 To compute the psuedo-random sequence, @code{rand} uses the Mersenne\n\ |
345 Twister with a period of 2^19937-1 (See M. Matsumoto and T. Nishimura,\n\ | |
346 ``Mersenne Twister: A 623-dimensionally\n\ | |
5730 | 347 equidistributed uniform pseudorandom number generator'', ACM Trans. on\n\ |
7001 | 348 Modeling and Computer Simulation Vol. 8, No. 1, January pp.3-30 1998,\n\ |
7171 | 349 @url{http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html}).\n\ |
6547 | 350 Do @strong{not} use for cryptography without securely hashing\n\ |
351 several returned values together, otherwise the generator state\n\ | |
352 can be learned after reading 624 consecutive values.\n\ | |
5730 | 353 \n\ |
7096 | 354 Older versions of Octave used a different random number generator.\n\ |
355 The new generator is used by default\n\ | |
5730 | 356 as it is significantly faster than the old generator, and produces\n\ |
5798 | 357 random numbers with a significantly longer cycle time. However, in\n\ |
358 some circumstances it might be desirable to obtain the same random\n\ | |
5730 | 359 sequences as used by the old generators. To do this the keyword\n\ |
360 \"seed\" is used to specify that the old generators should be use,\n\ | |
361 as in\n\ | |
2928 | 362 \n\ |
3369 | 363 @example\n\ |
5730 | 364 rand (\"seed\", val)\n\ |
3369 | 365 @end example\n\ |
366 \n\ | |
5730 | 367 which sets the seed of the generator to @var{val}. The seed of the\n\ |
368 generator can be queried with\n\ | |
369 \n\ | |
370 @example\n\ | |
371 s = rand (\"seed\")\n\ | |
372 @end example\n\ | |
373 \n\ | |
374 However, it should be noted that querying the seed will not cause\n\ | |
375 @code{rand} to use the old generators, only setting the seed will.\n\ | |
376 To cause @code{rand} to once again use the new generators, the\n\ | |
377 keyword \"state\" should be used to reset the state of the @code{rand}.\n\ | |
5798 | 378 @seealso{randn, rande, randg, randp}\n\ |
3369 | 379 @end deftypefn") |
2928 | 380 { |
4307 | 381 octave_value retval; |
2928 | 382 |
383 int nargin = args.length (); | |
384 | |
4543 | 385 retval = do_rand (args, nargin, "rand"); |
2928 | 386 |
387 return retval; | |
388 } | |
389 | |
5730 | 390 /* |
391 %!test # 'state' can be a scalar | |
392 %! rand('state',12); x = rand(1,4); | |
393 %! rand('state',12); y = rand(1,4); | |
394 %! assert(x,y); | |
395 %!test # 'state' can be a vector | |
396 %! rand('state',[12,13]); x=rand(1,4); | |
397 %! rand('state',[12;13]); y=rand(1,4); | |
398 %! assert(x,y); | |
399 %!test # querying 'state' doesn't disturb sequence | |
400 %! rand('state',12); rand(1,2); x=rand(1,2); | |
401 %! rand('state',12); rand(1,2); | |
402 %! s=rand('state'); y=rand(1,2); | |
403 %! assert(x,y); | |
404 %! rand('state',s); z=rand(1,2); | |
405 %! assert(x,z); | |
406 %!test # 'seed' must be a scalar | |
407 %! rand('seed',12); x = rand(1,4); | |
408 %! rand('seed',12); y = rand(1,4); | |
409 %! assert(x,y); | |
410 %!error(rand('seed',[12,13])) | |
411 %!test # querying 'seed' returns a value which can be used later | |
412 %! s=rand('seed'); x=rand(1,2); | |
413 %! rand('seed',s); y=rand(1,2); | |
414 %! assert(x,y); | |
415 %!test # querying 'seed' doesn't disturb sequence | |
416 %! rand('seed',12); rand(1,2); x=rand(1,2); | |
417 %! rand('seed',12); rand(1,2); | |
418 %! s=rand('seed'); y=rand(1,2); | |
419 %! assert(x,y); | |
420 %! rand('seed',s); z=rand(1,2); | |
421 %! assert(x,z); | |
422 */ | |
423 | |
424 /* | |
425 %!test | |
6437 | 426 %! % Test fixed state |
427 %! rand("state",1); | |
6443 | 428 %! assert (rand(1,6), [0.1343642441124013 0.8474337369372327 0.763774618976614 0.2550690257394218 0.495435087091941 0.4494910647887382],1e-6); |
6437 | 429 %!test |
6443 | 430 %! % Test fixed seed |
6437 | 431 %! rand("seed",1); |
6443 | 432 %! assert (rand(1,6), [0.8668024251237512 0.9126510815694928 0.09366085007786751 0.1664607301354408 0.7408077004365623 0.7615650338120759],1e-6); |
5730 | 433 %!test |
6437 | 434 %! if (__random_statistical_tests__) |
435 %! % statistical tests may fail occasionally. | |
436 %! rand("state",12); | |
437 %! x = rand(100000,1); | |
438 %! assert(max(x)<1.); %*** Please report this!!! *** | |
439 %! assert(min(x)>0.); %*** Please report this!!! *** | |
440 %! assert(mean(x),0.5,0.0024); | |
441 %! assert(var(x),1/48,0.0632); | |
442 %! assert(skewness(x),0,0.012); | |
443 %! assert(kurtosis(x),-6/5,0.0094); | |
444 %! endif | |
445 %!test | |
446 %! if (__random_statistical_tests__) | |
447 %! % statistical tests may fail occasionally. | |
448 %! rand("seed",12); | |
449 %! x = rand(100000,1); | |
450 %! assert(max(x)<1.); %*** Please report this!!! *** | |
451 %! assert(min(x)>0.); %*** Please report this!!! *** | |
452 %! assert(mean(x),0.5,0.0024); | |
453 %! assert(var(x),1/48,0.0632); | |
454 %! assert(skewness(x),0,0.012); | |
455 %! assert(kurtosis(x),-6/5,0.0094); | |
456 %! endif | |
5730 | 457 */ |
458 | |
459 | |
4307 | 460 static std::string current_distribution = octave_rand::distribution (); |
461 | |
2928 | 462 static void |
463 reset_rand_generator (void *) | |
464 { | |
4307 | 465 octave_rand::distribution (current_distribution); |
2928 | 466 } |
467 | |
4665 | 468 DEFUN_DLD (randn, args, , |
3369 | 469 "-*- texinfo -*-\n\ |
470 @deftypefn {Loadable Function} {} randn (@var{x})\n\ | |
471 @deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\ | |
5730 | 472 @deftypefnx {Loadable Function} {} randn (\"state\", @var{x})\n\ |
473 @deftypefnx {Loadable Function} {} randn (\"seed\", @var{x})\n\ | |
474 Return a matrix with normally distributed random elements. The\n\ | |
475 arguments are handled the same as the arguments for @code{rand}.\n\ | |
3369 | 476 \n\ |
5730 | 477 By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ |
478 transform from a uniform to a normal distribution. (G. Marsaglia and\n\ | |
479 W.W. Tsang, 'Ziggurat method for generating random variables',\n\ | |
480 J. Statistical Software, vol 5, 2000,\n\ | |
481 @url{http://www.jstatsoft.org/v05/i08/})\n\ | |
2928 | 482 \n\ |
6547 | 483 @seealso{rand, rande, randg, randp}\n\ |
3369 | 484 @end deftypefn") |
2928 | 485 { |
4307 | 486 octave_value retval; |
2928 | 487 |
488 int nargin = args.length (); | |
489 | |
4543 | 490 unwind_protect::begin_frame ("randn"); |
2928 | 491 |
4543 | 492 // This relies on the fact that elements are popped from the unwind |
493 // stack in the reverse of the order they are pushed | |
494 // (i.e. current_distribution will be reset before calling | |
495 // reset_rand_generator()). | |
2928 | 496 |
4543 | 497 unwind_protect::add (reset_rand_generator, 0); |
498 unwind_protect_str (current_distribution); | |
2928 | 499 |
4543 | 500 current_distribution = "normal"; |
2928 | 501 |
4543 | 502 octave_rand::distribution (current_distribution); |
2928 | 503 |
4543 | 504 retval = do_rand (args, nargin, "randn"); |
2928 | 505 |
4543 | 506 unwind_protect::run_frame ("randn"); |
2928 | 507 |
508 return retval; | |
509 } | |
510 | |
511 /* | |
5730 | 512 %!test |
6437 | 513 %! % Test fixed state |
514 %! randn("state",1); | |
6443 | 515 %! assert (randn(1,6), [-2.666521678978671 -0.7381719971724564 1.507903992673601 0.6019427189162239 -0.450661261143348 -0.7054431351574116],1e-6); |
6437 | 516 %!test |
6443 | 517 %! % Test fixed seed |
6437 | 518 %! randn("seed",1); |
6443 | 519 %! assert (randn(1,6), [-1.039402365684509 -1.25938892364502 0.1968704611063004 0.3874166905879974 -0.5976632833480835 -0.6615074276924133],1e-6); |
5730 | 520 %!test |
6437 | 521 %! if (__random_statistical_tests__) |
522 %! % statistical tests may fail occasionally. | |
523 %! randn("state",12); | |
524 %! x = randn(100000,1); | |
525 %! assert(mean(x),0,0.01); | |
526 %! assert(var(x),1,0.02); | |
527 %! assert(skewness(x),0,0.02); | |
528 %! assert(kurtosis(x),0,0.04); | |
529 %! endif | |
530 %!test | |
531 %! if (__random_statistical_tests__) | |
532 %! % statistical tests may fail occasionally. | |
533 %! randn("seed",12); | |
534 %! x = randn(100000,1); | |
535 %! assert(mean(x),0,0.01); | |
536 %! assert(var(x),1,0.02); | |
537 %! assert(skewness(x),0,0.02); | |
538 %! assert(kurtosis(x),0,0.04); | |
539 %! endif | |
5730 | 540 */ |
541 | |
542 DEFUN_DLD (rande, args, , | |
543 "-*- texinfo -*-\n\ | |
544 @deftypefn {Loadable Function} {} rande (@var{x})\n\ | |
545 @deftypefnx {Loadable Function} {} rande (@var{n}, @var{m})\n\ | |
546 @deftypefnx {Loadable Function} {} rande (\"state\", @var{x})\n\ | |
547 @deftypefnx {Loadable Function} {} rande (\"seed\", @var{x})\n\ | |
548 Return a matrix with exponentially distributed random elements. The\n\ | |
549 arguments are handled the same as the arguments for @code{rand}.\n\ | |
550 \n\ | |
551 By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ | |
552 transform from a uniform to a exponential distribution. (G. Marsaglia and\n\ | |
553 W.W. Tsang, 'Ziggurat method for generating random variables',\n\ | |
554 J. Statistical Software, vol 5, 2000,\n\ | |
555 @url{http://www.jstatsoft.org/v05/i08/})\n\ | |
6547 | 556 @seealso{rand, randn, randg, randp}\n\ |
5730 | 557 @end deftypefn") |
558 { | |
559 octave_value retval; | |
560 | |
561 int nargin = args.length (); | |
562 | |
563 unwind_protect::begin_frame ("rande"); | |
564 | |
565 // This relies on the fact that elements are popped from the unwind | |
566 // stack in the reverse of the order they are pushed | |
567 // (i.e. current_distribution will be reset before calling | |
568 // reset_rand_generator()). | |
569 | |
570 unwind_protect::add (reset_rand_generator, 0); | |
571 unwind_protect_str (current_distribution); | |
572 | |
573 current_distribution = "exponential"; | |
574 | |
575 octave_rand::distribution (current_distribution); | |
576 | |
577 retval = do_rand (args, nargin, "rande"); | |
578 | |
579 unwind_protect::run_frame ("rande"); | |
580 | |
581 return retval; | |
582 } | |
583 | |
584 /* | |
585 %!test | |
6437 | 586 %! % Test fixed state |
587 %! rande("state",1); | |
6443 | 588 %! assert (rande(1,6), [3.602973885835625 0.1386190677555021 0.6743112889616958 0.4512830847258422 0.7255744741233175 0.3415969205292291],1e-6); |
6437 | 589 %!test |
6443 | 590 %! % Test fixed seed |
6437 | 591 %! rande("seed",1); |
6443 | 592 %! assert (rande(1,6), [0.06492075175653866 1.717980206012726 0.4816154008731246 0.5231300676241517 0.103910739364359 1.668931916356087],1e-6); |
5730 | 593 %!test |
6437 | 594 %! if (__random_statistical_tests__) |
595 %! % statistical tests may fail occasionally | |
596 %! rande("state",1); | |
597 %! x = rande(100000,1); | |
598 %! assert(min(x)>0); % *** Please report this!!! *** | |
599 %! assert(mean(x),1,0.01); | |
600 %! assert(var(x),1,0.03); | |
601 %! assert(skewness(x),2,0.06); | |
602 %! assert(kurtosis(x),6,0.7); | |
603 %! endif | |
604 %!test | |
605 %! if (__random_statistical_tests__) | |
606 %! % statistical tests may fail occasionally | |
607 %! rande("seed",1); | |
608 %! x = rande(100000,1); | |
609 %! assert(min(x)>0); % *** Please report this!!! *** | |
610 %! assert(mean(x),1,0.01); | |
611 %! assert(var(x),1,0.03); | |
612 %! assert(skewness(x),2,0.06); | |
613 %! assert(kurtosis(x),6,0.7); | |
614 %! endif | |
5730 | 615 */ |
616 | |
617 DEFUN_DLD (randg, args, , | |
618 "-*- texinfo -*-\n\ | |
619 @deftypefn {Loadable Function} {} randg (@var{a}, @var{x})\n\ | |
620 @deftypefnx {Loadable Function} {} randg (@var{a}, @var{n}, @var{m})\n\ | |
621 @deftypefnx {Loadable Function} {} randg (\"state\", @var{x})\n\ | |
622 @deftypefnx {Loadable Function} {} randg (\"seed\", @var{x})\n\ | |
623 Return a matrix with @code{gamma(@var{a},1)} distributed random elements.\n\ | |
624 The arguments are handled the same as the arguments for @code{rand},\n\ | |
625 except for the argument @var{a}.\n\ | |
626 \n\ | |
627 This can be used to generate many distributions:\n\ | |
628 \n\ | |
629 @table @asis\n\ | |
6547 | 630 @item @code{gamma (a, b)} for @code{a > -1}, @code{b > 0}\n\ |
5730 | 631 @example\n\ |
6547 | 632 r = b * randg (a)\n\ |
5730 | 633 @end example\n\ |
6547 | 634 @item @code{beta (a, b)} for @code{a > -1}, @code{b > -1}\n\ |
5730 | 635 @example\n\ |
6547 | 636 r1 = randg (a, 1)\n\ |
637 r = r1 / (r1 + randg (b, 1))\n\ | |
5730 | 638 @end example\n\ |
6547 | 639 @item @code{Erlang (a, n)}\n\ |
5730 | 640 @example\n\ |
6547 | 641 r = a * randg (n)\n\ |
5730 | 642 @end example\n\ |
6547 | 643 @item @code{chisq (df)} for @code{df > 0}\n\ |
5730 | 644 @example\n\ |
6547 | 645 r = 2 * randg (df / 2)\n\ |
5730 | 646 @end example\n\ |
647 @item @code{t(df)} for @code{0 < df < inf} (use randn if df is infinite)\n\ | |
648 @example\n\ | |
6547 | 649 r = randn () / sqrt (2 * randg (df / 2) / df)\n\ |
5730 | 650 @end example\n\ |
6547 | 651 @item @code{F (n1, n2)} for @code{0 < n1}, @code{0 < n2}\n\ |
5730 | 652 @example\n\ |
7096 | 653 @group\n\ |
654 ## r1 equals 1 if n1 is infinite\n\ | |
655 r1 = 2 * randg (n1 / 2) / n1\n\ | |
656 ## r2 equals 1 if n2 is infinite\n\ | |
657 r2 = 2 * randg (n2 / 2) / n2\n\ | |
5730 | 658 r = r1 / r2\n\n\ |
7096 | 659 @end group\n\ |
5730 | 660 @end example\n\ |
661 @item negative @code{binomial (n, p)} for @code{n > 0}, @code{0 < p <= 1}\n\ | |
662 @example\n\ | |
6547 | 663 r = randp ((1 - p) / p * randg (n))\n\ |
5730 | 664 @end example\n\ |
6547 | 665 @item non-central @code{chisq (df, L)}, for @code{df >= 0} and @code{L > 0}\n\ |
5730 | 666 (use chisq if @code{L = 0})\n\ |
667 @example\n\ | |
6547 | 668 r = randp (L / 2)\n\ |
669 r(r > 0) = 2 * randg (r(r > 0))\n\ | |
670 r(df > 0) += 2 * randg (df(df > 0)/2)\n\ | |
5730 | 671 @end example\n\ |
6547 | 672 @item @code{Dirichlet (a1, ..., ak)}\n\ |
5730 | 673 @example\n\ |
6547 | 674 r = (randg (a1), ..., randg (ak))\n\ |
675 r = r / sum (r)\n\ | |
5730 | 676 @end example\n\ |
677 @end table\n\ | |
6547 | 678 @seealso{rand, randn, rande, randp}\n\ |
5730 | 679 @end deftypefn") |
680 { | |
681 octave_value retval; | |
682 | |
683 int nargin = args.length (); | |
684 | |
685 if (nargin < 1) | |
686 error ("randg: insufficient arguments"); | |
687 else | |
688 { | |
689 unwind_protect::begin_frame ("randg"); | |
690 | |
691 // This relies on the fact that elements are popped from the unwind | |
692 // stack in the reverse of the order they are pushed | |
693 // (i.e. current_distribution will be reset before calling | |
694 // reset_rand_generator()). | |
695 | |
696 unwind_protect::add (reset_rand_generator, 0); | |
697 unwind_protect_str (current_distribution); | |
698 | |
699 current_distribution = "gamma"; | |
700 | |
701 octave_rand::distribution (current_distribution); | |
702 | |
703 retval = do_rand (args, nargin, "randg", true); | |
704 | |
705 unwind_protect::run_frame ("randg"); | |
706 } | |
707 | |
708 return retval; | |
709 } | |
710 | |
711 /* | |
712 %!test | |
6437 | 713 %! randg("state",12) |
714 %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report | |
715 | |
716 | |
717 %!test | |
718 %! % Test fixed state | |
719 %! randg("state",1); | |
6443 | 720 %! assert (randg(0.1,1,6), [0.0103951513331241 8.335671459898252e-05 0.00138691397249762 0.000587308416993855 0.495590518784736 2.3921917414795e-12],1e-6); |
6437 | 721 %!test |
722 %! % Test fixed state | |
723 %! randg("state",1); | |
6443 | 724 %! assert (randg(0.95,1,6), [3.099382433255327 0.3974529788871218 0.644367450750855 1.143261091802246 1.964111762696822 0.04011915547957939],1e-6); |
6437 | 725 %!test |
726 %! % Test fixed state | |
727 %! randg("state",1); | |
6443 | 728 %! assert (randg(1,1,6), [0.2273389379645993 1.288822625058359 0.2406335209340746 1.218869553370733 1.024649860162554 0.09631230343599533],1e-6); |
6437 | 729 %!test |
730 %! % Test fixed state | |
731 %! randg("state",1); | |
6443 | 732 %! assert (randg(10,1,6), [3.520369644331133 15.15369864472106 8.332112081991205 8.406211067432674 11.81193475187611 10.88792728177059],1e-5); |
6437 | 733 %!test |
734 %! % Test fixed state | |
735 %! randg("state",1); | |
6443 | 736 %! assert (randg(100,1,6), [75.34570255262264 115.4911985594699 95.23493031356388 95.48926019250911 106.2397448229803 103.4813150404118],1e-4); |
6437 | 737 %!test |
738 %! % Test fixed seed | |
739 %! randg("seed",1); | |
6443 | 740 %! assert (randg(0.1,1,6), [0.07144210487604141 0.460641473531723 0.4749028384685516 0.06823389977216721 0.000293838675133884 1.802567535340305e-12],1e-6); |
6437 | 741 %!test |
742 %! % Test fixed seed | |
743 %! randg("seed",1); | |
6443 | 744 %! assert (randg(0.95,1,6), [1.664905071258545 1.879976987838745 1.905677795410156 0.9948706030845642 0.5606933236122131 0.0766092911362648],1e-6); |
6437 | 745 %!test |
746 %! % Test fixed seed | |
747 %! randg("seed",1); | |
6443 | 748 %! assert (randg(1,1,6), [0.03512085229158401 0.6488978862762451 0.8114678859710693 0.1666885763406754 1.60791552066803 1.90356981754303],1e-6); |
6437 | 749 %!test |
750 %! % Test fixed seed | |
751 %! randg("seed",1); | |
6443 | 752 %! assert (randg(10,1,6), [6.566435813903809 10.11648464202881 10.73162078857422 7.747178077697754 6.278522491455078 6.240195751190186],1e-5); |
6437 | 753 %!test |
754 %! % Test fixed seed | |
755 %! randg("seed",1); | |
6443 | 756 %! assert (randg(100,1,6), [89.40208435058594 101.4734725952148 103.4020004272461 93.62763214111328 88.33104705810547 88.1871337890625],1e-4); |
6437 | 757 %!test |
758 %! if (__random_statistical_tests__) | |
759 %! % statistical tests may fail occasionally. | |
760 %! randg("state",12) | |
761 %! a=0.1; x = randg(a,100000,1); | |
762 %! assert(mean(x), a, 0.01); | |
763 %! assert(var(x), a, 0.01); | |
764 %! assert(skewness(x),2/sqrt(a), 1.); | |
765 %! assert(kurtosis(x),6/a, 50.); | |
766 %! endif | |
767 %!test | |
768 %! if (__random_statistical_tests__) | |
769 %! % statistical tests may fail occasionally. | |
770 %! randg("state",12) | |
771 %! a=0.95; x = randg(a,100000,1); | |
772 %! assert(mean(x), a, 0.01); | |
773 %! assert(var(x), a, 0.04); | |
774 %! assert(skewness(x),2/sqrt(a), 0.2); | |
775 %! assert(kurtosis(x),6/a, 2.); | |
776 %! endif | |
777 %!test | |
778 %! if (__random_statistical_tests__) | |
779 %! % statistical tests may fail occasionally. | |
780 %! randg("state",12) | |
781 %! a=1; x = randg(a,100000,1); | |
782 %! assert(mean(x),a, 0.01); | |
783 %! assert(var(x),a, 0.04); | |
784 %! assert(skewness(x),2/sqrt(a), 0.2); | |
785 %! assert(kurtosis(x),6/a, 2.); | |
786 %! endif | |
787 %!test | |
788 %! if (__random_statistical_tests__) | |
789 %! % statistical tests may fail occasionally. | |
790 %! randg("state",12) | |
791 %! a=10; x = randg(a,100000,1); | |
792 %! assert(mean(x), a, 0.1); | |
793 %! assert(var(x), a, 0.5); | |
794 %! assert(skewness(x),2/sqrt(a), 0.1); | |
795 %! assert(kurtosis(x),6/a, 0.5); | |
796 %! endif | |
797 %!test | |
798 %! if (__random_statistical_tests__) | |
799 %! % statistical tests may fail occasionally. | |
800 %! randg("state",12) | |
801 %! a=100; x = randg(a,100000,1); | |
802 %! assert(mean(x), a, 0.2); | |
803 %! assert(var(x), a, 2.); | |
804 %! assert(skewness(x),2/sqrt(a), 0.05); | |
805 %! assert(kurtosis(x),6/a, 0.2); | |
806 %! endif | |
807 %!test | |
808 %! randg("seed",12) | |
5730 | 809 %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report |
810 %!test | |
6437 | 811 %! if (__random_statistical_tests__) |
812 %! % statistical tests may fail occasionally. | |
813 %! randg("seed",12) | |
814 %! a=0.1; x = randg(a,100000,1); | |
815 %! assert(mean(x), a, 0.01); | |
816 %! assert(var(x), a, 0.01); | |
817 %! assert(skewness(x),2/sqrt(a), 1.); | |
818 %! assert(kurtosis(x),6/a, 50.); | |
819 %! endif | |
5730 | 820 %!test |
6437 | 821 %! if (__random_statistical_tests__) |
822 %! % statistical tests may fail occasionally. | |
823 %! randg("seed",12) | |
824 %! a=0.95; x = randg(a,100000,1); | |
825 %! assert(mean(x), a, 0.01); | |
826 %! assert(var(x), a, 0.04); | |
827 %! assert(skewness(x),2/sqrt(a), 0.2); | |
828 %! assert(kurtosis(x),6/a, 2.); | |
829 %! endif | |
5730 | 830 %!test |
6437 | 831 %! if (__random_statistical_tests__) |
832 %! % statistical tests may fail occasionally. | |
833 %! randg("seed",12) | |
834 %! a=1; x = randg(a,100000,1); | |
835 %! assert(mean(x),a, 0.01); | |
836 %! assert(var(x),a, 0.04); | |
837 %! assert(skewness(x),2/sqrt(a), 0.2); | |
838 %! assert(kurtosis(x),6/a, 2.); | |
839 %! endif | |
5730 | 840 %!test |
6437 | 841 %! if (__random_statistical_tests__) |
842 %! % statistical tests may fail occasionally. | |
843 %! randg("seed",12) | |
844 %! a=10; x = randg(a,100000,1); | |
845 %! assert(mean(x), a, 0.1); | |
846 %! assert(var(x), a, 0.5); | |
847 %! assert(skewness(x),2/sqrt(a), 0.1); | |
848 %! assert(kurtosis(x),6/a, 0.5); | |
849 %! endif | |
5730 | 850 %!test |
6437 | 851 %! if (__random_statistical_tests__) |
852 %! % statistical tests may fail occasionally. | |
853 %! randg("seed",12) | |
854 %! a=100; x = randg(a,100000,1); | |
855 %! assert(mean(x), a, 0.2); | |
856 %! assert(var(x), a, 2.); | |
857 %! assert(skewness(x),2/sqrt(a), 0.05); | |
858 %! assert(kurtosis(x),6/a, 0.2); | |
859 %! endif | |
5730 | 860 */ |
861 | |
862 | |
863 DEFUN_DLD (randp, args, , | |
864 "-*- texinfo -*-\n\ | |
865 @deftypefn {Loadable Function} {} randp (@var{l}, @var{x})\n\ | |
866 @deftypefnx {Loadable Function} {} randp (@var{l}, @var{n}, @var{m})\n\ | |
867 @deftypefnx {Loadable Function} {} randp (\"state\", @var{x})\n\ | |
868 @deftypefnx {Loadable Function} {} randp (\"seed\", @var{x})\n\ | |
869 Return a matrix with Poisson distributed random elements. The arguments\n\ | |
870 are handled the same as the arguments for @code{rand}, except for the\n\ | |
871 argument @var{l}.\n\ | |
872 \n\ | |
873 Five different algorithms are used depending on the range of @var{l}\n\ | |
874 and whether or not @var{l} is a scalar or a matrix.\n\ | |
875 \n\ | |
876 @table @asis\n\ | |
877 @item For scalar @var{l} <= 12, use direct method.\n\ | |
878 Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ | |
879 @item For scalar @var{l} > 12, use rejection method.[1]\n\ | |
880 Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ | |
881 @item For matrix @var{l} <= 10, use inversion method.[2]\n\ | |
882 Stadlober E., et al., WinRand source code, available via FTP.\n\ | |
883 @item For matrix @var{l} > 10, use patchwork rejection method.\n\ | |
884 Stadlober E., et al., WinRand source code, available via FTP, or\n\ | |
885 H. Zechner, 'Efficient sampling from continuous and discrete\n\ | |
886 unimodal distributions', Doctoral Dissertaion, 156pp., Technical\n\ | |
887 University Graz, Austria, 1994.\n\ | |
888 @item For @var{l} > 1e8, use normal approximation.\n\ | |
889 L. Montanet, et al., 'Review of Particle Properties', Physical Review\n\ | |
890 D 50 p1284, 1994\n\ | |
891 @end table\n\ | |
6547 | 892 @seealso{rand, randn, rande, randg}\n\ |
5730 | 893 @end deftypefn") |
894 { | |
895 octave_value retval; | |
896 | |
897 int nargin = args.length (); | |
898 | |
899 if (nargin < 1) | |
900 error ("randp: insufficient arguments"); | |
901 else | |
902 { | |
903 unwind_protect::begin_frame ("randp"); | |
904 | |
905 // This relies on the fact that elements are popped from the unwind | |
906 // stack in the reverse of the order they are pushed | |
907 // (i.e. current_distribution will be reset before calling | |
908 // reset_rand_generator()). | |
909 | |
910 unwind_protect::add (reset_rand_generator, 0); | |
911 unwind_protect_str (current_distribution); | |
912 | |
913 current_distribution = "poisson"; | |
914 | |
915 octave_rand::distribution (current_distribution); | |
916 | |
917 retval = do_rand (args, nargin, "randp", true); | |
918 | |
919 unwind_protect::run_frame ("randp"); | |
920 } | |
921 | |
922 return retval; | |
923 } | |
924 | |
925 /* | |
926 %!test | |
6437 | 927 %! randp("state",12) |
928 %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report | |
929 %!test | |
930 %! % Test fixed state | |
931 %! randp("state",1); | |
932 %! assert(randp(5,1,6),[5 5 3 7 7 3]) | |
933 %!test | |
934 %! % Test fixed state | |
935 %! randp("state",1); | |
936 %! assert(randp(15,1,6),[13 15 8 18 18 15]) | |
937 %!test | |
938 %! % Test fixed state | |
939 %! randp("state",1); | |
7421 | 940 %! assert(randp(1e9,1,6),[999915677 999976657 1000047684 1000019035 999985749 999977692],-1e-6) |
6437 | 941 %!test |
942 %! % Test fixed state | |
943 %! randp("seed",1); | |
6449 | 944 %! %%assert(randp(5,1,6),[8 2 3 6 6 8]) |
945 %! assert(randp(5,1,5),[8 2 3 6 6]) | |
6437 | 946 %!test |
947 %! % Test fixed state | |
948 %! randp("seed",1); | |
949 %! assert(randp(15,1,6),[15 16 12 10 10 12]) | |
950 %!test | |
951 %! % Test fixed state | |
952 %! randp("seed",1); | |
7421 | 953 %! assert(randp(1e9,1,6),[1000006208 1000012224 999981120 999963520 999963072 999981440],-1e-6) |
6437 | 954 %!test |
955 %! if (__random_statistical_tests__) | |
956 %! % statistical tests may fail occasionally. | |
957 %! randp("state",12) | |
958 %! for a=[5, 15, 1e9; 0.03, 0.03, -5e-3; 0.03, 0.03, 0.03] | |
959 %! x = randp(a(1),100000,1); | |
960 %! assert(min(x)>=0); % *** Please report this!!! *** | |
961 %! assert(mean(x),a(1),a(2)); | |
962 %! assert(var(x),a(1),0.02*a(1)); | |
963 %! assert(skewness(x),1/sqrt(a(1)),a(3)); | |
964 %! assert(kurtosis(x),1/a(1),3*a(3)); | |
965 %! endfor | |
966 %! endif | |
967 %!test | |
968 %! if (__random_statistical_tests__) | |
969 %! % statistical tests may fail occasionally. | |
970 %! randp("state",12) | |
971 %! for a=[5, 15, 1e9; 0.03, 0.03, -5e-3; 0.03, 0.03, 0.03] | |
972 %! x = randp(a(1)*ones(100000,1),100000,1); | |
973 %! assert(min(x)>=0); % *** Please report this!!! *** | |
974 %! assert(mean(x),a(1),a(2)); | |
975 %! assert(var(x),a(1),0.02*a(1)); | |
976 %! assert(skewness(x),1/sqrt(a(1)),a(3)); | |
977 %! assert(kurtosis(x),1/a(1),3*a(3)); | |
978 %! endfor | |
979 %! endif | |
980 %!test | |
981 %! randp("seed",12) | |
5730 | 982 %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report |
983 %!test | |
6449 | 984 %! if (__random_statistical_tests__) |
985 %! % statistical tests may fail occasionally. | |
986 %! randp("seed",12) | |
987 %! for a=[5, 15, 1e9; 0.03, 0.03, -5e-3; 0.03, 0.03, 0.03] | |
988 %! x = randp(a(1),100000,1); | |
989 %! assert(min(x)>=0); % *** Please report this!!! *** | |
990 %! assert(mean(x),a(1),a(2)); | |
991 %! assert(var(x),a(1),0.02*a(1)); | |
992 %! assert(skewness(x),1/sqrt(a(1)),a(3)); | |
993 %! assert(kurtosis(x),1/a(1),3*a(3)); | |
994 %! endfor | |
995 %! endif | |
5730 | 996 %!test |
6449 | 997 %! if (__random_statistical_tests__) |
998 %! % statistical tests may fail occasionally. | |
999 %! randp("seed",12) | |
1000 %! for a=[5, 15, 1e9; 0.03, 0.03, -5e-3; 0.03, 0.03, 0.03] | |
1001 %! x = randp(a(1)*ones(100000,1),100000,1); | |
1002 %! assert(min(x)>=0); % *** Please report this!!! *** | |
1003 %! assert(mean(x),a(1),a(2)); | |
1004 %! assert(var(x),a(1),0.02*a(1)); | |
1005 %! assert(skewness(x),1/sqrt(a(1)),a(3)); | |
1006 %! assert(kurtosis(x),1/a(1),3*a(3)); | |
1007 %! endfor | |
1008 %! endif | |
5730 | 1009 */ |
1010 | |
1011 /* | |
2928 | 1012 ;;; Local Variables: *** |
1013 ;;; mode: C++ *** | |
1014 ;;; End: *** | |
1015 */ |