1/* 2 * randrun - perform a run test on rand() 3 * 4 * Copyright (C) 1999 David I. Bell 5 * 6 * Calc is open software; you can redistribute it and/or modify it under 7 * the terms of the version 2.1 of the GNU Lesser General Public License 8 * as published by the Free Software Foundation. 9 * 10 * Calc is distributed in the hope that it will be useful, but WITHOUT 11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 12 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General 13 * Public License for more details. 14 * 15 * A copy of version 2.1 of the GNU Lesser General Public License is 16 * distributed with calc under the filename COPYING-LGPL. You should have 17 * received a copy with calc; if not, write to Free Software Foundation, Inc. 18 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 19 * 20 * Under source code control: 1995/02/12 20:00:06 21 * File existed as early as: 1995 22 * 23 * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ 24 */ 25 26/* 27 * If X(j) < X(j+1) < ... X(j+k) >= X(j+k+1), then we have a run of 'k'. 28 * We ignore the run breaker, X(j+k+1), and start with X(j+k+2) when 29 * considering a new run in order to make our runs chi independent. 30 * 31 * See Knuth's "Art of Computer Programming - 2nd edition", 32 * Volume 2 ("Seminumerical Algorithms"), Section 3.3.2. 33 * "G. Run test", pp. 65-68, 34 * "problem #14", pp. 74, 536. 35 * 36 * We use the suggestion in problem #14 to allow an application of the 37 * chi-square test and to make estimating the run length probs easy. 38 */ 39 40 41define randrun(run_cnt) 42{ 43 local i; /* index */ 44 local max_run; /* longest run */ 45 local long_run_cnt; /* number of runs longer than MAX_RUN */ 46 local run; /* current run length */ 47 local tally_sum; /* sum of all tally values */ 48 local last; /* last random number */ 49 local current; /* current random number */ 50 local MAX_RUN = 9; /* max run we will keep track of */ 51 local mat tally[1:MAX_RUN]; /* tally of length of a rise run of 'x' */ 52 local mat prob[1:MAX_RUN]; /* prob[x] = probability of 'x' length run */ 53 54 /* 55 * parse args 56 */ 57 if (param(0) == 0) { 58 run_cnt = 65536; 59 } 60 61 /* 62 * run setup 63 */ 64 max_run = 0; /* no runs yet */ 65 long_run_cnt = 0; /* no long runs set */ 66 current = rand(); /* our first number */ 67 run = 1; 68 69 /* 70 * compute the run length probabilities 71 * 72 * A run length of 'r' occurs with a probability of: 73 * 74 * 1/r! - 1/(r+1)! 75 */ 76 for (i=1; i <= MAX_RUN; ++i) { 77 prob[i] = 1.0/fact(i) - 1.0/fact(i+1); 78 } 79 80 /* 81 * look at a number of random number trials 82 */ 83 for (i=0; i < run_cnt; ++i) { 84 85 /* get our current number */ 86 last = current; 87 current = rand(); 88 89 /* look for a run break */ 90 if (current < last) { 91 92 /* record the stats */ 93 if (run > max_run) { 94 max_run = run; 95 } 96 if (run > MAX_RUN) { 97 ++long_run_cnt; 98 } else { 99 ++tally[run]; 100 } 101 102 /* start a new run */ 103 current = rand(); 104 run = 1; 105 106 /* note the continuing run */ 107 } else { 108 ++run; 109 } 110 } 111 /* determine the number of runs found */ 112 tally_sum = matsum(tally) + long_run_cnt; 113 114 /* 115 * print the stats 116 */ 117 printf("rand run test used %d values to produce %d runs\n", 118 run_cnt, tally_sum); 119 for (i=1; i <= MAX_RUN; ++i) { 120 printf("length=%d\tprob=%9.7f\texpect=%d \tcount=%d \terr=%9.7f\n", 121 i, prob[i], round(tally_sum*prob[i]), tally[i], 122 (tally[i] - round(tally_sum*prob[i]))/tally_sum); 123 } 124 printf("length>%d\t\t\t\t\tcount=%d\n", MAX_RUN, long_run_cnt); 125 printf("max length=%d\n", max_run); 126} 127 128if (config("resource_debug") & 3) { 129 print "randrun([run_length]) defined"; 130} 131