FW: bad rand generator: distribution/prevalences/periodicity

Discussion specific to NXT-G, NXC, NBC, RobotC, Lejos, and more.
HaWe
Posts: 2500
Joined: 04 Nov 2014, 19:00

Re: FW: extremely poor equal distribution of random prevalences

Post by HaWe »

new test:
how is the distribution of two consecutive numbers?
(all consecutive numbers should be nearly equally distributed, otherwise there must be clusters or local minima or maxima or hyperplane behaviors).
This is the test code:

Code: Select all

// Test on statistical equal distributionof (n, n+1)
// of Lego and customized randomization functions
// ver. 1.1

#define printf1( _x, _y, _format1, _value1) { \
  string sval1 = FormatNum(_format1, _value1); \
  TextOut(_x, _y, sval1); \
}
char LCDline[]={56,48,40,32,24,16,8,0};


//********************************************
// randomize functions
//********************************************
//-------------------------------------------------------
// Kernighan & Ritchie LCG (linear congruence generator)
//-------------------------------------------------------

// global variables and functions to plant a random seed and to randomize
// how to use:
// sKrand(x)                         // x = val. for repet. sequences
// myVar = Krand() % 100;           // assignes a random number


unsigned long _RAND_SEED_ = 1;      // 1= default value (program start)
unsigned long _OLD_SEED_  = 1;


//--------------------------------------------

unsigned long Krand()               // custom random function
{
  _RAND_SEED_ = _RAND_SEED_ * 1103515245 + 12345;
  return (_RAND_SEED_ % (RAND_MAX + 1));
}

//--------------------------------------------


void sKrand(unsigned long seed)      // seeds for a new random series
{

  if (seed==0)                      // 0: a "real" randomized random seed
  {
    seed = ( CurrentTick()*BatteryLevel() );
  }  // substitute to time(0)

  else
  if (seed==-1) {                   // -1: restore last random series
    seed = _OLD_SEED_;
  }

  _OLD_SEED_  = seed;
  _RAND_SEED_ = seed;               // patch for rand_ function
}


//---------------------------------------------------------
// Mersenne Twister by Makoto Matsumoto & Takuji Nishimura
//---------------------------------------------------------
const int N_MTW = 25;
const int M_MTW = 7;
const unsigned long A_MTW[] = {0, 0x8ebfd028};

unsigned long y_MTW[25];
int i_MTW;
char virgin_MTW=TRUE, new_MTW=FALSE;

void InitMTW() {
  i_MTW = N_MTW+1;
}

//--------------------------------------------

unsigned long sMrand(unsigned long x) {
  unsigned long r, s;
  int k;
  if (x==1) InitMTW();

  else {
    InitMTW();
    if (x==0) {x=Random(CurrentTick());}
    for (k=0 ; k<25 ; ++k) {
      InitMTW();
      r = 9+x;
      s = 3402+x;
      for (k=0 ; k<25 ; ++k) {
           r = 509845221 * r + 3;
           s *= s + 1;
           y_MTW[k] = s + (r >> 10);
      }
      new_MTW=TRUE;
    }
  }
}

//--------------------------------------------

unsigned long Mrand(void) {
  unsigned long r, s, e;
  int k;

  if (virgin_MTW) {InitMTW(); virgin_MTW=FALSE; }

  if (i_MTW >= N_MTW) {
      if ((i_MTW > N_MTW) && !new_MTW) {
         r = 9;
         s = 3402;
         for (k=0 ; k<N_MTW ; ++k) {
           r = 509845221 * r + 3;
           s *= s + 1;
           y_MTW[k] = s + (r >> 10);
         }
      }
      for (k=0 ; k<N_MTW-M_MTW ; ++k)
         y_MTW[k] = y_MTW[k+M_MTW] ^ (y_MTW[k] >> 1) ^ A_MTW[y_MTW[k] & 1];
      for (; k<N_MTW ; ++k)
         y_MTW[k] = y_MTW[k+(M_MTW-N_MTW)] ^ (y_MTW[k] >> 1) ^ A_MTW[y_MTW[k] & 1];
      i_MTW = 0;
  }

  new_MTW=false;
  e = y_MTW[i_MTW++];
  e ^= (e << 7)  & 0x2b5b2500;
  e ^= (e << 15) & 0xdb8b0000;
  e ^= (e >> 16);
  return e;
}


//--------------------------------------------


