hmm.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import copy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A Toy HMM\n",
    "Jiaxin Shi, ishijiaxin@126.com"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 0. Implementation\n",
    "\n",
    "This implementation has following algorithms for HMM.\n",
    "- Forward/Backward (sum-product)\n",
    "- Viterbi (max-product)\n",
    "- Baum-Welch (EM with exact E)\n",
    "- Gibbs (EM with approx E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#!/usr/bin/env python\n",
    "# -*- coding: utf-8 -*-\n",
    "\n",
    "import numpy as np\n",
    "import copy\n",
    "\n",
    "\n",
    "class HMM:\n",
    "    def __init__(self, states, transition, emission, init):\n",
    "        self.state_names = copy.copy(states)\n",
    "        self.n_states = len(states)\n",
    "        self.A = transition.copy()\n",
    "        self.B = emission.copy()\n",
    "        self.n_emissions = self.B.shape[1]\n",
    "        self.init = init\n",
    "\n",
    "    def generate(self, length):\n",
    "        state = self.init\n",
    "        states = []\n",
    "        ret = []\n",
    "        for i in xrange(1, length + 1):\n",
    "            state = np.random.choice(range(self.n_states), p=self.A[state])\n",
    "            states.append(state)\n",
    "            ret.append(\n",
    "                np.random.choice(range(self.n_emissions), p=self.B[state]))\n",
    "        print 'Generating by states:', ''.join(self.state_names[i] for i in states)\n",
    "        ret = ''.join([str(i) for i in ret])\n",
    "        return ret\n",
    "\n",
    "    def _forward(self, seq_arr):\n",
    "        T = len(seq_arr)\n",
    "        alpha = np.zeros((T + 1, self.n_states))\n",
    "        alpha[0, self.init] = 1\n",
    "        log_px = 0.\n",
    "        for t in xrange(1, T + 1):\n",
    "            alpha[t] = self.B[:, seq_arr[t - 1]] * \\\n",
    "                       np.dot(alpha[t - 1], self.A)\n",
    "            pt = alpha[t].sum()\n",
    "            alpha[t] /= pt\n",
    "            log_px += np.log(pt)\n",
    "        return alpha, log_px\n",
    "\n",
    "    def _backward(self, seq_arr):\n",
    "        T = len(seq_arr)\n",
    "        beta = np.zeros((T + 1, self.n_states))\n",
    "        beta[T, :] = 1\n",
    "        log_px = 0.\n",
    "        for t in xrange(T, 0, -1):\n",
    "            beta[t - 1] = np.dot(self.A, beta[t] * self.B[:, seq_arr[t - 1]])\n",
    "            pt = beta[t - 1].sum()\n",
    "            beta[t - 1] /= pt\n",
    "            log_px += np.log(pt)\n",
    "        log_px += np.log(beta[0, self.init])\n",
    "        return beta, log_px\n",
    "\n",
    "    def viterbi(self, seq):\n",
    "        # := max-product\n",
    "        seq_arr = np.array([int(i) for i in seq])\n",
    "        T = len(seq_arr)\n",
    "        T1 = np.zeros((self.n_states, T + 1))\n",
    "        T1[self.init, 0] = 1\n",
    "        T2 = np.zeros((self.n_states, T + 1), dtype='int')\n",
    "        states = np.zeros(T + 1, dtype='int')\n",
    "        for t in xrange(1, T + 1):\n",
    "            for j in xrange(self.n_states):\n",
    "                T1[j, t] = np.max(T1[:, t - 1] * self.A[:, j])\n",
    "                T1[j, t] *= self.B[j, seq_arr[t - 1]]\n",
    "                T2[j, t] = np.argmax(T1[:, t - 1] * self.A[:, j])\n",
    "        states[T] = np.argmax(T1[:, T])\n",
    "        for t in xrange(T, 1, -1):\n",
    "            states[t - 1] = T2[states[t], t - 1]\n",
    "        return ''.join([self.state_names[s] for s in states[1:]])\n",
    "\n",
    "    def baum_welch(self, seq):\n",
    "        # := EM\n",
    "        seq_arr = np.array([int(i) for i in seq])\n",
    "        T = len(seq_arr)\n",
    "        kesi = np.zeros((T + 1, self.n_states, self.n_states))\n",
    "        log_px = None\n",
    "        iter = 0\n",
    "        while True:\n",
    "            iter += 1\n",
    "            alpha, alpha_log_px = self._forward(seq_arr)\n",
    "            print \"Iter %d\" % iter, \"log p(x): %s\" % alpha_log_px\n",
    "            if log_px and (np.abs(\n",
    "                    log_px - alpha_log_px) < np.abs(1e-6 * log_px)):\n",
    "                print \"Converged.\"\n",
    "                break\n",
    "            beta, beta_log_px = self._backward(seq_arr)\n",
    "            try:\n",
    "                assert np.abs(\n",
    "                    alpha_log_px - beta_log_px) < np.abs(1e-6 * alpha_log_px)\n",
    "            except AssertionError as e:\n",
    "                print \"alpha_log_px:\", alpha_log_px\n",
    "                print \"beta_log_px:\", beta_log_px\n",
    "                raise e\n",
    "            log_px = alpha_log_px\n",
    "            gamma = alpha * beta\n",
    "            gamma /= np.sum(gamma, axis=1, keepdims=True)\n",
    "            for t in xrange(1, T):\n",
    "                kesi[t] = np.outer(\n",
    "                    alpha[t],\n",
    "                    beta[t + 1] * self.B[:, seq_arr[t + 1 - 1]]) * self.A\n",
    "            kesi[1:T] = kesi[1:T] / kesi[1:T].sum(axis=(1, 2), keepdims=True)\n",
    "            self.A = kesi[1:T].sum(axis=0) / \\\n",
    "                     gamma[1:T].sum(axis=0)[:, np.newaxis]\n",
    "            assert np.all(np.abs(1. - self.A.sum(axis=1)) < 1e-6)\n",
    "            obs = np.zeros((T + 1, self.n_emissions))\n",
    "            obs[range(1, T + 1), seq_arr] = 1\n",
    "            self.B = np.dot(gamma[1:].T, obs[1:]) / \\\n",
    "                     gamma[1:].sum(axis=0)[:, np.newaxis]\n",
    "        print \"Estimate A:\"\n",
    "        print np.array_str(self.A, precision=3)\n",
    "        print \"Estimate B:\"\n",
    "        print np.array_str(self.B, precision=3)\n",
    "        return log_px, self.A, self.B\n",
    "\n",
    "    def gibbs(self, seq, steps=1, burn_in=0, max_iters=None):\n",
    "        seq_arr = np.array([int(i) for i in seq])\n",
    "        T = len(seq_arr)\n",
    "        states = np.zeros(T + 1, dtype='int')\n",
    "        iter = 0\n",
    "        log_px = None\n",
    "        states[0] = self.init\n",
    "        while True:\n",
    "            iter += 1\n",
    "            alpha, alpha_log_px = self._forward(seq_arr)\n",
    "            print \"Iter %d\" % iter, \"log p(x): %s\" % alpha_log_px\n",
    "            if log_px and (np.abs(\n",
    "                    log_px - alpha_log_px) < np.abs(1e-6 * log_px)):\n",
    "                print \"Converged.\"\n",
    "                break\n",
    "            log_px = alpha_log_px\n",
    "            if max_iters and (iter >= max_iters):\n",
    "                break\n",
    "            A = np.zeros_like(self.A)\n",
    "            B = np.zeros_like(self.B)\n",
    "            for t in xrange(1, T + 1):\n",
    "                states[t] = np.random.choice(range(3))\n",
    "            for step in xrange(steps):\n",
    "                for t in xrange(1, T + 1):\n",
    "                    p_state_t = self.B[:, seq_arr[t - 1]] * \\\n",
    "                                self.A[states[t - 1]]\n",
    "                    if t < T:\n",
    "                        p_state_t *= self.A[:, states[t + 1]]\n",
    "                    p_state_t /= p_state_t.sum()\n",
    "                    states[t] = np.random.choice(range(3), p=p_state_t)\n",
    "                if step >= burn_in:\n",
    "                    for t in xrange(1, T + 1):\n",
    "                        if t < T:\n",
    "                            A[states[t], states[t + 1]] += 1\n",
    "                        B[states[t], seq_arr[t - 1]] += 1\n",
    "            A = np.maximum(1., A)\n",
    "            B = np.maximum(1., B)\n",
    "            self.A = A / A.sum(axis=1, keepdims=True)\n",
    "            self.B = B / B.sum(axis=1, keepdims=True)\n",
    "        print \"Estimate A:\"\n",
    "        print np.array_str(self.A, precision=3)\n",
    "        print \"Estimate B:\"\n",
    "        print np.array_str(self.B, precision=3)\n",
    "        return log_px"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 1.1 Generation\n",
    "We comment out long sequence (1000/10000) code here (visual burden). To run long sequence, following the comments in the code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1 Generation\n",
      "\n",
      "Generating by states: AAAAAABBBBCCCCCCCCCCCCCCAAAAAAAABBBBBBBCCAABBBCCCCCCCCCCAABBBBBCCCCCCCCCCCCCCCCABBBABBAAAAAAAAAAAAAA\n",
      "Inferred optimal state series:\n",
      "AAAAAABBBBBCCCCCCCCCCCCCCCAAAAAAAAAABBBBBBBBBBBCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCAAAAAAAAAAAAAAAAAA\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(1236)\n",
    "states = ['A', 'B', 'C']\n",
    "\n",
    "print \"1.1 Generation\\n\"\n",
    "transition = np.array([\n",
    "    [0.8, 0.2, 0.0],\n",
    "    [0.1, 0.7, 0.2],\n",
    "    [0.1, 0.0, 0.9]\n",
    "])\n",
    "emission = np.array([\n",
    "    [0.9, 0.1],\n",
    "    [0.5, 0.5],\n",
    "    [0.1, 0.9]\n",
    "])\n",
    "init = 0\n",
    "hmm = HMM(states, transition, emission, init)\n",
    "seqs = []\n",
    "for seq_len in [100, 1000, 10000]:\n",
    "    seq = hmm.generate(seq_len)\n",
    "    seqs.append(seq)\n",
    "    print \"Inferred optimal state series:\"\n",
    "    print hmm.viterbi(seq)\n",
    "    # NOTE: To run chains with various length, REMOVE this break\n",
    "    break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 / 1.3 Baum Welch\n",
    "We comment out variance computation code here (time consuming). To run variance computation, following the comments in the code.\n",
    "So we put the variance computed here:\n",
    "\n",
    "#### Variance for chain with length 100\n",
    "```\n",
    "Var(A):\n",
    "[[ 0.09052057  0.06006296  0.04894882]\n",
    " [ 0.06632638  0.1313508   0.05201452]\n",
    " [ 0.06691535  0.0444691   0.08084979]]\n",
    "Var(B):\n",
    "[[ 0.06158714  0.06158714]\n",
    " [ 0.14604622  0.14604622]\n",
    " [ 0.0887852   0.0887852 ]]\n",
    "```\n",
    "\n",
    "#### Variance for chain with length 1000\n",
    "```\n",
    "Var(A):\n",
    "[[ 0.08107341  0.0389215   0.02910226]\n",
    " [ 0.07240883  0.10532198  0.07765103]\n",
    " [ 0.04485569  0.07224511  0.08158469]]\n",
    "Var(B):\n",
    "[[ 0.08268158  0.08268158]\n",
    " [ 0.1226752   0.1226752 ]\n",
    " [ 0.10824732  0.10824732]]\n",
    "```\n",
    "\n",
    "#### Variance for chain with length 10000\n",
    "```\n",
    "Var(A):\n",
    "[[ 0.02264952  0.01980305  0.02054278]\n",
    " [ 0.02993416  0.06203205  0.02176658]\n",
    " [ 0.03476353  0.03295623  0.07207851]]\n",
    "Var(B):\n",
    "[[ 0.13777976  0.13777976]\n",
    " [ 0.10384683  0.10384683]\n",
    " [ 0.12860774  0.12860774]]\n",
    "```\n",
    "\n",
    "The conclusion is that the longer chain, the less variance in estimation of transition matrix (A). No similar results for observation matrix (B)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "1.2/1.3 Baum Welch\n",
      "\n",
      "Sequence length: 100\n",
      "Run 0\n",
      "Init transition:\n",
      "[[ 0.41141391  0.43933319  0.1492529 ]\n",
      " [ 0.64224114  0.283969    0.07378985]\n",
      " [ 0.21495536  0.31890642  0.46613823]]\n",
      "Init emission:\n",
      "[[ 0.63570068  0.36429932]\n",
      " [ 0.09348741  0.90651259]\n",
      " [ 0.32872192  0.67127808]]\n",
      "Iter 1 log p(x): -72.1931863341\n",
      "Iter 2 log p(x): -69.6917131442\n",
      "Iter 3 log p(x): -68.7798199905\n",
      "Iter 4 log p(x): -67.6678265601\n",
      "Iter 5 log p(x): -65.8220378798\n",
      "Iter 6 log p(x): -62.7960291425\n",
      "Iter 7 log p(x): -59.4489171376\n",
      "Iter 8 log p(x): -57.5687741293\n",
      "Iter 9 log p(x): -56.9082650051\n",
      "Iter 10 log p(x): -56.642186168\n",
      "Iter 11 log p(x): -56.5049763982\n",
      "Iter 12 log p(x): -56.4228378285\n",
      "Iter 13 log p(x): -56.3650342775\n",
      "Iter 14 log p(x): -56.3170377132\n",
      "Iter 15 log p(x): -56.2717475601\n",
      "Iter 16 log p(x): -56.2252960643\n",
      "Iter 17 log p(x): -56.1751210366\n",
      "Iter 18 log p(x): -56.1191914118\n",
      "Iter 19 log p(x): -56.0558369205\n",
      "Iter 20 log p(x): -55.983959736\n",
      "Iter 21 log p(x): -55.9035071474\n",
      "Iter 22 log p(x): -55.8159897143\n",
      "Iter 23 log p(x): -55.7246439012\n",
      "Iter 24 log p(x): -55.6338528549\n",
      "Iter 25 log p(x): -55.5479203057\n",
      "Iter 26 log p(x): -55.4699130604\n",
      "Iter 27 log p(x): -55.4012618033\n",
      "Iter 28 log p(x): -55.3421033039\n",
      "Iter 29 log p(x): -55.2918501335\n",
      "Iter 30 log p(x): -55.2496079669\n",
      "Iter 31 log p(x): -55.2143891796\n",
      "Iter 32 log p(x): -55.1852110474\n",
      "Iter 33 log p(x): -55.1611474394\n",
      "Iter 34 log p(x): -55.141359324\n",
      "Iter 35 log p(x): -55.1251101402\n",
      "Iter 36 log p(x): -55.1117690744\n",
      "Iter 37 log p(x): -55.1008059433\n",
      "Iter 38 log p(x): -55.0917813187\n",
      "Iter 39 log p(x): -55.0843346731\n",
      "Iter 40 log p(x): -55.0781723298\n",
      "Iter 41 log p(x): -55.0730562179\n",
      "Iter 42 log p(x): -55.0687939036\n",
      "Iter 43 log p(x): -55.0652300485\n",
      "Iter 44 log p(x): -55.0622392614\n",
      "Iter 45 log p(x): -55.0597202196\n",
      "Iter 46 log p(x): -55.0575908933\n",
      "Iter 47 log p(x): -55.0557847019\n",
      "Iter 48 log p(x): -55.0542474378\n",
      "Iter 49 log p(x): -55.052934817\n",
      "Iter 50 log p(x): -55.0518105295\n",
      "Iter 51 log p(x): -55.0508446915\n",
      "Iter 52 log p(x): -55.0500126137\n",
      "Iter 53 log p(x): -55.0492938206\n",
      "Iter 54 log p(x): -55.048671268\n",
      "Iter 55 log p(x): -55.0481307161\n",
      "Iter 56 log p(x): -55.0476602262\n",
      "Iter 57 log p(x): -55.0472497548\n",
      "Iter 58 log p(x): -55.0468908257\n",
      "Iter 59 log p(x): -55.0465762632\n",
      "Iter 60 log p(x): -55.0462999745\n",
      "Iter 61 log p(x): -55.0460567725\n",
      "Iter 62 log p(x): -55.0458422287\n",
      "Iter 63 log p(x): -55.0456525536\n",
      "Iter 64 log p(x): -55.0454844964\n",
      "Iter 65 log p(x): -55.0453352625\n",
      "Iter 66 log p(x): -55.0452024442\n",
      "Iter 67 log p(x): -55.0450839633\n",
      "Iter 68 log p(x): -55.0449780225\n",
      "Iter 69 log p(x): -55.0448830649\n",
      "Iter 70 log p(x): -55.0447977394\n",
      "Iter 71 log p(x): -55.0447208721\n",
      "Iter 72 log p(x): -55.0446514409\n",
      "Iter 73 log p(x): -55.0445885555\n",
      "Iter 74 log p(x): -55.0445314385\n",
      "Iter 75 log p(x): -55.0444794112\n",
      "Converged.\n",
      "Estimate A:\n",
      "[[  9.222e-01   7.773e-02   5.672e-05]\n",
      " [  1.400e-01   5.146e-01   3.454e-01]\n",
      " [  1.348e-02   1.557e-01   8.308e-01]]\n",
      "Estimate B:\n",
      "[[  8.004e-01   1.996e-01]\n",
      " [  6.336e-01   3.664e-01]\n",
      " [  5.178e-06   1.000e+00]]\n",
      "Final log p(x): -55.0445314385\n",
      "Optimal state series: AAAAAAAAABBCCCCBBCCCCCCCCCBAAAAAAAAAAAAAAAAAABBCCCCCCCCCCCBBCCCCCCCCCCCCCBCCCCCCBBAAAAAAAAAAAAAAAAAA\n",
      "Variance of estimated:\n",
      "Var(A):\n",
      "[[ 0.  0.  0.]\n",
      " [ 0.  0.  0.]\n",
      " [ 0.  0.  0.]]\n",
      "Var(B):\n",
      "[[ 0.  0.]\n",
      " [ 0.  0.]\n",
      " [ 0.  0.]]\n"
     ]
    }
   ],
   "source": [
    "print \"\\n1.2/1.3 Baum Welch\"\n",
    "for seq in seqs:\n",
    "    print \"\\nSequence length:\", len(seq)\n",
    "    As = []\n",
    "    Bs = []\n",
    "    for run in xrange(10):\n",
    "        print \"Run\", run\n",
    "        transition2 = np.random.random((3, 3))\n",
    "        transition2 /= transition2.sum(axis=1, keepdims=True)\n",
    "        emission2 = np.random.random((3, 2))\n",
    "        emission2 /= emission2.sum(axis=1, keepdims=True)\n",
    "        print \"Init transition:\"\n",
    "        print transition2\n",
    "        print \"Init emission:\"\n",
    "        print emission2\n",
    "        hmm2 = HMM(states, transition2, emission2, init)\n",
    "        log_px, A, B = hmm2.baum_welch(seq)\n",
    "        As.append(A)\n",
    "        Bs.append(B)\n",
    "        print \"Final log p(x):\", log_px\n",
    "        print \"Optimal state series:\", hmm2.viterbi(seq)\n",
    "        # NOTE: To calculate variance, REMOVE this break\n",
    "        break\n",
    "    print \"Variance of estimated:\"\n",
    "    print \"Var(A):\"\n",
    "    print np.var(As, axis=0)\n",
    "    print \"Var(B):\"\n",
    "    print np.var(Bs, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare with true parameters/states\n",
    "Find estimated paramters/states in above results (only one run presented here).\n",
    "The estimated parameters in most run is very similar as the true parameter. But some run also got stuck in local optimal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Gibbs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "2.1 Gibbs\n",
      "Init transition:\n",
      "[[ 0.34110724  0.39255971  0.26633305]\n",
      " [ 0.24355523  0.37788669  0.37855808]\n",
      " [ 0.32565081  0.21240733  0.46194186]]\n",
      "Init emission:\n",
      "[[ 0.507358    0.492642  ]\n",
      " [ 0.37109127  0.62890873]\n",
      " [ 0.506693    0.493307  ]]\n",
      "Sequence length: 100\n",
      "Iter 1 log p(x): -68.9635658108\n",
      "Iter 2 log p(x): -68.9547155068\n",
      "Iter 3 log p(x): -68.9419872001\n",
      "Iter 4 log p(x): -68.9193644246\n",
      "Iter 5 log p(x): -68.8234956815\n",
      "Iter 6 log p(x): -68.7058509947\n",
      "Iter 7 log p(x): -68.5225882733\n",
      "Iter 8 log p(x): -68.2710770794\n",
      "Iter 9 log p(x): -67.955095736\n",
      "Iter 10 log p(x): -67.147018193\n",
      "Iter 11 log p(x): -65.453949634\n",
      "Iter 12 log p(x): -62.5769123357\n",
      "Iter 13 log p(x): -59.1442762918\n",
      "Iter 14 log p(x): -57.2899042789\n",
      "Iter 15 log p(x): -56.7861824909\n",
      "Iter 16 log p(x): -56.7136273509\n",
      "Iter 17 log p(x): -56.530892125\n",
      "Iter 18 log p(x): -56.4306142307\n",
      "Iter 19 log p(x): -56.3663842661\n",
      "Iter 20 log p(x): -56.3058989836\n",
      "Iter 21 log p(x): -56.2953490606\n",
      "Iter 22 log p(x): -56.2905485485\n",
      "Iter 23 log p(x): -56.3115588474\n",
      "Iter 24 log p(x): -56.2196551288\n",
      "Iter 25 log p(x): -56.2098604644\n",
      "Iter 26 log p(x): -56.2004136265\n",
      "Iter 27 log p(x): -56.3049885647\n",
      "Iter 28 log p(x): -56.2782506498\n",
      "Iter 29 log p(x): -56.3411704107\n",
      "Iter 30 log p(x): -56.3078423493\n",
      "Iter 31 log p(x): -56.3014630553\n",
      "Iter 32 log p(x): -56.2430228084\n",
      "Iter 33 log p(x): -56.2490560555\n",
      "Iter 34 log p(x): -56.2679411542\n",
      "Iter 35 log p(x): -56.2669972038\n",
      "Iter 36 log p(x): -56.2757044094\n",
      "Iter 37 log p(x): -56.2486786484\n",
      "Iter 38 log p(x): -56.1886326357\n",
      "Iter 39 log p(x): -56.179938699\n",
      "Iter 40 log p(x): -56.1710654908\n",
      "Iter 41 log p(x): -56.1782044197\n",
      "Iter 42 log p(x): -56.1531795046\n",
      "Iter 43 log p(x): -56.1639667295\n",
      "Iter 44 log p(x): -56.1699437791\n",
      "Iter 45 log p(x): -56.1656138302\n",
      "Iter 46 log p(x): -56.1425009422\n",
      "Iter 47 log p(x): -56.1680697995\n",
      "Iter 48 log p(x): -56.1537080292\n",
      "Iter 49 log p(x): -56.3718503968\n",
      "Iter 50 log p(x): -56.2160107945\n",
      "Iter 51 log p(x): -56.1660595442\n",
      "Iter 52 log p(x): -56.1733735962\n",
      "Iter 53 log p(x): -56.1715064641\n",
      "Iter 54 log p(x): -56.1541761832\n",
      "Iter 55 log p(x): -56.1760397114\n",
      "Iter 56 log p(x): -56.1590691818\n",
      "Iter 57 log p(x): -56.1948070584\n",
      "Iter 58 log p(x): -56.2151919709\n",
      "Iter 59 log p(x): -56.1837578861\n",
      "Iter 60 log p(x): -56.1926783361\n",
      "Iter 61 log p(x): -56.158465974\n",
      "Iter 62 log p(x): -56.2346167582\n",
      "Iter 63 log p(x): -56.179525632\n",
      "Iter 64 log p(x): -56.163651655\n",
      "Iter 65 log p(x): -56.1697728998\n",
      "Iter 66 log p(x): -56.1760555686\n",
      "Iter 67 log p(x): -56.1505635554\n",
      "Iter 68 log p(x): -56.1563390631\n",
      "Iter 69 log p(x): -56.1431079816\n",
      "Iter 70 log p(x): -56.1663123173\n",
      "Iter 71 log p(x): -56.1768555173\n",
      "Iter 72 log p(x): -56.173168148\n",
      "Iter 73 log p(x): -56.1728545061\n",
      "Iter 74 log p(x): -56.1449063307\n",
      "Iter 75 log p(x): -56.1415962865\n",
      "Iter 76 log p(x): -56.134820936\n",
      "Iter 77 log p(x): -56.1325251204\n",
      "Iter 78 log p(x): -56.1785449242\n",
      "Iter 79 log p(x): -56.1364104534\n",
      "Iter 80 log p(x): -56.1263359023\n",
      "Iter 81 log p(x): -56.1570526201\n",
      "Iter 82 log p(x): -56.1273921788\n",
      "Iter 83 log p(x): -56.1277285391\n",
      "Iter 84 log p(x): -56.1307451583\n",
      "Iter 85 log p(x): -56.1450480118\n",
      "Iter 86 log p(x): -56.1501074343\n",
      "Iter 87 log p(x): -56.1847383628\n",
      "Iter 88 log p(x): -56.1471443121\n",
      "Iter 89 log p(x): -56.1204377714\n",
      "Iter 90 log p(x): -56.1327416705\n",
      "Iter 91 log p(x): -56.1232768176\n",
      "Iter 92 log p(x): -56.1492637428\n",
      "Iter 93 log p(x): -56.1675329128\n",
      "Iter 94 log p(x): -56.1816697985\n",
      "Iter 95 log p(x): -56.1599693329\n",
      "Iter 96 log p(x): -56.1781762459\n",
      "Iter 97 log p(x): -56.1377558876\n",
      "Iter 98 log p(x): -56.1517287984\n",
      "Iter 99 log p(x): -56.2254692928\n",
      "Iter 100 log p(x): -56.1966604144\n",
      "Estimate A:\n",
      "[[ 0.301  0.115  0.584]\n",
      " [ 0.027  0.941  0.031]\n",
      " [ 0.418  0.027  0.555]]\n",
      "Estimate B:\n",
      "[[ 0.82   0.18 ]\n",
      " [ 0.127  0.873]\n",
      " [ 0.793  0.207]]\n",
      "Final log p(x): -56.1966604144\n",
      "Optimal state series: ACCCCCCCABBBBBBBBBBBBBBBBBCCCCCCCCCCCCCCCCCCABBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBCCCCCCCCCCCCC\n"
     ]
    }
   ],
   "source": [
    "print \"\\n2.1 Gibbs\"\n",
    "transition3 = np.random.random((3, 3))\n",
    "transition3 /= transition3.sum(axis=1, keepdims=True)\n",
    "emission3 = np.random.random((3, 2))\n",
    "emission3 /= emission3.sum(axis=1, keepdims=True)\n",
    "print \"Init transition:\"\n",
    "print transition3\n",
    "print \"Init emission:\"\n",
    "print emission3\n",
    "print \"Sequence length:\", len(seqs[0])\n",
    "hmm3 = HMM(states, transition3, emission3, init)\n",
    "log_px = hmm3.gibbs(seqs[0], steps=200, burn_in=100, max_iters=100)\n",
    "print \"Final log p(x):\", log_px\n",
    "print \"Optimal state series:\", hmm3.viterbi(seqs[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare with true parameters/states\n",
    "Find estimated paramters/states in above results (only one run presented here).\n",
    "The estimated parameters in most run is very similar as the true parameter (get similar likelihood with Baum Welch). But some run also got stuck in local optimal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 Compare results (Gibbs vs. Baum-Welch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Baum-Welch is an EM algorithm with exact inference for latent states. It is more likely to get stuck in local optimal. And due to exact inference, there is no chance for Baum-Welch to jump out of the local optimal.\n",
    "- EM with Gibbs inference for latent states is also an EM, but with E-step inference by gibbs sampling, which is approximating posterior by samples drawn. \n",
    "    - With short chain: This method has stochastic nature and can jump out of some local optimals caused by EM and get better results.\n",
    "    - With long chain: This method is approximating the posterior as well as exact inference by Baum Welch. So their results are similar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 MEMM: Gibbs pesudo code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "E: Gibbs sampling the latent states\n",
    "y[t] ~ p(y[t]|y[t-1], x[t])p(y[t+1]|y[t], x[t+1])\n",
    "M: Re-estimate parameters\n",
    "Restimate parameters of p(y[t]|y[t-1], x[t]) using sampled y[t] and y[t-1].\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12+"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}