Generate Random Samples from a Potts Model Using the Checkerboard Idea

Description

Generate random samples from a Potts model by Gibbs Sampling that takes advantage of conditional independence.

Usage

  BlocksGibbs(n, nvertex, ncolor, neighbors, blocks,
              weights=1, spatialMat=NULL, beta)

Arguments

Details

We use the Gibbs algorithm that takes advantage of conditional independence to speed up the generation of random samples from a Potts model. The idea is that if we can divide variables that need to be updated into different blocks and given the variables in other blocks, all the variables within the same block are conditionally independent, then we can update all blocks iteratively with the variables within the same block being updated simultaneously.

The spatialMat is the argument used to specify the relationship among vertices in neighbor. See rPotts1 for more information on the Potts model and spatialMat.

Value

The output is a nvertex by n matrix with the kth column being the kth sample.

References

Dai Feng (2008) Bayesian Hidden Markov Normal Mixture Models with Application to MRI Tissue Classification Ph. D. Dissertation, The University of Iowa

See Also

Wolff, SW

Examples

  #Example 1: Generate 100 samples from a repulsion Potts model with the
  #           neighborhood structure corresponding to a first-order
  #           Markov random field defined on a 3*3 2D graph.
  #           The number of colors is 3 and beta=0.1,a_1=2,a_2=1,a_3=0.
  #           All weights are equal to 1.
 
  neighbors <- getNeighbors(mask=matrix(1, 3, 3), neiStruc=c(2,2,0,0))
  blocks <- getBlocks(mask=matrix(1, 3, 3), nblock=2)
  spatialMat <- matrix(c(2,1,0, 1,2,1,0,1,2), ncol=3)
  BlocksGibbs(n=100, nvertex=9, ncolor=3, neighbors=neighbors, blocks=blocks,
              spatialMat=spatialMat, beta=0.1)

Generate Random Samples from a Compound Potts Model by the Swendsen-Wang Algorithm

Description

Generate random samples from a compound Potts model using the Swendsen-Wang algorithm.

Usage

  SW(n, nvertex, ncolor, edges, weights, beta)

Arguments

Details

We use the Swendsen-Wang algorithm to generate random samples from a compound Potts model. See rPotts1 for more information on the compound Potts model.

Value

The output is a nvertex by n matrix with the kth column being the kth sample.

References

Robert H. Swendsen and Jian-Sheng Wang (1987) Nonuniversal Critical Dynamics in Monte Carlo Simulations Physical Review Letters vol. 58, no. 2, 86-88

Dai Feng (2008) Bayesian Hidden Markov Normal Mixture Models with Application to MRI Tissue Classification Ph. D. Dissertation, The University of Iowa

See Also

Wolff, BlocksGibbs

Examples

  #Example 1: Generate 100 samples from a Potts model with the
  #           neighborhood structure corresponding to a
  #           second-order Markov random field defined on a 
  #           3*3 2D graph. The number of colors is 2.
  #           beta=0.1. All weights are equal to 1.
 
  edges <- getEdges(mask=matrix(1, 2, 2), neiStruc=rep(2,4))
  set.seed(100)
  SW(n=500, nvertex=4, ncolor=2, edges, beta=0.8)

Generate Random Samples from a Compound Potts Model by the Wolff Algorithm

Description

Generate random samples from a compound Potts model using the Wolff Algorithm.

Usage

  Wolff(n, nvertex, ncolor, neighbors, weights, beta)

Arguments

Details

We use the Wolff algorithm to generate random samples from a compound Potts model. See rPotts1 for more information on the compound Potts model.

Value

A nvertex by n matrix with the kth column being the kth sample.

References

Ulli Wolff (1989) Collective Monte Carlo Updating for Spin Systems Physical Review Letters vol. 62, no. 4, 361-364

Dai Feng (2008) Bayesian Hidden Markov Normal Mixture Models with Application to MRI Tissue Classification Ph. D. Dissertation, The University of Iowa

See Also

SW, BlocksGibbs

Examples

  #Example 1: Generate 100 samples from a Potts model with the
  #           neighborhood structure corresponding to a
  #           second-order Markov random field defined on a 
  #           3*3 2D graph. The number of colors is 2.
  #           beta=0.7. All weights are equal to 1.

  neighbors <- getNeighbors(mask=matrix(1, 3, 3), neiStruc=rep(2,4))
  Wolff(n=100, nvertex=9, ncolor=2, neighbors, beta=0.7)

