Difference between Matlab and Octave in random seeding

Hi,

In octave (3.0.0)
rand ('state', 1); randn
will always produce the same number. I.e. the seed on rand affects randn as
well.

In  Matlab seeding rand does not seed randn.




--
Mark van Rossum
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[John W. Eaton]

On 15-Feb-2008, Mark van Rossum wrote:

| In octave (3.0.0)
| rand ('state', 1); randn
| will always produce the same number. I.e. the seed on rand affects randn as
| well.
|
| In  Matlab seeding rand does not seed randn.

Please try the following patch.

Thanks,

jwe



# HG changeset patch
# User John W. Eaton <jwe@...>
# Date 1204021739 18000
# Node ID ff52243af934e339f4bef0d04e68779731be66d2
# Parent  493bb0de319977ab89626f0752e7c63ea7b24da1
save state separately for each MT random number generator

diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,4 +1,23 @@ 2008-02-26  John W. Eaton  <jwe@...
 2008-02-26  John W. Eaton  <jwe@...>
+
+ * oct-rand.cc (rand_states): New static variable.
+ (initialize_rand_states, get_dist_id, get_internal_state,
+ set_internal_state, switch_to_generator, save_state): New functions.
+ (octave_rand::state): New arg to specify distribution.
+ Save state in rand_states instead of setting internal state.
+ Return named state.  Use set_internal_state to generate proper
+ state vector from user supplied state.  Save and restore current
+ state if specified and current distributions are different.
+ (octave_rand::distribution (void)): Use switch rather than if/else.
+ (octave_rand::distribution (const std::string&)): Likewise.
+ (octave_rand::uniform_distribution,
+ octave_rand::normal_distribution,
+ octave_rand::exponential_distribution,
+ octave_rand::poisson_distribution,
+ octave_rand::gamma_distribution): Call switch_to_generator.
+ (octave_rand::state, maybe_initialize): For new_generators, just
+ call initialize_rand_states if not already initialized.
+ (octave_rand::scalar, fill_rand): Save state after generating value.
 
  * dMatrix.cc (Matrix::lssolve): Avoid another dgelsd lwork query bug.
  * CMatrix.cc (ComplexMatrix::lssolve): Likewise, for zgelsd
diff --git a/liboctave/oct-rand.cc b/liboctave/oct-rand.cc
--- a/liboctave/oct-rand.cc
+++ b/liboctave/oct-rand.cc
@@ -23,6 +23,8 @@ along with Octave; see the file COPYING.
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
+
+#include <map>
 #include <vector>
 
 #include "f77-fcn.h"
@@ -52,6 +54,8 @@ static bool old_initialized = false;
 static bool old_initialized = false;
 static bool new_initialized = false;
 static bool use_old_generators = false;
+
+std::map<int, ColumnVector> rand_states;
 
 extern "C"
 {
@@ -126,6 +130,46 @@ do_old_initialization (void)
   old_initialized = true;
 }
 
+static ColumnVector
+get_internal_state (void)
+{
+  ColumnVector s (MT_N + 1);
+
+  OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1);
+
+  oct_get_state (tmp);
+
+  for (octave_idx_type i = 0; i <= MT_N; i++)
+    s.elem (i) = static_cast<double> (tmp [i]);
+
+  return s;
+}
+
+static inline void
+save_state (void)
+{
+  rand_states[current_distribution] = get_internal_state ();;
+}
+
+static void
+initialize_rand_states (void)
+{
+  if (! new_initialized)
+    {
+      oct_init_by_entropy ();
+
+      ColumnVector s = get_internal_state ();
+
+      rand_states[uniform_dist] = s;
+      rand_states[normal_dist] = s;
+      rand_states[expon_dist] = s;
+      rand_states[poisson_dist] = s;
+      rand_states[gamma_dist] = s;
+
+      new_initialized = true;
+    }
+}
+
 static inline void
 maybe_initialize (void)
 {
@@ -137,10 +181,56 @@ maybe_initialize (void)
   else
     {
       if (! new_initialized)
- {
-  oct_init_by_entropy ();
-  new_initialized = true;
- }
+ initialize_rand_states ();
+    }
+}
+
+static int
+get_dist_id (const std::string& d)
+{
+  int retval;
+
+  if (d == "uniform" || d == "rand")
+    retval = uniform_dist;
+  else if (d == "normal" || d == "randn")
+    retval = normal_dist;
+  else if (d == "exponential" || d == "rande")
+    retval = expon_dist;
+  else if (d == "poisson" || d == "randp")
+    retval = poisson_dist;
+  else if (d == "gamma" || d == "rangd")
+    retval = gamma_dist;
+  else
+    (*current_liboctave_error_handler) ("rand: invalid distribution");
+
+  return retval;
+}
+
+static void
+set_internal_state (const ColumnVector& s)
+{
+  octave_idx_type len = s.length ();
+  octave_idx_type n = len < MT_N + 1 ? len : MT_N + 1;
+
+  OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1);
+
+  for (octave_idx_type i = 0; i < n; i++)
+    tmp[i] = static_cast<uint32_t> (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);
+}
+
+static inline void
+switch_to_generator (int dist)
+{
+  if (dist != current_distribution)
+    {
+      current_distribution = dist;
+
+      set_internal_state (rand_states[dist]);
     }
 }
 
