Changeset 36371 in project


Ignore:
Timestamp:
08/24/18 20:39:36 (3 months ago)
Author:
iraikov
Message:

statistics: yasos collections branch

Location:
release/5/statistics/branches
Files:
1 added
1 edited
1 copied

Legend:

Unmodified
Added
Removed
  • release/5/statistics/branches/collections/statistics.scm

    r36303 r36371  
    178178    )
    179179
    180   (import scheme (chicken base) (chicken foreign) (chicken format) (chicken keyword) (chicken sort)
    181           srfi-1 (only srfi-13 string<) srfi-25 srfi-69 vector-lib)
     180  (import scheme (chicken base) (chicken foreign) (chicken format)
     181          (chicken keyword) (prefix (only srfi-1 fold iota reverse) list.)
     182          (only srfi-13 string<) srfi-25 srfi-69 vector-lib
     183          yasos yasos-collections data-series)
    182184
    183185  ;; ---------------------------------------------------------------------------
     
    271273  (define PI 3.1415926535897932385)
    272274
     275  (define (position item sequence)
     276    (let ((index -1))
     277      (do-items-while (lambda (x)
     278                        (let ((res (equal? item (cadr x))))
     279                          (if res (set! index (car x)))
     280                          (not res)))
     281                      sequences)
     282      (if (> index -1)
     283        index
     284        (error "Position: item not in sequence"))))
     285
    273286  (define (array-shape arr)
    274287    (let ((r (array-rank arr)))
     
    295308  ;; belong in each.
    296309  (define (bin-and-count sequence n)
    297     (let ((smallest  (fold min (car sequence) (cdr sequence)))
    298           (increment (/ (range sequence)
    299                         n))
     310    (let ((smallest  (reduce* min sequence))
     311          (increment (/ (range sequence) n))
    300312          (bins (make-vector n 0)))
    301313      (let loop ((bin 0))
     
    381393    (* 0.5 (log (/ (+ 1 r) (- 1 r)))))
    382394
    383   ;; taken from TSPL 3
    384   (define (list-sort comparison items)
    385     (define (dosort pred? ls n)
    386       (if (= n 1)
    387         (list (car ls))
    388         (let ((i (quotient n 2)))
    389           (domerge pred?
    390                    (dosort pred? ls i)
    391                    (dosort pred? (list-tail ls i) (- n i))))))
    392     (define (domerge pred? l1 l2)
    393       (cond ((null? l1) l2)
    394             ((null? l2) l1)
    395             ((pred? (car l2) (car l1))
    396              (cons (car l2) (domerge pred? l1 (cdr l2))))
    397             (else
    398               (cons (car l1) (domerge pred? (cdr l1) l2)))))
    399     (if (null? items)
    400       items
    401       (dosort comparison items (length items))))
    402395
    403396  ;; PERMUTATIONS
     
    405398  ;; Rosner 88
    406399  (define (permutations n k)
    407     (fold * 1 (iota k n -1)))
    408 
    409   (define (position item items)
    410     (let ((index (list-index (lambda (x) (equal? item x))
    411                              items)))
    412       (if (number? index)
    413         index
    414         (error "Position: item not in list"))))
     400    (list.fold * 1 (list.iota k n -1)))
    415401
    416402  (define (cumsum sequence)
    417     (if (null? sequence)
     403    (if (empty? sequence)
    418404        0
    419         (reverse (fold (lambda (x ax) (cons (+ x (car ax)) ax)) (list (car sequence)) (cdr sequence)))))
     405        (list.reverse (reduce* (lambda (x ax) (cons (+ x (car ax)) ax))
     406                               sequence))))
    420407
    421408  (define (sign x)
     
    432419    x)
    433420
    434   (define (vector-sort comparison items)
    435     (list->vector (list-sort comparison (vector->list items))))
    436 
    437421  ;; RANDOM-NORMAL
    438422  ;; returns a random number with mean and standard-distribution as specified.
     
    440424    (+ mean (* sd (/ (random-uniform-int 1000000) 1000000))))
    441425
    442   ;; RANDOM-PICK
    443   ;; Random selection from list
    444   (define (random-pick items)
    445     (if (and (list? items) (not (null? items)))
    446       (list-ref items (random-uniform-int (length items)))
    447       #f))
    448426
    449427  ;; RANDOM-SAMPLE
     
    451429  ;; If N is equal to or greater than the length of the sequence, return
    452430  ;; the entire sequence.
    453   (define (random-sample n items)
    454     (cond ((<= n 0)
    455            '())
    456           ((>= n (length items))
    457            items)
    458           (else
    459             (let loop ((remaining items)
    460                        (kept '()))
    461               (if (= (length kept) n)
    462                 kept
    463                 (let ((one (random-pick remaining)))
    464                   (loop (delete one remaining)
    465                         (cons one kept))))))))
    466 
     431  (define (random-sample n sequence)
     432    (let ((result (make-vector n)))
     433      (cond ((<= n 0)
     434             (list->vector '()))
     435            (else
     436             (let loop ((remaining (gen-elts sequence))
     437                        (l 0))
     438               (let ((item (remaining)))
     439                 (if (not item)
     440                     result
     441                     (begin
     442                       (if (< l n)
     443                           (vector-set! result l item)
     444                           (let ((i (random-uniform-int n)))
     445                             (vector-set! result i item)))
     446                       (loop remaining (+ 1 l))
     447                       ))
     448                 ))
     449             ))
     450      ))
     451
     452 
    467453  ;; RANDOM-WEIGHTED-SAMPLE Return a random sample of size M from
    468454  ;; sequence of length N, without replacement, where each element has
    469455  ;; a defined probability of selection (weight) W.  If M is equal to
    470456  ;; or greater to N, return the entire sequence.
    471   (define (random-weighted-sample m items weights)
    472     (let ((n (length items)))
     457  (define (random-weighted-sample m sequence weights)
     458    (let ((n (size sequence)))
    473459      (cond ((<= m 0) '())
    474             ((>= m n) items)
     460            ((>= m n) sequence)
    475461            (else
    476              (let* ((keys (map (lambda (w) (expt (random-uniform) (/ 1 w))) weights))
    477                     (sorted-items (list-sort (lambda (x y) (> (car x) (car y))) (zip keys items))))
    478                (map cadr (take sorted-items m)))
     462             (let* ((keys (map-elts (lambda (w) (expt (random-uniform) (/ 1 w))) weights))
     463                    (sorted-items (sort (lambda (x y) (> (car (cadr x)) (car (cadr y))))
     464                                        (zip keys sequence))))
     465               (elt-take sorted-items m))
    479466             ))
    480467      ))
     
    487474  ;; Rosner 10
    488475  (define (mean sequence)
    489     (if (null? sequence)
     476    (if (empty? sequence)
    490477      0
    491       (/ (fold + 0 sequence) (length sequence))))
     478      (/ (reduce (lambda (x ax) (+ (cadr x) ax)) 0 sequence) (size sequence))))
    492479
    493480  ;; MEDIAN
     
    500487  ;; Returns two values: a list of the modes and the number of times they occur
    501488  (define (mode sequence)
    502     (if (null? sequence)
     489    (if (empty? sequence)
    503490      (error "Mode: Sequence must not be null")
    504491      (let ((count-table (make-hash-table eqv?))
    505492            (modes '())
    506493            (mode-count 0))
    507         (for-each (lambda (item)
    508                     (hash-table-set! count-table
    509                                      item
    510                                      (+ 1 (hash-table-ref count-table item (lambda () 0)))))
    511                   sequence)
    512         (for-each (lambda (key)
    513                     (let ((val (hash-table-ref count-table key (lambda () #f))))
    514                       (cond ((> val mode-count) ; keep mode
    515                              (set! modes (list key))
    516                              (set! mode-count val))
    517                             ((= val mode-count) ; store multiple modes
    518                              (set! modes (cons key modes))))))
    519                   (hash-table-keys count-table))
    520         (cond ((every number? modes) (set! modes (list-sort < modes)))
    521               ((every string? modes) (set! modes (list-sort string< modes)))
     494        (for-each-elt
     495         (lambda (item)
     496           (hash-table-set! count-table
     497                            item
     498                            (+ 1 (hash-table-ref count-table item (lambda () 0)))))
     499         sequence)
     500        (for-each
     501         (lambda (key)
     502           (let ((val (hash-table-ref count-table key (lambda () #f))))
     503             (cond ((> val mode-count) ; keep mode
     504                    (set! modes (list key))
     505                    (set! mode-count val))
     506                   ((= val mode-count) ; store multiple modes
     507                    (set! modes (cons key modes))))))
     508         (hash-table-keys count-table))
     509        (cond ((every number? modes) (set! modes (sort < modes)))
     510              ((every string? modes) (set! modes (sort string< modes)))
    522511              )
    523512        (values modes mode-count))))
     
    527516  (define (geometric-mean sequence . args)
    528517    (let ((base (if (null? args) 10 (car args))))
    529       (expt base (mean (map (lambda (x) (/ (log x)
    530                                            (log base)))
    531                             sequence)))))
     518      (expt base (mean (map-elts (lambda (x) (/ (log x)
     519                                                (log base)))
     520                                 sequence)))))
    532521
    533522  ;; RANGE
    534523  ;; Rosner 18
    535524  (define (range sequence)
    536     (if (null? sequence)
     525    (if (empty? sequence)
    537526      0
    538       (- (fold max (car sequence) (cdr sequence))
    539          (fold min (car sequence) (cdr sequence)))))
     527      (- (reduce* max sequence)
     528         (reduce* min sequence))))
    540529
    541530  ;; PERCENTILE
    542531  ;; Rosner 19
    543532  (define (percentile sequence percent)
    544     (if (null? sequence)
     533    (if (empty? sequence)
    545534      (error "Percentile: Sequence must not be null")
    546       (let* ((sorted-vec (vector-sort < (apply vector sequence)))
    547              (n (vector-length sorted-vec))
     535      (let* ((sorted-items (sort (lambda (i vi j vj) (< vi vj)) sequence))
     536             (n (size sorted-items))
    548537             (k (* n (/ percent 100)))
    549538             (floor-k (floor k)))
    550539        (if (= k floor-k)
    551           (/ (+ (vector-ref sorted-vec k)
    552                 (vector-ref sorted-vec (- k 1)))
     540          (/ (+ (elt-ref sorted-items k)
     541                (elt-ref sorted-items (- k 1)))
    553542             2)
    554           (vector-ref sorted-vec floor-k)))))
     543          (elt-ref sorted-items floor-k)))))
    555544
    556545  ;; VARIANCE
    557546  ;; Rosner 21
    558547  (define (variance sequence)
    559     (if (< (length sequence) 2)
     548    (if (< (size sequence) 2)
    560549      (error "variance: sequence must contain at least two elements")
    561550      (let ((mean1 (mean sequence)))
    562         (/ (fold + 0 (map (lambda (x) (square (- mean1 x))) sequence))
    563            (- (length sequence) 1)))))
     551        (/ (reduce (lambda (x ax) (+ (cadr x) ax)) 0
     552                   (map (lambda (x) (square (- mean1 x))) sequence))
     553           (- (size sequence) 1)))))
    564554
    565555  ;; STANDARD-DEVIATION
     
    11701160    (let ((tails (get-keyword #:tails args (lambda () ':both))))
    11711161      (let* ((nonzero-differences (filter (lambda (n) (not (zero? n))) differences))
    1172              (sorted-list (list-sort (lambda (x y) (< (car x) (car y)))
    1173                                      (map (lambda (dif)
    1174                                             (list (abs dif)
    1175                                                   (sign dif)))
    1176                                           nonzero-differences)))
     1162             (sorted-items (sort (lambda (x y) (< (car x) (car y)))
     1163                                 (map-items (lambda (dif)
     1164                                              (list (abs dif)
     1165                                                    (sign dif)))
     1166                                            nonzero-differences)))
    11771167             (distinct-values (delete-duplicates (map car sorted-list)))
    11781168             (ties '()))
    1179         (when (< (length nonzero-differences) 16)
     1169        (when (< (size nonzero-differences) 16)
    11801170          (error
    11811171            "This Wilcoxon Signed-Rank Test (normal approximation method) requires nonzero N > 15"))
     
    14701460  ;; the significance of the difference of the slope from 0.
    14711461  (define (linear-regression points)
    1472     (unless (> (length points) 2)
     1462    (unless (> (size points) 2)
    14731463      (error "Requires at least three points"))
    1474     (let ((xs (map car points))
    1475           (ys (map cadr points)))
     1464    (let ((xs (elt-map car points))
     1465          (ys (elt-map cadr points)))
    14761466      (let* ((x-bar (mean xs))
    14771467             (y-bar (mean ys))
    1478              (n (length points))
    1479              (Lxx (fold + 0
     1468             (n (size points))
     1469             (Lxx (reduce + 0
    14801470                        (map (lambda (xi) (square (- xi x-bar))) xs)))
    1481              (Lyy (fold + 0
     1471             (Lyy (reduce + 0
    14821472                        (map (lambda (yi) (square (- yi y-bar))) ys)))
    1483              (Lxy (fold + 0
    1484                         (map (lambda (xi yi) (* (- xi x-bar) (- yi y-bar)))
    1485                              xs ys)))
     1473             (Lxy (reduce + 0
     1474                          (elt-map (lambda (point) (let ((xi (car point))
     1475                                                         (yi (car point)))
     1476                                                     (* (- xi x-bar) (- yi y-bar))))
     1477                             points)))
    14861478             (b (if (zero? Lxx) 0 (/ Lxy Lxx)))
    14871479             (a (- y-bar (* b x-bar)))
     
    14991491  ;; Also called Pearson Correlation
    15001492  (define (correlation-coefficient points)
    1501     (let* ((xs (map car points))
    1502            (ys (map cadr points))
     1493    (let* ((xs (elt-map car points))
     1494           (ys (elt-map cadr points))
    15031495           (x-bar (mean xs))
    15041496           (y-bar (mean ys)))
    1505       (/ (fold + 0 (map (lambda (xi yi) (* (- xi x-bar)
    1506                                            (- yi y-bar)))
    1507                         xs
    1508                         ys))
    1509          (sqrt (* (fold + 0 (map (lambda (xi) (square (- xi x-bar)))
    1510                                  xs))
    1511                   (fold + 0 (map (lambda (yi) (square (- yi y-bar)))
    1512                                  ys)))))))
     1497      (/ (reduce + 0 (elt-map (lambda (point)
     1498                                (let ((xi (car point))
     1499                                      (yi (cadr point)))
     1500                                  (* (- xi x-bar) (- yi y-bar))))
     1501                              points))
     1502         (sqrt (* (reduce + 0 (elt-map (lambda (xi) (square (- xi x-bar)))
     1503                                       xs))
     1504                  (reduce + 0 (elt-map (lambda (yi) (square (- yi y-bar)))
     1505                                       ys)))))))
    15131506
    15141507
     
    15451538  ;; and its significance.
    15461539  (define (spearman-rank-correlation points)
    1547     (let ((xis (map car points))
    1548           (yis (map cadr points)))
    1549       (let* ((n (length points))
    1550              (sorted-xis (list-sort < xis))
    1551              (sorted-yis (list-sort < yis))
     1540    (let ((xis (elt-map car points))
     1541          (yis (elt-map cadr points)))
     1542      (let* ((n (size points))
     1543             (sorted-xis (sort (lambda (xi x yi y) (< x y)) xis))
     1544             (sorted-yis (sort (lambda (xi x yi y) (< x y)) yis))
    15521545             (average-x-ranks (map (lambda (x) (average-rank x sorted-xis)) xis))
    15531546             (average-y-ranks (map (lambda (y) (average-rank y sorted-yis)) yis))
    15541547             (mean-x-rank (mean average-x-ranks))
    15551548             (mean-y-rank (mean average-y-ranks))
    1556              (Lxx (fold + 0
    1557                         (map (lambda (xi-rank) (square (- xi-rank mean-x-rank)))
    1558                              average-x-ranks)))
    1559              (Lyy (fold + 0
    1560                         (map (lambda (yi-rank) (square (- yi-rank mean-y-rank)))
    1561                              average-y-ranks)))
    1562              (Lxy (fold + 0
    1563                         (map (lambda (xi-rank yi-rank)
    1564                                (* (- xi-rank mean-x-rank)
    1565                                   (- yi-rank mean-y-rank)))
    1566                              average-x-ranks average-y-ranks)))
     1549             (Lxx (reduce + 0
     1550                        (elt-map (lambda (xi-rank) (square (- xi-rank mean-x-rank)))
     1551                                 average-x-ranks)))
     1552             (Lyy (reduce + 0
     1553                          (elt-map (lambda (yi-rank) (square (- yi-rank mean-y-rank)))
     1554                                   average-y-ranks)))
     1555             (Lxy (reduce + 0
     1556                          (elt-map (lambda (xi-rank yi-rank)
     1557                                     (* (- xi-rank mean-x-rank)
     1558                                        (- yi-rank mean-y-rank)))
     1559                                   average-x-ranks average-y-ranks)))
    15671560             (rs (if (> (* Lxx Lyy) 0) (/ Lxy (sqrt (* Lxx Lyy))) 0)) ; TODO: think about rs = 1
    15681561             (ts (if (= 1 rs) 1 (/ (* rs (sqrt (- n 2))) (sqrt (- 1 (square rs))))))
Note: See TracChangeset for help on using the changeset viewer.