Re: generating random numbers -- what seeds to use?

[Paul Kienzle]

David Doolin wrote:
> On Sun, 27 Apr 2003, Paul Kienzle wrote:
> > winrand has a fast Poisson generator which I have
> > implemented as an oct-file.  The author says the
> > code is public domain so no reason not to redistribute
> > it other than it is tied up in the same file as the ziggurat
> > code.  The winrand code has other generators as well,
> > but I haven't needed them.  Let me know if you are
> > interested and I will send you what I've got.
> >   
>
> I am interested, send on along!  (I am 
> really picky about license issues, see 
> below.)
Here is a zip-file of  Ernst Stadlober's random number
generator package:

http://random.mat.sbg.ac.at/ftp/pub/software/winrand/

I'm using two sets of poisson generators: one
from numerical recipes for numbers with independent
lambda, and the winrand ones for numbers with
shared lambda.   I've reimplemented the winrand
code for small lambda, and I'm using pprsc.c directly
for large lambda.  Looking again at winrand, I see that 
there is more than one set of poisson generators, but
I didn't try them all.

My code assumes an RNOR macro which returns
the next normally distributed number.  I'm not including
rnorexp.c which defines this since I'm unclear on the
license, but it is available directly from JStatSoft
(not TOMS as I said earlier):

http://www.jstatsoft.org/v05/i08/rnorrexp.c

The article is there too if someone wants to make a
free implementation:

http://www.jstatsoft.org/v05/i08/ziggurat.pdf

My code was originally based on John Eaton's 
rand.cc.  My changes (pretty much everything
but the copyright by now, I think ;-) are public
domain.

Paul Kienzle
address@hidden


/*

Copyright (C) 1996, 1997 John W. Eaton

This file is part of Octave.

Octave is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.

Octave is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with Octave; see the file COPYING.  If not, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.


*/


// Octave modules for the Mersenne Twister (MT199337) Random Number Generator
// using Richard J. Wagner's C++ implementaion MersenneTwister.h
//
// This file provides two Octave functions:
//    rand              for Uniform random number
//    randn             for Normal random number
//
// Based on John Eaton's rand.cc and Dirk Eddelbuettel's randmt.cc
// Copyright (C) 1996, 1997 John W. Eaton
// Copyright (C) 1998, 1999 Dirk Eddelbuettel <address@hidden>
//
// 2001-02-10 Paul Kienzle
//   * use Richard J. Wagner's MersenneTwister.h
//   * use John Eaton's rand.cc for similar interface
//   * add rand("state") to obtain complete state
//   * renamed snorm to randn, and removed the `static' keyword from 
//     the variables in randu after thoroughly checking that they were
//     reset each time the function is entered.
// 2003-04-28 Paul Kienzle <address@hidden>
//   * includes poisson generator
// $Id: rand.cc,v 1.1.1.1 2001/10/10 19:54:49 pkienzle Exp $

#include <octave/oct.h>
#include <octave/lo-mappers.h>
#include <octave/lo-specfun.h>
#include <octave/lo-ieee.h>
#include "MersenneTwister.h"

// XXX FIXME XXX rnorrexp is not using proper seed initializer
// XXX FIXME XXX rnorrexp is not using Mersenne twister
#include "rnorrexp.c"

static MTRand randu;

// Octave interface starts here

static octave_value 
do_seed (octave_value_list args)
{
  octave_value retval;

  // Check if they said the magic words
  std::string s_arg = args(0).string_value ();
  if (s_arg == "seed")
    {
      // If they ask for the current "seed", then reseed with the next
      // available random number
      MTRand::uint32 a = randu.randInt();
      randu.seed(a);
      retval = (double)a;
    }
  else if (s_arg == "state")
    {
      MTRand::uint32 state[randu.SAVE];
      randu.save(state);
      RowVector a(randu.SAVE);
      for (int i=0; i < randu.SAVE; i++)
        a(i) = state[i];
      retval = a;
    }
  else
    {
      error ("rand: unrecognized string argument");
      return retval;
    }

  // Check if just getting state
  if (args.length() == 1)
    return retval;

  // Set the state from either a scalar or a previously returned state vector
  octave_value tmp = args(1);
  if (tmp.is_scalar_type ())
    {
      MTRand::uint32 n = MTRand::uint32(tmp.double_value());
      if (! error_state)
        randu.seed(n);
    }
  else if (tmp.is_matrix_type () 
           && tmp.rows() == randu.SAVE && tmp.columns() == 1)
    {
      Array<double> a(tmp.vector_value ());
      if (! error_state)
        {
          MTRand::uint32 state[randu.SAVE];
          for (int i = 0; i < randu.SAVE; i++)
            state[i] = MTRand::uint32(a(i));
          randu.load(state);
        }
    }
  else
    error ("rand: not a state vector");
  
  return retval;
}

