Changeset 36371 in project

Timestamp:

08/24/18 20:39:36 (11 months ago)

Author:

iraikov

Message:

statistics: yasos collections branch

Location:

release/5/statistics/branches

Files:

• . (added)
• collections (copied) (copied from release/5/statistics/trunk)
• collections/statistics.scm (modified) (15 diffs)

release/5/statistics/branches/collections/statistics.scm

@@ -178 +178 @@
     )
 
-  (import scheme (chicken base) (chicken foreign) (chicken format) (chicken keyword) (chicken sort)
-          srfi-1 (only srfi-13 string<) srfi-25 srfi-69 vector-lib)
+  (import scheme (chicken base) (chicken foreign) (chicken format)
+          (chicken keyword) (prefix (only srfi-1 fold iota reverse) list.)
+          (only srfi-13 string<) srfi-25 srfi-69 vector-lib
+          yasos yasos-collections data-series)
 
   ;; ---------------------------------------------------------------------------
@@ -271 +273 @@
   (define PI 3.1415926535897932385)
 
+  (define (position item sequence)
+    (let ((index -1))
+      (do-items-while (lambda (x)
+                        (let ((res (equal? item (cadr x))))
+                          (if res (set! index (car x)))
+                          (not res)))
+                      sequences)
+      (if (> index -1)
+        index
+        (error "Position: item not in sequence"))))
+
   (define (array-shape arr)
     (let ((r (array-rank arr)))
@@ -295 +308 @@
   ;; belong in each.
   (define (bin-and-count sequence n)
-    (let ((smallest  (fold min (car sequence) (cdr sequence)))
-          (increment (/ (range sequence)
-                        n))
+    (let ((smallest  (reduce* min sequence))
+          (increment (/ (range sequence) n))
           (bins (make-vector n 0)))
       (let loop ((bin 0))
@@ -381 +393 @@
     (* 0.5 (log (/ (+ 1 r) (- 1 r)))))
 
-  ;; taken from TSPL 3
-  (define (list-sort comparison items)
-    (define (dosort pred? ls n)
-      (if (= n 1)
-        (list (car ls))
-        (let ((i (quotient n 2)))
-          (domerge pred?
-                   (dosort pred? ls i)
-                   (dosort pred? (list-tail ls i) (- n i))))))
-    (define (domerge pred? l1 l2)
-      (cond ((null? l1) l2)
-            ((null? l2) l1)
-            ((pred? (car l2) (car l1))
-             (cons (car l2) (domerge pred? l1 (cdr l2))))
-            (else
-              (cons (car l1) (domerge pred? (cdr l1) l2)))))
-    (if (null? items)
-      items
-      (dosort comparison items (length items))))
 
   ;; PERMUTATIONS
@@ -405 +398 @@
   ;; Rosner 88
   (define (permutations n k)
-    (fold * 1 (iota k n -1)))
-
-  (define (position item items)
-    (let ((index (list-index (lambda (x) (equal? item x))
-                             items)))
-      (if (number? index)
-        index
-        (error "Position: item not in list"))))
+    (list.fold * 1 (list.iota k n -1)))
 
   (define (cumsum sequence)
-    (if (null? sequence)
+    (if (empty? sequence)
         0
-        (reverse (fold (lambda (x ax) (cons (+ x (car ax)) ax)) (list (car sequence)) (cdr sequence)))))
+        (list.reverse (reduce* (lambda (x ax) (cons (+ x (car ax)) ax))
+                               sequence))))
 
   (define (sign x)
@@ -432 +419 @@
     x)
 
-  (define (vector-sort comparison items)
-    (list->vector (list-sort comparison (vector->list items))))
-
   ;; RANDOM-NORMAL
   ;; returns a random number with mean and standard-distribution as specified.
@@ -440 +424 @@
     (+ mean (* sd (/ (random-uniform-int 1000000) 1000000))))
 
-  ;; RANDOM-PICK
-  ;; Random selection from list
-  (define (random-pick items)
-    (if (and (list? items) (not (null? items)))
-      (list-ref items (random-uniform-int (length items)))
-      #f))
 
   ;; RANDOM-SAMPLE
@@ -451 +429 @@
   ;; If N is equal to or greater than the length of the sequence, return 
   ;; the entire sequence.
-  (define (random-sample n items)
-    (cond ((<= n 0) 
-           '())
-          ((>= n (length items)) 
-           items)
-          (else
-            (let loop ((remaining items)
-                       (kept '()))
-              (if (= (length kept) n)
-                kept
-                (let ((one (random-pick remaining)))
-                  (loop (delete one remaining)
-                        (cons one kept))))))))
-
+  (define (random-sample n sequence)
+    (let ((result (make-vector n)))
+      (cond ((<= n 0) 
+             (list->vector '()))
+            (else
+             (let loop ((remaining (gen-elts sequence))
+                        (l 0))
+               (let ((item (remaining)))
+                 (if (not item)
+                     result
+                     (begin
+                       (if (< l n)
+                           (vector-set! result l item)
+                           (let ((i (random-uniform-int n)))
+                             (vector-set! result i item)))
+                       (loop remaining (+ 1 l))
+                       ))
+                 ))
+             ))
+      ))
+
+  
   ;; RANDOM-WEIGHTED-SAMPLE Return a random sample of size M from
   ;; sequence of length N, without replacement, where each element has
   ;; a defined probability of selection (weight) W.  If M is equal to
   ;; or greater to N, return the entire sequence.
-  (define (random-weighted-sample m items weights)
-    (let ((n (length items)))
+  (define (random-weighted-sample m sequence weights)
+    (let ((n (size sequence)))
       (cond ((<= m 0) '())
-            ((>= m n) items)
+            ((>= m n) sequence)
             (else
-             (let* ((keys (map (lambda (w) (expt (random-uniform) (/ 1 w))) weights))
-                    (sorted-items (list-sort (lambda (x y) (> (car x) (car y))) (zip keys items))))
-               (map cadr (take sorted-items m)))
+             (let* ((keys (map-elts (lambda (w) (expt (random-uniform) (/ 1 w))) weights))
+                    (sorted-items (sort (lambda (x y) (> (car (cadr x)) (car (cadr y))))
+                                        (zip keys sequence))))
+               (elt-take sorted-items m))
              ))
       ))
@@ -487 +474 @@
   ;; Rosner 10
   (define (mean sequence)
-    (if (null? sequence)
+    (if (empty? sequence)
       0
-      (/ (fold + 0 sequence) (length sequence))))
+      (/ (reduce (lambda (x ax) (+ (cadr x) ax)) 0 sequence) (size sequence))))
 
   ;; MEDIAN
@@ -500 +487 @@
   ;; Returns two values: a list of the modes and the number of times they occur
   (define (mode sequence)
-    (if (null? sequence)
+    (if (empty? sequence)
       (error "Mode: Sequence must not be null")
       (let ((count-table (make-hash-table eqv?))
             (modes '())
             (mode-count 0))
-        (for-each (lambda (item) 
-                    (hash-table-set! count-table
-                                     item
-                                     (+ 1 (hash-table-ref count-table item (lambda () 0)))))
-                  sequence)
-        (for-each (lambda (key) 
-                    (let ((val (hash-table-ref count-table key (lambda () #f))))
-                      (cond ((> val mode-count) ; keep mode
-                             (set! modes (list key))
-                             (set! mode-count val))
-                            ((= val mode-count) ; store multiple modes
-                             (set! modes (cons key modes))))))
-                  (hash-table-keys count-table))
-        (cond ((every number? modes) (set! modes (list-sort < modes)))
-              ((every string? modes) (set! modes (list-sort string< modes)))
+        (for-each-elt
+         (lambda (item) 
+           (hash-table-set! count-table
+                            item
+                            (+ 1 (hash-table-ref count-table item (lambda () 0)))))
+         sequence)
+        (for-each
+         (lambda (key) 
+           (let ((val (hash-table-ref count-table key (lambda () #f))))
+             (cond ((> val mode-count) ; keep mode
+                    (set! modes (list key))
+                    (set! mode-count val))
+                   ((= val mode-count) ; store multiple modes
+                    (set! modes (cons key modes))))))
+         (hash-table-keys count-table))
+        (cond ((every number? modes) (set! modes (sort < modes)))
+              ((every string? modes) (set! modes (sort string< modes)))
               )
         (values modes mode-count))))
@@ -527 +516 @@
   (define (geometric-mean sequence . args)
     (let ((base (if (null? args) 10 (car args))))
-      (expt base (mean (map (lambda (x) (/ (log x)
-                                           (log base))) 
-                            sequence)))))
+      (expt base (mean (map-elts (lambda (x) (/ (log x)
+                                                (log base))) 
+                                 sequence)))))
 
   ;; RANGE
   ;; Rosner 18
   (define (range sequence)
-    (if (null? sequence)
+    (if (empty? sequence)
       0
-      (- (fold max (car sequence) (cdr sequence))
-         (fold min (car sequence) (cdr sequence)))))
+      (- (reduce* max sequence)
+         (reduce* min sequence))))
 
   ;; PERCENTILE
   ;; Rosner 19
   (define (percentile sequence percent)
-    (if (null? sequence)
+    (if (empty? sequence)
       (error "Percentile: Sequence must not be null")
-      (let* ((sorted-vec (vector-sort < (apply vector sequence)))
-             (n (vector-length sorted-vec))
+      (let* ((sorted-items (sort (lambda (i vi j vj) (< vi vj)) sequence))
+             (n (size sorted-items))
              (k (* n (/ percent 100)))
              (floor-k (floor k)))
         (if (= k floor-k)
-          (/ (+ (vector-ref sorted-vec k)
-                (vector-ref sorted-vec (- k 1)))
+          (/ (+ (elt-ref sorted-items k)
+                (elt-ref sorted-items (- k 1)))
              2)
-          (vector-ref sorted-vec floor-k)))))
+          (elt-ref sorted-items floor-k)))))
 
   ;; VARIANCE
   ;; Rosner 21
   (define (variance sequence)
-    (if (< (length sequence) 2)
+    (if (< (size sequence) 2)
       (error "variance: sequence must contain at least two elements")
       (let ((mean1 (mean sequence)))
-        (/ (fold + 0 (map (lambda (x) (square (- mean1 x))) sequence))
-           (- (length sequence) 1)))))
+        (/ (reduce (lambda (x ax) (+ (cadr x) ax)) 0
+                   (map (lambda (x) (square (- mean1 x))) sequence))
+           (- (size sequence) 1)))))
 
   ;; STANDARD-DEVIATION
@@ -1170 +1160 @@
     (let ((tails (get-keyword #:tails args (lambda () ':both))))
       (let* ((nonzero-differences (filter (lambda (n) (not (zero? n))) differences))
-             (sorted-list (list-sort (lambda (x y) (< (car x) (car y)))
-                                     (map (lambda (dif)
-                                            (list (abs dif)
-                                                  (sign dif)))
-                                          nonzero-differences)))
+             (sorted-items (sort (lambda (x y) (< (car x) (car y)))
+                                 (map-items (lambda (dif)
+                                              (list (abs dif)
+                                                    (sign dif)))
+                                            nonzero-differences)))
              (distinct-values (delete-duplicates (map car sorted-list)))
              (ties '()))
-        (when (< (length nonzero-differences) 16)
+        (when (< (size nonzero-differences) 16)
           (error 
             "This Wilcoxon Signed-Rank Test (normal approximation method) requires nonzero N > 15"))
@@ -1470 +1460 @@
   ;; the significance of the difference of the slope from 0.
   (define (linear-regression points)
-    (unless (> (length points) 2)
+    (unless (> (size points) 2)
       (error "Requires at least three points"))
-    (let ((xs (map car points))
-          (ys (map cadr points)))
+    (let ((xs (elt-map car points))
+          (ys (elt-map cadr points)))
       (let* ((x-bar (mean xs))
              (y-bar (mean ys))
-             (n (length points))
-             (Lxx (fold + 0
+             (n (size points))
+             (Lxx (reduce + 0
                         (map (lambda (xi) (square (- xi x-bar))) xs)))
-             (Lyy (fold + 0
+             (Lyy (reduce + 0
                         (map (lambda (yi) (square (- yi y-bar))) ys)))
-             (Lxy (fold + 0
-                        (map (lambda (xi yi) (* (- xi x-bar) (- yi y-bar)))
-                             xs ys)))
+             (Lxy (reduce + 0
+                          (elt-map (lambda (point) (let ((xi (car point))
+                                                         (yi (car point)))
+                                                     (* (- xi x-bar) (- yi y-bar))))
+                             points)))
              (b (if (zero? Lxx) 0 (/ Lxy Lxx)))
              (a (- y-bar (* b x-bar)))
@@ -1499 +1491 @@
   ;; Also called Pearson Correlation
   (define (correlation-coefficient points)
-    (let* ((xs (map car points))
-           (ys (map cadr points))
+    (let* ((xs (elt-map car points))
+           (ys (elt-map cadr points))
            (x-bar (mean xs))
            (y-bar (mean ys)))
-      (/ (fold + 0 (map (lambda (xi yi) (* (- xi x-bar)
-                                           (- yi y-bar)))
-                        xs
-                        ys))
-         (sqrt (* (fold + 0 (map (lambda (xi) (square (- xi x-bar)))
-                                 xs))
-                  (fold + 0 (map (lambda (yi) (square (- yi y-bar)))
-                                 ys)))))))
+      (/ (reduce + 0 (elt-map (lambda (point)
+                                (let ((xi (car point))
+                                      (yi (cadr point)))
+                                  (* (- xi x-bar) (- yi y-bar))))
+                              points))
+         (sqrt (* (reduce + 0 (elt-map (lambda (xi) (square (- xi x-bar)))
+                                       xs))
+                  (reduce + 0 (elt-map (lambda (yi) (square (- yi y-bar)))
+                                       ys)))))))
 
 
@@ -1545 +1538 @@
   ;; and its significance. 
   (define (spearman-rank-correlation points)
-    (let ((xis (map car points))
-          (yis (map cadr points)))
-      (let* ((n (length points))
-             (sorted-xis (list-sort < xis))
-             (sorted-yis (list-sort < yis))
+    (let ((xis (elt-map car points))
+          (yis (elt-map cadr points)))
+      (let* ((n (size points))
+             (sorted-xis (sort (lambda (xi x yi y) (< x y)) xis))
+             (sorted-yis (sort (lambda (xi x yi y) (< x y))  yis))
              (average-x-ranks (map (lambda (x) (average-rank x sorted-xis)) xis))
              (average-y-ranks (map (lambda (y) (average-rank y sorted-yis)) yis))
              (mean-x-rank (mean average-x-ranks))
              (mean-y-rank (mean average-y-ranks))
-             (Lxx (fold + 0 
-                        (map (lambda (xi-rank) (square (- xi-rank mean-x-rank)))
-                             average-x-ranks)))
-             (Lyy (fold + 0
-                        (map (lambda (yi-rank) (square (- yi-rank mean-y-rank)))
-                             average-y-ranks)))
-             (Lxy (fold + 0
-                        (map (lambda (xi-rank yi-rank)
-                               (* (- xi-rank mean-x-rank)
-                                  (- yi-rank mean-y-rank)))
-                             average-x-ranks average-y-ranks)))
+             (Lxx (reduce + 0 
+                        (elt-map (lambda (xi-rank) (square (- xi-rank mean-x-rank)))
+                                 average-x-ranks)))
+             (Lyy (reduce + 0
+                          (elt-map (lambda (yi-rank) (square (- yi-rank mean-y-rank)))
+                                   average-y-ranks)))
+             (Lxy (reduce + 0
+                          (elt-map (lambda (xi-rank yi-rank)
+                                     (* (- xi-rank mean-x-rank)
+                                        (- yi-rank mean-y-rank)))
+                                   average-x-ranks average-y-ranks)))
              (rs (if (> (* Lxx Lyy) 0) (/ Lxy (sqrt (* Lxx Lyy))) 0)) ; TODO: think about rs = 1
              (ts (if (= 1 rs) 1 (/ (* rs (sqrt (- n 2))) (sqrt (- 1 (square rs))))))