Changeset 37172 in project


Ignore:
Timestamp:
02/01/19 01:28:45 (7 months ago)
Author:
Ivan Raikov
Message:

statistics: refactoring correlation and regression functions to take two dataset arguments instead of a single dataset of list elements

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

Legend:

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

    r36443 r37172  
    561561      (let ((mean1 (mean sequence)))
    562562        (/ (reduce + 0
    563                    (map (lambda (x) (square (- mean1 x))) sequence))
     563                   (map-elts (lambda (x) (square (- mean1 x))) sequence))
    564564           (- (size sequence) 1)))))
    565565
     
    14681468  ;; and reports the intercept, slope, correlation coefficient r, R^2, and
    14691469  ;; the significance of the difference of the slope from 0.
    1470   (define (linear-regression points)
    1471     (unless (> (size points) 2)
     1470  (define (linear-regression xs ys)
     1471    (unless (> (size xs) 2)
    14721472      (error "Requires at least three points"))
    1473     (let ((xs (map-elts car points))
    1474           (ys (map-elts cadr points)))
    1475       (let* ((x-bar (mean xs))
    1476              (y-bar (mean ys))
    1477              (n (size points))
    1478              (Lxx (reduce + 0.0
     1473    (let* ((x-bar (mean xs))
     1474           (y-bar (mean ys))
     1475           (n (size xs))
     1476           (Lxx (reduce + 0.0
    14791477                        (map-elts (lambda (xi) (square (- xi x-bar))) xs)))
    1480              (Lyy (reduce + 0.0
     1478           (Lyy (reduce + 0.0
    14811479                        (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)))
    1487              (b (if (zero? Lxx) 0 (/ Lxy Lxx)))
    1488              (a (- y-bar (* b x-bar)))
    1489              (reg-ss (* b Lxy))
    1490              (res-ms (/ (- Lyy reg-ss) (- n 2)))
    1491              (r (/ Lxy (sqrt (* Lxx Lyy))))
    1492              (r2 (/ reg-ss Lyy))
    1493              (t-test (/ b (sqrt (/ res-ms Lxx))))
    1494              (t-sign (t-significance t-test (- n 2) #:tails ':both)))
    1495         (format #t "~%Intercept = ~A, slope = ~A, r = ~A, R^2 = ~A, p = ~A~%"
    1496                 a b r r2 t-sign)
    1497         (values a b r r2 t-sign))))
     1480           (Lxy (reduce + 0.0
     1481                        (map-elts (lambda (xi yi) (* (- xi x-bar) (- yi y-bar)))
     1482                                  xs ys)))
     1483           (b (if (zero? Lxx) 0 (/ Lxy Lxx)))
     1484           (a (- y-bar (* b x-bar)))
     1485           (reg-ss (* b Lxy))
     1486           (res-ms (/ (- Lyy reg-ss) (- n 2)))
     1487           (r (/ Lxy (sqrt (* Lxx Lyy))))
     1488           (r2 (/ reg-ss Lyy))
     1489           (t-test (/ b (sqrt (/ res-ms Lxx))))
     1490           (t-sign (t-significance t-test (- n 2) #:tails ':both)))
     1491      (format #t "~%Intercept = ~A, slope = ~A, r = ~A, R^2 = ~A, p = ~A~%"
     1492              a b r r2 t-sign)
     1493      (values a b r r2 t-sign)))
    14981494
    14991495  ;; CORRELATION-COEFFICIENT
    15001496  ;; Also called Pearson Correlation
    1501   (define (correlation-coefficient points)
    1502     (let* ((xs (map-elts car points))
    1503            (ys (map-elts cadr points))
    1504            (x-bar (mean xs))
     1497  (define (correlation-coefficient xs ys)
     1498    (let* ((x-bar (mean xs))
    15051499           (y-bar (mean 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))
     1500      (/ (reduce + 0 (map-elts (lambda (xi yi)
     1501                                 (* (- xi x-bar) (- yi y-bar)))
     1502                              xs ys))
    15111503         (sqrt (* (reduce + 0 (map-elts (lambda (xi) (square (- xi x-bar)))
    15121504                                       xs))
     
    15461538  ;; linear-regression) and returns the spearman rank correlation coefficient
    15471539  ;; and its significance.
    1548   (define (spearman-rank-correlation points)
    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))
    1556              (mean-x-rank (mean average-x-ranks))
    1557              (mean-y-rank (mean average-y-ranks))
    1558              (Lxx (reduce + 0
     1540  (define (spearman-rank-correlation xis yis)
     1541    (let* ((n (size xis))
     1542           (sorted-xis (sort (lambda (xi x yi y) (< x y)) xis))
     1543           (sorted-yis (sort (lambda (xi x yi y) (< x y))  yis))
     1544           (average-x-ranks (map-elts (lambda (x) (average-rank x sorted-xis)) xis))
     1545           (average-y-ranks (map-elts (lambda (y) (average-rank y sorted-yis)) yis))
     1546           (mean-x-rank (mean average-x-ranks))
     1547           (mean-y-rank (mean average-y-ranks))
     1548           (Lxx (reduce + 0
    15591549                        (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)))
    1569              (rs (if (> (* Lxx Lyy) 0) (/ Lxy (sqrt (* Lxx Lyy))) 0)) ; TODO: think about rs = 1
    1570              (ts (if (= 1 rs) 1 (/ (* rs (sqrt (- n 2))) (sqrt (- 1 (square rs))))))
    1571              (p (t-significance ts (- n 2) #:tails ':both)))
    1572         (format #t "~%Spearman correlation coefficient ~A, p = ~A~%" rs p)
    1573         (values rs p))))
     1550                                  average-x-ranks)))
     1551           (Lyy (reduce + 0
     1552                        (map-elts (lambda (yi-rank) (square (- yi-rank mean-y-rank)))
     1553                                  average-y-ranks)))
     1554           (Lxy (reduce + 0
     1555                        (map-elts (lambda (xi-rank yi-rank)
     1556                                    (* (- xi-rank mean-x-rank)
     1557                                       (- yi-rank mean-y-rank)))
     1558                                  average-x-ranks average-y-ranks)))
     1559           (rs (if (> (* Lxx Lyy) 0) (/ Lxy (sqrt (* Lxx Lyy))) 0)) ; TODO: think about rs = 1
     1560           (ts (if (= 1 rs) 1 (/ (* rs (sqrt (- n 2))) (sqrt (- 1 (square rs))))))
     1561           (p (t-significance ts (- n 2) #:tails ':both)))
     1562      (format #t "~%Spearman correlation coefficient ~A, p = ~A~%" rs p)
     1563      (values rs p)))
    15741564
    15751565  ;; ---------------------------------------------------------------------------
  • release/5/statistics/trunk/tests/run.scm

    r36303 r37172  
    279279
    280280
    281 (let-values (((s p) (spearman-rank-correlation '((4 5) (10 8) (3 6) (1 2) (9 10) (2 3) (6 9) (7 4) (8 7) (5 1)))))
     281(let-values (((s p) (spearman-rank-correlation '(4 10 3 1 9 2 6 7 8 5)
     282                                               '(5 8 6 2 10 3 9 4 7 1))))
    282283  (test-assert "spearman-rank-correlation" (=5 s (/ 113 165)))
    283284  (test-assert "spearman-rank-correlation" (=5 p 0.0288827975067328)))
     
    296297(test-group "correlation and regression"
    297298
    298             (let-values (((a b r r2 p) (linear-regression '((1.0 0.1) (2.0 0.3) (3.0 0.8)))))
     299            (let-values (((a b r r2 p) (linear-regression '(1.0 2.0 3.0) '(0.1 0.3 0.8))))
    299300              (test "linear-regression"  -0.3 a)
    300301              (test-assert "linear-regression" (=5 b 0.35))
     
    305306
    306307            (test-assert "correlation-coefficient"
    307                          (=5 0.970725343394151 (correlation-coefficient '((1.0 0.1) (2.0 0.3) (3.0 0.8)))))
     308                         (=5 0.970725343394151 (correlation-coefficient  '(1.0 2.0 3.0) '(0.1 0.3 0.8))))
    308309            )
    309310
Note: See TracChangeset for help on using the changeset viewer.