AutoLISP STDLIB mail archive [part 2 ], 30.Jun.98 - 31.Jul.98
[Date Prev][Date Next] [Thread Prev][Thread Next] [Date Index] [Thread Index] [Author Index]

RANDTST (monotony test and max of 5 test)

Hello Reini,

Send you the last of empirical Knuth tests. Remain spectral test (it
was used e.g. for calculating the number of planes) but there are some
problems with its realization (we should doing operations like
multiplication and addition with accuracy  ~ 80 - 90 bits).

;;;Monotony interval test (test G from Knuth)
;;;Calculates monotony intervals lengths r of sequences
;;;of integers from 0 to (d-1).
(defun rnd-monotony-test(/ n tmax r rr
                           estval obsval DOF chi2)
  ;;n and tmax should be chosen in such a way that estimated numbers
  ;;for every interval length >= 5
  ;;pr = 1/r! - 1/(r+1)!, ptmax = 1/tmax!
  (setq n 2000 tmax 5 DOF (1- tmax))
  ;;calculates observed values
  (setq obsval (calc-monotony-interval n tmax))
  ;;calculate estimated value
  (setq r tmax rr (/ n (float (std-faculty r)))
          estval (list rr))
  (while (> (setq estval (cons (- (* r rr) rr) estval)
            rr (* r rr)
            r (1- r)) 1)
  (setq chi2 (calc-chi2 estval obsval))
  (list chi2 DOF)

;;;Calculates monotony interval statistics.
;;;Parameters: n - interval numbers
(defun calc-monotony-interval(n tmax / count r u uj)
  ;;Initial setting
  (setq count (std-make-list tmax 0))
  (repeat n
    (setq r 0 uj (std-random nil))
    (while (< uj (setq u (std-random nil)))
      (setq uj u r (1+ r)))
    ;;Register length
    (if (>= r tmax) (setq r (1- tmax)))
    (setq count (std-rplace count r (1+ (nth r count))))
    ;;Make dummy call to provide interval independence
    (std-random nil)

(setq r (std-rseq 0.1 1.0 10))
;;;Test "Maximum of 5" (test H in Knuth)
;;;Realization is different (using descrete intervals and chi2
;;;instead Kolmogorov-Smirnov criteria
;;;Returns chi-squared value with DOF = d-1,
;;;where d - group number.
(defun rnd-maxof5-test( / n d umax pr pr5 r
                          estval obsval DOF chi2
  (defun rand() (std-random nil))
  ;;ajdust counters
  (setq n 500 d 10 DOF (1- d))
  (setq obsval (std-make-list d 0)
        pr (std-rseq (/ 1.0 d) 1.0 d)
        pr5 (mapcar '(lambda(p) (expt p 0.2)) pr))
  ;;Calculate observed values
  (repeat n
    (setq umax (apply 'max (std-make-list 5 'rand)))
    (setq r (std-count-if '(lambda(x) (< x umax)) pr5))
    (setq obsval (std-rplace obsval r (1+ (nth r obsval))))
  ;;Calculate estimated values
  (setq estval (std-make-list d (/ (float n) d)))
  (setq chi2 (calc-chi2 estval obsval))
  (list chi2 DOF)

Best regards,