Changeset 36443 in project


Ignore:
Timestamp:
08/26/18 06:03:42 (4 weeks ago)
Author:
iraikov
Message:

statistics: merging yasos-based version into trunk

Location:
release/5/statistics/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • release/5/statistics/trunk/statistics.egg

    r36303 r36443  
    66 (category math)
    77 (build-dependencies compile-file srfi-13 )
    8  (dependencies srfi-1 srfi-25 srfi-69 vector-lib)
     8 (dependencies srfi-1 srfi-25 srfi-69 vector-lib yasos)
    99 (test-dependencies test srfi-1)
    1010 (components (extension statistics (custom-build "build.sh"))))
  • release/5/statistics/trunk/statistics.scm

    r36391 r36443  
    8080    permutations
    8181    random-normal
    82     random-pick
    8382    random-sample
    8483    random-weighted-sample
     
    178177    )
    179178
    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)
     179  (import scheme (chicken base) (chicken foreign) (chicken format)
     180          (chicken keyword) (prefix (chicken sort) list.)
     181          (prefix (only srfi-1 fold iota filter find delete-duplicates every first last) list.)
     182          (only srfi-13 string<) srfi-25 srfi-69 vector-lib
     183          yasos yasos-collections)
    182184
    183185  ;; ---------------------------------------------------------------------------
     
    271273  (define PI 3.1415926535897932385)
    272274
     275 
     276  (define (position item sequence)
     277    (let ((found-item
     278           (g-find
     279            (lambda (x) (equal? item (cadr x)))
     280            (g-map list
     281                   (gen-keys sequence)
     282                   (gen-elts sequence)))))
     283      (if (not (eof-object? found-item))
     284        (car found-item)
     285        (error "Position: item not in sequence"))))
     286
     287  (define (positions item sequence)
     288    (let ((found-items
     289           (g-filter
     290            (lambda (x) (equal? item (cadr x)))
     291            (g-map list
     292                   (gen-keys sequence)
     293                   (gen-elts sequence)))))
     294      (map car found-items)))
     295
    273296  (define (array-shape arr)
    274297    (let ((r (array-rank arr)))
     
    283306  ;; but lisp is 0 based, so add 1!
    284307  (define (average-rank value sorted-values)
    285     (let ((first (position value sorted-values))
    286           (last (- (- (length sorted-values) 1)
    287                    (position value (reverse sorted-values)))))
    288       (+ 1 (if (= first last)
    289              first
    290              (/ (+ first last) 2)))))
    291   (define (average-rank value sorted-values)
    292     (let ((first (position value sorted-values))
    293           (last (- (- (length sorted-values) 1)
    294                    (position value (reverse sorted-values)))))
    295       (let ((result (+ 1 (if (= first last)
    296                              first
    297                              (/ (+ first last) 2)))))
    298         (print "average-rank: value = " value " first = "  first " last = " last " result = " result)
     308    (let* ((idxs (positions value sorted-values)))
     309      (let ((result (+ 1 (mean idxs))))
    299310        result)))
    300311
     
    304315  ;; belong in each.
    305316  (define (bin-and-count sequence n)
    306     (let ((smallest  (fold min (car sequence) (cdr sequence)))
    307           (increment (/ (range sequence)
    308                         n))
     317    (let ((smallest  (reduce* min sequence))
     318          (increment (/ (range sequence) n))
    309319          (bins (make-vector n 0)))
    310320      (let loop ((bin 0))
     
    316326            (vector-set! bins
    317327                         bin
    318                          (length
    319                            (filter
    320                              (lambda (x) (and (>= x (+ smallest (* bin increment)))
    321                                               (< x  (+ smallest (* (+ 1 bin) increment)))))
    322                              sequence)))
     328                         (reduce
     329                          (lambda (x ax)
     330                            (if (and (>= x (+ smallest (* bin increment)))
     331                                     (< x  (+ smallest (* (+ 1 bin) increment))))
     332                                (+ 1 ax) ax)) 0
     333                                sequence))
    323334            (loop (+ 1 bin)))))))
    324335
     
    390401    (* 0.5 (log (/ (+ 1 r) (- 1 r)))))
    391402
    392   ;; taken from TSPL 3
    393   (define (list-sort comparison items)
    394     (define (dosort pred? ls n)
    395       (if (= n 1)
    396         (list (car ls))
    397         (let ((i (quotient n 2)))
    398           (domerge pred?
    399                    (dosort pred? ls i)
    400                    (dosort pred? (list-tail ls i) (- n i))))))
    401     (define (domerge pred? l1 l2)
    402       (cond ((null? l1) l2)
    403             ((null? l2) l1)
    404             ((pred? (car l2) (car l1))
    405              (cons (car l2) (domerge pred? l1 (cdr l2))))
    406             (else
    407               (cons (car l1) (domerge pred? (cdr l1) l2)))))
    408     (if (null? items)
    409       items
    410       (dosort comparison items (length items))))
    411403
    412404  ;; PERMUTATIONS
     
    414406  ;; Rosner 88
    415407  (define (permutations n k)
    416     (fold * 1 (iota k n -1)))
    417 
    418   (define (position item items)
    419     (let ((index (list-index (lambda (x) (equal? item x))
    420                              items)))
    421       (if (number? index)
    422         index
    423         (error "Position: item not in list"))))
     408    (list.fold * 1 (list.iota k n -1)))
    424409
    425410  (define (cumsum sequence)
    426     (if (null? sequence)
     411    (if (empty? sequence)
    427412        0
    428         (reverse (fold (lambda (x ax) (cons (+ x (car ax)) ax)) (list (car sequence)) (cdr sequence)))))
     413        (let ((g (gen-elts sequence)))
     414          (reverse (g-reduce (lambda (x ax) (cons (+ x (car ax)) ax))
     415                             (list (g)) g))
     416          )
     417        ))
    429418
    430419  (define (sign x)
     
    441430    x)
    442431
    443   (define (vector-sort comparison items)
    444     (list->vector (list-sort comparison (vector->list items))))
    445 
    446432  ;; RANDOM-NORMAL
    447433  ;; returns a random number with mean and standard-distribution as specified.
     
    449435    (+ mean (* sd (/ (random-uniform-int 1000000) 1000000))))
    450436
    451   ;; RANDOM-PICK
    452   ;; Random selection from list
    453   (define (random-pick items)
    454     (if (and (list? items) (not (null? items)))
    455       (list-ref items (random-uniform-int (length items)))
    456       #f))
    457437
    458438  ;; RANDOM-SAMPLE
     
    460440  ;; If N is equal to or greater than the length of the sequence, return
    461441  ;; the entire sequence.
    462   (define (random-sample n items)
    463     (cond ((<= n 0)
    464            '())
    465           ((>= n (length items))
    466            items)
    467           (else
    468             (let loop ((remaining items)
    469                        (kept '()))
    470               (if (= (length kept) n)
    471                 kept
    472                 (let ((one (random-pick remaining)))
    473                   (loop (delete one remaining)
    474                         (cons one kept))))))))
    475 
     442  (define (random-sample n sequence)
     443    (let ((result (make-vector n)))
     444      (cond ((<= n 0)
     445             (list->vector '()))
     446            (else
     447             (let loop ((remaining (gen-elts sequence))
     448                        (l 0))
     449               (let ((item (remaining)))
     450                 (if (not item)
     451                     result
     452                     (begin
     453                       (if (< l n)
     454                           (vector-set! result l item)
     455                           (let ((i (random-uniform-int n)))
     456                             (vector-set! result i item)))
     457                       (loop remaining (+ 1 l))
     458                       ))
     459                 ))
     460             ))
     461      ))
     462
     463 
    476464  ;; RANDOM-WEIGHTED-SAMPLE Return a random sample of size M from
    477465  ;; sequence of length N, without replacement, where each element has
    478466  ;; a defined probability of selection (weight) W.  If M is equal to
    479467  ;; or greater to N, return the entire sequence.
    480   (define (random-weighted-sample m items weights)
    481     (let ((n (length items)))
     468  (define (random-weighted-sample m sequence weights)
     469    (let ((n (size sequence)))
    482470      (cond ((<= m 0) '())
    483             ((>= m n) items)
     471            ((>= m n) sequence)
    484472            (else
    485              (let* ((keys (map (lambda (w) (expt (random-uniform) (/ 1 w))) weights))
    486                     (sorted-items (list-sort (lambda (x y) (> (car x) (car y))) (zip keys items))))
    487                (map cadr (take sorted-items m)))
     473             (let* ((keys (map-elts (lambda (w) (expt (random-uniform) (/ 1 w))) weights))
     474                    (sorted-items (sort (lambda (x y) (> (car x) (car y)))
     475                                        (g-map list (gen-elts keys) (gen-elts sequence)))))
     476               (elt-take sorted-items m))
    488477             ))
    489478      ))
     
    496485  ;; Rosner 10
    497486  (define (mean sequence)
    498     (if (null? sequence)
     487    (if (empty? sequence)
    499488      0
    500       (/ (fold + 0 sequence) (length sequence))))
     489      (/ (reduce + 0 sequence) (size sequence))))
    501490
    502491  ;; MEDIAN
     
    509498  ;; Returns two values: a list of the modes and the number of times they occur
    510499  (define (mode sequence)
    511     (if (null? sequence)
    512       (error "Mode: Sequence must not be null")
     500    (if (empty? sequence)
     501      (error "Mode: Sequence must not be empty")
    513502      (let ((count-table (make-hash-table eqv?))
    514             (modes '())
    515             (mode-count 0))
    516         (for-each (lambda (item)
    517                     (hash-table-set! count-table
    518                                      item
    519                                      (+ 1 (hash-table-ref count-table item (lambda () 0)))))
    520                   sequence)
    521         (for-each (lambda (key)
    522                     (let ((val (hash-table-ref count-table key (lambda () #f))))
    523                       (cond ((> val mode-count) ; keep mode
    524                              (set! modes (list key))
    525                              (set! mode-count val))
    526                             ((= val mode-count) ; store multiple modes
    527                              (set! modes (cons key modes))))))
    528                   (hash-table-keys count-table))
    529         (cond ((every number? modes) (set! modes (list-sort < modes)))
    530               ((every string? modes) (set! modes (list-sort string< modes)))
     503            (modes (make-parameter '()))
     504            (mode-count (make-parameter 0)))
     505        (for-each-elt
     506         (lambda (item)
     507           (hash-table-set! count-table
     508                            item
     509                            (+ 1 (hash-table-ref/default count-table item  0))))
     510         sequence)
     511        (for-each
     512         (lambda (key)
     513           (let ((val (hash-table-ref/default count-table key #f)))
     514             (cond ((> val (mode-count)) ; keep mode
     515                    (modes (list key))
     516                    (mode-count val))
     517                   ((= val (mode-count)) ; store multiple modes
     518                    (modes (cons key (modes)))))))
     519         (hash-table-keys count-table))
     520        (cond ((list.every number? (modes)) (modes (list.sort (modes) <)))
     521              ((list.every string? (modes)) (modes (list.sort (modes) string< )))
    531522              )
    532         (values modes mode-count))))
     523        (values (modes) (mode-count)))))
    533524
    534525  ;; GEOMETRIC-MEAN
     
    536527  (define (geometric-mean sequence . args)
    537528    (let ((base (if (null? args) 10 (car args))))
    538       (expt base (mean (map (lambda (x) (/ (log x)
    539                                            (log base)))
    540                             sequence)))))
     529      (expt base (mean (map-elts (lambda (x) (/ (log x)
     530                                                (log base)))
     531                                 sequence)))))
    541532
    542533  ;; RANGE
    543534  ;; Rosner 18
    544535  (define (range sequence)
    545     (if (null? sequence)
     536    (if (empty? sequence)
    546537      0
    547       (- (fold max (car sequence) (cdr sequence))
    548          (fold min (car sequence) (cdr sequence)))))
     538      (- (reduce* max sequence)
     539         (reduce* min sequence))))
    549540
    550541  ;; PERCENTILE
    551542  ;; Rosner 19
    552543  (define (percentile sequence percent)
    553     (if (null? sequence)
     544    (if (empty? sequence)
    554545      (error "Percentile: Sequence must not be null")
    555       (let* ((sorted-vec (vector-sort < (apply vector sequence)))
    556              (n (vector-length sorted-vec))
     546      (let* ((sorted-items (sort (lambda (i vi j vj) (< vi vj)) sequence))
     547             (n (size sorted-items))
    557548             (k (* n (/ percent 100)))
    558549             (floor-k (floor k)))
    559550        (if (= k floor-k)
    560           (/ (+ (vector-ref sorted-vec k)
    561                 (vector-ref sorted-vec (- k 1)))
     551          (/ (+ (elt-ref sorted-items k)
     552                (elt-ref sorted-items (- k 1)))
    562553             2)
    563           (vector-ref sorted-vec floor-k)))))
     554          (elt-ref sorted-items floor-k)))))
    564555
    565556  ;; VARIANCE
    566557  ;; Rosner 21
    567558  (define (variance sequence)
    568     (if (< (length sequence) 2)
     559    (if (< (size sequence) 2)
    569560      (error "variance: sequence must contain at least two elements")
    570561      (let ((mean1 (mean sequence)))
    571         (/ (fold + 0 (map (lambda (x) (square (- mean1 x))) sequence))
    572            (- (length sequence) 1)))))
     562        (/ (reduce + 0
     563                   (map (lambda (x) (square (- mean1 x))) sequence))
     564           (- (size sequence) 1)))))
    573565
    574566  ;; STANDARD-DEVIATION
     
    588580  (define (standard-error-of-the-mean sequence)
    589581    (/ (standard-deviation sequence)
    590        (sqrt (length sequence))))
     582       (sqrt (size sequence))))
    591583
    592584  ;; MEAN-SD-N
     
    596588    (values (mean sequence)
    597589            (standard-deviation sequence)
    598             (length sequence)))
     590            (size sequence)))
    599591
    600592  ;; ---------------------------------------------------------------------------
     
    814806  (define (normal-variance-ci-on-sequence sequence alpha)
    815807    (let ((variance (variance sequence))
    816           (n (length sequence)))
     808          (n (size sequence)))
    817809      (normal-variance-ci variance n alpha)))
    818810
     
    827819  (define (normal-sd-ci-on-sequence sequence alpha)
    828820    (let ((sd (standard-deviation sequence))
    829           (n (length sequence)))
     821          (n (size sequence)))
    830822      (normal-sd-ci sd n alpha)))
    831823
     
    864856          (tails (get-keyword #:tails args (lambda () ':both))))
    865857      (let ((x-bar (mean sequence))
    866             (n (length sequence)))
     858            (n (size sequence)))
    867859        (z-test x-bar n #:mu mu #:sigma sigma #:tails tails))))
    868860
     
    11621154    (let ((exact (get-keyword #:exact? args (lambda () #f)))
    11631155          (tails (get-keyword #:tails args (lambda () ':both))))
    1164       (let* ((differences (map - sequence1 sequence2))
    1165              (plus-count (length (filter positive? differences)))
    1166              (minus-count (length (filter negative? differences))))
     1156      (let* ((differences (map-elts - sequence1 sequence2))
     1157             (plus-count (reduce (lambda (x ax) (if (positive? x) (+ 1 ax) ax)) 0 differences))
     1158             (minus-count (reduce (lambda (x ax) (if (negative? x) (+ 1 ax) ax)) 0 differences)))
    11671159        (sign-test plus-count minus-count :exact? exact :tails tails))))
    11681160
     
    11781170  (define (wilcoxon-signed-rank-test differences . args)
    11791171    (let ((tails (get-keyword #:tails args (lambda () ':both))))
    1180       (let* ((nonzero-differences (filter (lambda (n) (not (zero? n))) differences))
    1181              (sorted-list (list-sort (lambda (x y) (< (car x) (car y)))
    1182                                      (map (lambda (dif)
    1183                                             (list (abs dif)
    1184                                                   (sign dif)))
    1185                                           nonzero-differences)))
    1186              (distinct-values (delete-duplicates (map car sorted-list)))
     1172      (let* ((nonzero-differences (list.filter (lambda (n) (not (zero? n))) differences))
     1173             (sorted-list (list.sort (map (lambda (dif) (list (abs dif) (sign dif)))
     1174                                          nonzero-differences)
     1175                                     (lambda (x y) (< (car x) (car y)))))
     1176             (distinct-values (list.delete-duplicates (map car sorted-list)))
    11871177             (ties '()))
    1188         (when (< (length nonzero-differences) 16)
     1178        (when (< (size nonzero-differences) 16)
    11891179          (error
    11901180            "This Wilcoxon Signed-Rank Test (normal approximation method) requires nonzero N > 15"))
     
    11951185                          (last (position value (reverse (map car sorted-list)))))
    11961186                      (if (= first last)
    1197                         (append (find (lambda (item) (= (car item) value))
    1198                                       sorted-list)
    1199                                 (list (+ 1 first)))
    1200                         (let ((number-tied (+ 1 (- last first)))
    1201                               (avg-rank (+ 1 (/ (+ first last) 2)))) ; + 1 since 0 based
    1202                           (set! ties (cons number-tied ties))
    1203                           (let loop ((i 0)
    1204                                      (result '()))
    1205                             (if (= i number-tied)
    1206                               (reverse result)
    1207                               (loop (+ 1 i)
    1208                                     (cons (cons (list-ref sorted-list (+ first i))
    1209                                                 (list avg-rank))
    1210                                           result))))))))
     1187                          (append (list.find (lambda (item) (= (car item) value))
     1188                                             sorted-list)
     1189                                  (list (+ 1 first)))
     1190                          (let ((number-tied (+ 1 (- last first)))
     1191                                (avg-rank (+ 1 (/ (+ first last) 2)))) ; + 1 since 0 based
     1192                            (set! ties (cons number-tied ties))
     1193                            (let loop ((i 0)
     1194                                       (result '()))
     1195                              (if (= i number-tied)
     1196                                  (reverse result)
     1197                                  (loop (+ 1 i)
     1198                                        (cons (cons (list-ref sorted-list (+ first i))
     1199                                                    (list avg-rank))
     1200                                              result))))))))
    12111201                  distinct-values)
    12121202        (set! ties (reverse ties))
    12131203        (let* ((direction (if (eq? tails ':negative) -1 1))
    1214                (r1 (fold + 0
    1215                          (map (lambda (entry)
    1216                                 (if (= (cadr entry) direction)
    1217                                   (caddr entry)
    1218                                   0))
    1219                               sorted-list)))
     1204               (r1 (list.fold + 0
     1205                              (map (lambda (entry)
     1206                                     (if (= (cadr entry) direction)
     1207                                         (caddr entry)
     1208                                         0))
     1209                                   sorted-list)))
    12201210               (n (length nonzero-differences))
    12211211               (expected-r1 (/ (* n (+ 1 n)) 4))
    12221212               (ties-factor (if ties
    1223                               (/ (fold + 0
    1224                                        (map (lambda (ti) (- (* ti ti ti) ti))
    1225                                             ties))
     1213                              (/ (list.fold + 0
     1214                                            (map (lambda (ti) (- (* ti ti ti) ti))
     1215                                                 ties))
    12261216                                 48)
    12271217                              0))
     
    13111301      (let* ((ns (map + row1-counts row2-counts))
    13121302             (p-hats (map / row1-counts ns))
    1313              (n (fold + 0 ns))
    1314              (p-bar (/ (fold + 0 row1-counts) n))
     1303             (n (list.fold + 0 ns))
     1304             (p-bar (/ (list.fold + 0 row1-counts) n))
    13151305             (q-bar (- 1 p-bar))
    13161306             (s-bar (mean scores))
    1317              (a (fold + 0.0
    1318                       (map (lambda (p-hat ni s)
    1319                              (* ni (- p-hat p-bar) (- s s-bar)))
    1320                            p-hats ns scores)))
    1321              (b (* 1.0 p-bar q-bar (- (fold + 0 (map (lambda (ni s) (* ni (square s)))
    1322                                                  ns scores))
    1323                                   (/ (square (fold + 0 (map (lambda (ni s) (* ni s))
    1324                                                             ns scores)))
    1325                                      n))))
     1307             (a (list.fold + 0.0
     1308                           (map (lambda (p-hat ni s)
     1309                                  (* ni (- p-hat p-bar) (- s s-bar)))
     1310                                p-hats ns scores)))
     1311             (b (* 1.0 p-bar q-bar (- (list.fold + 0 (map (lambda (ni s) (* ni (square s)))
     1312                                                          ns scores))
     1313                                      (/ (square (list.fold + 0 (map (lambda (ni s) (* ni s))
     1314                                                                     ns scores)))
     1315                                         n))))
    13261316             (x2 (/ (square a) b))
    13271317             (significance (- 1 (chi-square-cdf x2 1))))
     
    14791469  ;; the significance of the difference of the slope from 0.
    14801470  (define (linear-regression points)
    1481     (unless (> (length points) 2)
     1471    (unless (> (size points) 2)
    14821472      (error "Requires at least three points"))
    1483     (let ((xs (map car points))
    1484           (ys (map cadr points)))
     1473    (let ((xs (map-elts car points))
     1474          (ys (map-elts cadr points)))
    14851475      (let* ((x-bar (mean xs))
    14861476             (y-bar (mean ys))
    1487              (n (length points))
    1488              (Lxx (fold + 0
    1489                         (map (lambda (xi) (square (- xi x-bar))) xs)))
    1490              (Lyy (fold + 0
    1491                         (map (lambda (yi) (square (- yi y-bar))) ys)))
    1492              (Lxy (fold + 0
    1493                         (map (lambda (xi yi) (* (- xi x-bar) (- yi y-bar)))
    1494                              xs ys)))
     1477             (n (size points))
     1478             (Lxx (reduce + 0.0
     1479                        (map-elts (lambda (xi) (square (- xi x-bar))) xs)))
     1480             (Lyy (reduce + 0.0
     1481                        (map-elts (lambda (yi) (square (- yi y-bar))) ys)))
     1482             (Lxy (reduce + 0.0
     1483                          (map-elts (lambda (point) (let ((xi (car point))
     1484                                                          (yi (cadr point)))
     1485                                                      (* (- xi x-bar) (- yi y-bar))))
     1486                             points)))
    14951487             (b (if (zero? Lxx) 0 (/ Lxy Lxx)))
    14961488             (a (- y-bar (* b x-bar)))
     
    15081500  ;; Also called Pearson Correlation
    15091501  (define (correlation-coefficient points)
    1510     (let* ((xs (map car points))
    1511            (ys (map cadr points))
     1502    (let* ((xs (map-elts car points))
     1503           (ys (map-elts cadr points))
    15121504           (x-bar (mean xs))
    15131505           (y-bar (mean ys)))
    1514       (/ (fold + 0 (map (lambda (xi yi) (* (- xi x-bar)
    1515                                            (- yi y-bar)))
    1516                         xs
    1517                         ys))
    1518          (sqrt (* (fold + 0 (map (lambda (xi) (square (- xi x-bar)))
    1519                                  xs))
    1520                   (fold + 0 (map (lambda (yi) (square (- yi y-bar)))
    1521                                  ys)))))))
     1506      (/ (reduce + 0 (map-elts (lambda (point)
     1507                                (let ((xi (car point))
     1508                                      (yi (cadr point)))
     1509                                  (* (- xi x-bar) (- yi y-bar))))
     1510                              points))
     1511         (sqrt (* (reduce + 0 (map-elts (lambda (xi) (square (- xi x-bar)))
     1512                                       xs))
     1513                  (reduce + 0 (map-elts (lambda (yi) (square (- yi y-bar)))
     1514                                       ys)))))))
    15221515
    15231516
     
    15541547  ;; and its significance.
    15551548  (define (spearman-rank-correlation points)
    1556     (let ((xis (map car points))
    1557           (yis (map cadr points)))
    1558       (let* ((n (length points))
    1559              (sorted-xis (list-sort < xis))
    1560              (sorted-yis (list-sort < yis))
    1561              (average-x-ranks (map (lambda (x) (average-rank x sorted-xis)) xis))
    1562              (average-y-ranks (map (lambda (y) (average-rank y sorted-yis)) yis))
     1549    (let ((xis (map-elts car points))
     1550          (yis (map-elts cadr points)))
     1551      (let* ((n (size points))
     1552             (sorted-xis (sort (lambda (xi x yi y) (< x y)) xis))
     1553             (sorted-yis (sort (lambda (xi x yi y) (< x y)) yis))
     1554             (average-x-ranks (map-elts (lambda (x) (average-rank x sorted-xis)) xis))
     1555             (average-y-ranks (map-elts (lambda (y) (average-rank y sorted-yis)) yis))
    15631556             (mean-x-rank (mean average-x-ranks))
    15641557             (mean-y-rank (mean average-y-ranks))
    1565              (Lxx (fold + 0
    1566                         (map (lambda (xi-rank) (square (- xi-rank mean-x-rank)))
    1567                              average-x-ranks)))
    1568              (Lyy (fold + 0
    1569                         (map (lambda (yi-rank) (square (- yi-rank mean-y-rank)))
    1570                              average-y-ranks)))
    1571              (Lxy (fold + 0
    1572                         (map (lambda (xi-rank yi-rank)
    1573                                (* (- xi-rank mean-x-rank)
    1574                                   (- yi-rank mean-y-rank)))
    1575                              average-x-ranks average-y-ranks)))
     1558             (Lxx (reduce + 0
     1559                        (map-elts (lambda (xi-rank) (square (- xi-rank mean-x-rank)))
     1560                                 average-x-ranks)))
     1561             (Lyy (reduce + 0
     1562                          (map-elts (lambda (yi-rank) (square (- yi-rank mean-y-rank)))
     1563                                   average-y-ranks)))
     1564             (Lxy (reduce + 0
     1565                          (map-elts (lambda (xi-rank yi-rank)
     1566                                     (* (- xi-rank mean-x-rank)
     1567                                        (- yi-rank mean-y-rank)))
     1568                                   average-x-ranks average-y-ranks)))
    15761569             (rs (if (> (* Lxx Lyy) 0) (/ Lxy (sqrt (* Lxx Lyy))) 0)) ; TODO: think about rs = 1
    15771570             (ts (if (= 1 rs) 1 (/ (* rs (sqrt (- n 2))) (sqrt (- 1 (square rs))))))
Note: See TracChangeset for help on using the changeset viewer.