task main(){
  long i, j, old;
  long H[16];            // prevalence counter array
  float x;

  sKrand(2);  sMrand(2);

                         // print output matrix
  for (j=0; j<16; ++j) {
    printf1(0+(j%2)*50, LCDline[j/2], "%2d:", j);
  }

  for (i=-1; i<= (1024 << 7); ++i) {
  
    // 128 k  randomizations of numbers between 0...15
  
    j=rand()%16;        // std Lego FW rand
    //j=Random(16);       // NXC John Hansen rand
    //j=Krand()%16;       // Kernighan and Ritchie LCG rand
    //j=Mrand()%16;       // Mesenne Twister rand

    if ((j==old+1) || ((j==0 && old==15)))
    if (i>=0) H[j]+=1;   // add +1 to the prevalence to the related counter

    old=j;
                         // print prevalence test distribution
    printf1(16+(j%2)*50,LCDline[j/2], "%5d", H[j]);
  }
  
  PlaySound(SOUND_FAST_UP);
  
  while(!ButtonPressed(BTNCENTER, true));

  
  ClearScreen();
                               printf1(0,56, "count=%8d", i+1); //  !!
  j=ArrayMin (H, NA, NA);      printf1(0,48, "min  =%8d", j);
  j=ArrayMax (H, NA, NA);      printf1(0,40, "max  =%8d", j);
  j=ArrayMean(H, NA, NA);      printf1(0,32, "mean =%8d", j);
  x=ArrayStd (H, NA, NA);      printf1(0,24, "stand=%8.2f", x);
  
  while(1) {}

}
This is for lego rand()%16:
worst case:<br />Lego fw rand() % 16
worst case:
Lego fw rand() % 16
rand(%16)x_x+1.jpg (20.66 KiB) Viewed 7781 times
worst case:<br />Lego fw rand() % 16 statistics
worst case:
Lego fw rand() % 16 statistics
rand(%16)x_x+1_stat.jpg (11.3 KiB) Viewed 7781 times
Last edited by HaWe on 25 Jun 2011, 21:18, edited 8 times in total.
HaWe
Posts: 2500
Joined: 04 Nov 2014, 19:00

Re: FW: extremely poor equal distribution of random prevalences

Post by HaWe »

this is with John's NXC Random (16) - much better!
better: <br />John's NXC Random(16)
better:
John's NXC Random(16)
random(16)x_x+1.jpg (21.32 KiB) Viewed 7781 times
better: <br />John's NXC Random(16) statistics
better:
John's NXC Random(16) statistics
random(16)x_x+1_stat.jpg (11.14 KiB) Viewed 7781 times
Last edited by HaWe on 25 Jun 2011, 14:37, edited 1 time in total.
HaWe
Posts: 2500
Joined: 04 Nov 2014, 19:00

Re: FW: extremely poor equal distribution of random prevalences

Post by HaWe »

this is the result with Kernighan + Richie's LCG Krand()%16:
random seed sKrand(1);