DEFUN_DLD (rand, args, nargout, 
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {} rand (@var{x})\n\
@deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\
@deftypefnx {Loadable Function} address@hidden =} rand (\"state\", @var{x})\n\
@deftypefnx {Loadable Function} address@hidden =} rand (\"seed\", @var{x})\n\
Return a matrix with random elements uniformly distributed on the\n\
semi-open interval [0, 1).  The arguments are handled the same as the\n\
arguments for @code{eye}.\n\
\n\
You can reset the state of the random number generator using the\n\
form\n\
\n\
@example\n\
v = rand (\"state\", x)\n\
@end example\n\
\n\
where @var{x} is a scalar value.  This returns the current state\n\
of the random number generator in the column vector @var{v}.  If\n\
@var{x} is not given, then the state is returned but not changed.\n\
Later, you can restore the random number generator to the state @var{v}\n\
using the form\n\
\n\
@example\n\
u = rand (\"state\", v)\n\
@end example\n\
\n\
@noindent\n\
If instead of \"state\" you use \"seed\" to query the random\n\
number generator, then the state will be collapsed from roughly\n\
20000 bits down to 32 bits.\n\
\n\
@code{rand} uses the Mersenne Twister, with a period of 2^19937-1.\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\
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\
\n\
http://www.math.keio.ac.jp/~matumoto/emt.html\n\
@end deftypefn\n\
@seealso{randn, rande, randp}\n")
{
  octave_value_list retval;     // list of return values

  int nargin = args.length ();  // number of arguments supplied
  if (nargin > 2) 
    print_usage("rand");

  else if (nargin > 0 && args(0).is_string())
    retval(0) = do_seed (args);

  else
    {
      int nr=0, nc=0;
      switch (nargin) {
      case 0: return RUNI;
      case 1: get_dimensions(args(0), "rand", nr, nc); break;
      case 2: get_dimensions(args(0), args(1), "rand", nr, nc); break;
      }

      if (! error_state)
        {
          Matrix X(nr, nc);
          double *pX = X.fortran_vec();

          for (int i=nr*nc-1; i >= 0; i--)
            pX[i] = RUNI;

          retval(0) = X;
        }
    }

  return retval;
}

DEFUN_DLD (rand, args, nargout, 
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {} rand (@var{x})\n\
@deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\
@deftypefnx {Loadable Function} address@hidden =} rand (\"state\", @var{x})\n\
@deftypefnx {Loadable Function} address@hidden =} rand (\"seed\", @var{x})\n\
Return a matrix with random elements uniformly distributed on the\n\
semi-open interval [0, 1).  The arguments are handled the same as the\n\
arguments for @code{eye}.\n\
\n\
You can reset the state of the random number generator using the\n\
form\n\
\n\
@example\n\
v = rand (\"state\", x)\n\
@end example\n\
\n\
where @var{x} is a scalar value.  This returns the current state\n\
of the random number generator in the column vector @var{v}.  If\n\
@var{x} is not given, then the state is returned but not changed.\n\
Later, you can restore the random number generator to the state @var{v}\n\
using the form\n\
\n\
@example\n\
u = rand (\"state\", v)\n\
@end example\n\
\n\
@noindent\n\
If instead of \"state\" you use \"seed\" to query the random\n\
number generator, then the state will be collapsed from roughly\n\
20000 bits down to 32 bits.\n\
\n\
@code{rand} uses the Mersenne Twister, with a period of 2^19937-1.\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\
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\
\n\
http://www.math.keio.ac.jp/~matumoto/emt.html\n\
@end deftypefn\n\
@seealso{randn, rande, randp}\n")
{
  octave_value_list retval;     // list of return values

  int nargin = args.length ();  // number of arguments supplied
  if (nargin > 2) 
    print_usage("rand");

  else if (nargin > 0 && args(0).is_string())
    retval(0) = do_seed (args);

  else
    {
      int nr=0, nc=0;
      switch (nargin) {
      case 0: return RUNI;
      case 1: get_dimensions(args(0), "rand", nr, nc); break;
      case 2: get_dimensions(args(0), args(1), "rand", nr, nc); break;
      }

      if (! error_state)
        {
          Matrix X(nr, nc);
          double *pX = X.fortran_vec();

          for (int i=nr*nc-1; i >= 0; i--)
            pX[i] = RUNI;

          retval(0) = X;
        }
    }

  return retval;
}

