advanced/practicum/week-3/lab-3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Networks and HMM\n",
    "\n",
    "$\n",
    "\\newcommand{\\x}{\\mathbf{x}}\n",
    "\\newcommand{\\y}{\\mathbf{y}}\n",
    "$\n",
    "\n",
    "| | |\n",
    "|--|--|\n",
    "| **Names** | *Thomas Brus & Jan Ubbo van Baardewijk* |\n",
    "| **Group** | *ML_HMI_01* |\n",
    "\n",
    "\n",
    "In this lab, you are given a basic implementation of the components of Factor graphs and the sum-product algorithm. With this, you can compute marginal probabilities and, by incorporating extra factors, joint probabilities of nodes and observations. Familiarise yourself with the code, and then answer the questions below.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "import scipy.ndimage as ndim\n",
    "import matplotlib.pyplot as plt\n",
    "import operator\n",
    "import pylab\n",
    "\n",
    "pylab.rcParams['figure.figsize'] = (10.0, 8.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Factor:\n",
    "    \"\"\"Implement a factor node in a Bayesian Network\"\"\"\n",
    "    def __init__(self,name,factor):\n",
    "        \"\"\"Initialisation with the name of the node, and the implemented function of the connected variables\n",
    "        \n",
    "        The function `factor` must take a tuple as argument, with as many elements in the tuple as there \n",
    "        are connected variables to the factor node and the elements of the argument being the indices of the \n",
    "        possible values in the variable node. It returns the factor value for that setting of variables.\n",
    "        \n",
    "        Example: a factor implementing p(A), with A being a boolean variable with possible values {0: True, 1: False},\n",
    "        factor((1,)) should return p(~a).\n",
    "        \"\"\"\n",
    "        self.name = name\n",
    "        self.fn = factor\n",
    "        self.connections = [] # Pointer to the nodes this factor is connected to\n",
    "        self.msgRcvd = []     # Stores messages received from the connected nodes\n",
    "        self.msgSent = []     # Has a message been sent to the connected nodes\n",
    "        self.nMsgsRcvd = 0    # number of messages received\n",
    "    \n",
    "    def addConnection(self, node):\n",
    "        \"\"\"Register a connection to a variable. The connections must be added to the factor node\n",
    "        in the same order as the indices are expected by the factor function\n",
    "        \n",
    "        This automatically registers the inverse connection from variable to the factor.\"\"\"\n",
    "        self.connections.append(node)\n",
    "        self.msgSent.append(False)\n",
    "        self.msgRcvd.append(None)\n",
    "        node.addConnection(self)\n",
    "        \n",
    "    def __str__(self):\n",
    "        return self.name\n",
    "    \n",
    "    def __repr__(self):\n",
    "        return \"Factor \"+self.name\n",
    "        \n",
    "    def recPrint(self, parent=None, indent=\"\"):\n",
    "        \"\"\"Recursively print the tree\"\"\"\n",
    "        print indent,self.__str__()\n",
    "        for i in range(len(self.connections)):\n",
    "            if self.connections[i]!=parent:\n",
    "                self.connections[i].recPrint(self,indent+\"  \")\n",
    "            \n",
    "    def recvMsg(self, sender, msg):\n",
    "        \"\"\"Receive a message `msg` from a connected `sender` variable node\"\"\"\n",
    "#         print self.name,\"Receiving message from\",sender,\":\",msg\n",
    "        for i in range(len(self.connections)):\n",
    "            if self.connections[i] == sender:\n",
    "                if self.msgRcvd[i]!=None:\n",
    "                    print \"... ignoring\"\n",
    "                    return\n",
    "                self.msgRcvd[i] = msg\n",
    "                self.nMsgsRcvd += 1\n",
    "                return\n",
    "        print \"ERROR:\",self.__str__(),\"could not find the corresponding variable\",sender.__str__()\n",
    "   \n",
    "    def resetMsgs(self):\n",
    "        self.nMsgsRcvd=0\n",
    "        self.msgRcvd = [ None for m in self.msgRcvd ]\n",
    "        self.msgSent = [ False for m in self.msgSent ]\n",
    "\n",
    "    def recursiveEnumerate(self, outIndex, valIndices, values, res):\n",
    "        \"\"\"Stores, in list `res`, for each possible value of the node connected to connection `outindex`,\n",
    "            the sum over all possible combinations of values of the incoming connections of the product \n",
    "            of the factor of these values with the messages with these values. (Aaaargh.) \"\"\"\n",
    "        if len(valIndices) == len(self.connections): # end of recursion\n",
    "            prod = self.fn(tuple(valIndices))\n",
    "            for i in range(len(valIndices)):\n",
    "                if i != outIndex:\n",
    "                    prod *= self.msgRcvd[i][valIndices[i]]\n",
    "            \n",
    "            res[valIndices[outIndex]] += prod \n",
    "            return\n",
    "        \n",
    "        for i,v in enumerate(self.connections[len(valIndices)].values):\n",
    "            self.recursiveEnumerate(outIndex,valIndices+[i],values+[v],res)\n",
    "        \n",
    "        \n",
    "    def sendMessages(self):\n",
    "        \"\"\"Send a message to all connections on which no message has been sent yet, if all\n",
    "        other connections have received a message already\"\"\"        \n",
    "        \n",
    "        # Special case: single connection\n",
    "        if not self.msgSent[0] and len(self.connections) == 1:\n",
    "            self.msgSent[0] = True    # Record that a message has been sent on that connection\n",
    "            self.connections[0].recvMsg(self,np.array([self.fn((i,)) for i,v in enumerate(self.connections[0].values) ]))\n",
    "            return True\n",
    "        \n",
    "        if self.nMsgsRcvd+1<len(self.connections): # More than one connection without incoming message\n",
    "            return False                           # Can't send anything\n",
    "        \n",
    "        # If one connection has not received it's message, send a msg out on that connection\n",
    "        if self.nMsgsRcvd+1 == len(self.connections): \n",
    "            for i in range(len(self.connections)):    # For connections on which we could send a message                \n",
    "                if self.msgRcvd[i] is None:           # Find the one that hasn't received a message yet\n",
    "                    if not self.msgSent[i]:           # If no message has been sent on that connection yet\n",
    "                        msg = np.array([ 0. for c in self.connections[i].values ]) # Create a message with an entry for each possible value of the outgoing node\n",
    "                        self.recursiveEnumerate(i,[],[],msg) # Compute the message value\n",
    "                        self.connections[i].recvMsg(self,msg) \n",
    "                        self.msgSent[i] = True\n",
    "                        return True\n",
    "            \n",
    "        if self.nMsgsRcvd == len(self.connections): # message received on all connections, send wherever not done yet.\n",
    "            sentAny = False\n",
    "            for i in range(len(self.connections)):       # For connections on which we could send a message\n",
    "                if self.msgSent[i]:                      # ignore if handled already\n",
    "                    continue\n",
    "                msg = np.array([ 0. for c in self.connections[i].values ])\n",
    "                self.recursiveEnumerate(i,[],[],msg)\n",
    "                self.connections[i].recvMsg(self,msg)\n",
    "                self.msgSent[i] = True\n",
    "                sentAny = True\n",
    "            return sentAny\n",
    "                            \n",
    "        return False\n",
    "            \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class Variable:\n",
    "    \"\"\"Implement a Variable node with categorical values\"\"\"\n",
    "    def __init__(self,name,values):\n",
    "        \"\"\"Initialise with the name of the variable and a list of possible values for that variable\"\"\"\n",
    "        self.name = name\n",
    "        self.connections = []\n",
    "        self.values = values\n",
    "        self.msgRcvd = []     # Stores messages received from the connected nodes\n",
    "        self.msgSent = []     # Has a message been sent to the connected nodes\n",
    "        self.nMsgsRcvd = 0    # number of messages received\n",
    "    \n",
    "    def addConnection(self, node):\n",
    "        \"\"\"Add a connection to a factor node. This should never be called explicitly, except by\n",
    "        the addConnection method of the Factor class\"\"\"\n",
    "        self.connections.append(node)\n",
    "        self.msgSent.append(False)\n",
    "        self.msgRcvd.append(None)\n",
    "        \n",
    "    def __str__(self):\n",
    "        return self.name\n",
    "    def __repr__(self):\n",
    "        return \"Var \"+self.name\n",
    "    \n",
    "    \n",
    "    def recPrint(self, parent=None, indent=\"\"):\n",
    "        \"\"\"Recursively print the tree\"\"\"\n",
    "        print indent,self.__str__()\n",
    "        for i in range(len(self.connections)):\n",
    "            if self.connections[i]!=parent:\n",
    "                self.connections[i].recPrint(self,indent+\"  \")        \n",
    "    \n",
    "    def recvMsg(self, sender, msg):\n",
    "        \"\"\"Receive a message `msg` from Factor node `sender`\"\"\"\n",
    "#         print self.name,\"Receiving message from\",sender,\":\",msg\n",
    "        for i in range(len(self.connections)):\n",
    "            if self.connections[i] == sender:\n",
    "                if self.msgRcvd[i]!=None:\n",
    "                    print \"WARNING: ignoring a message\"\n",
    "                    return\n",
    "                self.msgRcvd[i] = msg\n",
    "                self.nMsgsRcvd += 1\n",
    "                return\n",
    "        print \"ERROR:\",self.__str__(),\"could not find the corresponding factor\",sender.__str__()\n",
    "        \n",
    "    def resetMsgs(self):\n",
    "        self.nMsgsRcvd=0\n",
    "        self.msgRcvd = [ None for m in self.msgRcvd ]\n",
    "        self.msgSent = [ False for m in self.msgSent ]\n",
    "        \n",
    "    \n",
    "    def sendMessages(self):\n",
    "        \"\"\"Send a message to all connections on which no message has been sent yet, if all\n",
    "        other connections have received a message already\n",
    "        \"\"\"    \n",
    "        # Special case: single connection\n",
    "        if not self.msgSent[0] and len(self.connections) == 1:\n",
    "            self.msgSent[0] = True\n",
    "            self.connections[0].recvMsg(self,np.array([1. for v in self.values]))\n",
    "            return True\n",
    "        \n",
    "        # If one connection has not received it's message, send a msg out on that connection\n",
    "        if self.nMsgsRcvd+1 == len(self.connections): \n",
    "            con = 0\n",
    "            for i in range(len(self.connections)):    # For connections on which we could send a message                \n",
    "                if self.msgRcvd[i] == None:           # Find the one that hasn't received a message yet\n",
    "                    con = i\n",
    "                    break\n",
    "                    \n",
    "            if self.msgSent[con]:                     # Don't send it again\n",
    "                return False\n",
    "            \n",
    "            msg = np.ones(len(self.values))           # The message is the product of the \n",
    "            for i in range(len(self.connections)):    # incoming messages on all connections,\n",
    "                if i == con:                          # except the outgoing one\n",
    "                    continue\n",
    "                msg *= self.msgRcvd[i]\n",
    "            self.connections[con].recvMsg(self,msg)   # Send the message\n",
    "            self.msgSent[con] = True\n",
    "            return True\n",
    "            \n",
    "        if self.nMsgsRcvd == len(self.connections): # message received on all connections, send wherever not done yet.\n",
    "            sentAny = False\n",
    "            for con in range(len(self.connections)):       # For connections on which we could send a message\n",
    "                if self.msgSent[con]:                      # ignore if handled already\n",
    "                    continue\n",
    "                msg = np.ones(len(self.values))            # Initialise the message\n",
    "                for i in range(len(self.connections)):     # multiply with all relevant incoming messages\n",
    "                    if i == con:\n",
    "                        continue\n",
    "                    msg *= self.msgRcvd[i]\n",
    "                self.connections[con].recvMsg(self,msg)    # and send it out\n",
    "                self.msgSent[con] = True                \n",
    "                sentAny = True\n",
    "            return sentAny\n",
    "\n",
    "        return False\n",
    "    \n",
    "    def getMarginals(self):\n",
    "        \"\"\"Once all messages have been received, the marginal probability \n",
    "        (or joint probability with specific values of other nodes, if such \n",
    "        factors have been implemented), is the product of all incoming messages\"\"\"\n",
    "        if self.nMsgsRcvd != len(self.connections):        # Check that all messages have been received\n",
    "            print \"ERROR printing marginals: messages missing\"\n",
    "        prod = self.msgRcvd[0]                             # Compute the product\n",
    "        for m in self.msgRcvd[1:]:\n",
    "            prod *= m\n",
    "        return prod                              "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A [ 0.77777778  0.22222222]\n",
      "C [ 0.33333333  0.66666667]\n",
      "B [ 0.08333333  0.91666667]\n",
      "D [ 0.16666667  0.83333333]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/thomasbrus/Code/machine-learning/lib/python2.7/site-packages/ipykernel/__main__.py:65: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.\n"
     ]
    }
   ],
   "source": [
    "class Graph:\n",
    "    def __init__(self):\n",
    "        self.var = {}\n",
    "        self.fact = {}\n",
    "    def addVariables(self,values, *names):\n",
    "        for n in names:\n",
    "            self.var[n] = Variable(n,values)\n",
    "    def addFactor(self,name,fTable, *connections):\n",
    "        self.fact[name] = Factor(name,lambda(v) : fTable[v])\n",
    "        for c in connections:\n",
    "            self.fact[name].addConnection(self.var[c])\n",
    "    def propagateMessages(self):\n",
    "        hasSent = True\n",
    "        while hasSent:\n",
    "            hasSent = False\n",
    "            for k,f in self.fact.iteritems():\n",
    "                hasSent |= f.sendMessages()\n",
    "            for k,v in self.var.iteritems():\n",
    "                hasSent |= v.sendMessages()\n",
    "    def printMarginals(self):\n",
    "        # Print out the marginals\n",
    "        for k,v in self.var.iteritems():\n",
    "            print v,v.getMarginals()\n",
    "            \n",
    "    def getMarginals(self, n):\n",
    "        return self.var[n].getMarginals()\n",
    "    \n",
    "    def printTree(self):\n",
    "        for k,v in self.var.iteritems():\n",
    "            print v.recPrint();\n",
    "\n",
    "    def reset(self):\n",
    "        for k,f in self.fact.iteritems():\n",
    "            f.resetMsgs()\n",
    "        for k,v in self.var.iteritems():\n",
    "            v.resetMsgs()\n",
    "\n",
    "\n",
    "pB = np.array([1./12,11./12])\n",
    "pC = np.array([1./3,2./3])\n",
    "pA_BC = np.array([[[0, 0.25],[.5,1.]],[[1,.75],[.5,0.]]])\n",
    "pD_A = np.array([[0,.75],[1,.25]])\n",
    "\n",
    "g = Graph()\n",
    "g.addVariables([\"True\",\"False\"], \"A\",\"B\",\"C\",\"D\")\n",
    "g.addFactor(\"p(B)\", pB, \"B\")\n",
    "g.addFactor(\"p(C)\", pC, \"C\")\n",
    "g.addFactor(\"p(A|B,C)\", pA_BC, \"A\",\"B\",\"C\")\n",
    "g.addFactor(\"P(D|A)\",pD_A,\"D\",\"A\")\n",
    "\n",
    "g.propagateMessages()\n",
    "g.printMarginals()\n",
    "g.reset()\n",
    "\n",
    "# g.addFactor(\"obs B\", np.array([1.,0.]),\"B\")\n",
    "# g.propagateMessages()\n",
    "# g.printMarginals()\n",
    "\n",
    "# g.printTree();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Question 1 [10 marks]** Using the above code, implement the computation of the following probabilities:\n",
    "\n",
    "* $p(B|\\neg A)$\n",
    "* $p(C|\\neg A)$\n",
    "* $p(C|\\neg A,\\neg B)$\n",
    "* $p(\\neg A, \\neg B, C, \\neg D)$\n",
    "* $p(B,C|D)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Q1.1:\n",
      "p(B|-A) =  0.0694444444444\n",
      "\n",
      "Q1,2:\n",
      "p(C|-A) =  0.180555555556\n",
      "\n",
      "Q1.3:\n",
      "p(C|-A,-B) =  0.152777777778\n",
      "\n",
      "Q1.4:\n",
      "P(C|-A,-B) = (see Q1.3) = 0.152777777778\n",
      "P(-B|-A) = (see Q1.1) = 0.025462962963\n",
      "P(-A) = (see marginal probabilities of initial graph) = 0.22222222\n",
      "P(-D|-A) = (given) = 0.25\n",
      "P(-A,-B,C,-D) = P(C|-A,-B) * P(B|-A) * P(-A) * P(-D|-A) =  0.000216120825457\n",
      "\n",
      "Q1.5:\n",
      "p(B|A) =  0.0138888888889\n",
      "p(C|A) =  0.152777777778\n",
      "p(D|A) =  0.0\n",
      "P(A) = 0.77777778\n",
      "\n",
      "P(B|-A) = 0.0694444444444\n",
      "P(C|-A) = 0.180555555556\n",
      "P(D|-A) = 0.166666666667\n",
      "P(-A) = 0.22222222\n",
      "\n",
      "P(B,C,D) = P(B|A) * P(C|A) * P(D|A) * P(A) + P(B|-A) * P(C|-A) * P(D|-A) * P(-A) =  0.000464391856355\n",
      "P(D) = (see marginial probabilities of initial graph) = 0.16666667\n",
      "\n",
      "P(B,C|D) = P(B,C,D) / P(D) =  0.0027863510824\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/thomasbrus/Code/machine-learning/lib/python2.7/site-packages/ipykernel/__main__.py:65: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.\n"
     ]
    }
   ],
   "source": [
    "# Answer Q1\n",
    "\n",
    "# https://www.cl.cam.ac.uk/teaching/1011/L101/ml4lp-lect3.pdf\n",
    "\n",
    "g.addFactor(\"obs -A\", np.array([0.,1.]) ,\"A\")\n",
    "g.propagateMessages()\n",
    "\n",
    "# Answer Q1.1\n",
    "print \"Q1.1:\"\n",
    "print \"p(B|-A) = \", g.getMarginals(\"B\")[0]\n",
    "# print \"p(-B|-A) = \", g.getMarginals(\"B\")[1]\n",
    "\n",
    "# Answer Q1.2\n",
    "print \"\\nQ1,2:\"\n",
    "print \"p(C|-A) = \", g.getMarginals(\"C\")[0]\n",
    "\n",
    "# print \"p(D|-A) = \", g.getMarginals(\"D\")[0]\n",
    "\n",
    "g.reset()\n",
    "\n",
    "g.addFactor(\"obs -B\", np.array([0.,1.]), \"B\")\n",
    "g.propagateMessages()\n",
    "\n",
    "# Answer Q1.3\n",
    "print \"\\nQ1.3:\"\n",
    "print \"p(C|-A,-B) = \", g.getMarginals(\"C\")[0]\n",
    "\n",
    "g.reset()\n",
    "\n",
    "# Answer Q1.4\n",
    "print \"\\nQ1.4:\"\n",
    "print \"P(C|-A,-B) = (see Q1.3) = 0.152777777778\"\n",
    "print \"P(-B|-A) = (see Q1.1) = 0.025462962963\"\n",
    "print \"P(-A) = (see marginal probabilities of initial graph) = 0.22222222\"\n",
    "print \"P(-D|-A) = (given) = 0.25\"\n",
    "\n",
    "print \"P(-A,-B,C,-D) = P(C|-A,-B) * P(B|-A) * P(-A) * P(-D|-A) = \", 0.152777777778 * 0.025462962963 * 0.22222222 * 0.25\n",
    "\n",
    "# Answer Q1.5\n",
    "g2 = Graph()\n",
    "g2.addVariables([\"True\",\"False\"], \"A\",\"B\",\"C\",\"D\")\n",
    "g2.addFactor(\"p(B)\", pB, \"B\")\n",
    "g2.addFactor(\"p(C)\", pC, \"C\")\n",
    "g2.addFactor(\"p(A|B,C)\", pA_BC, \"A\",\"B\",\"C\")\n",
    "g2.addFactor(\"P(D|A)\",pD_A,\"D\",\"A\")\n",
    "\n",
    "g2.addFactor(\"obs A\", np.array([1.,0.]), \"A\")\n",
    "g2.propagateMessages()\n",
    "\n",
    "print \"\\nQ1.5:\"\n",
    "print \"p(B|A) = \", g2.getMarginals(\"B\")[0]\n",
    "print \"p(C|A) = \", g2.getMarginals(\"C\")[0]\n",
    "print \"p(D|A) = \", g2.getMarginals(\"D\")[0]\n",
    "print \"P(A) = 0.77777778\"\n",
    "\n",
    "print \"\\nP(B|-A) = 0.0694444444444\"\n",
    "print \"P(C|-A) = 0.180555555556\"\n",
    "print \"P(D|-A) = 0.166666666667\"\n",
    "print \"P(-A) = 0.22222222\"\n",
    "\n",
    "print \"\\nP(B,C,D) = P(B|A) * P(C|A) * P(D|A) * P(A) + P(B|-A) * P(C|-A) * P(D|-A) * P(-A) = \", 0.0138888888889 * 0.152777777778 * 0.0 * 0.77777778 + 0.0694444444444 * 0.180555555556 * 0.166666666667 * 0.22222222\n",
    "print \"P(D) = (see marginial probabilities of initial graph) = 0.16666667\"\n",
    "\n",
    "print \"\\nP(B,C|D) = P(B,C,D) / P(D) = \", 0.000464391856355 / 0.16666667"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "\n",
    "Somebody is cheating at playing heads and tails: he's using a biased coin (probability of heads: 0.7). However, to avoid suspicion, he also uses a fair coin (probability of heads: 0.5) sometimes, and switches between the two coins occasionally. If he has just thrown a fair coin, he will switch to the biased coin with probability 0.2 (and keep the fair coin with probability 0.8). If he's using the biased coin, he'll switch to fair with probability 0.5. When he starts playing, he chooses a coin at random (p=0.5)\n",
    "\n",
    "You're playing with this person and observe the following sequence of outcomes:\n",
    "\n",
    "$$\n",
    "HHTHHHTTTHTH\n",
    "$$\n",
    "\n",
    "**Question 2 [20 marks]** Using the above code, implement a graph corresponding to the above problem. This is an example of a Hidden Markov Model. Given the probabilities listed above, what is the probability of this sequence? For each throw, what is the probability that a fair coin was used? \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Answer Q2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": []
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}