also much better (mean better but admittedly stand_dev a little worse than John's code)^^ - but that may depend on modulus ;)
also much better:<br />Kernighan + Richie's LCG Krand()%16
also much better:
Kernighan + Richie's LCG Krand()%16
krand(%16)x_x+1.jpg (14.41 KiB) Viewed 7781 times
also much better:<br />Kernighan + Richie's LCG Krand()%16 statistics
also much better:
Kernighan + Richie's LCG Krand()%16 statistics
krand(%16)x_x+1_stat.jpg (5.92 KiB) Viewed 7781 times
Last edited by HaWe on 25 Jun 2011, 17:40, edited 4 times in total.
HaWe
Posts: 2500
Joined: 04 Nov 2014, 19:00

Re: FW: extremely poor equal distribution of random prevalences

Post by HaWe »

best:
Mersenne Twister!

Mrand%16, random seed sMrand(1);
best:<br />Mersenne Twister Mrand()%16
best:
Mersenne Twister Mrand()%16
Mrand(%16)x_x+1.jpg (10.64 KiB) Viewed 7780 times
best:<br />Mersenne Twister Mrand()%16 statistics
best:
Mersenne Twister Mrand()%16 statistics
Mrand(%16)x_x+1_stat.jpg (10.93 KiB) Viewed 7780 times
HaWe
Posts: 2500
Joined: 04 Nov 2014, 19:00

Re: FW: extremely poor equal distribution of random prevalences

Post by HaWe »

next test:
test on periodicity length ("Infinite Monkey Test"):
the first 9 randomized numbers are stored and then a repetitive randomization is started until the start sequence is reached again.
By a random(100) the probability of an accidentally re-drawn start pattern is 1 : 100^9 (nearly inifinte improbability ^^).
So the longer the minimum periodicity length is until a start sequence is re-generated again - : the better is the random number generator.

This is the test source code (test algorithm: a finite state machine (FSM), inspired by a basic idea of hergipotter - I'm really really curious if it works and what the results are like...):

Code: Select all

// Test on periodicity length ("Infinite Monkey Test"):
// of Lego and customized randomization functions
// ver. 1.6

#define printf1( _x, _y, _format1, _value1) { \
  string sval1 = FormatNum(_format1, _value1); \
  TextOut(_x, _y, sval1); \
}


char LCDline[]={56,48,40,32,24,16,8,0};


//********************************************
// randomize functions
//********************************************
//-------------------------------------------------------
// Kernighan & Ritchie LCG (linear congruence generator)
//-------------------------------------------------------

unsigned long _RAND_SEED_ = 1;      // 1= default value (program start)
unsigned long _OLD_SEED_  = 1;


//--------------------------------------------

unsigned long Krand()               // custom random function
{
  _RAND_SEED_ = _RAND_SEED_ * 1103515245 + 12345;
  return (_RAND_SEED_ % (RAND_MAX + 1));
}

//--------------------------------------------


void sKrand(unsigned long seed)      // seeds for a new random series
{

  if (seed==0)                      // 0: a "real" randomized random seed
  {
    seed = ( CurrentTick()*BatteryLevel() );
  }  // substitute to time(0)

  else
  if (seed==-1) {                   // -1: restore last random series
    seed = _OLD_SEED_;
  }

  _OLD_SEED_  = seed;
  _RAND_SEED_ = seed;               // patch for rand_ function
}


//---------------------------------------------------------
// Mersenne Twister by Makoto Matsumoto & Takuji Nishimura
//---------------------------------------------------------
const int N_MTW = 25;
const int M_MTW = 7;
const unsigned long A_MTW[] = {0, 0x8ebfd028};

unsigned long y_MTW[25];
int i_MTW;
char virgin_MTW=TRUE, new_MTW=FALSE;

void InitMTW() {
  i_MTW = N_MTW+1;
}

//--------------------------------------------

unsigned long sMrand(unsigned long x) {
  unsigned long r, s;
  int k;
  if (x==1) InitMTW();

  else {
    InitMTW();
    if (x==0) {x=Random(CurrentTick());}
    for (k=0 ; k<25 ; ++k) {
      InitMTW();
      r = 9+x;
      s = 3402+x;
      for (k=0 ; k<25 ; ++k) {
           r = 509845221 * r + 3;
           s *= s + 1;
           y_MTW[k] = s + (r >> 10);
      }
      new_MTW=TRUE;
    }
  }
}

//--------------------------------------------

unsigned long Mrand(void) {
  unsigned long r, s, e;
  int k;

  if (virgin_MTW) {InitMTW(); virgin_MTW=FALSE; }

  if (i_MTW >= N_MTW) {
      if ((i_MTW > N_MTW) && !new_MTW) {
         r = 9;
         s = 3402;
         for (k=0 ; k<N_MTW ; ++k) {
           r = 509845221 * r + 3;
           s *= s + 1;
           y_MTW[k] = s + (r >> 10);
         }
      }
      for (k=0 ; k<N_MTW-M_MTW ; ++k)
         y_MTW[k] = y_MTW[k+M_MTW] ^ (y_MTW[k] >> 1) ^ A_MTW[y_MTW[k] & 1];
      for (; k<N_MTW ; ++k)
         y_MTW[k] = y_MTW[k+(M_MTW-N_MTW)] ^ (y_MTW[k] >> 1) ^ A_MTW[y_MTW[k] & 1];
      i_MTW = 0;
  }

  new_MTW=false;
  e = y_MTW[i_MTW++];
  e ^= (e << 7)  & 0x2b5b2500;
  e ^= (e << 15) & 0xdb8b0000;
  e ^= (e >> 16);
  return e;
}


//--------------------------------------------

unsigned long i;       // counter
unsigned long H[10];   // for start sequence
int maxZ;              // maximum level for subsequences

void showWhat(){
  float V;
  int r, c;
  
  TextOut( 0,56, "start rand sequence:");
  for (c=0; c<3; ++c) {
     r=H[c]; printf1(c*30,48,"%3u",r);
  }
  for (c=0; c<3; ++c) {
     r=H[c+3]; printf1(c*30,40,"%3u",r);
  }
  for (c=0; c<3; ++c) {
     r=H[c+6]; printf1(c*30,32,"%3u",r);
  }
  printf1( 0,24,"maxLevel=%2u", maxZ);   // max. reached level
  if (maxZ>0) {                          // number of max. reached level
    printf1(60,16,"sub%3u", H[maxZ-1]);
  }

  printf1(0, 8,"count=%d", i);
  V=BatteryLevel()/1000.0;
  printf1(0, 0,"Volt=%6.3f",V);
}


task main(){
  unsigned long j;

  float x;
  int level;

  sKrand(1);  sMrand(1);

  
  for (i=0; i<10; ++i) {  // start store 9 generated numbers (1st dropped)
  
    j=rand()%100;        // std Lego FW rand
    //j=Random(100);       // NXC John Hansen rand
    //j=Krand()%100;       // Kernighan and Ritchie LCG rand
    //j=Mrand()%100;       // Mesenne Twister rand
    if (i>0)
      H[i-1]=j;            // end store the first 9 generated numbers
  }
  
  showWhat();
  Wait(1000);
  level=0;
  
  while(1) {
  
      j=rand()%100;        // std Lego FW rand
      //j=Random(100);       // NXC John Hansen rand
      //j=Krand()%100;       // Kernighan and Ritchie LCG rand
      //j=Mrand()%100;       // Mesenne Twister rand


      if (j==H[level])  {
        level+=1;
        if (level>maxZ)   {  maxZ=level;  }
      }
      else
      {
        if (level>0)    {
          level=0;
        }
      }  // if
      if (level>maxZ)   {
        maxZ=level;
      }
      i+=1;

      while (ButtonPressed(BTNCENTER, false)) {
        showWhat();
        printf1( 0,16,"rand=%3u", j);          // current rand number

      }
      if (ButtonPressed(BTNRIGHT,  false))  {
        showWhat();
        printf1( 0,16,"rand=%3u", j);          // current rand number

        Wait(200);
      }

      if ((i%2000==0)|| (i<1000)) {
        showWhat();
        printf1( 0,16,"rand=%3u", j);          // current rand number

      }

      if (level>1) {
        if (level>3) PlayTone(200*level, 100);
        showWhat();
         printf1( 0,16,"rand=%3u", j);          // current rand number

        Wait(200); }

      if (maxZ>8)
      {
        PlaySound(SOUND_FAST_UP);
        ClearScreen();

        showWhat();
        printf1(0,16, "strtseq.reach%3u",j); // current number of final level
        while(1);
      }
  }
  
}
run (1): periodicity=113484 iterations
Lego fw rand()%100 <br />run (1): periodicity=113484 iterations
Lego fw rand()%100
run (1): periodicity=113484 iterations
Monkey9rand%16.jpg (15.8 KiB) Viewed 7754 times
run (2): periodicity=2134868 iterations
Lego fw rand()%100 <br />run (2): periodicity=2134868 iterations
Lego fw rand()%100
run (2): periodicity=2134868 iterations
Monkey9rand%16(2).jpg (18.23 KiB) Viewed 7746 times
run (3): periodicity=4180450
Lego fw rand()%100 <br />run (3): periodicity=4180450
Lego fw rand()%100
run (3): periodicity=4180450
Monkey9rand%16(3).jpg (18.89 KiB) Viewed 7740 times
Last edited by HaWe on 28 Jun 2011, 07:44, edited 2 times in total.
HaWe
Posts: 2500
Joined: 04 Nov 2014, 19:00

Re: FW: extremely poor equal distribution of random prevalences

Post by HaWe »

again Lego fw rand()%100.

that's really odd:
in run (4): repetitive sequence already after 50808 iterations!!
Lego fw rand()%100 <br />run (4): periodicity=  50808
Lego fw rand()%100
run (4): periodicity= 50808
Monkey9rand%100(4).jpg (17.88 KiB) Viewed 7738 times
run (5): periodicity = 1542188 iterations
Lego fw rand()%100 <br />run (5): periodicity= 1542188 iterations
Lego fw rand()%100
run (5): periodicity= 1542188 iterations
Monkey9rand%100(5).jpg (18.52 KiB) Viewed 7737 times
run (6): periodicity = 996500 iterations
Lego fw rand()%100 <br />run (6): periodicity= 996500  iterations
Lego fw rand()%100
run (6): periodicity= 996500 iterations
Monkey9rand%100(6).jpg (18.59 KiB) Viewed 7736 times
HaWe
Posts: 2500
Joined: 04 Nov 2014, 19:00

FW: extremely poor equal distribution of random prevalences

Post by HaWe »

now
NXC Random(100):
run (1): 111176 iterations
NXC Random(100)
NXC Random(100)
Monkey9random(100)(1).jpg (18.27 KiB) Viewed 7736 times
run (2): 3022238 iterations
NXC Random(100)
NXC Random(100)
Monkey9random(100)(2).jpg (18.73 KiB) Viewed 7729 times
run (3): 112760 iterations
NXC Random(100)
NXC Random(100)
Monkey9random(100)(3).jpg (18.26 KiB) Viewed 7729 times
HaWe
Posts: 2500
Joined: 04 Nov 2014, 19:00

FW: bad rand generator: distribution/prevalences/periodicity

Post by HaWe »

interim report:

3 NXT bricks using Krand() modulo 100 (Krand: the Kernighan & Richie LCG ) have been running simultaneously since Sunday evening 7 pm (GMST+1) , still running.
They had been initialized with different seeds (sKrand(1), sKrand(2), sKrand(3) ).

Meanwhile there have been > 300 million random iterations each, maximum subsequence levels are (as expected) 4 at all 3 NXTs , no complete level-9 start sequences have been found / re-generated so far:

(prematurely user interrupted because of other purposes !!)
Kernighan & Richie LCG Krand() % 100 seeded by sKrand(1)
Monkey9Krand%100(1i).jpg
Monkey9Krand%100(1i).jpg (18.13 KiB) Viewed 7620 times
Kernighan & Richie LCG Krand() % 100 seeded by sKrand(2)
Monkey9Krand%100(2i).jpg
Monkey9Krand%100(2i).jpg (18.63 KiB) Viewed 7620 times
Kernighan & Richie LCG Krand() % 100 seeded by sKrand(3)
Monkey9Krand%100(3i).jpg
Monkey9Krand%100(3i).jpg (18.87 KiB) Viewed 7647 times
to remember: As shown above, the Lego function rand() has regenerated that level-9 start sequence already again after (mostly) > 100,000 up to < 5 millions of iterations ("periodicity length") - by using rand()%100 as well as by using Random(100)!
(statistical probability = 1: 100 9 = 1: 1018 = 1: 1 quintillion ).
Again: does anybody know what algorithm exactly the Lego FW rand() uses? My suspect: it uses sth like RANDU :shock:

ps
I had to prematurely stop 1 of the 3 running NXTs because I need it for other purposes...
Last edited by HaWe on 02 Jul 2011, 11:13, edited 2 times in total.
HaWe
Posts: 2500
Joined: 04 Nov 2014, 19:00

Re: FW: bad rand generator: distribution/prevalences/periodicity

Post by HaWe »

This is the intermediate summary of the previous test results (tests yet running/pending)

Code: Select all



               equality of      equality of      random series  seed by
               prevalences of   prevalences of   periodicity    constant or
               single numbers   (n, n+1) tuples  length         random seeds
=============================================================================
rand() % n       - - - - -        - - - - -       - - - - -       - - - 

Random(n)           + +              + +          - - - - -       - - - 

K&R LCD
Krand() % n         + +              + +             + +          + + + 

Mersenne                                            + + + 
Twister            + + +            + + +           10^240         + + 
Mrand() % n                                      (not tested) 
The results of the Mersenne twister periodicity length (in theory:> 10240) have not been finished.

I want to stress the results of the periodicity length - what actually is the meaning of a very short periodicity length of e.g. 1 million?
It means, that every 1 million iterations all randomized patterns repeat themselves. (And for this point of view it actually does not matter if it's now 5 millions or just 50 thousand.)

That means, that (if randomized out of a set of 100) almost
1018 different possible combinations of 9 or
1016 different possible combinations of 8 or
1014 different possible combinations of 7 or
1012 different possible combinations of 6 or
1010 different possible combinations of 5 consecutive numbers
WILL NEVER BE GENERATED.

As we know, even every LCG processes sorts of hyperplane behaviours.
But for not-critical, fast number generation he is fast and good enough and it could be the second best choice.
But for a high quality needed for stochastical filters such as Monte Carlo filters the slower but very high-quality Mersenne Twister should be preferred. Nevertheless, some results of the Mersenne Twister are still pending.
Last edited by HaWe on 02 Jul 2011, 11:05, edited 1 time in total.
afanofosc
Site Admin
Posts: 1256
Joined: 26 Sep 2010, 19:36
Location: Nashville, TN
Contact:

Re: FW: bad rand generator: distribution/prevalences/periodicity

Post by afanofosc »

What sort of real-world LEGO NXT device are you going to run with 1 million iterations of random numbers? A rover? A humanoid? A robotic arm? A maze solver? A balancing bot?

And why do you persist in using a widely known to be bad algorithm for generating random numbers within a small numeric range?

John Hansen
Multi-platform LEGO MINDSTORMS programming
http://bricxcc.sourceforge.net/
Post Reply

Who is online

Users browsing this forum: No registered users and 5 guests