gistfile1.r

library(CRF)

# there are 10 nodes.  Each node can take on one of 3 states.
nNodes <- 10
nStates <- 3

# make the adjacency matrices for each graph

# a fully-connected graph of 10 nodes.  Loopiest possible.
# we expect a lot of error: the inferred node beliefs should
# be pretty different from the exact result.
adj_complete <- matrix(1, nrow=nNodes, ncol=nNodes)

# a chain graph, but with a few extra loops added in (draw the graph!)
# we expect a little error relative to the exact answer
adj_loopy <- matrix(0, nrow=nNodes, ncol=nNodes)
for (i in 1:(nNodes-3))
{
  adj_loopy[i,i+1] <- 1
  adj_loopy[i+1,i] <- 1
  adj_loopy[i,i+2] <- 1
  adj_loopy[i+2,i] <- 1
  adj_loopy[ i,i+3] <- 1
  adj_loopy[i+3,i] <- 1
}

# a chain graph.  this is a special case of a tree graph, and so 
# belief prop should be exact.
adj_chain <- matrix(0, nrow=nNodes, ncol=nNodes)
for (i in 1:(nNodes-1))
{
  adj_chain[i,i+1] <- 1
  adj_chain[i+1,i] <- 1
}

# now we have to actually build our conditional random fields
# from the adjacency matrices.

# first, we build the crf object
crf_chain <- make.crf(adj_chain, nStates)

# now we set the node potentials.  I'm 90% sure these are the 
# \phi_i( x_i ) from equation (12) in 
# http://www.merl.com/publications/docs/TR2001-22.pdf

crf_chain$node.pot[1,] <- c(1, 3,4)
crf_chain$node.pot[2,] <- c(9, 1,3)
crf_chain$node.pot[3,] <- c(1, 3,5)
crf_chain$node.pot[4,] <- c(9, 1,1)

# next, we set the edge potentials.  I'm 90% sure these are the
# \psi_{ij}( x_i, x_j ) from the paper, 
# http://www.merl.com/publications/docs/TR2001-22.pdf
#  that are used in equation (12).  The x_i are the values
# of the states (1,2, or 3), and the i, j are the indices
# for which edge of the adj matrix you're talking about.  For
# each edge, there's a (nStates, nStates ) matrix, since the x_i
# and x_j are discrete random variables with nStates values. We're
# only setting 2 of the 3 rows, leaving the last implicitly all 0s.
# Wei had a good question: should the edge potential matrix actually 
# be symmetric, since the graph is undirected?

for (i in 1:crf_chain$n.edges)
{
  crf_chain$edge.pot[[i]][1,] <- c(2, 1,1)
  crf_chain$edge.pot[[i]][2,] <- c(1, 2,1)
}

# do the same for the loopy graph
crf_loopy<- make.crf(adj_loopy, nStates)
crf_loopy$node.pot[1,] <- c(1, 3,4)
crf_loopy$node.pot[2,] <- c(9, 1,3)
crf_loopy$node.pot[3,] <- c(1, 3,5)
crf_loopy$node.pot[4,] <- c(9, 1,1)

for (i in 1:crf_loopy$n.edges)
{
  crf_loopy$edge.pot[[i]][1,] <- c(2, 1,1)
  crf_loopy$edge.pot[[i]][2,] <- c(1, 2,1)
}

# and once more for the complete graph
crf_complete <- make.crf(adj_complete, nStates)
crf_complete$node.pot[1,] <- c(1, 3,4)
crf_complete$node.pot[2,] <- c(9, 1,3)
crf_complete$node.pot[3,] <- c(1, 3,5)
crf_complete$node.pot[4,] <- c(9, 1,1)

for (i in 1:crf_complete$n.edges)
{
  crf_complete$edge.pot[[i]][1,] <- c(2, 1,1)
  crf_complete$edge.pot[[i]][2,] <- c(1, 2,1)
}

# now, we do the inference!  

# this prints the beliefs for the exact answer
i <- infer.exact(crf_chain)
i$node.bel
# and this prints the beliefs from loopy belief propagation
i <- infer.lbp( crf_chain )
i$node.bel
# since that was the chain graph, the answers should be very
# close: BP is exact for tree graphs.

# now, we'll get more error in the loopy graph:
i <- infer.exact(crf_loopy)
i$node.bel
i <- infer.lbp( crf_loopy )
i$node.bel

# and finally, the complete graph should be way off.
i <- infer.exact(crf_complete)
i$node.bel
i <- infer.lbp( crf_complete )
i$node.bel