@@ -172,6 +262,7 @@ octave_rand::seed (double s)
 octave_rand::seed (double s)
 {
   use_old_generators = true;
+
   maybe_initialize ();
 
   int i0, i1;
@@ -197,77 +288,104 @@ octave_rand::seed (double s)
 }
 
 ColumnVector
-octave_rand::state (void)
+octave_rand::state (const std::string& d)
 {
-  ColumnVector s (MT_N + 1);
   if (! new_initialized)
-    {
-      oct_init_by_entropy ();
-      new_initialized = true;
-    }
+    initialize_rand_states ();
 
-  OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1);
-  oct_get_state (tmp);
-  for (octave_idx_type i = 0; i <= MT_N; i++)
-    s.elem (i) = static_cast<double>(tmp [i]);
-  return s;
+  return rand_states[d.empty () ? current_distribution : get_dist_id (d)];
 }
 
 void
-octave_rand::state (const ColumnVector &s)
+octave_rand::state (const ColumnVector& s, const std::string& d)
 {
   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 (uint32_t, tmp, MT_N + 1);
-  for (octave_idx_type i = 0; i < n; i++)
-    tmp[i] = static_cast<uint32_t> (s.elem(i));
+  int old_dist = current_distribution;
 
-  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);
+  int new_dist = d.empty () ? current_distribution : get_dist_id (d);
+
+  ColumnVector saved_state;
+
+  if (old_dist != new_dist)
+    saved_state = get_internal_state ();
+
+  set_internal_state (s);
+
+  rand_states[new_dist] = get_internal_state ();
+
+  if (old_dist != new_dist)
+    rand_states[old_dist] = saved_state;
 }
 
 std::string
 octave_rand::distribution (void)
 {
+  std::string retval;
+
   maybe_initialize ();
 
-  if (current_distribution == uniform_dist)
-    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
+  switch (current_distribution)
     {
-      abort ();
-      return "";
+    case uniform_dist:
+      retval = "uniform";
+      break;
+
+    case normal_dist:
+      retval = "normal";
+      break;
+
+    case expon_dist:
+      retval = "exponential";
+      break;
+
+    case poisson_dist:
+      retval = "poisson";
+      break;
+
+    case gamma_dist:
+      retval = "gamma";
+      break;
+
+    default:
+      (*current_liboctave_error_handler) ("rand: invalid distribution");
+      break;
     }
+
+  return retval;
 }
 
 void
 octave_rand::distribution (const std::string& d)
 {
-  if (d == "uniform")
-    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");
+  switch (get_dist_id (d))
+    {
+    case uniform_dist:
+      octave_rand::uniform_distribution ();
+      break;
+
+    case normal_dist:
+      octave_rand::normal_distribution ();
+      break;
+
+    case expon_dist:
+      octave_rand::exponential_distribution ();
+      break;
+
+    case poisson_dist:
+      octave_rand::poisson_distribution ();
+      break;
+
+    case gamma_dist:
+      octave_rand::gamma_distribution ();
+      break;
+
+    default:
+      (*current_liboctave_error_handler) ("rand: invalid distribution");
+      break;
+    }
 }
 
 void
@@ -275,7 +393,7 @@ octave_rand::uniform_distribution (void)
 {
   maybe_initialize ();
 
-  current_distribution = uniform_dist;
+  switch_to_generator (uniform_dist);
 
   F77_FUNC (setcgn, SETCGN) (uniform_dist);
 }
@@ -285,7 +403,7 @@ octave_rand::normal_distribution (void)
 {
   maybe_initialize ();
 
-  current_distribution = normal_dist;
+  switch_to_generator (normal_dist);
 
   F77_FUNC (setcgn, SETCGN) (normal_dist);
 }
@@ -295,7 +413,7 @@ octave_rand::exponential_distribution (v
 {
   maybe_initialize ();
 
-  current_distribution = expon_dist;
+  switch_to_generator (expon_dist);
 
   F77_FUNC (setcgn, SETCGN) (expon_dist);
 }
@@ -305,7 +423,7 @@ octave_rand::poisson_distribution (void)
 {
   maybe_initialize ();
 
-  current_distribution = poisson_dist;
+  switch_to_generator (poisson_dist);
 
   F77_FUNC (setcgn, SETCGN) (poisson_dist);
 }
@@ -315,7 +433,7 @@ octave_rand::gamma_distribution (void)
 {
   maybe_initialize ();
 
-  current_distribution = gamma_dist;
+  switch_to_generator (gamma_dist);
 
   F77_FUNC (setcgn, SETCGN) (gamma_dist);
 }
@@ -363,7 +481,7 @@ octave_rand::scalar (double a)
   break;
 
  default:
-  abort ();
+  (*current_liboctave_error_handler) ("rand: invalid distribution");
   break;
  }
     }
