home *** CD-ROM | disk | FTP | other *** search
/ Garbo / Garbo.cdr / mac / science / xlspstrd.sit / glim / bradleyterry.lsp next >
Lisp/Scheme  |  1991-01-08  |  5KB  |  143 lines
  1. ;;;;
  2. ;;;;                  Fitting a Bradley-Terry Model with the
  3. ;;;;                          Lisp-Stat GLIM System
  4.  
  5. (require "glim")
  6.  
  7. ;;;
  8. ;;; Function for fitting the model as a binomial regression model
  9. ;;;
  10.  
  11. (defun bradley-terry-model (counts &key labels)
  12.   (let* ((n (round (sqrt (length counts))))
  13.          (i (iseq 1 (- n 1)))
  14.          (low-i (apply #'append (+ (* n i) (mapcar #'iseq i))))
  15.          (p-names (if labels
  16.                       (rest labels) 
  17.                       (level-names (iseq n) :prefix "Choice"))))
  18.     (labels ((tr (x)
  19.                (apply #'append (transpose (split-list (coerce x 'list) n))))
  20.              (lower (x) (select x low-i))
  21.              (low-indicators (x) (mapcar #'lower (indicators x))))
  22.       (let ((wins (lower counts))
  23.             (losses (lower (tr counts)))
  24.             (rows (low-indicators (repeat (iseq n) (repeat n n))))
  25.             (cols (low-indicators (repeat (iseq n) n))))
  26.         (binomialreg-model (- rows cols)
  27.                            wins 
  28.                            (+ wins losses)
  29.                            :intercept nil
  30.                            :predictor-names p-names)))))
  31.  
  32. ;;;;
  33. ;;;;
  34. ;;;; Bradley-Terry Model Prototype
  35. ;;;;
  36. ;;;;
  37.  
  38. (defproto bradley-terry-proto () () binomialreg-proto)
  39.  
  40. ;;;
  41. ;;; Constructor Function
  42. ;;;
  43.  
  44. (defun bradley-terry-model (counts &key labels)
  45.   (let* ((n (round (sqrt (length counts))))
  46.          (i (iseq 1 (- n 1)))
  47.          (low-i (apply #'append (+ (* n i) (mapcar #'iseq i))))
  48.          (p-names (if labels
  49.                       (rest labels) 
  50.                       (level-names (iseq n) :prefix "Choice"))))
  51.     (labels ((tr (x)
  52.                (apply #'append (transpose (split-list (coerce x 'list) n))))
  53.              (lower (x) (select x low-i))
  54.              (low-indicators (x) (mapcar #'lower (indicators x))))
  55.       (let ((wins (lower counts))
  56.             (losses (lower (tr counts)))
  57.             (rows (low-indicators (repeat (iseq n) (repeat n n))))
  58.             (cols (low-indicators (repeat (iseq n) n))))
  59.         (send bradley-terry-proto :new
  60.               :x (- rows cols)
  61.               :y wins
  62.               :trials (+ wins losses)
  63.               :intercept nil
  64.               :predictor-names p-names)))))
  65.  
  66. ;;;;
  67. ;;;; Methods for Estimating Success Probabilities
  68. ;;;;
  69.  
  70. ;;;
  71. ;;; Simplest version of :SUCCESS-PROB
  72. ;;;
  73. (defmeth bradley-terry-proto :success-prob (i j)
  74.   (let* ((phi (cons 0 (send self :coef-estimates)))
  75.          (exp-logit (exp (- (select phi i) (select phi j)))))
  76.     (/ exp-logit (+ 1 exp-logit))))
  77.  
  78. ;;;
  79. ;;; Version with standard errors based on :SUCCESS-LOGIT
  80. ;;;
  81. (defmeth bradley-terry-proto :success-prob (i j &optional stdev)
  82.   (let* ((success-logit (send self :success-logit i j stdev))
  83.          (exp-logit (exp (if stdev (first success-logit) success-logit)))
  84.          (p (/ exp-logit (+ 1 exp-logit)))
  85.          (s (if stdev (* p (- 1 p) (second success-logit)))))
  86.     (if stdev (list p s) p)))
  87.  
  88. ;;;
  89. ;;; Non-vectorized version of :SUCCESS-LOGIT
  90. ;;;
  91. (defmeth bradley-terry-proto :success-logit (i j &optional stdev)
  92.   (let ((coefs (send self :coef-estimates)))
  93.     (flet ((lincomb (i j)
  94.              (let ((v (repeat 0 (length coefs))))
  95.                (if (/= 0 i) (setf (select v (- i 1)) 1))
  96.                (if (/= 0 j) (setf (select v (- j 1)) -1))
  97.                v)))
  98.       (let* ((v (lincomb i j))
  99.              (logit (inner-product v coefs))
  100.              (var (if stdev (matmult v (send self :xtxinv) v))))
  101.         (if stdev (list logit (sqrt var)) logit)))))
  102.  
  103. ;;;
  104. ;;; Vectorized version of :SUCCESS-LOGIT
  105. ;;;
  106. (defmeth bradley-terry-proto :success-logit (i j &optional stdev)
  107.   (let ((coefs (send self :coef-estimates)))
  108.     (flet ((lincomb (i j)
  109.              (let ((v (repeat 0 (length coefs))))
  110.                (if (/= 0 i) (setf (select v (- i 1)) 1))
  111.                (if (/= 0 j) (setf (select v (- j 1)) -1))
  112.                v))
  113.            (matrix-or-list (x)
  114.              (let ((x (coerce x 'list)))
  115.                (if (numberp (first x)) x (apply #'bind-rows x))))
  116.            (apply-sum (x)
  117.              (if (matrixp x) (mapcar #'sum (row-list x)) (apply #'sum x))))
  118.       (let* ((v (matrix-or-list (map-elements #'lincomb i j)))
  119.              (xtxinv (if stdev (send self :xtxinv)))
  120.              (logit (matmult v coefs))
  121.              (var (if stdev (apply-sum (* (matmult v xtxinv) v)))))
  122.         (if stdev (list logit (sqrt var)) logit)))))
  123.  
  124. ;;;
  125. ;;; Baseball Data Example from Agresti
  126. ;;;
  127.  
  128. #|
  129. (def wins-losses '( -  7  9  7  7  9 11
  130.                     6  -  7  5 11  9  9
  131.                     4  6  -  7  7  8 12
  132.                     6  8  6  -  6  7 10
  133.                     6  2  6  7  -  7 12
  134.                     4  4  5  6  6  -  6
  135.                     2  4  1  3  1  7  -))
  136.  
  137. (def teams '("Milwaukee" "Detroit" "Toronto" "New York"
  138.              "Boston" "Cleveland" "Baltimore"))
  139.  
  140. (def wl (bradley-terry-model wins-losses :labels teams))
  141. (send wl :success-prob 4 3)
  142. |#
  143.