Get Blocks of a Graph

Description

Obtain blocks of vertices of a 1D, 2D, or 3D graph, in order to use the conditional independence to speed up the simulation (checkerboard idea).

Usage

   getBlocks(mask, nblock)

Arguments

Details

The vertices within each block are mutually independent given the vertices in other blocks. Some blocks could be empty.

Value

A list with the number of components equal to nblock. Each component consists of vertices within the same block.

References

Darren J. Wilkinson Parallel Bayesian Computation Handbook of Parallel Computing and Statistics 481-512 Marcel Dekker/CRC Press 2005

Examples

  #Example 1: split a line into 2 blocks
  getBlocks(mask=c(1,1,1,1,0,0,1,1,0), nblock=2)

  #Example 2: split a 4*4 2D graph into 4 blocks in order
  #           to use the checkerboard idea for a neighborhood structure
  #           corresponding to the second-order Markov random field.
  getBlocks(mask=matrix(1, nrow=4, ncol=4), nblock=4)

  #Example 3: split a 3*3*3 3D graph into 8 blocks
  #           in order to use the checkerboard idea for a neighborhood
  #           structure based on the 18 neighbors definition, where the
  #           neighbors of a vertex comprise its available
  #           adjacencies sharing the same edges or faces.
  mask <- array(1, dim=rep(3,3))
  getBlocks(mask, nblock=8)

Generate Configurations of a Graph

Description

Using recursive method to generate all possible configurations of a graph.

Usage

   getConfs(nvertex, ncolor)

Arguments

Details

Suppose there are n vertices and each can take values from 1,2, \ldots, ncolor. This function generates all possible configurations. For example, if there are two vertices and each can be either 1 or 2, then the possible configurations are (1,1), (1,2), (2,1) and (2,2).

Value

A matrix of all possible configurations. Each column corresponds to one configuration.

Examples

  #Example 1: There are two vertices and each is either of
  #           color 1 or 2.
  getConfs(2,2)

Get Edges of a Graph

Description

Obtain edges of a 1D, 2D, or 3D graph based on the neighborhood structure.

Usage

   getEdges(mask, neiStruc)

Arguments

Details

There could be more than one way to define the same 3D neighborhood structure for a graph (see Example 4 for illustration).

Value

A matrix of two columns with one edge per row. The edges connecting vertices and their corresponding first neighbors are listed first, and then those corresponding to the second neighbors, and so on and so forth. The order of neighbors is the same as in getNeighbors.

References

Gerhard Winkler (1995) Image Analysis, Random Fields and Dynamic Monte Carlo Methods Springer-Verlag

Dai Feng (2008) Bayesian Hidden Markov Normal Mixture Models with Application to MRI Tissue Classification Ph. D. Dissertation, The University of Iowa

Examples

  #Example 1: get all edges of a 1D graph. 
  mask <- c(0,0,rep(1,4),0,1,1,0,0)
  getEdges(mask, neiStruc=2)

  #Example 2: get all edges of a 2D graph based on neighborhood structure
  #           corresponding to the first-order Markov random field.
  mask <- matrix(1 ,nrow=2, ncol=3)
  getEdges(mask, neiStruc=c(2,2,0,0))

  #Example 3: get all edges of a 2D graph based on neighborhood structure
  #           corresponding to the second-order Markov random field.
  mask <- matrix(1 ,nrow=3, ncol=3)
  getEdges(mask, neiStruc=c(2,2,2,2))

  #Example 4: get all edges of a 3D graph based on 6 neighbors structure
  #           where the neighbors of a vertex comprise its available
  #           N,S,E,W, upper and lower adjacencies. To achieve it, there
  #           are several ways, including the two below.
  mask <- array(1, dim=rep(3,3))
  n61 <- matrix(c(2,2,0,0,
                  0,2,0,0,
                  0,0,0,0), nrow=3, byrow=TRUE)
  n62 <- matrix(c(2,0,0,0,
                  0,2,0,0,
                  2,0,0,0), nrow=3, byrow=TRUE)
  e1 <- getEdges(mask, neiStruc=n61)
  e2 <- getEdges(mask, neiStruc=n62)
  e1 <- e1[order(e1[,1], e1[,2]),]
  e2 <- e2[order(e2[,1], e2[,2]),]
  all(e1==e2)

  #Example 5: get all edges of a 3D graph based on 18 neighbors structure
  #           where the neighbors of a vertex comprise its available
  #           adjacencies sharing the same edges or faces.
  #           To achieve it, there are several ways, including the one below.

  n18 <- matrix(c(2,2,2,2,
                  0,2,2,2,
                  0,0,2,2), nrow=3, byrow=TRUE)  
  mask <- array(1, dim=rep(3,3))
  getEdges(mask, neiStruc=n18)
  
  

