Changeset 37172 in project
 Timestamp:
 02/01/19 01:28:45 (2 weeks ago)
 Location:
 release/5/statistics/trunk
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

release/5/statistics/trunk/statistics.scm
r36443 r37172 561 561 (let ((mean1 (mean sequence))) 562 562 (/ (reduce + 0 563 (map (lambda (x) (square ( mean1 x))) sequence))563 (mapelts (lambda (x) (square ( mean1 x))) sequence)) 564 564 ( (size sequence) 1))))) 565 565 … … 1468 1468 ;; and reports the intercept, slope, correlation coefficient r, R^2, and 1469 1469 ;; the significance of the difference of the slope from 0. 1470 (define (linearregression points)1471 (unless (> (size points) 2)1470 (define (linearregression xs ys) 1471 (unless (> (size xs) 2) 1472 1472 (error "Requires at least three points")) 1473 (let ((xs (mapelts car points)) 1474 (ys (mapelts cadr points))) 1475 (let* ((xbar (mean xs)) 1476 (ybar (mean ys)) 1477 (n (size points)) 1478 (Lxx (reduce + 0.0 1473 (let* ((xbar (mean xs)) 1474 (ybar (mean ys)) 1475 (n (size xs)) 1476 (Lxx (reduce + 0.0 1479 1477 (mapelts (lambda (xi) (square ( xi xbar))) xs))) 1480 1478 (Lyy (reduce + 0.0 1481 1479 (mapelts (lambda (yi) (square ( yi ybar))) ys))) 1482 (Lxy (reduce + 0.0 1483 (mapelts (lambda (point) (let ((xi (car point)) 1484 (yi (cadr point))) 1485 (* ( xi xbar) ( yi ybar)))) 1486 points))) 1487 (b (if (zero? Lxx) 0 (/ Lxy Lxx))) 1488 (a ( ybar (* b xbar))) 1489 (regss (* b Lxy)) 1490 (resms (/ ( Lyy regss) ( n 2))) 1491 (r (/ Lxy (sqrt (* Lxx Lyy)))) 1492 (r2 (/ regss Lyy)) 1493 (ttest (/ b (sqrt (/ resms Lxx)))) 1494 (tsign (tsignificance ttest ( n 2) #:tails ':both))) 1495 (format #t "~%Intercept = ~A, slope = ~A, r = ~A, R^2 = ~A, p = ~A~%" 1496 a b r r2 tsign) 1497 (values a b r r2 tsign)))) 1480 (Lxy (reduce + 0.0 1481 (mapelts (lambda (xi yi) (* ( xi xbar) ( yi ybar))) 1482 xs ys))) 1483 (b (if (zero? Lxx) 0 (/ Lxy Lxx))) 1484 (a ( ybar (* b xbar))) 1485 (regss (* b Lxy)) 1486 (resms (/ ( Lyy regss) ( n 2))) 1487 (r (/ Lxy (sqrt (* Lxx Lyy)))) 1488 (r2 (/ regss Lyy)) 1489 (ttest (/ b (sqrt (/ resms Lxx)))) 1490 (tsign (tsignificance ttest ( n 2) #:tails ':both))) 1491 (format #t "~%Intercept = ~A, slope = ~A, r = ~A, R^2 = ~A, p = ~A~%" 1492 a b r r2 tsign) 1493 (values a b r r2 tsign))) 1498 1494 1499 1495 ;; CORRELATIONCOEFFICIENT 1500 1496 ;; Also called Pearson Correlation 1501 (define (correlationcoefficient points) 1502 (let* ((xs (mapelts car points)) 1503 (ys (mapelts cadr points)) 1504 (xbar (mean xs)) 1497 (define (correlationcoefficient xs ys) 1498 (let* ((xbar (mean xs)) 1505 1499 (ybar (mean ys))) 1506 (/ (reduce + 0 (mapelts (lambda (point) 1507 (let ((xi (car point)) 1508 (yi (cadr point))) 1509 (* ( xi xbar) ( yi ybar)))) 1510 points)) 1500 (/ (reduce + 0 (mapelts (lambda (xi yi) 1501 (* ( xi xbar) ( yi ybar))) 1502 xs ys)) 1511 1503 (sqrt (* (reduce + 0 (mapelts (lambda (xi) (square ( xi xbar))) 1512 1504 xs)) … … 1546 1538 ;; linearregression) and returns the spearman rank correlation coefficient 1547 1539 ;; and its significance. 1548 (define (spearmanrankcorrelation points) 1549 (let ((xis (mapelts car points)) 1550 (yis (mapelts cadr points))) 1551 (let* ((n (size points)) 1552 (sortedxis (sort (lambda (xi x yi y) (< x y)) xis)) 1553 (sortedyis (sort (lambda (xi x yi y) (< x y)) yis)) 1554 (averagexranks (mapelts (lambda (x) (averagerank x sortedxis)) xis)) 1555 (averageyranks (mapelts (lambda (y) (averagerank y sortedyis)) yis)) 1556 (meanxrank (mean averagexranks)) 1557 (meanyrank (mean averageyranks)) 1558 (Lxx (reduce + 0 1540 (define (spearmanrankcorrelation xis yis) 1541 (let* ((n (size xis)) 1542 (sortedxis (sort (lambda (xi x yi y) (< x y)) xis)) 1543 (sortedyis (sort (lambda (xi x yi y) (< x y)) yis)) 1544 (averagexranks (mapelts (lambda (x) (averagerank x sortedxis)) xis)) 1545 (averageyranks (mapelts (lambda (y) (averagerank y sortedyis)) yis)) 1546 (meanxrank (mean averagexranks)) 1547 (meanyrank (mean averageyranks)) 1548 (Lxx (reduce + 0 1559 1549 (mapelts (lambda (xirank) (square ( xirank meanxrank))) 1560 averagexranks)))1561 1562 1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 (values rs p))))1550 averagexranks))) 1551 (Lyy (reduce + 0 1552 (mapelts (lambda (yirank) (square ( yirank meanyrank))) 1553 averageyranks))) 1554 (Lxy (reduce + 0 1555 (mapelts (lambda (xirank yirank) 1556 (* ( xirank meanxrank) 1557 ( yirank meanyrank))) 1558 averagexranks averageyranks))) 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 (tsignificance ts ( n 2) #:tails ':both))) 1562 (format #t "~%Spearman correlation coefficient ~A, p = ~A~%" rs p) 1563 (values rs p))) 1574 1564 1575 1565 ;;  
release/5/statistics/trunk/tests/run.scm
r36303 r37172 279 279 280 280 281 (letvalues (((s p) (spearmanrankcorrelation '((4 5) (10 8) (3 6) (1 2) (9 10) (2 3) (6 9) (7 4) (8 7) (5 1))))) 281 (letvalues (((s p) (spearmanrankcorrelation '(4 10 3 1 9 2 6 7 8 5) 282 '(5 8 6 2 10 3 9 4 7 1)))) 282 283 (testassert "spearmanrankcorrelation" (=5 s (/ 113 165))) 283 284 (testassert "spearmanrankcorrelation" (=5 p 0.0288827975067328))) … … 296 297 (testgroup "correlation and regression" 297 298 298 (letvalues (((a b r r2 p) (linearregression '( (1.0 0.1) (2.0 0.3) (3.0 0.8)))))299 (letvalues (((a b r r2 p) (linearregression '(1.0 2.0 3.0) '(0.1 0.3 0.8)))) 299 300 (test "linearregression" 0.3 a) 300 301 (testassert "linearregression" (=5 b 0.35)) … … 305 306 306 307 (testassert "correlationcoefficient" 307 (=5 0.970725343394151 (correlationcoefficient '((1.0 0.1) (2.0 0.3) (3.0 0.8)))))308 (=5 0.970725343394151 (correlationcoefficient '(1.0 2.0 3.0) '(0.1 0.3 0.8)))) 308 309 ) 309 310
Note: See TracChangeset
for help on using the changeset viewer.