# HG changeset patch # User jwe # Date 1144091906 0 # Node ID 109fdf7b3dcb3109eacd1b381edcbb05afbe987b # Parent e065f7c18bdc7c16377f2ae6bf81c1c4f95fb623 [project @ 2006-04-03 19:18:26 by jwe] diff --git a/doc/ChangeLog b/doc/ChangeLog --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,3 +1,8 @@ +2006-04-03 David Bateman + + * interpreter/matrix.txi: Add rande, randp, randg and update + for different random generator behavior. + 2006-03-28 John W. Eaton * texinfo.tex: Update FSF address. diff --git a/doc/interpreter/matrix.txi b/doc/interpreter/matrix.txi --- a/doc/interpreter/matrix.txi +++ b/doc/interpreter/matrix.txi @@ -143,8 +143,15 @@ @DOCSTRING(randn) -The @code{rand} and @code{randn} functions use separate generators. -This ensures that +@DOCSTRING(rande) + +@DOCSTRING(randp) + +@DOCSTRING(randg) + +The new random generators all use a common Mersenne Twister generator, +and so the state of only one of the generators needs to be reset. +The old generator function use separate generators. This ensures that @example @group diff --git a/libcruft/ChangeLog b/libcruft/ChangeLog --- a/libcruft/ChangeLog +++ b/libcruft/ChangeLog @@ -1,3 +1,7 @@ +2006-04-03 David Bateman + + * ranlib/wrap.f (dgenexp, dgengam, dignpoi): New functions. + 2006-03-21 John W. Eaton * misc/f77-fcn.h (F77_XFCN): Save octave_interrupt_immediately and diff --git a/libcruft/ranlib/wrap.f b/libcruft/ranlib/wrap.f --- a/libcruft/ranlib/wrap.f +++ b/libcruft/ranlib/wrap.f @@ -8,3 +8,18 @@ result = genunf (real (low), real (high)) return end + subroutine dgenexp (av, result) + double precision av, result + result = genexp (real (av)) + return + end + subroutine dgengam (a, r, result) + double precision a, r, result + result = gengam (real (a), real (r)) + return + end + subroutine dignpoi (mu, result) + double precision mu, result + result = ignpoi (real (mu)) + return + end diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog --- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,32 @@ +2006-04-03 David Bateman + + * Makefile.in (INCLUDES): Add randgamma.h, randpoisson.h and + randmtzig.h to the list. + (LIBOCTAVE_C_SOURCES): Add randgamma.c, randpoisson.c and + randmtzig.c to the list. + * oct-rand.cc (do_old_initialization): Rename from do_initialization. + (use_old_generators): New variable. + (old_initialized): Rename from initialized. + (new_initialized): New variable. + (oct_init_by_entropy): New function. + (maybe_initialize): Initialize new or old generator depending on + value of use_old_generators. + (octave_rand::state): New functions. + (octave_rand::distribution): Add gamma, exponential and poisson + distributions. + (octave_rand::exponential_distribution, + octave_rand::poisson_distribution, + octave_rand::gamma_distribution): New methods to select + exponential, poisson or gamma distribution. + (octave_rand::scalar, octave_rand::matrix, octave_rand::nd_array, + octave_rand::vector): Add new distributions. + * oct-rand.h: Provide decls for new functions. + (octave_rand::matrix, octave_rand::scalar, octave_rand:: + (octave_rand::scalar, octave_rand::matrix, octave_rand::nd_array, + octave_rand::vector): New arg A, for gamma and poisson distributions. + * randpoisson.c, randpoisson.h, randgamma.c, randmtzig.c, + randmtzig.h: New files. + 2006-03-24 John W. Eaton * dSparse.cc (SparseMatrix::bsolve): Integer work vector is diff --git a/liboctave/Makefile.in b/liboctave/Makefile.in --- a/liboctave/Makefile.in +++ b/liboctave/Makefile.in @@ -69,8 +69,9 @@ oct-passwd.h oct-rand.h oct-rl-edit.h oct-rl-hist.h \ oct-shlib.h oct-sort.h oct-spparms.h oct-syscalls.h \ oct-sparse.h oct-time.h oct-types.h oct-uname.h \ - pathlen.h pathsearch.h \ - prog-args.h so-array.h sparse-sort.h statdefs.h str-vec.h \ + pathlen.h pathsearch.h prog-args.h \ + randgamma.h randmtzig.h randpoisson.h \ + so-array.h sparse-sort.h statdefs.h str-vec.h \ sparse-util.h sun-utils.h sysdir.h systime.h syswait.h \ $(OPTS_INC) \ $(MATRIX_INC) \ @@ -128,7 +129,8 @@ $(SPARSE_MX_OP_SRC) LIBOCTAVE_C_SOURCES := f2c-main.c filemode.c getopt.c getopt1.c \ - lo-cieee.c lo-cutils.c mkdir.c oct-getopt.c rename.c \ + lo-cieee.c lo-cutils.c mkdir.c oct-getopt.c \ + randgamma.c randmtzig.c randpoisson.c rename.c \ rmdir.c strftime.c strptime.c tempname.c tempnam.c LIBOCTAVE_SOURCES := $(LIBOCTAVE_CXX_SOURCES) $(LIBOCTAVE_C_SOURCES) diff --git a/liboctave/oct-rand.cc b/liboctave/oct-rand.cc --- a/liboctave/oct-rand.cc +++ b/liboctave/oct-rand.cc @@ -24,22 +24,34 @@ #ifdef HAVE_CONFIG_H #include #endif +#include #include "f77-fcn.h" +#include "lo-ieee.h" #include "lo-error.h" +#include "lo-mappers.h" #include "oct-rand.h" #include "oct-time.h" +#include "data-conv.h" +#include "randmtzig.h" +#include "randpoisson.h" +#include "randgamma.h" // Possible distributions of random numbers. This was handled with an // enum, but unwind_protecting that doesn't work so well. #define uniform_dist 1 #define normal_dist 2 +#define expon_dist 3 +#define poisson_dist 4 +#define gamma_dist 5 // Current distribution of random numbers. static int current_distribution = uniform_dist; -// Has the seed been set yet? -static bool initialized = false; +// Has the seed/state been set yet? +static bool old_initialized = false; +static bool new_initialized = false; +static bool use_old_generators = false; extern "C" { @@ -50,6 +62,15 @@ F77_FUNC (dgenunf, DGENUNF) (const double&, const double&, double&); F77_RET_T + F77_FUNC (dgenexp, DGENEXP) (const double&, double&); + + F77_RET_T + F77_FUNC (dignpoi, DIGNPOI) (const double&, double&); + + F77_RET_T + F77_FUNC (dgengam, DGENGAM) (const double&, const double&, double&); + + F77_RET_T F77_FUNC (setall, SETALL) (const octave_idx_type&, const octave_idx_type&); F77_RET_T @@ -83,7 +104,7 @@ // work ok to give fairly different seeds each time Octave starts. static void -do_initialization (void) +do_old_initialization (void) { octave_localtime tm; @@ -99,20 +120,32 @@ F77_FUNC (setall, SETALL) (s0, s1); - initialized = true; + old_initialized = true; } static inline void maybe_initialize (void) { - if (! initialized) - do_initialization (); + if (use_old_generators) + { + if (! old_initialized) + do_old_initialization (); + } + else + { + if (! new_initialized) + { + oct_init_by_entropy (); + new_initialized = true; + } + } } double octave_rand::seed (void) { - maybe_initialize (); + if (! old_initialized) + do_old_initialization (); union d2i { double d; octave_idx_type i[2]; }; union d2i u; @@ -123,6 +156,7 @@ void octave_rand::seed (double s) { + use_old_generators = true; maybe_initialize (); union d2i { double d; octave_idx_type i[2]; }; @@ -133,6 +167,41 @@ F77_FUNC (setsd, SETSD) (i0, i1); } +ColumnVector +octave_rand::state (void) +{ + ColumnVector s (MT_N + 1); + if (! new_initialized) + { + oct_init_by_entropy (); + new_initialized = true; + } + + OCTAVE_LOCAL_BUFFER (unsigned FOUR_BYTE_INT, tmp, MT_N + 1); + oct_get_state (tmp); + for (octave_idx_type i = 0; i <= MT_N; i++) + s.elem (i) = static_cast(tmp [i]); + return s; +} + +void +octave_rand::state (const ColumnVector &s) +{ + use_old_generators = false; + maybe_initialize (); + + octave_idx_type len = s.length(); + octave_idx_type n = len < MT_N + 1 ? len : MT_N + 1; + OCTAVE_LOCAL_BUFFER (unsigned FOUR_BYTE_INT, tmp, MT_N + 1); + for (octave_idx_type i = 0; i < n; i++) + tmp[i] = static_cast (s.elem(i)); + + if (len == MT_N + 1 && tmp[MT_N] <= MT_N && tmp[MT_N] > 0) + oct_set_state (tmp); + else + oct_init_by_array (tmp, len); +} + std::string octave_rand::distribution (void) { @@ -142,6 +211,12 @@ return "uniform"; else if (current_distribution == normal_dist) return "normal"; + else if (current_distribution == expon_dist) + return "exponential"; + else if (current_distribution == poisson_dist) + return "poisson"; + else if (current_distribution == gamma_dist) + return "gamma"; else { abort (); @@ -156,6 +231,12 @@ octave_rand::uniform_distribution (); else if (d == "normal") octave_rand::normal_distribution (); + else if (d == "exponential") + octave_rand::exponential_distribution (); + else if (d == "poisson") + octave_rand::poisson_distribution (); + else if (d == "gamma") + octave_rand::gamma_distribution (); else (*current_liboctave_error_handler) ("rand: invalid distribution"); } @@ -180,114 +261,105 @@ F77_FUNC (setcgn, SETCGN) (normal_dist); } +void +octave_rand::exponential_distribution (void) +{ + maybe_initialize (); + + current_distribution = expon_dist; + + F77_FUNC (setcgn, SETCGN) (expon_dist); +} + +void +octave_rand::poisson_distribution (void) +{ + maybe_initialize (); + + current_distribution = poisson_dist; + + F77_FUNC (setcgn, SETCGN) (poisson_dist); +} + +void +octave_rand::gamma_distribution (void) +{ + maybe_initialize (); + + current_distribution = gamma_dist; + + F77_FUNC (setcgn, SETCGN) (gamma_dist); +} + + double -octave_rand::scalar (void) +octave_rand::scalar (double a) { maybe_initialize (); double retval = 0.0; - switch (current_distribution) + if (use_old_generators) { - case uniform_dist: - F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval); - break; + switch (current_distribution) + { + case uniform_dist: + F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval); + break; - case normal_dist: - F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval); - break; + case normal_dist: + F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval); + break; - default: - abort (); - break; - } - - return retval; -} + case expon_dist: + F77_FUNC (dgenexp, DGENEXP) (1.0, retval); + break; -#define MAKE_RAND_MAT(mat, nr, nc, f, F) \ - do \ - { \ - double val; \ - for (volatile octave_idx_type j = 0; j < nc; j++) \ - for (volatile octave_idx_type i = 0; i < nr; i++) \ - { \ - OCTAVE_QUIT; \ - F77_FUNC (f, F) (0.0, 1.0, val); \ - mat(i,j) = val; \ - } \ - } \ - while (0) - -Matrix -octave_rand::matrix (octave_idx_type n, octave_idx_type m) -{ - maybe_initialize (); - - Matrix retval; + case poisson_dist: + if (a < 0.0 || xisnan(a) || xisinf(a)) + retval = octave_NaN; + else + { + // workaround bug in ignpoi, by calling with different Mu + F77_FUNC (dignpoi, DIGNPOI) (a + 1, retval); + F77_FUNC (dignpoi, DIGNPOI) (a, retval); + } + break; - if (n >= 0 && m >= 0) - { - retval.resize (n, m); + case gamma_dist: + if (a <= 0.0 || xisnan(a) || xisinf(a)) + retval = octave_NaN; + else + F77_FUNC (dgengam, DGENGAM) (1.0, a, retval); + break; - if (n > 0 && m > 0) - { - switch (current_distribution) - { - case uniform_dist: - MAKE_RAND_MAT (retval, n, m, dgenunf, DGENUNF); - break; - - case normal_dist: - MAKE_RAND_MAT (retval, n, m, dgennor, DGENNOR); - break; - - default: - abort (); - break; - } + default: + abort (); + break; } } else - (*current_liboctave_error_handler) ("rand: invalid negative argument"); - - return retval; -} - -#define MAKE_RAND_ND_ARRAY(mat, len, f, F) \ - do \ - { \ - double val; \ - for (volatile octave_idx_type i = 0; i < len; i++) \ - { \ - OCTAVE_QUIT; \ - F77_FUNC (f, F) (0.0, 1.0, val); \ - mat(i) = val; \ - } \ - } \ - while (0) - -NDArray -octave_rand::nd_array (const dim_vector& dims) -{ - maybe_initialize (); - - NDArray retval; - - if (! dims.all_zero ()) { - retval.resize (dims); - - octave_idx_type len = retval.length (); - switch (current_distribution) { case uniform_dist: - MAKE_RAND_ND_ARRAY (retval, len, dgenunf, DGENUNF); + retval = oct_randu(); break; case normal_dist: - MAKE_RAND_ND_ARRAY (retval, len, dgennor, DGENNOR); + retval = oct_randn(); + break; + + case expon_dist: + retval = oct_rande(); + break; + + case poisson_dist: + retval = oct_randp(a); + break; + + case gamma_dist: + retval = oct_randg(a); break; default: @@ -299,21 +371,142 @@ return retval; } -#define MAKE_RAND_ARRAY(array, n, f, F) \ +#define MAKE_RAND(len) \ do \ { \ double val; \ - for (volatile octave_idx_type i = 0; i < n; i++) \ + for (volatile octave_idx_type i = 0; i < len; i++) \ { \ OCTAVE_QUIT; \ - F77_FUNC (f, F) (0.0, 1.0, val); \ - array(i) = val; \ + RAND_FUNC (val); \ + v[i] = val; \ } \ } \ while (0) +static void +fill_rand (octave_idx_type len, double *v, double a) +{ + maybe_initialize (); + + if (len < 1) + return; + + switch (current_distribution) + { + case uniform_dist: + if (use_old_generators) + { +#define RAND_FUNC(x) F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, x) + MAKE_RAND (len); +#undef RAND_FUNC + } + else + oct_fill_randu (len, v); + break; + + case normal_dist: + if (use_old_generators) + { +#define RAND_FUNC(x) F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, x) + MAKE_RAND (len); +#undef RAND_FUNC + } + else + oct_fill_randn (len, v); + break; + + case expon_dist: + if (use_old_generators) + { +#define RAND_FUNC(x) F77_FUNC (dgenexp, DGENEXP) (1.0, x) + MAKE_RAND (len); +#undef RAND_FUNC + } + else + oct_fill_rande (len, v); + break; + + case poisson_dist: + if (use_old_generators) + { + if (a < 0.0 || xisnan(a) || xisinf(a)) +#define RAND_FUNC(x) x = octave_NaN; + MAKE_RAND (len); +#undef RAND_FUNC + else + { + // workaround bug in ignpoi, by calling with different Mu + double tmp; + F77_FUNC (dignpoi, DIGNPOI) (a + 1, tmp); +#define RAND_FUNC(x) F77_FUNC (dignpoi, DIGNPOI) (a, x) + MAKE_RAND (len); +#undef RAND_FUNC + } + } + else + oct_fill_randp (a, len, v); + break; + + case gamma_dist: + if (use_old_generators) + { + if (a <= 0.0 || xisnan(a) || xisinf(a)) +#define RAND_FUNC(x) x = octave_NaN; + MAKE_RAND (len); +#undef RAND_FUNC + else +#define RAND_FUNC(x) F77_FUNC (dgengam, DGENGAM) (1.0, a, x) + MAKE_RAND (len); +#undef RAND_FUNC + } + else + oct_fill_randg (a, len, v); + break; + + default: + abort (); + break; + } + + return; +} + +Matrix +octave_rand::matrix (octave_idx_type n, octave_idx_type m, double a) +{ + Matrix retval; + + if (n >= 0 && m >= 0) + { + retval.resize (n, m); + + if (n > 0 && m > 0) + fill_rand (retval.capacity(), retval.fortran_vec(), a); + } + else + (*current_liboctave_error_handler) ("rand: invalid negative argument"); + + return retval; +} + +NDArray +octave_rand::nd_array (const dim_vector& dims, double a) +{ + NDArray retval; + + if (! dims.all_zero ()) + { + retval.resize (dims); + + fill_rand (retval.capacity(), retval.fortran_vec(), a); + } + + return retval; +} + Array -octave_rand::vector (octave_idx_type n) +octave_rand::vector (octave_idx_type n, double a) { maybe_initialize (); @@ -323,20 +516,7 @@ { retval.resize (n); - switch (current_distribution) - { - case uniform_dist: - MAKE_RAND_ARRAY (retval, n, dgenunf, DGENUNF); - break; - - case normal_dist: - MAKE_RAND_ARRAY (retval, n, dgennor, DGENNOR); - break; - - default: - abort (); - break; - } + fill_rand (retval.capacity(), retval.fortran_vec(), a); } else if (n < 0) (*current_liboctave_error_handler) ("rand: invalid negative argument"); diff --git a/liboctave/oct-rand.h b/liboctave/oct-rand.h --- a/liboctave/oct-rand.h +++ b/liboctave/oct-rand.h @@ -26,6 +26,7 @@ #include +#include "dColVector.h" #include "dMatrix.h" #include "dNDArray.h" @@ -38,6 +39,12 @@ // Set the seed. static void seed (double s); + // Return the current state. + static ColumnVector state (void); + + // Set the current state/ + static void state (const ColumnVector &s); + // Return the current distribution. static std::string distribution (void); @@ -49,19 +56,25 @@ static void normal_distribution (void); + static void exponential_distribution (void); + + static void poisson_distribution (void); + + static void gamma_distribution (void); + // Return the next number from the sequence. - static double scalar (void); + static double scalar (double a = 1.); // Return a matrix of numbers from the sequence, filled in column // major order. - static Matrix matrix (octave_idx_type r, octave_idx_type c); + static Matrix matrix (octave_idx_type r, octave_idx_type c, double a = 1.); // Return an N-dimensional array of numbers from the sequence, // filled in column major order. - static NDArray nd_array (const dim_vector& dims); + static NDArray nd_array (const dim_vector& dims, double a = 1.); // Return an array of numbers from the sequence. - static Array vector (octave_idx_type n); + static Array vector (octave_idx_type n, double a = 1.); }; #endif diff --git a/src/ChangeLog b/src/ChangeLog --- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,5 +1,12 @@ 2006-04-03 David Bateman + * DLD-FUNCTIONS/rand.cc (do_rand): Additional argument for + gamma and poisson distributions. Change "state" and "seed" + arguments so that they choose between generators. + Add, poisson, gamma and exponential generators. + (Frand, Frandn): Update docs for new generators, add tests. + (Frande, Frandp, Frandg): New generators, with test code. + * DLD-FUNCTIONS/daspk.cc (Fdaspk): Allow functions to be passed using function handles, inline functions, and cell arrays of strings, inline and function handles. diff --git a/src/DLD-FUNCTIONS/rand.cc b/src/DLD-FUNCTIONS/rand.cc --- a/src/DLD-FUNCTIONS/rand.cc +++ b/src/DLD-FUNCTIONS/rand.cc @@ -42,28 +42,56 @@ #include "utils.h" static octave_value -do_rand (const octave_value_list& args, int nargin, const char *fcn) +do_rand (const octave_value_list& args, int nargin, const char *fcn, + bool additional_arg = false) { octave_value retval; + NDArray a; + int idx = 0; + dim_vector dims; - dim_vector dims; + if (additional_arg) + { + if (nargin == 0) + { + error ("%s: expecting at least one argument", fcn); + goto done; + } + else if (args(0).is_string()) + additional_arg = false; + else + { + a = args(0).array_value (); + if (error_state) + { + error ("%s: expecting scalar or matrix arguments", fcn); + goto done; + } + idx++; + nargin--; + } + } switch (nargin) { case 0: { - dims.resize (2); + if (additional_arg) + dims = a.dims (); + else + { + dims.resize (2); - dims(0) = 1; - dims(1) = 1; - + dims(0) = 1; + dims(1) = 1; + } goto gen_matrix; } break; case 1: { - octave_value tmp = args(0); + octave_value tmp = args(idx); if (tmp.is_string ()) { @@ -73,10 +101,14 @@ { retval = octave_rand::distribution (); } - else if (s_arg == "seed" || s_arg == "state") + else if (s_arg == "seed") { retval = octave_rand::seed (); } + else if (s_arg == "state" || s_arg == "twister") + { + retval = octave_rand::state (); + } else if (s_arg == "uniform") { octave_rand::uniform_distribution (); @@ -85,6 +117,18 @@ { octave_rand::normal_distribution (); } + else if (s_arg == "exponential") + { + octave_rand::exponential_distribution (); + } + else if (s_arg == "poisson") + { + octave_rand::poisson_distribution (); + } + else if (s_arg == "gamma") + { + octave_rand::gamma_distribution (); + } else error ("%s: unrecognized string argument", fcn); } @@ -176,19 +220,27 @@ default: { - octave_value tmp = args(0); + octave_value tmp = args(idx); if (nargin == 2 && tmp.is_string ()) { std::string ts = tmp.string_value (); - if (ts == "seed" || ts == "state") + if (ts == "seed") { - double d = args(1).double_value (); + double d = args(idx+1).double_value (); if (! error_state) octave_rand::seed (d); } + else if (ts == "state" || ts == "twister") + { + ColumnVector s = + ColumnVector (args(idx+1).vector_value(false, true)); + + if (! error_state) + octave_rand::state (s); + } else error ("%s: unrecognized string argument", fcn); } @@ -198,7 +250,7 @@ for (int i = 0; i < nargin; i++) { - dims(i) = (octave_idx_type)args(i).int_value (); + dims(i) = (octave_idx_type)args(idx+i).int_value (); if (error_state) { @@ -221,33 +273,97 @@ dims.chop_trailing_singletons (); - return octave_rand::nd_array (dims); + if (additional_arg) + { + if (a.length() == 1) + return octave_rand::nd_array (dims, a(0)); + else + { + if (a.dims() != dims) + { + error ("%s: mismatch in argument size", fcn); + return retval; + } + octave_idx_type len = a.length (); + NDArray m (dims); + double *v = m.fortran_vec (); + for (octave_idx_type i = 0; i < len; i++) + v[i] = octave_rand::scalar (a(i)); + return m; + } + } + else + return octave_rand::nd_array (dims); } DEFUN_DLD (rand, args, , "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {} rand (@var{x})\n\ @deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\ -@deftypefnx {Loadable Function} {} rand (@code{\"seed\"}, @var{x})\n\ +@deftypefnx {Loadable Function} {} rand (\"state\", @var{x})\n\ +@deftypefnx {Loadable Function} {} rand (\"seed\", @var{x})\n\ Return a matrix with random elements uniformly distributed on the\n\ interval (0, 1). The arguments are handled the same as the arguments\n\ -for @code{eye}. In\n\ -addition, you can set the seed for the random number generator using the\n\ +for @code{eye}.\n\ +\n\ +You can query the state of the random number generator using the\n\ form\n\ \n\ @example\n\ -rand (\"seed\", @var{x})\n\ +v = rand (\"state\")\n\ +@end example\n\ +\n\ +This returns a column vector @var{v} of length 625. Later, you can\n\ +restore the random number generator to the state @var{v}\n\ +using the form\n\ +\n\ +@example\n\ +rand (\"state\", v)\n\ @end example\n\ \n\ @noindent\n\ -where @var{x} is a scalar value. If called as\n\ +You may also initialize the state vector from an arbitrary vector of\n\ +length <= 625 for @var{v}. This new state will be a hash based on the\n\ +the value of @var{v}, not @var{v} itself.\n\ +\n\ +By default, the generator is initialized from @code{/dev/urandom} if it is\n\ +available, otherwise from cpu time, wall clock time and the current\n\ +fraction of a second.\n\ +\n\ +@code{rand} uses the Mersenne Twister with a period of 2^19937-1\n\ +(See M. Matsumoto and T. Nishimura, ``Mersenne Twister: A 623-dimensionally\n\ +equidistributed uniform pseudorandom number generator'', ACM Trans. on\n\ +Modeling and Computer Simulation Vol. 8, No. 1, Januray pp.3-30 1998,\n\ +@url{http://www.math.keio.ac.jp/~matumoto/emt.html}).\n\ +Do NOT use for CRYPTOGRAPHY without securely hashing several returned\n\ +values together, otherwise the generator state can be learned after\n\ +reading 624 consecutive values.\n\ +\n\ +@code{rand} includes a second random number generator, that was the\n\ +previous generator used in octave. The new generator is used by default\n\ +as it is significantly faster than the old generator, and produces\n\ +random numebrs with a significantly longer cycle time. However, in\n\ +some circumstances it might be desireable to obtain the same random\n\ +sequences as used by the old generators. To do this the keyword\n\ +\"seed\" is used to specify that the old generators should be use,\n\ +as in\n\ \n\ @example\n\ -rand (\"seed\")\n\ +rand (\"seed\", val)\n\ @end example\n\ \n\ -@noindent\n\ -@code{rand} returns the current value of the seed.\n\ +which sets the seed of the generator to @var{val}. The seed of the\n\ +generator can be queried with\n\ +\n\ +@example\n\ +s = rand (\"seed\")\n\ +@end example\n\ +\n\ +However, it should be noted that querying the seed will not cause\n\ +@code{rand} to use the old generators, only setting the seed will.\n\ +To cause @code{rand} to once again use the new generators, the\n\ +keyword \"state\" should be used to reset the state of the @code{rand}.\n\ +@seealso{randn,rande,randg,randp}\n\ @end deftypefn") { octave_value retval; @@ -259,6 +375,64 @@ return retval; } +/* +%!test # 'state' can be a scalar +%! rand('state',12); x = rand(1,4); +%! rand('state',12); y = rand(1,4); +%! assert(x,y); +%!test # 'state' can be a vector +%! rand('state',[12,13]); x=rand(1,4); +%! rand('state',[12;13]); y=rand(1,4); +%! assert(x,y); +%!test # querying 'state' doesn't disturb sequence +%! rand('state',12); rand(1,2); x=rand(1,2); +%! rand('state',12); rand(1,2); +%! s=rand('state'); y=rand(1,2); +%! assert(x,y); +%! rand('state',s); z=rand(1,2); +%! assert(x,z); +%!test # 'seed' must be a scalar +%! rand('seed',12); x = rand(1,4); +%! rand('seed',12); y = rand(1,4); +%! assert(x,y); +%!error(rand('seed',[12,13])) +%!test # querying 'seed' returns a value which can be used later +%! s=rand('seed'); x=rand(1,2); +%! rand('seed',s); y=rand(1,2); +%! assert(x,y); +%!test # querying 'seed' doesn't disturb sequence +%! rand('seed',12); rand(1,2); x=rand(1,2); +%! rand('seed',12); rand(1,2); +%! s=rand('seed'); y=rand(1,2); +%! assert(x,y); +%! rand('seed',s); z=rand(1,2); +%! assert(x,z); +*/ + +/* +%!test +%! % statistical tests may fail occasionally. +%! rand("state",12); +%! x = rand(100000,1); +%! assert(max(x)<1.); %*** Please report this!!! *** +%! assert(min(x)>0.); %*** Please report this!!! *** +%! assert(mean(x),0.5,0.0024); +%! assert(var(x),1/48,0.0632); +%! assert(skewness(x),0,0.012); +%! assert(kurtosis(x),-6/5,0.0094); +%!test +%! % statistical tests may fail occasionally. +%! rand("seed",12); +%! x = rand(100000,1); +%! assert(max(x)<1.); %*** Please report this!!! *** +%! assert(min(x)>0.); %*** Please report this!!! *** +%! assert(mean(x),0.5,0.0024); +%! assert(var(x),1/48,0.0632); +%! assert(skewness(x),0,0.012); +%! assert(kurtosis(x),-6/5,0.0094); +*/ + + static std::string current_distribution = octave_rand::distribution (); static void @@ -271,25 +445,18 @@ "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {} randn (@var{x})\n\ @deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\ -@deftypefnx {Loadable Function} {} randn (@code{\"seed\"}, @var{x})\n\ -Return a matrix with normally distributed random elements. The\n\ -arguments are handled the same as the arguments for @code{eye}. In\n\ -addition, you can set the seed for the random number generator using the\n\ -form\n\ -\n\ -@example\n\ -randn (\"seed\", @var{x})\n\ -@end example\n\ +@deftypefnx {Loadable Function} {} randn (\"state\", @var{x})\n\ +@deftypefnx {Loadable Function} {} randn (\"seed\", @var{x})\n\ +Return a matrix with normally distributed random elements. The\n\ +arguments are handled the same as the arguments for @code{rand}.\n\ \n\ -@noindent\n\ -where @var{x} is a scalar value. If called as\n\ +By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ +transform from a uniform to a normal distribution. (G. Marsaglia and\n\ +W.W. Tsang, 'Ziggurat method for generating random variables',\n\ +J. Statistical Software, vol 5, 2000,\n\ +@url{http://www.jstatsoft.org/v05/i08/})\n\ \n\ -@example\n\ -randn (\"seed\")\n\ -@end example\n\ -\n\ -@noindent\n\ -@code{randn} returns the current value of the seed.\n\ +@seealso{rand,rande,randg,randp}\n\ @end deftypefn") { octave_value retval; @@ -318,6 +485,369 @@ } /* +%!test +%! % statistical tests may fail occasionally. +%! rand("state",12); +%! x = randn(100000,1); +%! assert(mean(x),0,0.01); +%! assert(var(x),1,0.02); +%! assert(skewness(x),0,0.02); +%! assert(kurtosis(x),0,0.04); +%!test +%! % statistical tests may fail occasionally. +%! rand("seed",12); +%! x = randn(100000,1); +%! assert(mean(x),0,0.01); +%! assert(var(x),1,0.02); +%! assert(skewness(x),0,0.02); +%! assert(kurtosis(x),0,0.04); +*/ + +DEFUN_DLD (rande, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} rande (@var{x})\n\ +@deftypefnx {Loadable Function} {} rande (@var{n}, @var{m})\n\ +@deftypefnx {Loadable Function} {} rande (\"state\", @var{x})\n\ +@deftypefnx {Loadable Function} {} rande (\"seed\", @var{x})\n\ +Return a matrix with exponentially distributed random elements. The\n\ +arguments are handled the same as the arguments for @code{rand}.\n\ +\n\ +By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ +transform from a uniform to a exponential distribution. (G. Marsaglia and\n\ +W.W. Tsang, 'Ziggurat method for generating random variables',\n\ +J. Statistical Software, vol 5, 2000,\n\ +@url{http://www.jstatsoft.org/v05/i08/})\n\ +@seealso{rand,randn,randg,randp}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + unwind_protect::begin_frame ("rande"); + + // This relies on the fact that elements are popped from the unwind + // stack in the reverse of the order they are pushed + // (i.e. current_distribution will be reset before calling + // reset_rand_generator()). + + unwind_protect::add (reset_rand_generator, 0); + unwind_protect_str (current_distribution); + + current_distribution = "exponential"; + + octave_rand::distribution (current_distribution); + + retval = do_rand (args, nargin, "rande"); + + unwind_protect::run_frame ("rande"); + + return retval; +} + +/* +%!test +%! % statistical tests may fail occasionally +%! rand("state",12); +%! x = rande(100000,1); +%! assert(min(x)>0); % *** Please report this!!! *** +%! assert(mean(x),1,0.01); +%! assert(var(x),1,0.03); +%! assert(skewness(x),2,0.06); +%! assert(kurtosis(x),6,0.7); +%!test +%! % statistical tests may fail occasionally +%! rand("seed",12); +%! x = rande(100000,1); +%! assert(min(x)>0); % *** Please report this!!! *** +%! assert(mean(x),1,0.01); +%! assert(var(x),1,0.03); +%! assert(skewness(x),2,0.06); +%! assert(kurtosis(x),6,0.7); +*/ + +DEFUN_DLD (randg, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} randg (@var{a}, @var{x})\n\ +@deftypefnx {Loadable Function} {} randg (@var{a}, @var{n}, @var{m})\n\ +@deftypefnx {Loadable Function} {} randg (\"state\", @var{x})\n\ +@deftypefnx {Loadable Function} {} randg (\"seed\", @var{x})\n\ +Return a matrix with @code{gamma(@var{a},1)} distributed random elements.\n\ +The arguments are handled the same as the arguments for @code{rand},\n\ +except for the argument @var{a}.\n\ +\n\ +This can be used to generate many distributions:\n\ +\n\ +@table @asis\n\ +@item @code{gamma (a,b)} for @code{a > -1}, @code{b > 0}\n\ +@example\n\ +r = b*randg(a)\n\ +@end example\n\ +@item @code{beta(a,b)} for @code{a > -1}, @code{b > -1}\n\ +@example\n\ +r1 = randg(a,1)\n\ +r = r1 / (r1 + randg(b,1))\n\ +@end example\n\ +@item @code{Erlang(a, n)}\n\ +@example\n\ +r = a*randg(n)\n\ +@end example\n\ +@item @code{chisq(df)} for @code{df > 0}\n\ +@example\n\ +r = 2*randg(df/2)\n\ +@end example\n\ +@item @code{t(df)} for @code{0 < df < inf} (use randn if df is infinite)\n\ +@example\n\ +r = randn() / sqrt(2*randg(df/2)/df)\n\ +@end example\n\ +@item @code{F(n1,n2)} for @code{0 < n1}, @code{0 < n2}\n\ +@example\n\ +r1 = 2*randg(n1/2)/n1 or 1 if n1 is infinite\n\ +r2 = 2*randg(n2/2)/n2 or 1 if n2 is infinite\n\ +r = r1 / r2\n\n\ +@end example\n\ +@item negative @code{binomial (n, p)} for @code{n > 0}, @code{0 < p <= 1}\n\ +@example\n\ +r = randp((1-p)/p * randg(n))\n\ +@end example\n\ +@item non-central @code{chisq(df,L)}, for @code{df >= 0} and @code{L > 0}\n\ +(use chisq if @code{L = 0})\n\ +@example\n\ +r = randp(L/2)\n\ +r(r > 0) = 2*randg(r(r > 0))\n\ +r(df > 0) += 2*randg(df(df > 0)/2)\n\ +@end example\n\ +@item @code{Dirichlet(a1,...,ak)}\n\ +@example\n\ +r = (randg(a1),...,randg(ak))\n\ +r = r / sum(r)\n\ +@end example\n\ +@end table\n\ +@seealso{rand,randn,rande,randp}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin < 1) + error ("randg: insufficient arguments"); + else + { + unwind_protect::begin_frame ("randg"); + + // This relies on the fact that elements are popped from the unwind + // stack in the reverse of the order they are pushed + // (i.e. current_distribution will be reset before calling + // reset_rand_generator()). + + unwind_protect::add (reset_rand_generator, 0); + unwind_protect_str (current_distribution); + + current_distribution = "gamma"; + + octave_rand::distribution (current_distribution); + + retval = do_rand (args, nargin, "randg", true); + + unwind_protect::run_frame ("randg"); + } + + return retval; +} + +/* +%!test +%! rand("state",12) +%!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report +%!test +%! % statistical tests may fail occasionally. +%! a=0.1; x = randg(a,100000,1); +%! assert(mean(x), a, 0.01); +%! assert(var(x), a, 0.01); +%! assert(skewness(x),2/sqrt(a), 1.); +%! assert(kurtosis(x),6/a, 50.); +%!test +%! % statistical tests may fail occasionally. +%! a=0.95; x = randg(a,100000,1); +%! assert(mean(x), a, 0.01); +%! assert(var(x), a, 0.04); +%! assert(skewness(x),2/sqrt(a), 0.2); +%! assert(kurtosis(x),6/a, 2.); +%!test +%! % statistical tests may fail occasionally. +%! a=1; x = randg(a,100000,1); +%! assert(mean(x),a, 0.01); +%! assert(var(x),a, 0.04); +%! assert(skewness(x),2/sqrt(a), 0.2); +%! assert(kurtosis(x),6/a, 2.); +%!test +%! % statistical tests may fail occasionally. +%! a=10; x = randg(a,100000,1); +%! assert(mean(x), a, 0.1); +%! assert(var(x), a, 0.5); +%! assert(skewness(x),2/sqrt(a), 0.1); +%! assert(kurtosis(x),6/a, 0.5); +%!test +%! % statistical tests may fail occasionally. +%! a=100; x = randg(a,100000,1); +%! assert(mean(x), a, 0.2); +%! assert(var(x), a, 2.); +%! assert(skewness(x),2/sqrt(a), 0.05); +%! assert(kurtosis(x),6/a, 0.2); +%!test +%! rand("seed",12) +%!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report +%!test +%! % statistical tests may fail occasionally. +%! a=0.1; x = randg(a,100000,1); +%! assert(mean(x), a, 0.01); +%! assert(var(x), a, 0.01); +%! assert(skewness(x),2/sqrt(a), 1.); +%! assert(kurtosis(x),6/a, 50.); +%!test +%! % statistical tests may fail occasionally. +%! a=0.95; x = randg(a,100000,1); +%! assert(mean(x), a, 0.01); +%! assert(var(x), a, 0.04); +%! assert(skewness(x),2/sqrt(a), 0.2); +%! assert(kurtosis(x),6/a, 2.); +%!test +%! % statistical tests may fail occasionally. +%! a=1; x = randg(a,100000,1); +%! assert(mean(x),a, 0.01); +%! assert(var(x),a, 0.04); +%! assert(skewness(x),2/sqrt(a), 0.2); +%! assert(kurtosis(x),6/a, 2.); +%!test +%! % statistical tests may fail occasionally. +%! a=10; x = randg(a,100000,1); +%! assert(mean(x), a, 0.1); +%! assert(var(x), a, 0.5); +%! assert(skewness(x),2/sqrt(a), 0.1); +%! assert(kurtosis(x),6/a, 0.5); +%!test +%! % statistical tests may fail occasionally. +%! a=100; x = randg(a,100000,1); +%! assert(mean(x), a, 0.2); +%! assert(var(x), a, 2.); +%! assert(skewness(x),2/sqrt(a), 0.05); +%! assert(kurtosis(x),6/a, 0.2); +*/ + + +DEFUN_DLD (randp, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} randp (@var{l}, @var{x})\n\ +@deftypefnx {Loadable Function} {} randp (@var{l}, @var{n}, @var{m})\n\ +@deftypefnx {Loadable Function} {} randp (\"state\", @var{x})\n\ +@deftypefnx {Loadable Function} {} randp (\"seed\", @var{x})\n\ +Return a matrix with Poisson distributed random elements. The arguments\n\ +are handled the same as the arguments for @code{rand}, except for the\n\ +argument @var{l}.\n\ +\n\ +Five different algorithms are used depending on the range of @var{l}\n\ +and whether or not @var{l} is a scalar or a matrix.\n\ +\n\ +@table @asis\n\ +@item For scalar @var{l} <= 12, use direct method.\n\ +Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ +@item For scalar @var{l} > 12, use rejection method.[1]\n\ +Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ +@item For matrix @var{l} <= 10, use inversion method.[2]\n\ +Stadlober E., et al., WinRand source code, available via FTP.\n\ +@item For matrix @var{l} > 10, use patchwork rejection method.\n\ +Stadlober E., et al., WinRand source code, available via FTP, or\n\ +H. Zechner, 'Efficient sampling from continuous and discrete\n\ +unimodal distributions', Doctoral Dissertaion, 156pp., Technical\n\ +University Graz, Austria, 1994.\n\ +@item For @var{l} > 1e8, use normal approximation.\n\ +L. Montanet, et al., 'Review of Particle Properties', Physical Review\n\ +D 50 p1284, 1994\n\ +@end table\n\ +@seealso{rand,randn,rande,randg}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin < 1) + error ("randp: insufficient arguments"); + else + { + unwind_protect::begin_frame ("randp"); + + // This relies on the fact that elements are popped from the unwind + // stack in the reverse of the order they are pushed + // (i.e. current_distribution will be reset before calling + // reset_rand_generator()). + + unwind_protect::add (reset_rand_generator, 0); + unwind_protect_str (current_distribution); + + current_distribution = "poisson"; + + octave_rand::distribution (current_distribution); + + retval = do_rand (args, nargin, "randp", true); + + unwind_protect::run_frame ("randp"); + } + + return retval; +} + +/* +%!test +%! rand("state",12) +%!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report +%!test +%! % statistical tests may fail occasionally. +%! for a=[5 15] +%! x = randp(a,100000,1); +%! assert(min(x)>=0); % *** Please report this!!! *** +%! assert(mean(x),a,0.03); +%! assert(var(x),a,0.2); +%! assert(skewness(x),1/sqrt(a),0.03); +%! assert(kurtosis(x),1/a,0.08); +%! end +%!test +%! % statistical tests may fail occasionally. +%! for a=[5 15] +%! x = randp(a*ones(100000,1),100000,1); +%! assert(min(x)>=0); % *** Please report this!!! *** +%! assert(mean(x),a,0.03); +%! assert(var(x),a,0.2); +%! assert(skewness(x),1/sqrt(a),0.03); +%! assert(kurtosis(x),1/a,0.08); +%! end +%!test +%! rand("seed",12) +%!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report +%!test +%! % statistical tests may fail occasionally. +%! for a=[5 15] +%! x = randp(a,100000,1); +%! assert(min(x)>=0); % *** Please report this!!! *** +%! assert(mean(x),a,0.03); +%! assert(var(x),a,0.2); +%! assert(skewness(x),1/sqrt(a),0.03); +%! assert(kurtosis(x),1/a,0.08); +%! end +%!test +%! % statistical tests may fail occasionally. +%! for a=[5 15] +%! x = randp(a*ones(100000,1),100000,1); +%! assert(min(x)>=0); % *** Please report this!!! *** +%! assert(mean(x),a,0.03); +%! assert(var(x),a,0.2); +%! assert(skewness(x),1/sqrt(a),0.03); +%! assert(kurtosis(x),1/a,0.08); +%! end +*/ + +/* ;;; Local Variables: *** ;;; mode: C++ *** ;;; End: ***