Mercurial > octave
comparison src/DLD-FUNCTIONS/rand.cc @ 5730:109fdf7b3dcb
[project @ 2006-04-03 19:18:26 by jwe]
author | jwe |
---|---|
date | Mon, 03 Apr 2006 19:18:26 +0000 |
parents | cf44c749ba52 |
children | 8d7162924bd3 |
comparison
equal
deleted
inserted
replaced
5729:e065f7c18bdc | 5730:109fdf7b3dcb |
---|---|
40 #include "oct-obj.h" | 40 #include "oct-obj.h" |
41 #include "unwind-prot.h" | 41 #include "unwind-prot.h" |
42 #include "utils.h" | 42 #include "utils.h" |
43 | 43 |
44 static octave_value | 44 static octave_value |
45 do_rand (const octave_value_list& args, int nargin, const char *fcn) | 45 do_rand (const octave_value_list& args, int nargin, const char *fcn, |
46 bool additional_arg = false) | |
46 { | 47 { |
47 octave_value retval; | 48 octave_value retval; |
48 | 49 NDArray a; |
50 int idx = 0; | |
49 dim_vector dims; | 51 dim_vector dims; |
52 | |
53 if (additional_arg) | |
54 { | |
55 if (nargin == 0) | |
56 { | |
57 error ("%s: expecting at least one argument", fcn); | |
58 goto done; | |
59 } | |
60 else if (args(0).is_string()) | |
61 additional_arg = false; | |
62 else | |
63 { | |
64 a = args(0).array_value (); | |
65 if (error_state) | |
66 { | |
67 error ("%s: expecting scalar or matrix arguments", fcn); | |
68 goto done; | |
69 } | |
70 idx++; | |
71 nargin--; | |
72 } | |
73 } | |
50 | 74 |
51 switch (nargin) | 75 switch (nargin) |
52 { | 76 { |
53 case 0: | 77 case 0: |
54 { | 78 { |
55 dims.resize (2); | 79 if (additional_arg) |
56 | 80 dims = a.dims (); |
57 dims(0) = 1; | 81 else |
58 dims(1) = 1; | 82 { |
59 | 83 dims.resize (2); |
84 | |
85 dims(0) = 1; | |
86 dims(1) = 1; | |
87 } | |
60 goto gen_matrix; | 88 goto gen_matrix; |
61 } | 89 } |
62 break; | 90 break; |
63 | 91 |
64 case 1: | 92 case 1: |
65 { | 93 { |
66 octave_value tmp = args(0); | 94 octave_value tmp = args(idx); |
67 | 95 |
68 if (tmp.is_string ()) | 96 if (tmp.is_string ()) |
69 { | 97 { |
70 std::string s_arg = tmp.string_value (); | 98 std::string s_arg = tmp.string_value (); |
71 | 99 |
72 if (s_arg == "dist") | 100 if (s_arg == "dist") |
73 { | 101 { |
74 retval = octave_rand::distribution (); | 102 retval = octave_rand::distribution (); |
75 } | 103 } |
76 else if (s_arg == "seed" || s_arg == "state") | 104 else if (s_arg == "seed") |
77 { | 105 { |
78 retval = octave_rand::seed (); | 106 retval = octave_rand::seed (); |
79 } | 107 } |
108 else if (s_arg == "state" || s_arg == "twister") | |
109 { | |
110 retval = octave_rand::state (); | |
111 } | |
80 else if (s_arg == "uniform") | 112 else if (s_arg == "uniform") |
81 { | 113 { |
82 octave_rand::uniform_distribution (); | 114 octave_rand::uniform_distribution (); |
83 } | 115 } |
84 else if (s_arg == "normal") | 116 else if (s_arg == "normal") |
85 { | 117 { |
86 octave_rand::normal_distribution (); | 118 octave_rand::normal_distribution (); |
119 } | |
120 else if (s_arg == "exponential") | |
121 { | |
122 octave_rand::exponential_distribution (); | |
123 } | |
124 else if (s_arg == "poisson") | |
125 { | |
126 octave_rand::poisson_distribution (); | |
127 } | |
128 else if (s_arg == "gamma") | |
129 { | |
130 octave_rand::gamma_distribution (); | |
87 } | 131 } |
88 else | 132 else |
89 error ("%s: unrecognized string argument", fcn); | 133 error ("%s: unrecognized string argument", fcn); |
90 } | 134 } |
91 else if (tmp.is_scalar_type ()) | 135 else if (tmp.is_scalar_type ()) |
174 } | 218 } |
175 break; | 219 break; |
176 | 220 |
177 default: | 221 default: |
178 { | 222 { |
179 octave_value tmp = args(0); | 223 octave_value tmp = args(idx); |
180 | 224 |
181 if (nargin == 2 && tmp.is_string ()) | 225 if (nargin == 2 && tmp.is_string ()) |
182 { | 226 { |
183 std::string ts = tmp.string_value (); | 227 std::string ts = tmp.string_value (); |
184 | 228 |
185 if (ts == "seed" || ts == "state") | 229 if (ts == "seed") |
186 { | 230 { |
187 double d = args(1).double_value (); | 231 double d = args(idx+1).double_value (); |
188 | 232 |
189 if (! error_state) | 233 if (! error_state) |
190 octave_rand::seed (d); | 234 octave_rand::seed (d); |
235 } | |
236 else if (ts == "state" || ts == "twister") | |
237 { | |
238 ColumnVector s = | |
239 ColumnVector (args(idx+1).vector_value(false, true)); | |
240 | |
241 if (! error_state) | |
242 octave_rand::state (s); | |
191 } | 243 } |
192 else | 244 else |
193 error ("%s: unrecognized string argument", fcn); | 245 error ("%s: unrecognized string argument", fcn); |
194 } | 246 } |
195 else | 247 else |
196 { | 248 { |
197 dims.resize (nargin); | 249 dims.resize (nargin); |
198 | 250 |
199 for (int i = 0; i < nargin; i++) | 251 for (int i = 0; i < nargin; i++) |
200 { | 252 { |
201 dims(i) = (octave_idx_type)args(i).int_value (); | 253 dims(i) = (octave_idx_type)args(idx+i).int_value (); |
202 | 254 |
203 if (error_state) | 255 if (error_state) |
204 { | 256 { |
205 error ("%s: expecting integer arguments", fcn); | 257 error ("%s: expecting integer arguments", fcn); |
206 goto done; | 258 goto done; |
219 | 271 |
220 gen_matrix: | 272 gen_matrix: |
221 | 273 |
222 dims.chop_trailing_singletons (); | 274 dims.chop_trailing_singletons (); |
223 | 275 |
224 return octave_rand::nd_array (dims); | 276 if (additional_arg) |
277 { | |
278 if (a.length() == 1) | |
279 return octave_rand::nd_array (dims, a(0)); | |
280 else | |
281 { | |
282 if (a.dims() != dims) | |
283 { | |
284 error ("%s: mismatch in argument size", fcn); | |
285 return retval; | |
286 } | |
287 octave_idx_type len = a.length (); | |
288 NDArray m (dims); | |
289 double *v = m.fortran_vec (); | |
290 for (octave_idx_type i = 0; i < len; i++) | |
291 v[i] = octave_rand::scalar (a(i)); | |
292 return m; | |
293 } | |
294 } | |
295 else | |
296 return octave_rand::nd_array (dims); | |
225 } | 297 } |
226 | 298 |
227 DEFUN_DLD (rand, args, , | 299 DEFUN_DLD (rand, args, , |
228 "-*- texinfo -*-\n\ | 300 "-*- texinfo -*-\n\ |
229 @deftypefn {Loadable Function} {} rand (@var{x})\n\ | 301 @deftypefn {Loadable Function} {} rand (@var{x})\n\ |
230 @deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\ | 302 @deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\ |
231 @deftypefnx {Loadable Function} {} rand (@code{\"seed\"}, @var{x})\n\ | 303 @deftypefnx {Loadable Function} {} rand (\"state\", @var{x})\n\ |
304 @deftypefnx {Loadable Function} {} rand (\"seed\", @var{x})\n\ | |
232 Return a matrix with random elements uniformly distributed on the\n\ | 305 Return a matrix with random elements uniformly distributed on the\n\ |
233 interval (0, 1). The arguments are handled the same as the arguments\n\ | 306 interval (0, 1). The arguments are handled the same as the arguments\n\ |
234 for @code{eye}. In\n\ | 307 for @code{eye}.\n\ |
235 addition, you can set the seed for the random number generator using the\n\ | 308 \n\ |
309 You can query the state of the random number generator using the\n\ | |
236 form\n\ | 310 form\n\ |
237 \n\ | 311 \n\ |
238 @example\n\ | 312 @example\n\ |
239 rand (\"seed\", @var{x})\n\ | 313 v = rand (\"state\")\n\ |
314 @end example\n\ | |
315 \n\ | |
316 This returns a column vector @var{v} of length 625. Later, you can\n\ | |
317 restore the random number generator to the state @var{v}\n\ | |
318 using the form\n\ | |
319 \n\ | |
320 @example\n\ | |
321 rand (\"state\", v)\n\ | |
240 @end example\n\ | 322 @end example\n\ |
241 \n\ | 323 \n\ |
242 @noindent\n\ | 324 @noindent\n\ |
243 where @var{x} is a scalar value. If called as\n\ | 325 You may also initialize the state vector from an arbitrary vector of\n\ |
244 \n\ | 326 length <= 625 for @var{v}. This new state will be a hash based on the\n\ |
245 @example\n\ | 327 the value of @var{v}, not @var{v} itself.\n\ |
246 rand (\"seed\")\n\ | 328 \n\ |
247 @end example\n\ | 329 By default, the generator is initialized from @code{/dev/urandom} if it is\n\ |
248 \n\ | 330 available, otherwise from cpu time, wall clock time and the current\n\ |
249 @noindent\n\ | 331 fraction of a second.\n\ |
250 @code{rand} returns the current value of the seed.\n\ | 332 \n\ |
333 @code{rand} uses the Mersenne Twister with a period of 2^19937-1\n\ | |
334 (See M. Matsumoto and T. Nishimura, ``Mersenne Twister: A 623-dimensionally\n\ | |
335 equidistributed uniform pseudorandom number generator'', ACM Trans. on\n\ | |
336 Modeling and Computer Simulation Vol. 8, No. 1, Januray pp.3-30 1998,\n\ | |
337 @url{http://www.math.keio.ac.jp/~matumoto/emt.html}).\n\ | |
338 Do NOT use for CRYPTOGRAPHY without securely hashing several returned\n\ | |
339 values together, otherwise the generator state can be learned after\n\ | |
340 reading 624 consecutive values.\n\ | |
341 \n\ | |
342 @code{rand} includes a second random number generator, that was the\n\ | |
343 previous generator used in octave. The new generator is used by default\n\ | |
344 as it is significantly faster than the old generator, and produces\n\ | |
345 random numebrs with a significantly longer cycle time. However, in\n\ | |
346 some circumstances it might be desireable to obtain the same random\n\ | |
347 sequences as used by the old generators. To do this the keyword\n\ | |
348 \"seed\" is used to specify that the old generators should be use,\n\ | |
349 as in\n\ | |
350 \n\ | |
351 @example\n\ | |
352 rand (\"seed\", val)\n\ | |
353 @end example\n\ | |
354 \n\ | |
355 which sets the seed of the generator to @var{val}. The seed of the\n\ | |
356 generator can be queried with\n\ | |
357 \n\ | |
358 @example\n\ | |
359 s = rand (\"seed\")\n\ | |
360 @end example\n\ | |
361 \n\ | |
362 However, it should be noted that querying the seed will not cause\n\ | |
363 @code{rand} to use the old generators, only setting the seed will.\n\ | |
364 To cause @code{rand} to once again use the new generators, the\n\ | |
365 keyword \"state\" should be used to reset the state of the @code{rand}.\n\ | |
366 @seealso{randn,rande,randg,randp}\n\ | |
251 @end deftypefn") | 367 @end deftypefn") |
252 { | 368 { |
253 octave_value retval; | 369 octave_value retval; |
254 | 370 |
255 int nargin = args.length (); | 371 int nargin = args.length (); |
256 | 372 |
257 retval = do_rand (args, nargin, "rand"); | 373 retval = do_rand (args, nargin, "rand"); |
258 | 374 |
259 return retval; | 375 return retval; |
260 } | 376 } |
377 | |
378 /* | |
379 %!test # 'state' can be a scalar | |
380 %! rand('state',12); x = rand(1,4); | |
381 %! rand('state',12); y = rand(1,4); | |
382 %! assert(x,y); | |
383 %!test # 'state' can be a vector | |
384 %! rand('state',[12,13]); x=rand(1,4); | |
385 %! rand('state',[12;13]); y=rand(1,4); | |
386 %! assert(x,y); | |
387 %!test # querying 'state' doesn't disturb sequence | |
388 %! rand('state',12); rand(1,2); x=rand(1,2); | |
389 %! rand('state',12); rand(1,2); | |
390 %! s=rand('state'); y=rand(1,2); | |
391 %! assert(x,y); | |
392 %! rand('state',s); z=rand(1,2); | |
393 %! assert(x,z); | |
394 %!test # 'seed' must be a scalar | |
395 %! rand('seed',12); x = rand(1,4); | |
396 %! rand('seed',12); y = rand(1,4); | |
397 %! assert(x,y); | |
398 %!error(rand('seed',[12,13])) | |
399 %!test # querying 'seed' returns a value which can be used later | |
400 %! s=rand('seed'); x=rand(1,2); | |
401 %! rand('seed',s); y=rand(1,2); | |
402 %! assert(x,y); | |
403 %!test # querying 'seed' doesn't disturb sequence | |
404 %! rand('seed',12); rand(1,2); x=rand(1,2); | |
405 %! rand('seed',12); rand(1,2); | |
406 %! s=rand('seed'); y=rand(1,2); | |
407 %! assert(x,y); | |
408 %! rand('seed',s); z=rand(1,2); | |
409 %! assert(x,z); | |
410 */ | |
411 | |
412 /* | |
413 %!test | |
414 %! % statistical tests may fail occasionally. | |
415 %! rand("state",12); | |
416 %! x = rand(100000,1); | |
417 %! assert(max(x)<1.); %*** Please report this!!! *** | |
418 %! assert(min(x)>0.); %*** Please report this!!! *** | |
419 %! assert(mean(x),0.5,0.0024); | |
420 %! assert(var(x),1/48,0.0632); | |
421 %! assert(skewness(x),0,0.012); | |
422 %! assert(kurtosis(x),-6/5,0.0094); | |
423 %!test | |
424 %! % statistical tests may fail occasionally. | |
425 %! rand("seed",12); | |
426 %! x = rand(100000,1); | |
427 %! assert(max(x)<1.); %*** Please report this!!! *** | |
428 %! assert(min(x)>0.); %*** Please report this!!! *** | |
429 %! assert(mean(x),0.5,0.0024); | |
430 %! assert(var(x),1/48,0.0632); | |
431 %! assert(skewness(x),0,0.012); | |
432 %! assert(kurtosis(x),-6/5,0.0094); | |
433 */ | |
434 | |
261 | 435 |
262 static std::string current_distribution = octave_rand::distribution (); | 436 static std::string current_distribution = octave_rand::distribution (); |
263 | 437 |
264 static void | 438 static void |
265 reset_rand_generator (void *) | 439 reset_rand_generator (void *) |
269 | 443 |
270 DEFUN_DLD (randn, args, , | 444 DEFUN_DLD (randn, args, , |
271 "-*- texinfo -*-\n\ | 445 "-*- texinfo -*-\n\ |
272 @deftypefn {Loadable Function} {} randn (@var{x})\n\ | 446 @deftypefn {Loadable Function} {} randn (@var{x})\n\ |
273 @deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\ | 447 @deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\ |
274 @deftypefnx {Loadable Function} {} randn (@code{\"seed\"}, @var{x})\n\ | 448 @deftypefnx {Loadable Function} {} randn (\"state\", @var{x})\n\ |
275 Return a matrix with normally distributed random elements. The\n\ | 449 @deftypefnx {Loadable Function} {} randn (\"seed\", @var{x})\n\ |
276 arguments are handled the same as the arguments for @code{eye}. In\n\ | 450 Return a matrix with normally distributed random elements. The\n\ |
277 addition, you can set the seed for the random number generator using the\n\ | 451 arguments are handled the same as the arguments for @code{rand}.\n\ |
278 form\n\ | 452 \n\ |
279 \n\ | 453 By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ |
280 @example\n\ | 454 transform from a uniform to a normal distribution. (G. Marsaglia and\n\ |
281 randn (\"seed\", @var{x})\n\ | 455 W.W. Tsang, 'Ziggurat method for generating random variables',\n\ |
282 @end example\n\ | 456 J. Statistical Software, vol 5, 2000,\n\ |
283 \n\ | 457 @url{http://www.jstatsoft.org/v05/i08/})\n\ |
284 @noindent\n\ | 458 \n\ |
285 where @var{x} is a scalar value. If called as\n\ | 459 @seealso{rand,rande,randg,randp}\n\ |
286 \n\ | |
287 @example\n\ | |
288 randn (\"seed\")\n\ | |
289 @end example\n\ | |
290 \n\ | |
291 @noindent\n\ | |
292 @code{randn} returns the current value of the seed.\n\ | |
293 @end deftypefn") | 460 @end deftypefn") |
294 { | 461 { |
295 octave_value retval; | 462 octave_value retval; |
296 | 463 |
297 int nargin = args.length (); | 464 int nargin = args.length (); |
314 | 481 |
315 unwind_protect::run_frame ("randn"); | 482 unwind_protect::run_frame ("randn"); |
316 | 483 |
317 return retval; | 484 return retval; |
318 } | 485 } |
486 | |
487 /* | |
488 %!test | |
489 %! % statistical tests may fail occasionally. | |
490 %! rand("state",12); | |
491 %! x = randn(100000,1); | |
492 %! assert(mean(x),0,0.01); | |
493 %! assert(var(x),1,0.02); | |
494 %! assert(skewness(x),0,0.02); | |
495 %! assert(kurtosis(x),0,0.04); | |
496 %!test | |
497 %! % statistical tests may fail occasionally. | |
498 %! rand("seed",12); | |
499 %! x = randn(100000,1); | |
500 %! assert(mean(x),0,0.01); | |
501 %! assert(var(x),1,0.02); | |
502 %! assert(skewness(x),0,0.02); | |
503 %! assert(kurtosis(x),0,0.04); | |
504 */ | |
505 | |
506 DEFUN_DLD (rande, args, , | |
507 "-*- texinfo -*-\n\ | |
508 @deftypefn {Loadable Function} {} rande (@var{x})\n\ | |
509 @deftypefnx {Loadable Function} {} rande (@var{n}, @var{m})\n\ | |
510 @deftypefnx {Loadable Function} {} rande (\"state\", @var{x})\n\ | |
511 @deftypefnx {Loadable Function} {} rande (\"seed\", @var{x})\n\ | |
512 Return a matrix with exponentially distributed random elements. The\n\ | |
513 arguments are handled the same as the arguments for @code{rand}.\n\ | |
514 \n\ | |
515 By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ | |
516 transform from a uniform to a exponential distribution. (G. Marsaglia and\n\ | |
517 W.W. Tsang, 'Ziggurat method for generating random variables',\n\ | |
518 J. Statistical Software, vol 5, 2000,\n\ | |
519 @url{http://www.jstatsoft.org/v05/i08/})\n\ | |
520 @seealso{rand,randn,randg,randp}\n\ | |
521 @end deftypefn") | |
522 { | |
523 octave_value retval; | |
524 | |
525 int nargin = args.length (); | |
526 | |
527 unwind_protect::begin_frame ("rande"); | |
528 | |
529 // This relies on the fact that elements are popped from the unwind | |
530 // stack in the reverse of the order they are pushed | |
531 // (i.e. current_distribution will be reset before calling | |
532 // reset_rand_generator()). | |
533 | |
534 unwind_protect::add (reset_rand_generator, 0); | |
535 unwind_protect_str (current_distribution); | |
536 | |
537 current_distribution = "exponential"; | |
538 | |
539 octave_rand::distribution (current_distribution); | |
540 | |
541 retval = do_rand (args, nargin, "rande"); | |
542 | |
543 unwind_protect::run_frame ("rande"); | |
544 | |
545 return retval; | |
546 } | |
547 | |
548 /* | |
549 %!test | |
550 %! % statistical tests may fail occasionally | |
551 %! rand("state",12); | |
552 %! x = rande(100000,1); | |
553 %! assert(min(x)>0); % *** Please report this!!! *** | |
554 %! assert(mean(x),1,0.01); | |
555 %! assert(var(x),1,0.03); | |
556 %! assert(skewness(x),2,0.06); | |
557 %! assert(kurtosis(x),6,0.7); | |
558 %!test | |
559 %! % statistical tests may fail occasionally | |
560 %! rand("seed",12); | |
561 %! x = rande(100000,1); | |
562 %! assert(min(x)>0); % *** Please report this!!! *** | |
563 %! assert(mean(x),1,0.01); | |
564 %! assert(var(x),1,0.03); | |
565 %! assert(skewness(x),2,0.06); | |
566 %! assert(kurtosis(x),6,0.7); | |
567 */ | |
568 | |
569 DEFUN_DLD (randg, args, , | |
570 "-*- texinfo -*-\n\ | |
571 @deftypefn {Loadable Function} {} randg (@var{a}, @var{x})\n\ | |
572 @deftypefnx {Loadable Function} {} randg (@var{a}, @var{n}, @var{m})\n\ | |
573 @deftypefnx {Loadable Function} {} randg (\"state\", @var{x})\n\ | |
574 @deftypefnx {Loadable Function} {} randg (\"seed\", @var{x})\n\ | |
575 Return a matrix with @code{gamma(@var{a},1)} distributed random elements.\n\ | |
576 The arguments are handled the same as the arguments for @code{rand},\n\ | |
577 except for the argument @var{a}.\n\ | |
578 \n\ | |
579 This can be used to generate many distributions:\n\ | |
580 \n\ | |
581 @table @asis\n\ | |
582 @item @code{gamma (a,b)} for @code{a > -1}, @code{b > 0}\n\ | |
583 @example\n\ | |
584 r = b*randg(a)\n\ | |
585 @end example\n\ | |
586 @item @code{beta(a,b)} for @code{a > -1}, @code{b > -1}\n\ | |
587 @example\n\ | |
588 r1 = randg(a,1)\n\ | |
589 r = r1 / (r1 + randg(b,1))\n\ | |
590 @end example\n\ | |
591 @item @code{Erlang(a, n)}\n\ | |
592 @example\n\ | |
593 r = a*randg(n)\n\ | |
594 @end example\n\ | |
595 @item @code{chisq(df)} for @code{df > 0}\n\ | |
596 @example\n\ | |
597 r = 2*randg(df/2)\n\ | |
598 @end example\n\ | |
599 @item @code{t(df)} for @code{0 < df < inf} (use randn if df is infinite)\n\ | |
600 @example\n\ | |
601 r = randn() / sqrt(2*randg(df/2)/df)\n\ | |
602 @end example\n\ | |
603 @item @code{F(n1,n2)} for @code{0 < n1}, @code{0 < n2}\n\ | |
604 @example\n\ | |
605 r1 = 2*randg(n1/2)/n1 or 1 if n1 is infinite\n\ | |
606 r2 = 2*randg(n2/2)/n2 or 1 if n2 is infinite\n\ | |
607 r = r1 / r2\n\n\ | |
608 @end example\n\ | |
609 @item negative @code{binomial (n, p)} for @code{n > 0}, @code{0 < p <= 1}\n\ | |
610 @example\n\ | |
611 r = randp((1-p)/p * randg(n))\n\ | |
612 @end example\n\ | |
613 @item non-central @code{chisq(df,L)}, for @code{df >= 0} and @code{L > 0}\n\ | |
614 (use chisq if @code{L = 0})\n\ | |
615 @example\n\ | |
616 r = randp(L/2)\n\ | |
617 r(r > 0) = 2*randg(r(r > 0))\n\ | |
618 r(df > 0) += 2*randg(df(df > 0)/2)\n\ | |
619 @end example\n\ | |
620 @item @code{Dirichlet(a1,...,ak)}\n\ | |
621 @example\n\ | |
622 r = (randg(a1),...,randg(ak))\n\ | |
623 r = r / sum(r)\n\ | |
624 @end example\n\ | |
625 @end table\n\ | |
626 @seealso{rand,randn,rande,randp}\n\ | |
627 @end deftypefn") | |
628 { | |
629 octave_value retval; | |
630 | |
631 int nargin = args.length (); | |
632 | |
633 if (nargin < 1) | |
634 error ("randg: insufficient arguments"); | |
635 else | |
636 { | |
637 unwind_protect::begin_frame ("randg"); | |
638 | |
639 // This relies on the fact that elements are popped from the unwind | |
640 // stack in the reverse of the order they are pushed | |
641 // (i.e. current_distribution will be reset before calling | |
642 // reset_rand_generator()). | |
643 | |
644 unwind_protect::add (reset_rand_generator, 0); | |
645 unwind_protect_str (current_distribution); | |
646 | |
647 current_distribution = "gamma"; | |
648 | |
649 octave_rand::distribution (current_distribution); | |
650 | |
651 retval = do_rand (args, nargin, "randg", true); | |
652 | |
653 unwind_protect::run_frame ("randg"); | |
654 } | |
655 | |
656 return retval; | |
657 } | |
658 | |
659 /* | |
660 %!test | |
661 %! rand("state",12) | |
662 %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report | |
663 %!test | |
664 %! % statistical tests may fail occasionally. | |
665 %! a=0.1; x = randg(a,100000,1); | |
666 %! assert(mean(x), a, 0.01); | |
667 %! assert(var(x), a, 0.01); | |
668 %! assert(skewness(x),2/sqrt(a), 1.); | |
669 %! assert(kurtosis(x),6/a, 50.); | |
670 %!test | |
671 %! % statistical tests may fail occasionally. | |
672 %! a=0.95; x = randg(a,100000,1); | |
673 %! assert(mean(x), a, 0.01); | |
674 %! assert(var(x), a, 0.04); | |
675 %! assert(skewness(x),2/sqrt(a), 0.2); | |
676 %! assert(kurtosis(x),6/a, 2.); | |
677 %!test | |
678 %! % statistical tests may fail occasionally. | |
679 %! a=1; x = randg(a,100000,1); | |
680 %! assert(mean(x),a, 0.01); | |
681 %! assert(var(x),a, 0.04); | |
682 %! assert(skewness(x),2/sqrt(a), 0.2); | |
683 %! assert(kurtosis(x),6/a, 2.); | |
684 %!test | |
685 %! % statistical tests may fail occasionally. | |
686 %! a=10; x = randg(a,100000,1); | |
687 %! assert(mean(x), a, 0.1); | |
688 %! assert(var(x), a, 0.5); | |
689 %! assert(skewness(x),2/sqrt(a), 0.1); | |
690 %! assert(kurtosis(x),6/a, 0.5); | |
691 %!test | |
692 %! % statistical tests may fail occasionally. | |
693 %! a=100; x = randg(a,100000,1); | |
694 %! assert(mean(x), a, 0.2); | |
695 %! assert(var(x), a, 2.); | |
696 %! assert(skewness(x),2/sqrt(a), 0.05); | |
697 %! assert(kurtosis(x),6/a, 0.2); | |
698 %!test | |
699 %! rand("seed",12) | |
700 %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report | |
701 %!test | |
702 %! % statistical tests may fail occasionally. | |
703 %! a=0.1; x = randg(a,100000,1); | |
704 %! assert(mean(x), a, 0.01); | |
705 %! assert(var(x), a, 0.01); | |
706 %! assert(skewness(x),2/sqrt(a), 1.); | |
707 %! assert(kurtosis(x),6/a, 50.); | |
708 %!test | |
709 %! % statistical tests may fail occasionally. | |
710 %! a=0.95; x = randg(a,100000,1); | |
711 %! assert(mean(x), a, 0.01); | |
712 %! assert(var(x), a, 0.04); | |
713 %! assert(skewness(x),2/sqrt(a), 0.2); | |
714 %! assert(kurtosis(x),6/a, 2.); | |
715 %!test | |
716 %! % statistical tests may fail occasionally. | |
717 %! a=1; x = randg(a,100000,1); | |
718 %! assert(mean(x),a, 0.01); | |
719 %! assert(var(x),a, 0.04); | |
720 %! assert(skewness(x),2/sqrt(a), 0.2); | |
721 %! assert(kurtosis(x),6/a, 2.); | |
722 %!test | |
723 %! % statistical tests may fail occasionally. | |
724 %! a=10; x = randg(a,100000,1); | |
725 %! assert(mean(x), a, 0.1); | |
726 %! assert(var(x), a, 0.5); | |
727 %! assert(skewness(x),2/sqrt(a), 0.1); | |
728 %! assert(kurtosis(x),6/a, 0.5); | |
729 %!test | |
730 %! % statistical tests may fail occasionally. | |
731 %! a=100; x = randg(a,100000,1); | |
732 %! assert(mean(x), a, 0.2); | |
733 %! assert(var(x), a, 2.); | |
734 %! assert(skewness(x),2/sqrt(a), 0.05); | |
735 %! assert(kurtosis(x),6/a, 0.2); | |
736 */ | |
737 | |
738 | |
739 DEFUN_DLD (randp, args, , | |
740 "-*- texinfo -*-\n\ | |
741 @deftypefn {Loadable Function} {} randp (@var{l}, @var{x})\n\ | |
742 @deftypefnx {Loadable Function} {} randp (@var{l}, @var{n}, @var{m})\n\ | |
743 @deftypefnx {Loadable Function} {} randp (\"state\", @var{x})\n\ | |
744 @deftypefnx {Loadable Function} {} randp (\"seed\", @var{x})\n\ | |
745 Return a matrix with Poisson distributed random elements. The arguments\n\ | |
746 are handled the same as the arguments for @code{rand}, except for the\n\ | |
747 argument @var{l}.\n\ | |
748 \n\ | |
749 Five different algorithms are used depending on the range of @var{l}\n\ | |
750 and whether or not @var{l} is a scalar or a matrix.\n\ | |
751 \n\ | |
752 @table @asis\n\ | |
753 @item For scalar @var{l} <= 12, use direct method.\n\ | |
754 Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ | |
755 @item For scalar @var{l} > 12, use rejection method.[1]\n\ | |
756 Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ | |
757 @item For matrix @var{l} <= 10, use inversion method.[2]\n\ | |
758 Stadlober E., et al., WinRand source code, available via FTP.\n\ | |
759 @item For matrix @var{l} > 10, use patchwork rejection method.\n\ | |
760 Stadlober E., et al., WinRand source code, available via FTP, or\n\ | |
761 H. Zechner, 'Efficient sampling from continuous and discrete\n\ | |
762 unimodal distributions', Doctoral Dissertaion, 156pp., Technical\n\ | |
763 University Graz, Austria, 1994.\n\ | |
764 @item For @var{l} > 1e8, use normal approximation.\n\ | |
765 L. Montanet, et al., 'Review of Particle Properties', Physical Review\n\ | |
766 D 50 p1284, 1994\n\ | |
767 @end table\n\ | |
768 @seealso{rand,randn,rande,randg}\n\ | |
769 @end deftypefn") | |
770 { | |
771 octave_value retval; | |
772 | |
773 int nargin = args.length (); | |
774 | |
775 if (nargin < 1) | |
776 error ("randp: insufficient arguments"); | |
777 else | |
778 { | |
779 unwind_protect::begin_frame ("randp"); | |
780 | |
781 // This relies on the fact that elements are popped from the unwind | |
782 // stack in the reverse of the order they are pushed | |
783 // (i.e. current_distribution will be reset before calling | |
784 // reset_rand_generator()). | |
785 | |
786 unwind_protect::add (reset_rand_generator, 0); | |
787 unwind_protect_str (current_distribution); | |
788 | |
789 current_distribution = "poisson"; | |
790 | |
791 octave_rand::distribution (current_distribution); | |
792 | |
793 retval = do_rand (args, nargin, "randp", true); | |
794 | |
795 unwind_protect::run_frame ("randp"); | |
796 } | |
797 | |
798 return retval; | |
799 } | |
800 | |
801 /* | |
802 %!test | |
803 %! rand("state",12) | |
804 %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report | |
805 %!test | |
806 %! % statistical tests may fail occasionally. | |
807 %! for a=[5 15] | |
808 %! x = randp(a,100000,1); | |
809 %! assert(min(x)>=0); % *** Please report this!!! *** | |
810 %! assert(mean(x),a,0.03); | |
811 %! assert(var(x),a,0.2); | |
812 %! assert(skewness(x),1/sqrt(a),0.03); | |
813 %! assert(kurtosis(x),1/a,0.08); | |
814 %! end | |
815 %!test | |
816 %! % statistical tests may fail occasionally. | |
817 %! for a=[5 15] | |
818 %! x = randp(a*ones(100000,1),100000,1); | |
819 %! assert(min(x)>=0); % *** Please report this!!! *** | |
820 %! assert(mean(x),a,0.03); | |
821 %! assert(var(x),a,0.2); | |
822 %! assert(skewness(x),1/sqrt(a),0.03); | |
823 %! assert(kurtosis(x),1/a,0.08); | |
824 %! end | |
825 %!test | |
826 %! rand("seed",12) | |
827 %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report | |
828 %!test | |
829 %! % statistical tests may fail occasionally. | |
830 %! for a=[5 15] | |
831 %! x = randp(a,100000,1); | |
832 %! assert(min(x)>=0); % *** Please report this!!! *** | |
833 %! assert(mean(x),a,0.03); | |
834 %! assert(var(x),a,0.2); | |
835 %! assert(skewness(x),1/sqrt(a),0.03); | |
836 %! assert(kurtosis(x),1/a,0.08); | |
837 %! end | |
838 %!test | |
839 %! % statistical tests may fail occasionally. | |
840 %! for a=[5 15] | |
841 %! x = randp(a*ones(100000,1),100000,1); | |
842 %! assert(min(x)>=0); % *** Please report this!!! *** | |
843 %! assert(mean(x),a,0.03); | |
844 %! assert(var(x),a,0.2); | |
845 %! assert(skewness(x),1/sqrt(a),0.03); | |
846 %! assert(kurtosis(x),1/a,0.08); | |
847 %! end | |
848 */ | |
319 | 849 |
320 /* | 850 /* |
321 ;;; Local Variables: *** | 851 ;;; Local Variables: *** |
322 ;;; mode: C++ *** | 852 ;;; mode: C++ *** |
323 ;;; End: *** | 853 ;;; End: *** |