@@ -372,29 +490,31 @@ octave_rand::scalar (double a)
       switch (current_distribution)
  {
  case uniform_dist:
-  retval = oct_randu();
+  retval = oct_randu ();
   break;
 
  case normal_dist:
-  retval = oct_randn();
+  retval = oct_randn ();
   break;
 
  case expon_dist:
-  retval = oct_rande();
+  retval = oct_rande ();
   break;
 
  case poisson_dist:
-  retval = oct_randp(a);
+  retval = oct_randp (a);
   break;
 
  case gamma_dist:
-  retval = oct_randg(a);
+  retval = oct_randg (a);
   break;
 
  default:
-  abort ();
+  (*current_liboctave_error_handler) ("rand: invalid distribution");
   break;
  }
+
+      save_state ();
     }
 
   return retval;
@@ -494,9 +614,11 @@ fill_rand (octave_idx_type len, double *
       break;
 
     default:
-      abort ();
+      (*current_liboctave_error_handler) ("rand: invalid distribution");
       break;
     }
+
+  save_state ();
 
   return;
 }
diff --git a/liboctave/oct-rand.h b/liboctave/oct-rand.h
--- a/liboctave/oct-rand.h
+++ b/liboctave/oct-rand.h
@@ -40,16 +40,17 @@ octave_rand
   static void seed (double s);
 
   // Return the current state.
-  static ColumnVector state (void);
+  static ColumnVector state (const std::string& d = std::string ());
 
   // Set the current state/
-  static void state (const ColumnVector &s);
+  static void state (const ColumnVector &s,
+     const std::string& d = std::string ());
   
   // Return the current distribution.
   static std::string distribution (void);
 
   // Set the current distribution.  May be either "uniform" (the
-  // default) or "normal".
+  // default), "normal", "exponential", "poisson", or "gamma".
   static void distribution (const std::string& d);
 
   static void uniform_distribution (void);
diff --git a/liboctave/randmtzig.c b/liboctave/randmtzig.c
--- a/liboctave/randmtzig.c
+++ b/liboctave/randmtzig.c
@@ -203,7 +203,7 @@ oct_init_by_int (uint32_t s)
 /* init_key is the array for initializing keys */
 /* key_length is its length */
 void
-oct_init_by_array (uint32_t init_key[], int key_length)
+oct_init_by_array (uint32_t *init_key, int key_length)
 {
   int i, j, k;
   oct_init_by_int (19650218UL);
@@ -281,17 +281,17 @@ oct_init_by_entropy (void)
 }
 
 void
-oct_set_state (uint32_t save[])
+oct_set_state (uint32_t *save)
 {
   int i;
-  for (i=0; i < MT_N; i++)
+  for (i = 0; i < MT_N; i++)
     state[i] = save[i];
   left = save[MT_N];
   next = state + (MT_N - left + 1);
 }
 
 void
