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

(def vocabulary '(Call me Ishmael))

(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?)))))

(defn all-but-last [l]
  (cond (empty? l) (throw "bad thing")
        (empty? (rest l)) '()
        :else (cons (first l) (all-but-last (rest l)))))

The algorithm we examined for HMM parameter estimation, the Baum-Welch algorithm, is a special case of a very general and useful class of algorithms known as expectation maximization (EM) algorithms.

Expectation maximization algorithms are a general way to search for a maximum likelihood or maximum a posteriori point estimate of some parameters when there are hidden variables that must be marginalized away thus making an direct solution impossible. In other words, it is a method for approximately optimizing the objective function below, which contains a sum over some latent variables.

\[\DeclareMathOperator*{\argmax}{arg\,max} \underset {\theta }\argmax \log \Pr(D | \theta) = \underset {\theta }\argmax \log \sum_{H} \Pr(D, H \mid \theta)\]

It can also be used for searching for a point estimate of \(\theta\) that optimizes a Bayesian criterion \(\log \Pr(D | \theta)\Pr(\theta) = \log \sum_{H} \Pr(D, H \mid \theta) \Pr(\theta)\). Note that this is not full Bayesian inference, since it only searches for a single set of parameter values for the whole parameter space.

The basic idea behind EM is that instead of directly optimizing the intractable objective function above, we work with an alternative objective.

\[\DeclareMathOperator*{\argmax}{arg\,max} \underset {\theta }\argmax \sum_{H} P(H | D, \theta^{\mathrm{old}}) \log P(D,H \mid \theta)\]

This objective makes use of the posterior over our latent variables \(P(H | D, \theta^\mathrm{old})\) which must be computed before the parameters can be optimized. This leads to a two-step algorithm which alternates between computing this posterior in terms of our old parameters \(\theta^{\mathrm{old}}\) and then opimizing the parameters in terms of this posterior.

It is beyong the scope of this course, but it can be proven that this two stage process is equivalent to optimizing a lower bound on the incomplete data log likelihood. However, EM algorithms do not provide any general guarantees that the parameters they discover are globally optimal. They simply guarantee that the likelihood of the data will not be decreased (and will likely increase) on each step of the iteration. But they can get stuck in local maxima of the search space.

---

← 36 The Baum-Welch Algorithm 38 Gibbs Sampling →