source: project/release/4/random-test/trunk/random-test.scm @ 14480

Last change on this file since 14480 was 14480, checked in by Ivan Raikov, 10 years ago

ported random-test to Chicken 4

File size: 11.6 KB
Line 
1;;
2;;
3;; Perform some simple randomness tests on a sequence of numerical
4;; values.
5;;
6;; Based on the ENT program by John Walker
7;;
8;;
9;; Copyright 2007-2009 Ivan Raikov.
10;;
11;;
12;; This program is free software: you can redistribute it and/or modify
13;; it under the terms of the GNU General Public License as published by
14;; the Free Software Foundation, either version 3 of the License, or (at
15;; your option) any later version.
16;;
17;; This program is distributed in the hope that it will be useful, but
18;; WITHOUT ANY WARRANTY; without even the implied warranty of
19;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
20;; General Public License for more details.
21;;
22;; A full copy of the GPL license can be found at
23;; <http://www.gnu.org/licenses/>."))))
24;;
25;;
26
27(module random-test
28   
29  ( make-random-test format-random-stats )
30
31  (import scheme chicken foreign extras data-structures srfi-4)
32 
33  (require-extension srfi-4 easyffi)
34
35
36;;
37;;
38;; Chi-Square Test
39;;
40;;    In general, the Chi-Square distribution for an experiment with k
41;; possible outcomes, performed n times, in which Y1, Y2,... Yk are
42;; the number of experiments which resulted in each possible outcome,
43;; and probabilities of each outcome p1, p2,... pk, is given as:
44;;
45;;   \Chi^2 = \sum_{1 <= i <= k}{(Y_i - m * p_i)^2 / np_i}
46;;
47;; \Chi^2 will grow larger as the measured results diverge from those
48;; expected by pure chance. The probability Q that a Chi-Square value
49;; calculated for an experiment with d degrees of freedom (where
50;; d=k-1, one less the number of possible outcomes) is due to chance
51;; is:
52;;
53;;   Q(\Chi^2,d) =
54;;     [2^{d/2} * \Gamma(d/2)]^{-1} *
55;;     \int_{\Chi^2}^{\infty}(t)^{d/2-1} * exp(-t/2) * dt
56;;
57;; Where Gamma is the generalisation of the factorial function to real
58;; and complex arguments:
59;;
60;;   Gamma(x) = \int_{0}^{\infty} t^{x-1} * exp(-t) * dt
61;;
62;;
63;; There is no closed form solution for Q, so it must be evaluated
64;; numerically.  Note that the probability calculated from the \Chi^2
65;; is an approximation which is valid only for large values of n, and
66;; is therefore only meaningful when calculated from a large number of
67;; independent experiments.
68;;
69;; In this implementation, the Chi-Square distribution is calculated
70;; for the list of values given as argument to the random-test
71;; procedure and expressed as an absolute number and a percentage
72;; which indicates how frequently a truly random sequence would exceed
73;; the value calculated.
74;;
75;; The percentage can be interpreted as the degree to which the
76;; sequence tested is suspected of being non-random. If the percentage
77;; is greater than 99% or less than 1%, the sequence is almost
78;; certainly not random. If the percentage is between 99% and 95% or
79;; between 1% and 5%, the sequence is suspect. Percentages between 90%
80;; and 95% and 5% and 10% indicate the sequence is "almost
81;; suspect".
82;;
83;;
84;; Arithmetic Mean
85;;
86;;     This is simply the result of summing the all the values in the
87;; sequence and dividing by the sequence length. If the data are close
88;; to random, the mean should be about (2^b - 1)/2 where b is the
89;; number of bits used to represent a value. If the mean departs from
90;; this value, the values are consistently high or low.
91;;
92;;
93;; Monte Carlo Value for Pi
94;;
95;;     Each pair of two values in the input sequence is used as X and
96;; Y coordinates within a square with side N (the length of the input
97;; sequence). If the distance of the randomly-generated point is less
98;; than the radius of a circle inscribed within the square, the pair
99;; of values is considered a "hit". The percentage of hits can be used
100;; to calculate the value of Pi:
101;;
102;;
103;;    # points within circle     1/4 * pi * r^2
104;;    ----------------------  =  -------------- = 1/4 * pi
105;;    # points within square          r^2
106;;
107;;
108;; Therefore,
109;;                   # points within circle
110;;          pi = 4 * ----------------------
111;;                   # points within square
112;;
113;;
114;;  For very long sequences (this approximation converges very
115;;  slowly), the value will approach the correct value of Pi if the
116;;  sequence is close to random.
117;;
118;;
119;; Serial Correlation Coefficient
120;;
121;;     This statistic measures the extent to which each value in the
122;; sequence depends upon the previous value. For random sequences,
123;; this value (which can be positive or negative) must be close to
124;; zero.
125;;
126;;
127
128; Parse & embed
129#>!
130  double pochisq (double x, double df);
131  double log2 (double x);
132<#
133
134;; Number of bins for Chi-Square test
135(define N 1024)
136
137;; Circle radius for the Monte Carlo calculation of Pi
138(define R N)
139
140;; Number of coordinates for the Monte Carlo calculation of Pi
141(define M 2)
142
143
144(define PI 3.14159265358979323846)
145
146(define (make-random-test . rest)
147  (let-optionals rest ((scar car) (scdr cdr) (snull? null?))
148   (lambda (seq)
149
150     ;; Bins to count occurrences of values
151     (define bins (make-u32vector N 0))
152
153     ;; Monte Carlo coordinates
154     (define mcoord (make-u32vector M 0))
155     
156     (let ((nmin       0.0)   ;; Minimum and maximum values encountered
157           (nmax       0.0)
158           (mp          0)   ;; Monte Carlo accumulator pointer
159           (mcount      0)   ;; Monte Carlo tries
160           (incircle    0))  ;; Monte Carlo circle inside count
161       (let* ;; Serial correlation terms
162           ((sccn0     (and seq (scar seq)))
163            (sccn-1    #f)   
164            (scct1     0.0)   
165            (scct2     0.0)   
166            (scct3     0.0)) 
167         
168         (let loop ((total 0) (seq seq))
169           (if (snull? seq)
170               (random-stats total nmin nmax bins incircle mcount sccn0 scct1 scct2 scct3 sccn-1 )
171               (let ((n  (scar seq)))
172                 ;; Updated minimum and maximum
173                 ;; Update bin counters
174                 (let ((nbin (inexact->exact (modulo n N))))
175                   (u32vector-set! bins nbin (fx+ 1 (u32vector-ref bins nbin)))
176                   (let ((bi (u32vector-ref bins nbin)))
177                     (if (< bi nmin) (set! nmin bi))
178                     (if (> bi nmax) (set! nmax bi)))
179                   ;; Update inside / outside circle counts for Monte Carlo
180                   ;; calculation of Pi
181                   (u32vector-set! mcoord mp nbin)
182                   (set! mp (fx+ 1 mp))
183                   (if (fx<= M mp)
184                      (begin (set! mcount (fx+ mcount 1))
185                             (set! mp 0)
186                             (let ((x  (u32vector-ref mcoord 0))
187                                   (y  (u32vector-ref mcoord 1)))
188                               ;; determine if the random "coordinates" are
189                               ;; inside or outside the circle
190                               (let ((dist  (sqrt (+ (expt x 2) (expt y 2)))))
191                                 (if (<= dist R) (set! incircle (fx+ 1 incircle)))))))
192                  ;; Update calculation of serial correlation coefficient
193                   (if sccn-1 (set! scct1 (+ scct1  (* sccn-1 n))))
194                   (set! scct2 (+ scct2 n))
195                   (set! scct3 (+ scct3 (* n n)))
196                   (set! sccn-1 n)
197                   (loop (fx+ 1 total) (scdr seq)))))))))))
198
199(define (random-stats total nmin nmax bins incircle mcount sccn0 scct1 scct2 scct3 sccn)
200
201  ;; Probabilities per bin for entropy
202  (define probs (make-f64vector N 0))
203
204  ;; Compute serial correlation coefficient
205  (let* ((scct1     (+ scct1 (* sccn sccn0)))
206         (scct2     (* scct2 scct2))
207         (sccdenom  (- (* total scct3) scct2))
208         (scc       (if (zero? sccdenom)
209                        -10000.0
210                        (/ (- (* total scct1) scct2) sccdenom))))
211         
212    ;; Scan bins and calculate probability for each bin and
213    ;; Chi-Square distribution
214    (let ((cexp      (/ total N))) ;; Expected count per bin
215      (let-values (((chisq pochisq datasum) 
216          (let loop ((i 0) (chisq 0.0) (datasum 0.0))
217            (if (fx< i N)
218                (let ((bi (u32vector-ref bins i)))
219                  (let ((a (- bi cexp)) (p (/ bi total)))
220                    (f64vector-set! probs i p) 
221                    (loop (fx+ 1 i) 
222                          (+ chisq (/ (* a a) cexp))
223                          (+ datasum (* i bi)))))
224                (values chisq (pochisq chisq (- total 1)) datasum)))))
225
226       ;; Calculate entropy
227       (let ((ent  (let loop ((i 0) (ent 0))
228                     (if (fx< i N)
229                         (let ((p (f64vector-ref probs i)))
230                           (loop (fx+ 1 i) (if (< 0.0 p) (+ ent (* p (log2 (/ 1.0 p)))) ent)))
231                         ent))))
232
233         ;; Calculate Monte Carlo value for Pi from percentage of hits
234         ;; within the circle
235         (let ((montepi (* 4 (/ incircle mcount))))
236           `((total    . ,total)
237             (ent      . ,ent)
238             (chisq    . ,chisq)
239             (pochisq  . ,pochisq)
240             (mean     . ,(/ datasum total))
241             (min      . ,nmin)
242             (max      . ,nmax)
243             (montepi  . ,montepi)
244             (scc      . ,scc))))))))
245
246(define (format-random-stats out stats)
247  (fprintf out "Entropy is ~S.\n" (alist-ref 'ent stats))
248  (fprintf out "Chi-Square distribution for ~S samples is ~S (~S percent random).\n"
249           (alist-ref 'total stats) 
250           (alist-ref 'chisq stats) 
251           (* 100 (alist-ref 'pochisq stats)))
252  (fprintf out "Arithmetic mean value of the sequence is ~S (~S = random).\n"
253           (alist-ref 'mean stats)
254           (/ (alist-ref 'total stats) N))
255  (fprintf out "Monte-Carlo value for Pi is ~S (error ~S percent).\n"
256           (alist-ref 'montepi stats)
257           (* 100 (/ (abs (- PI (alist-ref 'montepi stats))) PI)))
258  (fprintf out "Serial correlation coefficient is ~S.\n"
259           (let ((scc (alist-ref 'scc stats)))
260             (if (>= scc -99999) scc "undefined (all values equal)"))))
261
262
263; Include into generated code, but don't parse:
264#>
265
266#include <math.h>
267double poz  (double z);
268
269#define LOG_SQRT_PI     0.5723649429247000870717135 /* log (sqrt (pi)) */
270#define I_SQRT_PI       0.5641895835477562869480795 /* 1 / sqrt (pi) */
271#define BIGX           20.0         /* max value to represent exp (x) */
272#define ex(x)          (((x) < -BIGX) ? 0.0 : exp (x))
273
274
275/* ALGORITHM Compute probability of chi square value.
276;;      Adapted from:
277;;              Hill, I. D. and Pike, M. C.  Algorithm 299
278;;              Collected Algorithms for the CACM 1967 p. 243
279;;      Updated for rounding errors based on remark in
280;;              ACM TOMS June 1985, page 185
281*/
282
283/* 
284;;  x: obtained chi-square value
285;;  df: degrees of freedom
286*/
287
288double pochisq (double x, double df)
289{
290
291   double a, y, s;
292   double e, c, z;
293   int  even;     /* true if df is an even number */
294       
295   if (x <= 0.0 || df < 1)
296     return (1.0);
297       
298     a = 0.5 * x;
299     even = (2*(df/2)) == df;
300     if (df > 1)
301     y = ex (-a);
302     s = (even ? y : (2.0 * poz (-sqrt (x))));
303     if (df > 2)
304     {
305      x = 0.5 * (df - 1.0);
306        z = (even ? 1.0 : 0.5);
307        if (a > BIGX)
308        {
309         e = (even ? 0.0 : LOG_SQRT_PI);
310           c = log (a);
311           while (z <= x)
312           {
313            e = log (z) + e;
314              s += ex (c*z-a-e);
315              z += 1.0;
316              }
317           return (s);
318           }
319        else
320        {
321         e = (even ? 1.0 : (I_SQRT_PI / sqrt (a)));
322           c = 0.0;
323           while (z <= x)
324           {
325            e = e * (a / z);
326              c = c + e;
327              z += 1.0;
328              }
329           return (c * y + s);
330           }
331        }
332     else
333     return (s);
334}
335
336/* Maximum meaningful z value */
337#define Z_MAX 6.0   
338
339
340/*
341;; FUNCTION poz: probability of normal z value
342;;
343;;      Adapted from a polynomial approximation in:
344;;              Ibbetson D, Algorithm 209
345;;              Collected Algorithms of the CACM 1963 p. 616
346;;      Note:
347;;
348;;              This routine has six digit accuracy, so it is only
349;;              useful for absolute z values < 6.  For z values >= to
350;;              6.0, poz() returns 0.0.
351;;
352*/
353double poz (double z)
354{
355
356 double y, x, w;
357       
358 if (z == 0.0)
359    x = 0.0;
360 else
361 {
362    y = 0.5 * fabs (z);
363    if (y >= (Z_MAX * 0.5))
364    {
365       x = 1.0;
366    }
367    else if (y < 1.0)
368         {
369            w = y * y;
370            x = ((((((((0.000124818987 * w
371                       -0.001075204047) * w +0.005198775019) * w
372                       -0.019198292004) * w +0.059054035642) * w
373                       -0.151968751364) * w +0.319152932694) * w
374                       -0.531923007300) * w +0.797884560593) * y * 2.0;
375          }
376          else
377          {
378             y -= 2.0;
379               x = (((((((((((((-0.000045255659 * y
380                                +0.000152529290) * y -0.000019538132) * y
381                                -0.000676904986) * y +0.001390604284) * y
382                                -0.000794620820) * y -0.002034254874) * y
383                                +0.006549791214) * y -0.010557625006) * y
384                                +0.011630447319) * y -0.009279453341) * y
385                                +0.005353579108) * y -0.002141268741) * y
386                                +0.000535310849) * y +0.999936657524;
387          }
388 }
389
390 return (z > 0.0 ? ((x + 1.0) * 0.5) : ((1.0 - x) * 0.5));
391}
392
393<#
394)
Note: See TracBrowser for help on using the repository browser.