-oct_get_state (uint32_t save[])
+oct_get_state (uint32_t *save)
 {
   int i;
   for (i = 0; i < MT_N; i++)
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,4 +1,7 @@ 2008-02-26  John W. Eaton  <jwe@...
 2008-02-26  John W. Eaton  <jwe@...>
+
+ * DLD-FUNCTIONS/rand.cc (do_rand): Pass name of calling function
+ to octave_rand::state.
 
  * variables.cc (bind_ans): Handle cs-lists recursively.
 
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
@@ -113,7 +113,7 @@ do_rand (const octave_value_list& args,
       }
     else if (s_arg == "state" || s_arg == "twister")
       {
- retval = octave_rand::state ();
+ retval = octave_rand::state (fcn);
       }
     else if (s_arg == "uniform")
       {
@@ -250,7 +250,7 @@ do_rand (const octave_value_list& args,
   ColumnVector (args(idx+1).vector_value(false, true));
 
  if (! error_state)
-  octave_rand::state (s);
+  octave_rand::state (s, fcn);
       }
     else
       error ("%s: unrecognized string argument", fcn);

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[dbateman]

John W. Eaton wrote:

On 15-Feb-2008, Mark van Rossum wrote:

| In octave (3.0.0)
| rand ('state', 1); randn
| will always produce the same number. I.e. the seed on rand affects randn as
| well.
|
| In  Matlab seeding rand does not seed randn.

Please try the following patch.

Thanks,

jwe

I asked at the time that the Mersenne twister code was included in Octave whether we wanted independent generator states per generator or not, and if I remember correctly it didn't seem important. Note also that there are a number of additional random generator functions that are derived from the base generator in the statistics package, and so the state of these generators will depend on the underlying generator functions. I don't have access to matlab and so can't check, but how do things like exprnd depend on the state of the rand and randn  generators in Matlab? I suspect that they rely on the same underlying generators. However as Octave has a ziggurat implementation of rande that is also used for exprnd, then this change will mean that Octave is now different than Matlab in a different manner. Should we care as this will only affect the users of the matlab statistics toolbox..

D.

[John W. Eaton]

On  1-Mar-2008, dbateman wrote:

| I asked at the time that the Mersenne twister code was included in Octave
| whether we wanted independent generator states per generator or not, and if
| I remember correctly it didn't seem important. Note also that there are a
| number of additional random generator functions that are derived from the
| base generator in the statistics package, and so the state of these
| generators will depend on the underlying generator functions. I don't have
| access to matlab and so can't check, but how do things like exprnd depend on
| the state of the rand and randn  generators in Matlab? I suspect that they
| rely on the same underlying generators. However as Octave has a ziggurat
| implementation of rande that is also used for exprnd, then this change will
| mean that Octave is now different than Matlab in a different manner. Should
| we care as this will only affect the users of the matlab statistics
| toolbox..

OK, I don't know what is best here.  Would someone (not David) please
check the Matlab behavior of exprnd, rand, and randn when setting the
state?  If you set the state for rand or randn, is the exprnd sequence
also reinitialized, or is it independent?

Thanks,

jwe
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[Jaroslav Hajek-2]

On Mon, Mar 3, 2008 at 8:55 AM, John W. Eaton <jwe@...> wrote:

> On  1-Mar-2008, dbateman wrote:
>
>  | I asked at the time that the Mersenne twister code was included in Octave
>  | whether we wanted independent generator states per generator or not, and if
>  | I remember correctly it didn't seem important. Note also that there are a
>  | number of additional random generator functions that are derived from the
>  | base generator in the statistics package, and so the state of these
>  | generators will depend on the underlying generator functions. I don't have
>  | access to matlab and so can't check, but how do things like exprnd depend on
>  | the state of the rand and randn  generators in Matlab? I suspect that they
>  | rely on the same underlying generators. However as Octave has a ziggurat
>  | implementation of rande that is also used for exprnd, then this change will
>  | mean that Octave is now different than Matlab in a different manner. Should
>  | we care as this will only affect the users of the matlab statistics
>  | toolbox..
>
>  OK, I don't know what is best here.  Would someone (not David) please
>  check the Matlab behavior of exprnd, rand, and randn when setting the
>  state?  If you set the state for rand or randn, is the exprnd sequence
>  also reinitialized, or is it independent?
>

In matlab, rand and randn seem to have independent internals, but the
rest like exprnd, gamrnd are implemented using calls to rand and/or
randn.

>  Thanks,
>
>  jwe
>
>
> _______________________________________________
>  Bug-octave mailing list
>  Bug-octave@...
>  https://www.cae.wisc.edu/mailman/listinfo/bug-octave
>

--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[John W. Eaton]

On  3-Mar-2008, Jaroslav Hajek wrote:

| On Mon, Mar 3, 2008 at 8:55 AM, John W. Eaton <jwe@...> wrote:
| > On  1-Mar-2008, dbateman wrote:
| >
| >  | I asked at the time that the Mersenne twister code was included in Octave
| >  | whether we wanted independent generator states per generator or not, and if
| >  | I remember correctly it didn't seem important. Note also that there are a
| >  | number of additional random generator functions that are derived from the
| >  | base generator in the statistics package, and so the state of these
| >  | generators will depend on the underlying generator functions. I don't have
| >  | access to matlab and so can't check, but how do things like exprnd depend on
| >  | the state of the rand and randn  generators in Matlab? I suspect that they
| >  | rely on the same underlying generators. However as Octave has a ziggurat
| >  | implementation of rande that is also used for exprnd, then this change will
| >  | mean that Octave is now different than Matlab in a different manner. Should
| >  | we care as this will only affect the users of the matlab statistics
| >  | toolbox..
| >
| >  OK, I don't know what is best here.  Would someone (not David) please
| >  check the Matlab behavior of exprnd, rand, and randn when setting the
| >  state?  If you set the state for rand or randn, is the exprnd sequence
| >  also reinitialized, or is it independent?
| >
|
| In matlab, rand and randn seem to have independent internals, but the
| rest like exprnd, gamrnd are implemented using calls to rand and/or
| randn.

OK, again, I'm not sure what is best.  We could make the states
different for rand and randn, and have the others share the state with
either rand or randn, but that doesn't solve the problem if rand and
randn are both used to generate some other rand sequence (I don't know
whether that ever happens).  In any case, isn't it better that the
states can be set independently?

jwe
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[Jaroslav Hajek-2]

On Mon, Mar 3, 2008 at 9:32 AM, John W. Eaton <jwe@...> wrote:

> On  3-Mar-2008, Jaroslav Hajek wrote:
>
>
> | On Mon, Mar 3, 2008 at 8:55 AM, John W. Eaton <jwe@...> wrote:
>  | > On  1-Mar-2008, dbateman wrote:
>  | >
>  | >  | I asked at the time that the Mersenne twister code was included in Octave
>  | >  | whether we wanted independent generator states per generator or not, and if
>  | >  | I remember correctly it didn't seem important. Note also that there are a
>  | >  | number of additional random generator functions that are derived from the
>  | >  | base generator in the statistics package, and so the state of these
>  | >  | generators will depend on the underlying generator functions. I don't have
>  | >  | access to matlab and so can't check, but how do things like exprnd depend on
>  | >  | the state of the rand and randn  generators in Matlab? I suspect that they
>  | >  | rely on the same underlying generators. However as Octave has a ziggurat
>  | >  | implementation of rande that is also used for exprnd, then this change will
>  | >  | mean that Octave is now different than Matlab in a different manner. Should
>  | >  | we care as this will only affect the users of the matlab statistics
>  | >  | toolbox..
>  | >
>  | >  OK, I don't know what is best here.  Would someone (not David) please
>  | >  check the Matlab behavior of exprnd, rand, and randn when setting the
>  | >  state?  If you set the state for rand or randn, is the exprnd sequence
>  | >  also reinitialized, or is it independent?
>  | >
>  |
>  | In matlab, rand and randn seem to have independent internals, but the
>  | rest like exprnd, gamrnd are implemented using calls to rand and/or
>  | randn.
>
>  OK, again, I'm not sure what is best.  We could make the states
>  different for rand and randn, and have the others share the state with
>  either rand or randn, but that doesn't solve the problem if rand and
>  randn are both used to generate some other rand sequence (I don't know
>  whether that ever happens).  In any case, isn't it better that the
>  states can be set independently?
>

I guess that depends. Personally, I'd prefer to have a single seed for
all random numbers, because that's what I would expect, and that's
also how I would write a new one - calling the existing ones. The
independence of rand and randn in Matlab just surprised me.

>  jwe
>



--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[David Bateman]

Jaroslav Hajek wrote:

> On Mon, Mar 3, 2008 at 9:32 AM, John W. Eaton <jwe@...> wrote:
>  
>> On  3-Mar-2008, Jaroslav Hajek wrote:
>>
>>
>> | On Mon, Mar 3, 2008 at 8:55 AM, John W. Eaton <jwe@...> wrote:
>>  | > On  1-Mar-2008, dbateman wrote:
>>  | >
>>  | >  | I asked at the time that the Mersenne twister code was included in Octave
>>  | >  | whether we wanted independent generator states per generator or not, and if
>>  | >  | I remember correctly it didn't seem important. Note also that there are a
>>  | >  | number of additional random generator functions that are derived from the
>>  | >  | base generator in the statistics package, and so the state of these
>>  | >  | generators will depend on the underlying generator functions. I don't have
>>  | >  | access to matlab and so can't check, but how do things like exprnd depend on
>>  | >  | the state of the rand and randn  generators in Matlab? I suspect that they
>>  | >  | rely on the same underlying generators. However as Octave has a ziggurat
>>  | >  | implementation of rande that is also used for exprnd, then this change will
>>  | >  | mean that Octave is now different than Matlab in a different manner. Should
>>  | >  | we care as this will only affect the users of the matlab statistics
>>  | >  | toolbox..
>>  | >
>>  | >  OK, I don't know what is best here.  Would someone (not David) please
>>  | >  check the Matlab behavior of exprnd, rand, and randn when setting the
>>  | >  state?  If you set the state for rand or randn, is the exprnd sequence
>>  | >  also reinitialized, or is it independent?
>>  | >
>>  |
>>  | In matlab, rand and randn seem to have independent internals, but the
>>  | rest like exprnd, gamrnd are implemented using calls to rand and/or
>>  | randn.
>>
>>  OK, again, I'm not sure what is best.  We could make the states
>>  different for rand and randn, and have the others share the state with
>>  either rand or randn, but that doesn't solve the problem if rand and
>>  randn are both used to generate some other rand sequence (I don't know
>>  whether that ever happens).  In any case, isn't it better that the
>>  states can be set independently?
>>
>>    
>
> I guess that depends. Personally, I'd prefer to have a single seed for
> all random numbers, because that's what I would expect, and that's
> also how I would write a new one - calling the existing ones. The
> independence of rand and randn in Matlab just surprised me.
>
>  

I don't think there is anyway to address this if we keep rande, randg
and randp functions in Octave. I'd hate to loose these as they are so
much faster than depending on a derivation from rand/randn. I'd say
leave it as is...

D.

--
David Bateman                                David.Bateman@...
Motorola Labs - Paris                        [phone] (Ph)
Parc Les Algorithmes, Commune de St Aubin    [phone] (Mob)
91193 Gif-Sur-Yvette FRANCE                  [phone] (Fax)

The information contained in this communication has been classified as:

[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[Michael Creel]

On this issue, I vote for having the seed set in a single spot for all RNGs. This facilitates greatly parallel use of Octave. Sometimes it's important to be able to reset the seed, and it's a lot less likely to screw up if only one needs to be reset.
Best, Michael

[John W. Eaton]

On  5-Mar-2008, Michael Creel wrote:

| On this issue, I vote for having the seed set in a single spot for all RNGs.
| This facilitates greatly parallel use of Octave. Sometimes it's important to
| be able to reset the seed, and it's a lot less likely to screw up if only
| one needs to be reset.

I guess I'm missing something because I don't see how having separate
states for rand and randn has anything to do with parallelism.

jwe
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[Michael Creel]

On Thu, Mar 6, 2008 at 6:49 PM, John W. Eaton <jwe@...> wrote:

On  5-Mar-2008, Michael Creel wrote:

| On this issue, I vote for having the seed set in a single spot for all RNGs.
| This facilitates greatly parallel use of Octave. Sometimes it's important to
| be able to reset the seed, and it's a lot less likely to screw up if only
| one needs to be reset.

I guess I'm missing something because I don't see how having separate
states for rand and randn has anything to do with parallelism.

jwe

It doesn't have any intrinsic effect, of course, but from the point of view of ease of programming and likelihood of making hard to detect errors, it is not neutral. If RNG has its own seed, then the programmer must control the value of that seed on each node. If two RNGs that originally share a seed become endowed with separate seeds in a new version of Octave, the programmer must be aware of this. If he or she is not, the code will continue to run, but with bogus results. If we're only talking about having one or two seeds for uniform and normal, it's not that big of a deal. But if RNGs for other distributions also go this route, it will be an inconvenience. Basically, I find that the current situation with one seed for all RNGs is ideal. Just my 2 euro cents.
Michael

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[David Bateman]

The change that was made results in separate seeds for the rand, rand, rande, randg and rande generators.. The other generators in the statistics toolbox rely on these five intrinsic generators and so will of cause now be seeded differently.

However, I'm not sure I see the need for seeding in a parallel situation being any different from the normal case. The programmer needs to be aware of the change in both cases.

BTW, the old generators had/have seperate seeds, and the change makes the behavior of seeding with 'state' or 'seed' similar.


-----Original Message-----
From: michael.creel@... on behalf of Michael Creel
Sent: Fri 07-Mar-08 9:27 AM
To: John W. Eaton
Cc: bug-octave@...
Subject: Re: Difference between Matlab and Octave in random seeding

On Thu, Mar 6, 2008 at 6:49 PM, John W. Eaton <jwe@...> wrote:

> On  5-Mar-2008, Michael Creel wrote:
>
> | On this issue, I vote for having the seed set in a single spot for all
> RNGs.
> | This facilitates greatly parallel use of Octave. Sometimes it's
> important to
> | be able to reset the seed, and it's a lot less likely to screw up if
> only
> | one needs to be reset.
>
> I guess I'm missing something because I don't see how having separate
> states for rand and randn has anything to do with parallelism.
>
> jwe


It doesn't have any intrinsic effect, of course, but from the point of view
of ease of programming and likelihood of making hard to detect errors, it is
not neutral. If RNG has its own seed, then the programmer must control the
value of that seed on each node. If two RNGs that originally share a seed
become endowed with separate seeds in a new version of Octave, the
programmer must be aware of this. If he or she is not, the code will
continue to run, but with bogus results. If we're only talking about having
one or two seeds for uniform and normal, it's not that big of a deal. But if
RNGs for other distributions also go this route, it will be an
inconvenience. Basically, I find that the current situation with one seed
for all RNGs is ideal. Just my 2 euro cents.
Michael

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[Michael Creel]

The need for controlling how seeding is done is not different, nor is the need to be aware of the change. However, the change makes life more difficult and error-prone, I think, which is why I registered a vote against it.
Cheers, Michael

On Fri, Mar 7, 2008 at 3:44 PM, Bateman David-ADB014 <David.Bateman@...> wrote:

The change that was made results in separate seeds for the rand, rand, rande, randg and rande generators.. The other generators in the statistics toolbox rely on these five intrinsic generators and so will of cause now be seeded differently.

However, I'm not sure I see the need for seeding in a parallel situation being any different from the normal case. The programmer needs to be aware of the change in both cases.

BTW, the old generators had/have seperate seeds, and the change makes the behavior of seeding with 'state' or 'seed' similar.


-----Original Message-----
From: michael.creel@... on behalf of Michael Creel
Sent: Fri 07-Mar-08 9:27 AM
To: John W. Eaton
Cc: bug-octave@...
Subject: Re: Difference between Matlab and Octave in random seeding

On Thu, Mar 6, 2008 at 6:49 PM, John W. Eaton <jwe@...> wrote:

> On  5-Mar-2008, Michael Creel wrote:
>
> | On this issue, I vote for having the seed set in a single spot for all
> RNGs.
> | This facilitates greatly parallel use of Octave. Sometimes it's
> important to
> | be able to reset the seed, and it's a lot less likely to screw up if
> only
> | one needs to be reset.
>
> I guess I'm missing something because I don't see how having separate
> states for rand and randn has anything to do with parallelism.
>
> jwe


It doesn't have any intrinsic effect, of course, but from the point of view
of ease of programming and likelihood of making hard to detect errors, it is
not neutral. If RNG has its own seed, then the programmer must control the
value of that seed on each node. If two RNGs that originally share a seed
become endowed with separate seeds in a new version of Octave, the
programmer must be aware of this. If he or she is not, the code will
continue to run, but with bogus results. If we're only talking about having
one or two seeds for uniform and normal, it's not that big of a deal. But if
RNGs for other distributions also go this route, it will be an
inconvenience. Basically, I find that the current situation with one seed
for all RNGs is ideal. Just my 2 euro cents.
Michael

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[pbellec]

My understanding is that any difference of behavior between matlab and octave should be regarded as a bug, so there should be different seeds for rand and randn.

However, I must confess I have a very hard time with the use of different seeds in matlab. I was not aware of this subtility when I started to perform Monte Carlo simulations that mixed rand and randn calls, so I seeded rand only. Because the Monte Carlo chain actually depend on both seeds, it took me weeks before I finally found the bug.

Suggestion of workaround :
make Octave behave as Matlab, with specific seeds for rand and randn (and rande and randg...), but create an octave-specific instruction rand('all_state',S) which would seed all RNGs (rand, randn, rande, randg). Adding this instruction in the help of rand* would force users to be aware of the subtility.

Pierre

Michael Creel wrote:

The need for controlling how seeding is done is not different, nor is the
need to be aware of the change. However, the change makes life more
difficult and error-prone, I think, which is why I registered a vote against
it.
Cheers, Michael

On Fri, Mar 7, 2008 at 3:44 PM, Bateman David-ADB014 <
[e-mail]> wrote:

>
> The change that was made results in separate seeds for the rand, rand,
> rande, randg and rande generators.. The other generators in the statistics
> toolbox rely on these five intrinsic generators and so will of cause now be
> seeded differently.
>
> However, I'm not sure I see the need for seeding in a parallel situation
> being any different from the normal case. The programmer needs to be aware
> of the change in both cases.
>
> BTW, the old generators had/have seperate seeds, and the change makes the
> behavior of seeding with 'state' or 'seed' similar.
>
>
> -----Original Message-----
> From: [e-mail] on behalf of Michael Creel
> Sent: Fri 07-Mar-08 9:27 AM
> To: John W. Eaton
> Cc: [e-mail]
> Subject: Re: Difference between Matlab and Octave in random seeding
>
> On Thu, Mar 6, 2008 at 6:49 PM, John W. Eaton <[e-mail]>
> wrote:
>
> > On  5-Mar-2008, Michael Creel wrote:
> >
> > | On this issue, I vote for having the seed set in a single spot for all
> > RNGs.
> > | This facilitates greatly parallel use of Octave. Sometimes it's
> > important to
> > | be able to reset the seed, and it's a lot less likely to screw up if
> > only
> > | one needs to be reset.
> >
> > I guess I'm missing something because I don't see how having separate
> > states for rand and randn has anything to do with parallelism.
> >
> > jwe
>
>
> It doesn't have any intrinsic effect, of course, but from the point of
> view
> of ease of programming and likelihood of making hard to detect errors, it
> is
> not neutral. If RNG has its own seed, then the programmer must control the
> value of that seed on each node. If two RNGs that originally share a seed
> become endowed with separate seeds in a new version of Octave, the
> programmer must be aware of this. If he or she is not, the code will
> continue to run, but with bogus results. If we're only talking about
> having
> one or two seeds for uniform and normal, it's not that big of a deal. But
> if
> RNGs for other distributions also go this route, it will be an
> inconvenience. Basically, I find that the current situation with one seed
> for all RNGs is ideal. Just my 2 euro cents.
> Michael
>
>

_______________________________________________
Bug-octave mailing list
[e-mail]
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[John W. Eaton]

On  7-Mar-2008, pbellec wrote:

| My understanding is that any difference of behavior between matlab and octave
| should be regarded as a bug,

That's false.

jwe
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[Ben Abbott]

On Mar 7, 2008, at 6:25 PM, John W. Eaton wrote:

> On  7-Mar-2008, pbellec wrote:
>
> | My understanding is that any difference of behavior between matlab  
> and octave
> | should be regarded as a bug,
>
> That's false.
>
> jwe

I have no vested interest in this particular discussion, but do have  
an interest in consistency with Matlab.

Does this "feature" cross that line?

Ben


_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[pbellec]

John W. Eaton wrote:

On  7-Mar-2008, pbellec wrote:

| My understanding is that any difference of behavior between matlab and octave
| should be regarded as a bug,

That's false.

jwe

Well, then it would certainly be better to have rand('state',S) reset the state of all RNGs as was suggested by Michael.

[John W. Eaton]

On  7-Mar-2008, pbellec wrote:

| Well, then it would certainly be better to have rand('state',S) reset the
| state of all RNGs as was suggested by Michael.

I don't see why.  I don't see anything that says the generators must
be based on the same algorithm, so it might not even be possible for
them to share the same state.

jwe
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[Jaroslav Hajek-2]

On Sat, Mar 8, 2008 at 12:04 AM, pbellec <pbellec@...> wrote:

>
>  My understanding is that any difference of behavior between matlab and octave
>  should be regarded as a bug, so there should be different seeds for rand and
>  randn.
>
>  However, I must confess I have a very hard time with the use of different
>  seeds in matlab. I was not aware of this subtility when I started to perform
>  Monte Carlo simulations that mixed rand and randn calls, so I seeded rand
>  only. Because the Monte Carlo chain actually depend on both seeds, it took
>  me weeks before I finally found the bug.
>
>  Suggestion of workaround :
>  make Octave behave as Matlab, with specific seeds for rand and randn (and
>  rande and randg...), but create an octave-specific instruction
>  rand('all_state',S) which would seed all RNGs (rand, randn, rande, randg).
>  Adding this instruction in the help of rand* would force users to be aware
>  of the subtility.

Personally, I think this difference is too subtle to worry about. The
Matlab behaviour is not clearly better, it's just different, and I can
imagine it changing in future Matlab versions. Adding a workaround to
code that is supposed to behave identically under Matlab and Octave is
easy enough.

>
>  Pierre
>
>
>
>
>  Michael Creel wrote:
>  >
>  > The need for controlling how seeding is done is not different, nor is the
>  > need to be aware of the change. However, the change makes life more
>  > difficult and error-prone, I think, which is why I registered a vote
>  > against
>  > it.
>  > Cheers, Michael
>  >
>  > On Fri, Mar 7, 2008 at 3:44 PM, Bateman David-ADB014 <
>  > David.Bateman@...> wrote:
>  >
>  >>
>  >> The change that was made results in separate seeds for the rand, rand,
>  >> rande, randg and rande generators.. The other generators in the
>  >> statistics
>  >> toolbox rely on these five intrinsic generators and so will of cause now
>  >> be
>  >> seeded differently.
>  >>
>  >> However, I'm not sure I see the need for seeding in a parallel situation
>  >> being any different from the normal case. The programmer needs to be
>  >> aware
>  >> of the change in both cases.
>  >>
>  >> BTW, the old generators had/have seperate seeds, and the change makes the
>  >> behavior of seeding with 'state' or 'seed' similar.
>  >>
>  >>
>  >> -----Original Message-----
>  >> From: michael.creel@... on behalf of Michael Creel
>  >> Sent: Fri 07-Mar-08 9:27 AM
>  >> To: John W. Eaton
>  >> Cc: bug-octave@...
>  >> Subject: Re: Difference between Matlab and Octave in random seeding
>  >>
>  >> On Thu, Mar 6, 2008 at 6:49 PM, John W. Eaton <jwe@...>
>  >> wrote:
>  >>
>  >> > On  5-Mar-2008, Michael Creel wrote:
>  >> >
>  >> > | On this issue, I vote for having the seed set in a single spot for
>  >> all
>  >> > RNGs.
>  >> > | This facilitates greatly parallel use of Octave. Sometimes it's
>  >> > important to
>  >> > | be able to reset the seed, and it's a lot less likely to screw up if
>  >> > only
>  >> > | one needs to be reset.
>  >> >
>  >> > I guess I'm missing something because I don't see how having separate
>  >> > states for rand and randn has anything to do with parallelism.
>  >> >
>  >> > jwe
>  >>
>  >>
>  >> It doesn't have any intrinsic effect, of course, but from the point of
>  >> view
>  >> of ease of programming and likelihood of making hard to detect errors, it
>  >> is
>  >> not neutral. If RNG has its own seed, then the programmer must control
>  >> the
>  >> value of that seed on each node. If two RNGs that originally share a seed
>  >> become endowed with separate seeds in a new version of Octave, the
>  >> programmer must be aware of this. If he or she is not, the code will
>  >> continue to run, but with bogus results. If we're only talking about
>  >> having
>  >> one or two seeds for uniform and normal, it's not that big of a deal. But
>  >> if
>  >> RNGs for other distributions also go this route, it will be an
>  >> inconvenience. Basically, I find that the current situation with one seed
>  >> for all RNGs is ideal. Just my 2 euro cents.
>  >> Michael
>  >>
>  >>
>  >
>
> > _______________________________________________
>  > Bug-octave mailing list
>  > Bug-octave@...
>  > https://www.cae.wisc.edu/mailman/listinfo/bug-octave
>  >
>  >
>
>  --
>  View this message in context: http://www.nabble.com/Difference-between-Matlab-and-Octave-in-random-seeding-tp15508096p15908425.html
>
> Sent from the Octave - Bugs mailing list archive at Nabble.com.
>
>  _______________________________________________
>
>
> Bug-octave mailing list
>  Bug-octave@...
>  https://www.cae.wisc.edu/mailman/listinfo/bug-octave
>

--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[pbellec]

John W. Eaton wrote:

On  7-Mar-2008, pbellec wrote:

| Well, then it would certainly be better to have rand('state',S) reset the
| state of all RNGs as was suggested by Michael.

I don't see why.  I don't see anything that says the generators must
be based on the same algorithm, so it might not even be possible for
them to share the same state.

jwe
_______________________________________________
Bug-octave mailing list
[e-mail]
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

I get your point. However, in matlab, help messages for rand* functions are confusing. The help of rand says something as "You can query the state of the random number generator using...", while the help of randn says "The sequence of numbers generated is determined by the state of the generator...". This can give the false idea that there is such a thing as a single underlying "random number generator" whose state would impact both randn and rand. My comment would simply be to make sure that the help of octave is a bit more explicit than the help of matlab regarding the state of generators, e.g.:
"You can query the state of the UNIFORM random number generator using..."
and
"The sequence of numbers generated is determined by the state of the NORMAL generator..."
and maybe add a warning :
"Please note that the states of the uniform, normal, exponential and gamma random number generators are independent."