30 #if defined (HAVE_UNORDERED_MAP)
31 #include <unordered_map>
32 #elif defined (HAVE_TR1_UNORDERED_MAP)
33 #include <tr1/unordered_map>
58 const std::string& distribution,
bool additional_arg =
false)
64 bool is_single =
false;
73 if (nargin > 0 && args(nargin-1).is_string ())
75 std::string s_arg = args(nargin-1).string_value ();
77 if (s_arg ==
"single")
82 else if (s_arg ==
"double")
90 error (
"%s: expecting at least one argument", fcn);
93 else if (args(0).is_string ())
94 additional_arg =
false;
100 error (
"%s: expecting scalar or matrix arguments", fcn);
137 else if (s_arg ==
"seed")
141 else if (s_arg ==
"state" || s_arg ==
"twister")
145 else if (s_arg ==
"uniform")
149 else if (s_arg ==
"normal")
153 else if (s_arg ==
"exponential")
157 else if (s_arg ==
"poisson")
161 else if (s_arg ==
"gamma")
166 error (
"%s: unrecognized string argument", fcn);
174 error (
"%s: NaN is invalid matrix dimension", fcn);
204 dims(i) = base >= 0 ? base : 0;
212 error (
"%s: all elements of range must be integers",
230 dims(i) = elt >=0 ? elt : 0;
236 error (
"%s: expecting integer vector", fcn);
256 if (args(idx+1).is_real_scalar ())
258 double d = args(idx+1).double_value ();
263 else if (args(idx+1).is_string ()
264 && args(idx+1).string_value () ==
"reset")
267 error (
"%s: seed must be a real scalar", fcn);
269 else if (ts ==
"state" || ts ==
"twister")
271 if (args(idx+1).is_string ()
272 && args(idx+1).string_value () ==
"reset")
284 error (
"%s: unrecognized string argument", fcn);
290 for (
int i = 0; i < nargin; i++)
295 error (
"%s: expecting integer arguments", fcn);
299 dims(i) = elt >= 0 ? elt : 0;
324 if (a.
dims () != dims)
326 error (
"%s: mismatch in argument size", fcn);
348 if (a.
dims () != dims)
350 error (
"%s: mismatch in argument size", fcn);
368 @deftypefn {Built-in Function} {} rand (@var{n})\n\
369 @deftypefnx {Built-in Function} {} rand (@var{m}, @var{n}, @dots{})\n\
370 @deftypefnx {Built-in Function} {} rand ([@var{m} @var{n} @dots{}])\n\
371 @deftypefnx {Built-in Function} {@var{v} =} rand (\"state\")\n\
372 @deftypefnx {Built-in Function} {} rand (\"state\", @var{v})\n\
373 @deftypefnx {Built-in Function} {} rand (\"state\", \"reset\")\n\
374 @deftypefnx {Built-in Function} {@var{v} =} rand (\"seed\")\n\
375 @deftypefnx {Built-in Function} {} rand (\"seed\", @var{v})\n\
376 @deftypefnx {Built-in Function} {} rand (\"seed\", \"reset\")\n\
377 @deftypefnx {Built-in Function} {} rand (@dots{}, \"single\")\n\
378 @deftypefnx {Built-in Function} {} rand (@dots{}, \"double\")\n\
379 Return a matrix with random elements uniformly distributed on the\n\
382 The arguments are handled the same as the arguments for @code{eye}.\n\
384 You can query the state of the random number generator using the form\n\
387 v = rand (\"state\")\n\
390 This returns a column vector @var{v} of length 625. Later, you can restore\n\
391 the random number generator to the state @var{v} using the form\n\
394 rand (\"state\", v)\n\
398 You may also initialize the state vector from an arbitrary vector of length\n\
399 @leq{} 625 for @var{v}. This new state will be a hash based on the value of\n\
400 @var{v}, not @var{v} itself.\n\
402 By default, the generator is initialized from @code{/dev/urandom} if it is\n\
403 available, otherwise from CPU time, wall clock time, and the current\n\
404 fraction of a second. Note that this differs from @sc{matlab}, which\n\
405 always initializes the state to the same state at startup. To obtain\n\
406 behavior comparable to @sc{matlab}, initialize with a deterministic state\n\
407 vector in Octave's startup files (@pxref{Startup Files}).\n\
409 To compute the pseudo-random sequence, @code{rand} uses the Mersenne\n\
410 Twister with a period of @math{2^{19937}-1}\n\
411 (See @nospell{M. Matsumoto and T. Nishimura},\n\
412 @cite{Mersenne Twister: A 623-dimensionally equidistributed uniform\n\
413 pseudorandom number generator},\n\
414 ACM Trans. on Modeling and Computer Simulation Vol. 8, No. 1, pp. 3--30,\n\
416 @url{http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html}).\n\
417 Do @strong{not} use for cryptography without securely hashing several\n\
418 returned values together, otherwise the generator state can be learned after\n\
419 reading 624 consecutive values.\n\
421 Older versions of Octave used a different random number generator.\n\
422 The new generator is used by default as it is significantly faster than the\n\
423 old generator, and produces random numbers with a significantly longer cycle\n\
424 time. However, in some circumstances it might be desirable to obtain the\n\
425 same random sequences as produced by the old generators. To do this the\n\
426 keyword @qcode{\"seed\"} is used to specify that the old generators should\n\
430 rand (\"seed\", val)\n\
434 which sets the seed of the generator to @var{val}. The seed of the\n\
435 generator can be queried with\n\
438 s = rand (\"seed\")\n\
441 However, it should be noted that querying the seed will not cause\n\
442 @code{rand} to use the old generators, only setting the seed will. To cause\n\
443 @code{rand} to once again use the new generators, the keyword\n\
444 @qcode{\"state\"} should be used to reset the state of the @code{rand}.\n\
446 The state or seed of the generator can be reset to a new random value using\n\
447 the @qcode{\"reset\"} keyword.\n\
449 The class of the value returned can be controlled by a trailing\n\
450 @qcode{\"double\"} or @qcode{\"single\"} argument. These are the only valid\n\
452 @seealso{randn, rande, randg, randp}\n\
457 int nargin = args.
length ();
459 retval =
do_rand (args, nargin,
"rand",
"uniform");
551 DEFUN (randn, args, ,
553 @deftypefn {Built-in Function} {} randn (@var{n})\n\
554 @deftypefnx {Built-in Function} {} randn (@var{m}, @var{n}, @dots{})\n\
555 @deftypefnx {Built-in Function} {} randn ([@var{m} @var{n} @dots{}])\n\
556 @deftypefnx {Built-in Function} {@var{v} =} randn (\"state\")\n\
557 @deftypefnx {Built-in Function} {} randn (\"state\", @var{v})\n\
558 @deftypefnx {Built-in Function} {} randn (\"state\", \"reset\")\n\
559 @deftypefnx {Built-in Function} {@var{v} =} randn (\"seed\")\n\
560 @deftypefnx {Built-in Function} {} randn (\"seed\", @var{v})\n\
561 @deftypefnx {Built-in Function} {} randn (\"seed\", \"reset\")\n\
562 @deftypefnx {Built-in Function} {} randn (@dots{}, \"single\")\n\
563 @deftypefnx {Built-in Function} {} randn (@dots{}, \"double\")\n\
564 Return a matrix with normally distributed random elements having zero mean\n\
567 The arguments are handled the same as the arguments for @code{rand}.\n\
569 By default, @code{randn} uses the @nospell{Marsaglia and Tsang}\n\
570 ``Ziggurat technique'' to transform from a uniform to a normal distribution.\n\
572 The class of the value returned can be controlled by a trailing\n\
573 @qcode{\"double\"} or @qcode{\"single\"} argument. These are the only valid\n\
576 Reference: @nospell{G. Marsaglia and W.W. Tsang},\n\
577 @cite{Ziggurat Method for Generating Random Variables},\n\
578 J. Statistical Software, vol 5, 2000,\n\
579 @url{http://www.jstatsoft.org/v05/i08/}\n\
581 @seealso{rand, rande, randg, randp}\n\
586 int nargin = args.
length ();
588 retval =
do_rand (args, nargin,
"randn",
"normal");
624 DEFUN (rande, args, ,
626 @deftypefn {Built-in Function} {} rande (@var{n})\n\
627 @deftypefnx {Built-in Function} {} rande (@var{m}, @var{n}, @dots{})\n\
628 @deftypefnx {Built-in Function} {} rande ([@var{m} @var{n} @dots{}])\n\
629 @deftypefnx {Built-in Function} {@var{v} =} rande (\"state\")\n\
630 @deftypefnx {Built-in Function} {} rande (\"state\", @var{v})\n\
631 @deftypefnx {Built-in Function} {} rande (\"state\", \"reset\")\n\
632 @deftypefnx {Built-in Function} {@var{v} =} rande (\"seed\")\n\
633 @deftypefnx {Built-in Function} {} rande (\"seed\", @var{v})\n\
634 @deftypefnx {Built-in Function} {} rande (\"seed\", \"reset\")\n\
635 @deftypefnx {Built-in Function} {} rande (@dots{}, \"single\")\n\
636 @deftypefnx {Built-in Function} {} rande (@dots{}, \"double\")\n\
637 Return a matrix with exponentially distributed random elements.\n\
639 The arguments are handled the same as the arguments for @code{rand}.\n\
641 By default, @code{randn} uses the @nospell{Marsaglia and Tsang}\n\
642 ``Ziggurat technique'' to transform from a uniform to a normal distribution.\n\
644 The class of the value returned can be controlled by a trailing\n\
645 @qcode{\"double\"} or @qcode{\"single\"} argument. These are the only valid\n\
648 Reference: @nospell{G. Marsaglia and W.W. Tsang},\n\
649 @cite{Ziggurat Method for Generating Random Variables},\n\
650 J. Statistical Software, vol 5, 2000,\n\
651 @url{http://www.jstatsoft.org/v05/i08/}\n\
653 @seealso{rand, randn, randg, randp}\n\
658 int nargin = args.
length ();
660 retval =
do_rand (args, nargin,
"rande",
"exponential");
698 DEFUN (randg, args, ,
700 @deftypefn {Built-in Function} {} randg (@var{n})\n\
701 @deftypefnx {Built-in Function} {} randg (@var{m}, @var{n}, @dots{})\n\
702 @deftypefnx {Built-in Function} {} randg ([@var{m} @var{n} @dots{}])\n\
703 @deftypefnx {Built-in Function} {@var{v} =} randg (\"state\")\n\
704 @deftypefnx {Built-in Function} {} randg (\"state\", @var{v})\n\
705 @deftypefnx {Built-in Function} {} randg (\"state\", \"reset\")\n\
706 @deftypefnx {Built-in Function} {@var{v} =} randg (\"seed\")\n\
707 @deftypefnx {Built-in Function} {} randg (\"seed\", @var{v})\n\
708 @deftypefnx {Built-in Function} {} randg (\"seed\", \"reset\")\n\
709 @deftypefnx {Built-in Function} {} randg (@dots{}, \"single\")\n\
710 @deftypefnx {Built-in Function} {} randg (@dots{}, \"double\")\n\
711 Return a matrix with @code{gamma (@var{a},1)} distributed random elements.\n\
713 The arguments are handled the same as the arguments for @code{rand}, except\n\
714 for the argument @var{a}.\n\
716 This can be used to generate many distributions:\n\
719 @item @code{gamma (a, b)} for @code{a > -1}, @code{b > 0}\n\
725 @item @code{beta (a, b)} for @code{a > -1}, @code{b > -1}\n\
730 r = r1 / (r1 + randg (b, 1))\n\
734 @item @code{Erlang (a, n)}\n\
740 @item @code{chisq (df)} for @code{df > 0}\n\
743 r = 2 * randg (df / 2)\n\
746 @item @code{t (df)} for @code{0 < df < inf} (use randn if df is infinite)\n\
749 r = randn () / sqrt (2 * randg (df / 2) / df)\n\
752 @item @code{F (n1, n2)} for @code{0 < n1}, @code{0 < n2}\n\
756 ## r1 equals 1 if n1 is infinite\n\
757 r1 = 2 * randg (n1 / 2) / n1\n\
758 ## r2 equals 1 if n2 is infinite\n\
759 r2 = 2 * randg (n2 / 2) / n2\n\
764 @item negative @code{binomial (n, p)} for @code{n > 0}, @code{0 < p <= 1}\n\
767 r = randp ((1 - p) / p * randg (n))\n\
770 @item non-central @code{chisq (df, L)}, for @code{df >= 0} and @code{L > 0}\n\
771 (use chisq if @code{L = 0})\n\
776 r(r > 0) = 2 * randg (r(r > 0))\n\
777 r(df > 0) += 2 * randg (df(df > 0)/2)\n\
781 @item @code{Dirichlet (a1, @dots{} ak)}\n\
785 r = (randg (a1), @dots{}, randg (ak))\n\
792 The class of the value returned can be controlled by a trailing\n\
793 @qcode{\"double\"} or @qcode{\"single\"} argument. These are the only valid\n\
795 @seealso{rand, randn, rande, randp}\n\
800 int nargin = args.
length ();
803 error (
"randg: insufficient arguments");
805 retval =
do_rand (args, nargin,
"randg",
"gamma",
true);
976 DEFUN (randp, args, ,
978 @deftypefn {Built-in Function} {} randp (@var{l}, @var{n})\n\
979 @deftypefnx {Built-in Function} {} randp (@var{l}, @var{m}, @var{n}, @dots{})\n\
980 @deftypefnx {Built-in Function} {} randp (@var{l}, [@var{m} @var{n} @dots{}])\n\
981 @deftypefnx {Built-in Function} {@var{v} =} randp (\"state\")\n\
982 @deftypefnx {Built-in Function} {} randp (\"state\", @var{v})\n\
983 @deftypefnx {Built-in Function} {} randp (\"state\", \"reset\")\n\
984 @deftypefnx {Built-in Function} {@var{v} =} randp (\"seed\")\n\
985 @deftypefnx {Built-in Function} {} randp (\"seed\", @var{v})\n\
986 @deftypefnx {Built-in Function} {} randp (\"seed\", \"reset\")\n\
987 @deftypefnx {Built-in Function} {} randp (@dots{}, \"single\")\n\
988 @deftypefnx {Built-in Function} {} randp (@dots{}, \"double\")\n\
989 Return a matrix with Poisson distributed random elements with mean value\n\
990 parameter given by the first argument, @var{l}.\n\
992 The arguments are handled the same as the arguments for @code{rand}, except\n\
993 for the argument @var{l}.\n\
995 Five different algorithms are used depending on the range of @var{l} and\n\
996 whether or not @var{l} is a scalar or a matrix.\n\
999 @item For scalar @var{l} @leq{} 12, use direct method.\n\
1000 W.H. Press, et al., @cite{Numerical Recipes in C},\n\
1001 Cambridge University Press, 1992.\n\
1003 @item For scalar @var{l} > 12, use rejection method.[1]\n\
1004 W.H. Press, et al., @cite{Numerical Recipes in C},\n\
1005 Cambridge University Press, 1992.\n\
1007 @item For matrix @var{l} @leq{} 10, use inversion method.[2]\n\
1008 @nospell{E. Stadlober, et al., WinRand source code}, available via FTP.\n\
1010 @item For matrix @var{l} > 10, use patchwork rejection method.\n\
1011 @nospell{E. Stadlober, et al., WinRand source code}, available via FTP, or\n\
1012 @nospell{H. Zechner}, @cite{Efficient sampling from continuous and discrete\n\
1013 unimodal distributions}, Doctoral Dissertation, 156pp., Technical\n\
1014 University @nospell{Graz}, Austria, 1994.\n\
1016 @item For @var{l} > 1e8, use normal approximation.\n\
1017 @nospell{L. Montanet}, et al., @cite{Review of Particle Properties},\n\
1018 Physical Review D 50 p1284, 1994.\n\
1021 The class of the value returned can be controlled by a trailing\n\
1022 @qcode{\"double\"} or @qcode{\"single\"} argument. These are the only valid\n\
1024 @seealso{rand, randn, rande, randg}\n\
1029 int nargin = args.
length ();
1032 error (
"randp: insufficient arguments");
1034 retval =
do_rand (args, nargin,
"randp",
"poisson",
true);
1125 DEFUN (randperm, args, ,
1127 @deftypefn {Built-in Function} {} randperm (@var{n})\n\
1128 @deftypefnx {Built-in Function} {} randperm (@var{n}, @var{m})\n\
1129 Return a row vector containing a random permutation of @code{1:@var{n}}.\n\
1131 If @var{m} is supplied, return @var{m} unique entries, sampled without\n\
1132 replacement from @code{1:@var{n}}.\n\
1134 The complexity is O(@var{n}) in memory and O(@var{m}) in time, unless\n\
1135 @var{m} < @var{n}/5, in which case O(@var{m}) memory is used as well. The\n\
1136 randomization is performed using rand(). All permutations are equally\n\
1142 #ifdef USE_UNORDERED_MAP_WITH_TR1
1143 using std::tr1::unordered_map;
1145 using std::unordered_map;
1148 int nargin = args.length ();
1151 if (nargin == 1 || nargin == 2)
1155 n = args(0).idx_type_value (
true);
1158 m = args(1).idx_type_value (
true);
1163 error (
"randperm: M and N must be non-negative");
1166 error (
"randperm: M must be less than or equal to N");
1170 bool short_shuffle = m < n/5;
1184 catch (std::bad_alloc)
1189 short_shuffle =
true;
1200 unordered_map<octave_idx_type, octave_idx_type> map (m);
1213 std::swap (ivec[i], ivec[k]);
1217 if (map.find (k) == map.end ())
1220 std::swap (ivec[i], map[k]);
1232 std::swap (ivec[i], ivec[k]);
1238 rvec[i] = ivec[i] + 1;
bool is_range(void) const
void gripe_wrong_type_arg(const char *name, const char *s, bool is_error)
static void exponential_distribution(void)
static octave_value do_rand(const octave_value_list &args, int nargin, const char *fcn, const std::string &distribution, bool additional_arg=false)
OCTINTERP_API void print_usage(void)
octave_idx_type nelem(void) const
bool is_scalar_type(void) const
void resize(int n, int fill_value=0)
#define DEFUN(name, args_name, nargout_name, doc)
void error(const char *fmt,...)
static float float_scalar(float a=1.0)
F77_RET_T const double const double double * d
const dim_vector & dims(void) const
Return a const-reference so that dims ()(i) works efficiently.
static std::string current_distribution
void add_fcn(void(*fcn)(void))
std::string string_value(bool force=false) const
static void gamma_distribution(void)
Range range_value(void) const
bool is_matrix_type(void) const
bool is_string(void) const
void resize(const dim_vector &dv, const T &rfv)
static void poisson_distribution(void)
octave_idx_type length(void) const
static ColumnVector state(const std::string &d=std::string())
static NDArray nd_array(const dim_vector &dims, double a=1.0)
Handles the reference counting for all the derived classes.
octave_idx_type length(void) const
Number of elements in the array.
static FloatNDArray float_nd_array(const dim_vector &dims, float a=1.0)
Array< octave_value > array_value(void) const
static void normal_distribution(void)
bool all_elements_are_ints(void) const
static std::string distribution(void)
Array< int > int_vector_value(bool req_int=false, bool frc_str_conv=false, bool frc_vec_conv=false) const
std::complex< T > floor(const std::complex< T > &x)
const T * fortran_vec(void) const
octave_idx_type NINTbig(double x)
static double scalar(double a=1.0)
double double_value(bool frc_str_conv=false) const
void chop_trailing_singletons(void)
static void uniform_distribution(void)