(defn logsumexp [& log-vals]
  (let [mx (apply max log-vals)]
    (+ mx
       (Math/log2
	     (apply +
	       (map (fn [z] (Math/pow 2 z))
		    (map (fn [x] (- x mx))
			log-vals)))))))
			
(defn flip [p]
  (if (< (rand 1) p)
    true
    false))

(defn sample-categorical [outcomes params]
  (if (flip (first params))
    (first outcomes)
    (sample-categorical (rest outcomes)
                        (normalize (rest params)))))

(defn score-categorical [outcome outcomes params]
  (if (empty? params)
    (throw "no matching outcome")
    (if (= outcome (first outcomes))
     (Math/log2 (first params))
      (score-categorical outcome (rest outcomes) (rest params)))))

(defn normalize [params]
  (let [sum (apply + params)]
  (map (fn [x] (/ x sum)) params)))

(defn sample-gamma [shape scale]
(apply + (repeatedly
  shape  (fn []
		     (- (Math/log2 (rand))))
		)))

(defn sample-dirichlet [pseudos]
  (let [gammas (map (fn [sh]
                     (sample-gamma sh 1))
                pseudos)]
				(normalize gammas)))

(defn update-context [order old-context new-symbol]
  (if (>= (count old-context) order)
    (throw "Context too long!")
    (if (= (count old-context) (- order 1))
      (concat  (rest old-context) (list new-symbol))
      (concat old-context (list new-symbol)))))

(defn hmm-unfold [transition observation order context current stop?]
  (if (stop? current)
    (list current)
    (let [new-context (update-context
	           order
		       context
		       current)
	        nxt (transition new-context)]
	  (cons [current (observation current)]
	  (hmm-unfold
	         transition
	         observation
	         n-gram-order 
			 new-context
			 nxt
			 stop?)))))

(def categories '(N V Adj Adv P stop))

(def vocabulary '(Call me Ishmael))

(def category->observation-model
     (memoize (fn [category]
		(sample-dirichlet
		(repeatedly
		 (count vocabulary)
		  (fn [] 1))))))

(def category->transition-model
     (memoize (fn [category]
		(sample-dirichlet
		(repeatedly
		 (count categories)
		  (fn [] 1))))))

(defn sample-category [preceding]
  (sample-categorical
   categories
   (category->transition-model  preceding)))

(defn sample-observation [category]
     (sample-categorical
      vocabulary
      (category->observation-model category)))

(defn sample-categories-words [prec]
  (let [cat (sample-category prec)]
    (if (= cat 'stop)
      '()
      (cons
       [cat (sample-observation cat)]
       (sample-categories-words cat)))))

(defn forward-probability [category sentence]
  (if (empty? sentence)
      0
      (+ 
       (apply logsumexp
	      (map (fn [cat]  
		    (+ (forward-probability cat (all-but-last sentence))
		       (score-categorical
			cat
			categories
			(category->transition-model cat))))
                  categories))
       (if (= category 'stop)
           0
           (score-categorical
	    (last sentence)
	    vocabulary
	    (category->observation-model category))))))
		


  

In the last unit, we saw how to use the conditional independence assumptions of HMMs to efficiently compute the marginal score of a string based on the decomposition of forward probabilities. Recall that a the forward probability is defined as the probability of generating the first \(i\) words of a string and stopping in some particular state \(c^{(i+1)}\), that is,

\[\mathbf{fw}_i(c^{(i+1)}) = \Pr(w^{(1)},\dots,w^{(i)},c^{(i+1)})\]

This marginal probability can be computed recursively by traversing the trellis from left to right, computing the forward values at each step, in turn. If the length of our string is \(|\vec{w}|\), then

\[\mathbf{fw}_{|\vec{w}|}(\ltimes) = \Pr(w^{(1)},\dots,w^{(|\vec{w}|)},C^{(|\vec{w}|+1)}=\ltimes)\]

is the marginal probability of our whole string under the HMM.

This, however, is not the only way to efficiently compute the marginal probability of a string under an HMM. Consider the following quantity:

\[\mathbf{bk}_{k}(c^{(k-1)}) = \Pr(w^{(k)}, \dots, w^{(|\vec{w}|)} \mid c^{(k-1)})\]

This quantity can be thought of as the “reverse” of the forward probability. It measures the probability that, given we start in some state, we generate the remainder of the string (and stop in an accepting state). This is known as the backwards probability, and it can be decomposed recursively in a manner that is analogous to the forward probability.

---

← 25 The Forward Algorithm 27 Finite-State Automata →