DEFUN_DLD (randn, args, nargout, 
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {} randn (@var{x})\n\
@deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\
@deftypefnx {Loadable Function} address@hidden =} randn (\"state\", @var{x})\n\
@deftypefnx {Loadable Function} address@hidden =} 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\
@end deftypefn\n\
@seealso{rand, rande, randp}\n")
{
  octave_value_list retval;     // list of return values
  static int first=1;
  if (first) zigset(123456);
  first = 0;

  int nargin = args.length ();  // number of arguments supplied
  if (nargin > 2) 
    print_usage("randn");

  else if (nargin > 0 && args(0).is_string())
    retval(0) = do_seed (args);

  else
    {
      int nr=0, nc=0;
      switch (nargin) {
      case 0: return RNOR;
      case 1: get_dimensions(args(0), "randn", nr, nc); break;
      case 2: get_dimensions(args(0), args(1), "randn", nr, nc); break;
      }

      if (! error_state)
        {
          Matrix X(nr, nc);
          double *pX = X.fortran_vec();

          for (int i=nr*nc-1; i >= 0; i--)
            pX[i] = RNOR;

          retval(0) = X;
        }
    }


  return retval;
}

DEFUN_DLD (rande, args, nargout, 
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {} rande (@var{L}, @var{x})\n\
@deftypefnx {Loadable Function} {} rande (@var{L}, @var{n}, @var{m})\n\
Return a matrix with exponentially distributed random elements.  The\n\
arguments are handled the same as the arguments for @code{rand}.\n\
@end deftypefn\n\
@seealso{rand, randn, randp}\n")
{
  octave_value_list retval;     // list of return values

  int nargin = args.length ();  // number of arguments supplied
  if (nargin > 3)
    {
      print_usage("rande");
      return retval;
    }

  if (nargin == 0) return REXP;

  Matrix lambda(args(0).matrix_value());
  if (error_state) return retval;

  int nr=0, nc=0;
  switch (nargin) {
  case 1: nr = lambda.rows(); nc = lambda.columns(); break;
  case 2: get_dimensions(args(1), "rande", nr, nc); break;
  case 3: get_dimensions(args(1), args(2), "rande", nr, nc); break;
  }

  if ( (nr != lambda.rows() && lambda.rows() != 1)
       || (nc != lambda.columns() && lambda.columns() != 1) )
    {
      error("rande: dimensions of lambda must match requested matrix size");
      return retval;
    }

  Matrix X(nr, nc);
  double *pX = X.fortran_vec();

  if (lambda.length() == 1)
    {
      const double L = lambda(0,0);
      if (L < 0.0 || xisinf(L)) 
        {
          for (int i=nr*nc-1; i >= 0; i--)
            pX[i] = octave_NaN;
        }
      else
        {
          for (int i=nr*nc-1; i >= 0; i--)
            pX[i] = REXP/L;
        }
    }
  else
    {
      const double *pL = lambda.data();
      for (int i=nr*nc-1; i >= 0; i--)
        {
          const double L = pL[i];
          if (L < 0.0 || xisinf(L)) 
            pX[i] = octave_NaN;
          else
            pX[i] = REXP/L;
        }
    }
  
  retval(0) = X;


  return retval;
}

DEFUN_DLD (rande2, args, nargout, 
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {} rande (@var{L}, @var{x})\n\
@deftypefnx {Loadable Function} {} rande (@var{L}, @var{n}, @var{m})\n\
Return a matrix with exponentially distributed random elements.  The\n\
arguments are handled the same as the arguments for @code{rand}.\n\
@end deftypefn\n\
@seealso{rand, randn, randp}\n")
{
  octave_value_list retval;     // list of return values

  int nargin = args.length ();  // number of arguments supplied
  if (nargin > 3)
    {
      print_usage("rande");
      return retval;
    }

  if (nargin == 0) return REXP;

  Matrix lambda(args(0).matrix_value());
  if (error_state) return retval;

  int nr=0, nc=0;
  switch (nargin) {
  case 1: nr = lambda.rows(); nc = lambda.columns(); break;
  case 2: get_dimensions(args(1), "rande", nr, nc); break;
  case 3: get_dimensions(args(1), args(2), "rande", nr, nc); break;
  }

  if ( (nr != lambda.rows() && lambda.rows() != 1)
       || (nc != lambda.columns() && lambda.columns() != 1) )
    {
      error("rande: dimensions of lambda must match requested matrix size");
      return retval;
    }

  Matrix X(nr, nc);
  double *pX = X.fortran_vec();

  if (lambda.length() == 1)
    {
      const double L = lambda(0,0);
      if (L < 0.0 || xisinf(L)) 
        {
          for (int i=nr*nc-1; i >= 0; i--)
            pX[i] = octave_NaN;
        }
      else
        {
          for (int i=nr*nc-1; i >= 0; i--)
            pX[i] = -log(RUNI)/L;
        }
    }
  else
    {
      const double *pL = lambda.data();
      for (int i=nr*nc-1; i >= 0; i--)
        {
          const double L = pL[i];
          if (L < 0.0 || xisinf(L)) 
            pX[i] = octave_NaN;
          else
            pX[i] = -log(RUNI)/L;
        }
    }
  
  retval(0) = X;


  return retval;
}

// **** POISSON ****

// XXX FIXME XXX check what happens with pprsc if U is [0,1] rather than (0,1)

#include "pprsc.c"

// Given uniform u, find x such that CDF(L,x)==u.  Return x.
static void 
poisson_cdf_lookup(double lambda, double *pX, int n)
{
  // Table size is predicated on the maximum value of lambda
  // we want to store in the table, and the maximum value of
  // returned by the uniform random number generator on [0,1).
  // With lambda==10 and u_max = 1 - 1/(2^32+1), we
  // have poisson_pdf(lambda,36) < 1-u_max.  If instead our
  // generator uses more bits of mantissa or returns a value
  // in the range [0,1], then for lambda==10 we need a table 
  // size of 46 instead.  For long doubles, the table size 
  // will need to be longer still.
#define TABLESIZE 46
  double t[TABLESIZE];
  
  // Precompute the table for the u up to and including 0.458.
  // We will almost certainly need it.
  int intlambda = (int)floor(lambda);
  double p;
  int tableidx;
  t[0] = p = exp(-lambda);
  for (tableidx = 1; tableidx <= intlambda; tableidx++) {
    p = p*lambda/(double)tableidx;
    t[tableidx] = t[tableidx-1] + p;
  }
  
  for (int i=n-1; i >= 0; i--) {
    double u = RUNI;
    
    // If u > 0.458 we know we can jump to floor(lambda) before
    // comparing (this observation is based on Stadlober's winrand
    // code). For lambda >= 1, this will be a win.  Lambda < 1
    // is already fast, so adding an extra comparison is not a
    // problem.
    int k = (u > 0.458 ? intlambda : 0);
    
    // We aren't using a for loop here because when we find the
    // right k we want to jump to the next iteration of the
    // outer loop, and the continue statement will only work for 
    // the inner loop.
  nextk:
    if ( u <= t[k] ) {
      pX[i] = (double) k;
      continue;
    }
    if (++k < tableidx) goto nextk;
    
    // We only need high values of the table very rarely so we 
    // don't automatically compute the entire table.
    while (tableidx < TABLESIZE) {
      p = p*lambda/(double)tableidx;
      t[tableidx] = t[tableidx-1] + p;
      // Make sure we converge to 1.0 just in case u is uniform
      // on [0,1] rather than [0,1).
      if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0;
      tableidx++;
      if (u <= t[tableidx-1]) break;
    }
    
    // We are assuming that the table size is big enough here.
    // This should be true even if RUNI is returning values in
    // the range [0,1] rather than [0,1).
    pX[i] = (double)(tableidx-1);
  }
}

// From Press, et al., Numerical Recipes
static void
poisson_rejection (double lambda, double *pX, int n)
{
  double sq = sqrt(2.0*lambda);
  double alxm = log(lambda);
  double g = lambda*alxm - xlgamma(lambda+1.0);
  for (int i=0; i < n; i++) 
    {
      double y, em, t;
      do {
        do {
          y = tan(M_PI*RUNI);
          em = sq * y + lambda;
        } while (em < 0.0);
        em = floor(em);
        t = 0.9*(1.0+y*y)*exp(em*alxm-flogfak(em)-g);
      } while (RUNI > t);
      pX[i] = em;
    }
}

DEFUN_DLD (logfak, args,,"return the log factorial of integer k")
{
  Matrix M(args(0).matrix_value());
  double *pM = M.fortran_vec();
  for (int i=M.rows()*M.columns()-1; i >= 0; i--)
    pM[i] = flogfak(pM[i]);
  return M;
}

DEFUN_DLD (randp, args, nargout, 
  "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {} randp (@var{l})\n\
@deftypefnx {Loadable Function} {} randp (@var{l}, address@hidden, @var{m}])\n\
@deftypefnx {Loadable Function} {} randp (@var{l}, @var{n}, @var{m})\n\
Return a matrix with Poisson distributed random elements.  The\n\
arguments are handled the same as the arguments for @code{rand}.\n\
\n\
@end deftypefn\n\
@seealso{rand, randn, rande}\n")
{
  octave_value_list retval;     // list of return values

  int nargin = args.length ();  // number of arguments supplied
  if (nargin > 3 || nargin < 1) 
    {
      print_usage("randp");
      return retval;
    }

  Matrix lambda(args(0).matrix_value());
  if (error_state) return retval;

  int nr=0, nc=0;
  switch (nargin) {
  case 1: nr = lambda.rows(); nc = lambda.columns(); break;
  case 2: get_dimensions(args(1), "randp", nr, nc); break;
  case 3: get_dimensions(args(1), args(2), "randp", nr, nc); break;
  }

  if (error_state) return retval;

  if ( (nr != lambda.rows() && lambda.rows() != 1)
       || (nc != lambda.columns() && lambda.columns() != 1) )
    {
      error("randp: dimensions of lambda must match requested matrix size");
      return retval;
    }

  Matrix X(nr, nc);
  double *pX = X.fortran_vec();

  // The cutoff of L <= 1e8 before using the normal approximation is based on:
  //   > L=1e8; x=floor(linspace(0,2*L,1000)); 
  //   > max(abs(normal_pdf(x,L,L)-poisson_pdf(x,L)))
  //   ans =  1.1376e-28
  // For L=1e7, the max is around 1e-9, which is within the step size of RUNI.
  // For L>1e10 the pprsc function breaks down, as I saw from the histogram
  // of a large sample, so 1e8 is both small enough and large enough.

  if (nr*nc > 1 && lambda.length()==1)
    {
      const double L = lambda(0,0);
      if (L < 0.0)
        {
          for (int i=nr*nc-1; i>=0; i--) pX[i] = octave_NaN;
        }
      else if (L <= 10.0)
        {
          poisson_cdf_lookup(L, pX, nr*nc);
        }
      else if (L <= 1e8)
        {
          for (int i=nr*nc-1; i>=0; i--) pX[i] = pprsc(L);
        }
      else 
        {
          /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
          const double sqrtL = sqrt(L);
          for (int i=nr*nc-1; i>=0; i--) {
            pX[i] = floor(RNOR*sqrtL + L + 0.5);
            if (pX[i] < 0.0) pX[i] = 0.0; /* will probably never happen */
          }
        }
    }
  else
    {
      const double *pL = lambda.data();
      for (int i=nr*nc-1; i >= 0; i--) 
        {
          double const L = pL[i];
          if (L < 0.0)
            {
              pX[i] = octave_NaN;
            }
          else if (L <= 12.0)
            {
              // From Press, et al. Numerical recipes
              double g = exp(-L);
              int em = -1;
              double t = 1.0;
              do {
                ++em;
                t *= RUNI;
              } while (t > g);
              pX[i] = em;
            }
          else if (L <= 1e8)
            {
              /* numerical recipes */
              poisson_rejection(L, pX+i, 1);
            }
          else
            {
              /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
              pX[i] = floor(RNOR*sqrt(L) + L + 0.5);
              if (pX[i] < 0.0) pX[i] = 0.0; /* will probably never happen */
            }
        }
    }

  retval(0) = X;
  return retval;
}

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/