markov chain calculator

markov.rkt

(require racket)
(require math/matrix)

(define (array-sum an-array)
  (for/sum ([x an-array]) x))



;; absorbing markov chain
;; transition matrix P
; canonical form
; Q R
; 0 I

(define Q
  (matrix [[0 0 0 0 0 0 0]
           [0 0 0 0.5 0 0 0]
           [0 0.5 0 0.5 0 0 0]
           [0 0 0 0 0 0 0]
           [0 0 0 0 0 0 0]
           [0 0 0 0 0 0 0]
           [0 0 0 0 0 0 0]]))
(define R
  (matrix [[0 0]
           [0 0.5]
           [0 0]
           [0.5 0.5]
           [0 0]
           [0 0]
           [0 0]
          ]))

(define (real->decimal matrix)
  (matrix-map string->number
              (matrix-map real->decimal-string matrix)))
(define (cal-N matQ) ; N = fundamental matrix
  (real->decimal
   (matrix-inverse
    (matrix-
     (identity-matrix (matrix-num-rows matQ))
     matQ))))
(define (cal-Nc matN) ; t = steps in each state
  (matrix* matN
           (make-matrix (matrix-num-rows matN) 1 1)))
(define (cal-NR matN matR) ; from each state, proba to end up in absorbed
  (real->decimal
   (matrix* matN
           matR)))

;; ergodic
; from every state can go to every state
;; regular
; not too many zeros
; some power has no zeros (all positive)
; one max eigen value = 1
; w: fixed vector after power n
; equilibrium state (time in each state or independent trial)

(define P (matrix
           [[1/2 1/4 1/4]
            [1/2 0 1/2]
            [1/4 1/4 1/2]]))

(define (equations matP)
  (let ([row (matrix-num-rows matP)])
    (list->matrix row row
                  (append
                   (make-list row 0)
                   (for*/list ([j (sub1 row)]
                               [i row])
                     (matrix-ref matP i j))))))

(define P*
  (matrix [[1 0 0]
           [0 0 0]
           [0 -1 0]]))

(define augmented
  (matrix [[1][1][0]]))

(define (cal-w* matP) ; stationary distribution
  (matrix->list
   (matrix-solve (matrix+ (equations matP) P*)
                augmented)))

(define (cal-w matP)
  (let* ([l (cal-w* matP)]
        [s (sum l)])
    (for/list ([i (length l)])
      (/ (list-ref l i) s))))

;; fundamental matrix
(define (cal-Z matP matW)
  (matrix-inverse
   (matrix+
    (matrix- (identity-matrix (matrix-num-rows matP)) matP)
    matW)))

(define T (matrix [[.5 .5 0]
                   [.25 .5 .25]
                   [0 1/3 2/3]]))

(define (cal-mij i j matZ matW)
  (/ (- (matrix-ref matZ j j) (matrix-ref matZ i j))
     (matrix-ref matW 0 j)))

(define (cal-M matZ matW)
  (let ([row (matrix-num-rows matZ)])
    (list->matrix row row
                  (for*/list ([i row]
                              [j row])
                    (cal-mij i j matZ matW)))))

(define (cal-R matW)
  (real->decimal
   (matrix-map (lambda (x) (/ 1 x))
               matW)))