Calculate the Normalizing Constant of a Simple Potts Model

Description

Use the thermodynamic integration approach to calculate the normalizing constant of a Simple Potts Model.

Usage

  getNC(beta, subbetas, nvertex, ncolor,
        edges, neighbors=NULL, blocks=NULL, 
        algorithm=c("SwendsenWang", "Gibbs", "Wolff"), n, burn)

Arguments

Details

Use the thermodynamic integration approach to calculate the normalizing constant from a simple Potts model. See rPotts1 for more information on the simple Potts model.

By the thermodynamic integration method,

\log{C(\beta)} = N\log{k} + \int_{0}^{\beta}E(U({\bf z})|\beta^{'}, k)d\beta^{'}

where N is the total number of vertices (nvertex), k is the number of colors (ncolor), and U({\bf z}) = \sum_{i \sim j}\textrm{I}(z_{i}=z_{j}). Calculate E(U({\bf z}) for subbetas based on samples, and then compute the integral by numerical integration.

Value

The corresponding normalizing constant.

References

Peter J. Green and Sylvia Richardson (2002) Hidden Markov Models and Disease Mapping Journal of the American Statistical Association vol. 97, no. 460, 1055-1070

See Also

BlocksGibbs, SW, Wolff

Examples

  ## Not run: 
  #Example 1: Calculate the normalizing constant of a simple Potts model
  #           with the neighborhood structure corresponding to a
  #           first-order Markov random field defined on a
  #           3*3 2D graph. The number of colors is 2 and beta=2.
  #           Use 11 subbetas evenly distributed between 0 and 2.
  #           The sampling algorithm is Swendsen-Wang with 10000
  #           iterations and 1000 burn-in. 
 
  edges <- getEdges(mask=matrix(1,3,3), neiStruc=c(2,2,0,0))
  getNC(beta=2, subbetas=seq(0,2,by=0.2), nvertex=3*3, ncolor=2,
        edges, algorithm="S", n=10000, burn=1000)
  
## End(Not run)

Get Neighbors of All Vertices of a Graph

Description

Obtain neighbors of vertices of a 1D, 2D, or 3D graph.

Usage

   getNeighbors(mask, neiStruc)

Arguments

Details

There could be more than one way to define the same 3D neighborhood structure for a graph (see Example 3 for illustration).

Value

A matrix with each row giving the neighbors of a vertex. The number of the rows is equal to the number of vertices within the graph and the number or columns is the number of neighbors of each vertex.

For a 1D graph, if each vertex has two neighbors, The first column are the neighbors on the left-hand side of corresponding vertices and the second column the right-hand side. For the vertices on boundaries, missing neighbors are represented by the number of vertices within a graph plus 1. When neiStruc is bigger than 2, The first two columns are the same as when neiStruc is equal to 2; the third column are the neighbors on the left-hand side of the vertices on the first column; the forth column are the neighbors on the right-hand side of the vertices on the second column, and so on and so forth. And again for the vertices on boundaries, their missing neighbors are represented by the number of vertices within a graph plus 1.

For a 2D graph, the index to vertices is column-wised. For each vertex, the order of neighbors are as follows. First are those on the vertical direction, second the horizontal direction, third the NW to SE diagonal direction, and forth the SW to NE diagonal direction. For each direction, the neighbors of every vertex are arranged in the same way as in a 1D graph.

For a 3D graph, the index to vertices is that the leftmost subscript of the array moves the fastest. For each vertex, the neighbors from the 1-2 perspective appear first and then the 1-3 perspective and finally the 2-3 perspective. For each perspective, the neighbors are arranged in the same way as in a 2D graph.

References

Gerhard Winkler (1995) Image Analysis, Random Fields and Dynamic Monte Carlo Methods Springer-Verlag

Dai Feng (2008) Bayesian Hidden Markov Normal Mixture Models with Application to MRI Tissue Classification Ph. D. Dissertation, The University of Iowa

Examples

  #Example 1: get all neighbors of a 1D graph.
  mask <- c(0,0,rep(1,4),0,1,1,0,0,1,1,1)
  getNeighbors(mask, neiStruc=2)

  #Example 2: get all neighbors of a 2D graph based on neighborhood structure
  #           corresponding to the second-order Markov random field.
  mask <- matrix(1, nrow=2, ncol=3)
  getNeighbors(mask, neiStruc=c(2,2,2,2))

  #Example 3: get all neighbors of a 3D graph based on 6 neighbors structure
  #           where the neighbors of a vertex comprise its available
  #           N,S,E,W, upper and lower adjacencies. To achieve it, there
  #           are several ways, including the two below.
  mask <- array(1, dim=rep(3,3))
  n61 <- matrix(c(2,2,0,0,
                  0,2,0,0,
                  0,0,0,0), nrow=3, byrow=TRUE)
  n62 <- matrix(c(2,0,0,0,
                  0,2,0,0,
                  2,0,0,0), nrow=3, byrow=TRUE)
  n1 <- getNeighbors(mask, neiStruc=n61)
  n2 <- getNeighbors(mask, neiStruc=n62)
  n1 <- apply(n1, 1, sort)
  n2 <- apply(n2, 1, sort)
  all(n1==n2)

  #Example 4: get all neighbors of a 3D graph based on 18 neighbors structure
  #           where the neighbors of a vertex comprise its available
  #           adjacencies sharing the same edges or faces.
  #           To achieve it, there are several ways, including the one below.

  n18 <- matrix(c(2,2,2,2,
                  0,2,2,2,
                  0,0,2,2), nrow=3, byrow=TRUE)  
  mask <- array(1, dim=rep(3,3))
  getNeighbors(mask, neiStruc=n18)
  
  

Get Patches of a Graph

Description

Obtain patches of a graph by Rem's algorithm.

Usage

   getPatches(bonds, nvertex)

Arguments

Details

Given all bonds and the number of vertices in a graph, this function provides all patches.

Value

A list comprises all patches in a graph. Each component of the list consists of vertices within one patch.

References

Edsger W. Dijkstra (1976) A Discipline of Programming Englewood Cliffs, New Jersey : Prentice-Hall, Inc

Examples

  #Example 1: Find patches of a 3*3 2D graph with 6 bonds.

  bonds <- matrix(c(1,2,2,5,5,6,3,6,5,8,7,8), ncol=2, byrow=TRUE)
  getPatches(bonds, 9)

Get All Weights of a Graph

Description

Obtain weights of edges of a 1D, 2D, or 3D graph based on the neighborhood structure.

Usage

   getWeights(mask, neiStruc, format=1)

Arguments

Details

The weights are equal to the reciprocals of the distance between neighboring vertices.

Value

A vector of weights, one component corresponding to an edge of a graph. Or a matrix of weights, one component corresponding to two vertices in neighbor.

Examples

  #Example 1: get all weights of a 2D graph based on neighborhood structure
  #           corresponding to the first-order Markov random field.
  mask <- matrix(1 ,nrow=2, ncol=3)
  getWeights(mask, neiStruc=c(2,2,0,0))

  #Example 2: get all weights of a 2D graph based on neighborhood structure
  #           corresponding to the second-order Markov random field.
  #           Put the weights in a matrix form corresponding to
  #           neighbors of vertices.
  mask <- matrix(1 ,nrow=3, ncol=3)
  getWeights(mask, neiStruc=c(2,2,2,2), format=2)

  #Example 3: get all weights of a 3D graph based on 6 neighbors structure
  #           where the neighbors of a vertex comprise its available
  #           N,S,E,W, upper and lower adjacencies. 
  mask <- array(1, dim=rep(3,3))
  n61 <- matrix(c(2,2,0,0,
                  0,2,0,0,
                  0,0,0,0), nrow=3, byrow=TRUE)
  getWeights(mask, neiStruc=n61)
 

Generate One Random Sample from a Potts Model

Description

Generate one random sample from a Potts model with external field by Gibbs Sampling that takes advantage of conditional independence, or the partial decoupling method.

Usage

  rPotts1(nvertex, ncolor, neighbors, blocks, edges=NULL, weights=1,
          spatialMat=NULL, beta, external, colors,
          algorithm=c("Gibbs", "PartialDecoupling"))

Arguments

Details

This function generates random samples from a Potts model as follows:

p({\bf z})=C(\beta)^{-1} \exp\{\sum_i \alpha_i(z_i) + \beta \sum_{i \sim j} w_{ij} f(z_{i},z_{j})\}

where C(\beta) is a normalizing constant and i \sim j indicates neighboring vertices. The parameter \beta is called the "inverse temperature", which determines the level of spatial homogeneity between neighboring vertices in the graph. We assume \beta > 0. The set {\bf z}=\{z_{1}, z_{2},\ldots,\} comprises the indices to the colors of all vertices. Function f(z_{i}, z_{j}) determines the relationship among vertices in neighbor. Parameter w_{ij} is the weight between vertex i and j. The term \sum_i \alpha_i(z_i) is called the "external field".

For the simple, the compound, and the simple repulsion Potts models, the external field is equal to 0. For the simple and the compound Potts model f(z_{i}, z_{j}) = I(z_{i}=z_{j}). Parameters w_{ij} are all equal for the simple Potts model but not so for the compound model.

For the repulsion Potts model f(z_i , z_j) = \beta_1 if z_i = z_j; f(z_i , z_j) = \beta_2 if |z_i - z_j| = 1; f(z_i , z_j) = \beta_3 otherwise.

The argument spatialMat is used to specify the relationship among vertices in neighbor. The default value is NULL corresponding to the simple or the compound Potts model. The component at the ith row and jth column defining the relationship when the color of a vertex is i and the color of its neighbors is j. Besides the default setup, for the simple and the compound Potts models spatailMat could be an identity matrix also. For the repulsion Potts model, it is

\left(\begin{array}{ccccc} a_1 & a_2 & a_3 & \ldots & a_3 \\ a_2 & a_1 & a_2 & \ldots & a_3 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ a_3 & a_3 & a_3 & \ldots & a_1 \end{array}\right)

Other relationships among neighboring vertices can be specified through it as well.

Gibbs sampling can be used to generate samples from all kinds of Potts models. We use the method that takes advantage of conditional independence to speed up the simulation. See BlocksGibbs for details.

The partial decoupling method could be used to generate samples from the simple Potts model plus the external field. The \delta_{ij}s are specified through the argument weights.

Value

The output is a vector with the kth component being the new color of vertex k.

References

Dai Feng (2008) Bayesian Hidden Markov Normal Mixture Models with Application to MRI Tissue Classification Ph. D. Dissertation, The University of Iowa

David M. Higdon (1998) Auxiliary variable methods for Markov Chain Monte Carlo with applications Journal of the American Statistical Association vol. 93 585-595

See Also

BlocksGibbs, Wolff SW

Examples

  ## Not run: 
  neighbors <- getNeighbors(matrix(1, 16, 16), c(2,2,0,0))
  blocks <- getBlocks(matrix(1, 16, 16), 2)
  spatialMat <- matrix(c(2, 0, -1, 0, 2, 0, -1, 0, 2), ncol=3)
  mu <- c(22, 70 ,102)
  sigma <- c(17, 16, 19)
  count <- c(40, 140, 76)
  y <- unlist(lapply(1:3, function(i) rnorm(count[i], mu[i], sigma[i])))
  external <- do.call(cbind,
                      lapply(1:3, function(i) dnorm(y, mu[i],sigma[i])))
  current.colors <- rep(1:3, count)
  rPotts1(nvertex=16^2, ncolor=3, neighbors=neighbors, blocks=blocks,  
          spatialMat=spatialMat, beta=0.3, external=external,
          colors=current.colors, algorithm="G")
  edges <- getEdges(matrix(1, 16, 16), c(2,2,0,0))
  rPotts1(nvertex=16^2, ncolor=3, edges=edges, beta=0.3,
          external=external, colors=current.colors, algorithm="P")
 
## End(Not run)