Title: | Markov Chain Monte Carlo for Potts Models |
---|---|
Description: | Do Markov chain Monte Carlo (MCMC) simulation of Potts models (Potts, 1952, <doi:10.1017/S0305004100027419>), which are the multi-color generalization of Ising models (so, as as special case, also simulates Ising models). Use the Swendsen-Wang algorithm (Swendsen and Wang, 1987, <doi:10.1103/PhysRevLett.58.86>) so MCMC is fast. Do maximum composite likelihood estimation of parameters (Besag, 1975, <doi:10.2307/2987782>, Lindsay, 1988, <doi:10.1090/conm/080>). |
Authors: | Charles J. Geyer <[email protected]> and Leif Johnson <[email protected]> |
Maintainer: | Charles J. Geyer <[email protected]> |
License: | GPL (>= 2) |
Version: | 0.5-11 |
Built: | 2025-03-01 06:07:15 UTC |
Source: | https://github.com/cran/potts |
Variety of functions to support caching of calculated canonical statistics for Potts Models. There is some attempt at being 'smart' with when to regenerate the statistics.
generate_t_cache(x, ncolor, t_stat, sizeA, npixel, f, fapply=lapply, gridcache=NULL) gengridcache(ncolor, sizeCA, ncol) gensingleton(ncolor) singleton(x, ncolor, a, idx, gridcache=NULL) gentwopixel(ncolor) twopixel(x, ncolor, a, idx, gridcache=NULL) twopixel.nonoverlap(x, ncolor, a, idx, gridcache=NULL) genfourpixel(ncolor) fourpixel(x, ncolor, a, idx, gridcache=NULL) fourpixel.nonoverlap(x, ncolor, a, idx, gridcache=NULL) genthreebythree(ncolor) ninepixel.nonoverlap(x, ncolor, a, idx, gridcache=NULL) genfourbyfour(ncolor) sixteenpixel.nonoverlap(x, ncolor, a, idx, gridcache=NULL)
generate_t_cache(x, ncolor, t_stat, sizeA, npixel, f, fapply=lapply, gridcache=NULL) gengridcache(ncolor, sizeCA, ncol) gensingleton(ncolor) singleton(x, ncolor, a, idx, gridcache=NULL) gentwopixel(ncolor) twopixel(x, ncolor, a, idx, gridcache=NULL) twopixel.nonoverlap(x, ncolor, a, idx, gridcache=NULL) genfourpixel(ncolor) fourpixel(x, ncolor, a, idx, gridcache=NULL) fourpixel.nonoverlap(x, ncolor, a, idx, gridcache=NULL) genthreebythree(ncolor) ninepixel.nonoverlap(x, ncolor, a, idx, gridcache=NULL) genfourbyfour(ncolor) sixteenpixel.nonoverlap(x, ncolor, a, idx, gridcache=NULL)
t_stat |
numerical vector of length |
sizeA |
numerical. The number of elements in |
sizeCA |
numerical. The number of elements in |
npixel |
numerical. The number of pixels in one element of
|
f |
function. Takes arguments |
fapply |
function. It should behave exactly as lapply does. You can use this argument to enable parallel computing. |
gridcache |
list. Optional. If non-null, it is a list of the
elements of |
x |
numeric, 2 dimensional matrix, elements in 1, ..., |
ncolor |
numeric. Number of colors in this Potts Model. |
ncol |
numeric. Gives the number of columns in a rectangular window. |
a |
numeric. Indicates which member of |
idx |
numeric. Indicates which element of |
For a description of notation and terminology, see
composite.ll
.
This set of functions is used to generate cached calculations of the
canonical statistic of a Potts model suitable for passing into
composite.ll
or gr.composite.ll
.
All of the calculations using composite.ll
and these
caching functions need one of the color components to be dropped for
the model to be identifiable. For simplicity, the first color is
dropped by generate_t_cache
. In computing the composite log
likelihood for a Potts model with ncolor
colors, we are
interested in many calculations across , the set of all
permutations of colors across a window. These functions facilitate
those calculations.
gridcache
is a list of .
generate_t_cache
is the main function, and the others are
intended to be used in conjunction with it. generate_t_cache
creates a list of arrays. Each array represents one window in the
image, and each row of the array contains the value of
(with one component dropped) found by replacing the
pixels in that window with one of the elements of
.
gengridcache
can generate the gridcache
for any rectangular
window, give the number of colors, size of , and number
of columns in the window.
gensingleton
, gentwopixel
,
genfourpixel
, genthreebythree
and genfourbyfour
are all just simple wrappers for gengridcache
.
singleton
, twopixel
, twopixel.nonoverlap
,
fourpixel
, fourpixel.nonoverlap
,
ninepixel.nonoverlap
and sixteenpixel.nonoverlap
are
intended to be passed to generate_t_cache
in the argument
f
. They are used to calculate for the
element of
.
Functions that have overlap
and nonoverlap
versions
generate a overlapping and nonoverlapping set of windows respectively.
singleton
is for a single pixel window (Besag or MPLE).
twopixel
does a two horizontal pixel window.
fourpixel
does a two by two pixel window.
ninepixel
does a three by three pixel window.
sixteenpixel
does a four by four pixel window.
Functions that start with gen
return a list
of the
elements of .
The other functions (e.g. twopixel
, fourpixel
, ...)
return the result of replacing the a
-th window of x
with
the idx
-th element of and calculating
calc_t_innergrid
for that window.
ncolor <- 4 beta <- log(1+sqrt(ncolor)) theta <- c(rep(0,ncolor), beta) nrow <- 32 ncol <- 32 x <- matrix(sample(ncolor, nrow*ncol, replace=TRUE), nrow=nrow, ncol=ncol) foo <- packPotts(x, ncolor) out <- potts(foo, theta, nbatch=10) x <- unpackPotts(out$final) t_stat <- calc_t(x, ncolor) t_cache_mple <- generate_t_cache(x, ncolor, t_stat, nrow*ncol, 1, singleton) ## Not run: # use multicore to speed things up. library(multicore) t_cache_mple <- generate_t_cache(x, ncolor, t_stat, nrow*ncol, 1, singleton, fapply=mclapply) ## End(Not run)
ncolor <- 4 beta <- log(1+sqrt(ncolor)) theta <- c(rep(0,ncolor), beta) nrow <- 32 ncol <- 32 x <- matrix(sample(ncolor, nrow*ncol, replace=TRUE), nrow=nrow, ncol=ncol) foo <- packPotts(x, ncolor) out <- potts(foo, theta, nbatch=10) x <- unpackPotts(out$final) t_stat <- calc_t(x, ncolor) t_cache_mple <- generate_t_cache(x, ncolor, t_stat, nrow*ncol, 1, singleton) ## Not run: # use multicore to speed things up. library(multicore) t_cache_mple <- generate_t_cache(x, ncolor, t_stat, nrow*ncol, 1, singleton, fapply=mclapply) ## End(Not run)
Calculate the canonical statistic 't' for a realization of a Potts Model
calc_t_full(x,ncolor) calc_t_innergrid(x, ncolor, grid, i, j) calc_t(x, ncolor, grid=NULL, i=NULL, j=NULL)
calc_t_full(x,ncolor) calc_t_innergrid(x, ncolor, grid, i, j) calc_t(x, ncolor, grid=NULL, i=NULL, j=NULL)
x |
2 dimensional matrix, elements in 1, ..., |
ncolor |
numeric. Number of colors in this Potts Model. |
grid |
numeric. 2 dimensional matrix, elements in 1, ...,
|
i |
numeric. Row to place the grid. |
j |
numeric. Column to place the grid. |
For a description of notation and terminology, see
composite.ll
.
Calculates the canonical statistics for a realized Potts Model.
calc_t
calls calc_t_full
if grid
is NULL
and calc_t_innergrid
otherwise.
calc_t_full
calculates the canonical statistics for the full image.
calc_t_innergrid
calculates the canonical statistics for the a
window of the image, but with that window replaced by grid
,
with the upper left corner of grid
located at x[i,j]
.
For a description of notation and terminology, see
composite.ll
.
All functions return a vector of length ncolor+1
. Elements
1,...,ncolor
contain the number of pixels of each color.
Element ncolor+1
contains the number of matching neighbor pairs
for the image.
calc_t_full
returns the values for the whole image.
calc_t_innergrid
returns the value for just the selected
window, but this includes the number of matching pairs from the
replaced window to it's neighbors. E.g. if X
is the full
image, and is the value of some window in the image
and we want to know the value of
this would be
calc_t_full(X,
ncolor) + calc_t_innergrid(X, ncolor, y, i, j) - calc_t_innergrid(X,
ncolor, A(a), i, j)
generate_t_cache
, composite.ll
.
ncolor <- 4 beta <- log(1+sqrt(ncolor)) theta <- c(rep(0,ncolor), beta) nrow <- 32 ncol <- 32 x <- matrix(sample(ncolor, nrow*ncol, replace=TRUE), nrow=nrow, ncol=ncol) foo <- packPotts(x, ncolor) out <- potts(foo, theta, nbatch=10) x <- unpackPotts(out$final) t_stat <- calc_t(x, ncolor) t_stat_inner <- calc_t(x, ncolor, matrix(1, nrow=2, ncol=2), 1, 1)
ncolor <- 4 beta <- log(1+sqrt(ncolor)) theta <- c(rep(0,ncolor), beta) nrow <- 32 ncol <- 32 x <- matrix(sample(ncolor, nrow*ncol, replace=TRUE), nrow=nrow, ncol=ncol) foo <- packPotts(x, ncolor) out <- potts(foo, theta, nbatch=10) x <- unpackPotts(out$final) t_stat <- calc_t(x, ncolor) t_stat_inner <- calc_t(x, ncolor, matrix(1, nrow=2, ncol=2), 1, 1)
Calculate Composite Log Likelihood (CLL) and the gradient of the CLL for Potts models.
composite.ll(theta, t_stat, t_cache=NULL, fapply=lapply) gr.composite.ll(theta, t_stat, t_cache=NULL, fapply=lapply)
composite.ll(theta, t_stat, t_cache=NULL, fapply=lapply) gr.composite.ll(theta, t_stat, t_cache=NULL, fapply=lapply)
theta |
numeric canonical parameter vector. The CLL will be evaluated at this point. It is assumed that the component corresponding to the first color has been dropped. |
t_stat |
numeric, canonical statistic vector. The value of the canonical statistic for the full image. |
t_cache |
list of arrays. |
fapply |
function. Expected to function as |
For the given value of theta
composite.ll
and
gr.composite.ll
calculate the CLL and the gradient of the CLL
respectively for a realized Potts model represented by t_stat
and t_cache
.
is the set of all windows to be used in
calculating the Composite Log Likelihood (CLL) for a Potts model. A
window is a collection of adjacent pixels on the lattice of the
Potts model.
is used to represent a generic window in
, the code in this package
expects that all the windows in
have the same
size and shape.
is used to denote the size, or number
of pixels in a window. Each pixel in a Potts takes on a value in
, the set of possible colors. For simplicity, this
implementation takes
. Elements of
will be referenced
using
with
.
is used to denote all
the permutations of
across the window
, and
is used to denote the size of
.
In an abuse of notation, we use
to refer to the
element of
. No ordinal or
numerical properties of
,
or
are used, only that each element in the sets are
referenced by one and only one indexing value.
composite.ll
returns CLL evaluated at theta
.
gr.composite.ll
returns a numeric vector of length
length(theta)
containing the gradient of the CLL at theta
.
ncolor <- 4 beta <- log(1+sqrt(ncolor)) theta <- c(rep(0,ncolor), beta) nrow <- 32 ncol <- 32 x <- matrix(sample(ncolor, nrow*ncol, replace=TRUE), nrow=nrow, ncol=ncol) foo <- packPotts(x, ncolor) out <- potts(foo, theta, nbatch=10) x <- unpackPotts(out$final) t_stat <- calc_t(x, ncolor) t_cache_mple <- generate_t_cache(x, ncolor, t_stat, nrow*ncol, 1, singleton) t_cache_two <- generate_t_cache(x, ncolor, t_stat, nrow*ncol/2, 2, twopixel.nonoverlap) composite.ll(theta[-1], t_stat, t_cache_mple) gr.composite.ll(theta[-1], t_stat, t_cache_mple) ## Not run: optim.mple <- optim(theta.initial, composite.ll, gr=gr.composite.ll, t_stat, t_cache_mple, method="BFGS", control=list(fnscale=-1)) optim.mple$par optim.two <- optim(theta.initial, composite.ll, gr=gr.composite.ll, t_stat, t_cache_two, method="BFGS", control=list(fnscale=-1)) optim.two$par ## End(Not run) ## Not run: # or use mclapply to speed things up. library(multicore) optim.two <- optim(theta.initial, composite.ll, gr=gr.composite.ll, t_stat, t_cache_two, mclapply, method="BFGS", control=list(fnscale=-1)) optim.two$par ## End(Not run)
ncolor <- 4 beta <- log(1+sqrt(ncolor)) theta <- c(rep(0,ncolor), beta) nrow <- 32 ncol <- 32 x <- matrix(sample(ncolor, nrow*ncol, replace=TRUE), nrow=nrow, ncol=ncol) foo <- packPotts(x, ncolor) out <- potts(foo, theta, nbatch=10) x <- unpackPotts(out$final) t_stat <- calc_t(x, ncolor) t_cache_mple <- generate_t_cache(x, ncolor, t_stat, nrow*ncol, 1, singleton) t_cache_two <- generate_t_cache(x, ncolor, t_stat, nrow*ncol/2, 2, twopixel.nonoverlap) composite.ll(theta[-1], t_stat, t_cache_mple) gr.composite.ll(theta[-1], t_stat, t_cache_mple) ## Not run: optim.mple <- optim(theta.initial, composite.ll, gr=gr.composite.ll, t_stat, t_cache_mple, method="BFGS", control=list(fnscale=-1)) optim.mple$par optim.two <- optim(theta.initial, composite.ll, gr=gr.composite.ll, t_stat, t_cache_two, method="BFGS", control=list(fnscale=-1)) optim.two$par ## End(Not run) ## Not run: # or use mclapply to speed things up. library(multicore) optim.two <- optim(theta.initial, composite.ll, gr=gr.composite.ll, t_stat, t_cache_two, mclapply, method="BFGS", control=list(fnscale=-1)) optim.two$par ## End(Not run)
plot Potts model data.
## S3 method for class 'raw' image(x, col = c("white", "red", "blue", "green", "black", "cyan", "yellow", "magenta"), ...)
## S3 method for class 'raw' image(x, col = c("white", "red", "blue", "green", "black", "cyan", "yellow", "magenta"), ...)
x |
an R vector of class |
col |
a vector of colors. Must be as many as number of colors of Potts model. |
... |
other arguments passed to |
Too slow for large images. Needs to be rewritten for efficient plotting.
transform Potts model data from integer matrix to raw vector and vice versa.
packPotts(x, ncolor) inspectPotts(raw) unpackPotts(raw)
packPotts(x, ncolor) inspectPotts(raw) unpackPotts(raw)
x |
integer matrix containing Potts model data. Colors are coded
from one to |
ncolor |
integer scalar, number of colors. |
raw |
vector of type |
for packPotts
a vector of type "raw"
.
for inspectPotts
a list containing components
ncolor
, nrow
, and ncol
.
for unpackPotts
an integer matrix.
x <- matrix(sample(4, 2 * 3, replace = TRUE), nrow = 2) x foo <- packPotts(x, ncolor = 4) foo inspectPotts(foo) unpackPotts(foo)
x <- matrix(sample(4, 2 * 3, replace = TRUE), nrow = 2) x foo <- packPotts(x, ncolor = 4) foo inspectPotts(foo) unpackPotts(foo)
Simulate Potts model using Swendsen-Wang algorithm.
potts(obj, param, nbatch, blen = 1, nspac = 1, boundary = c("torus", "free", "condition"), debug = FALSE, outfun = NULL, ...)
potts(obj, param, nbatch, blen = 1, nspac = 1, boundary = c("torus", "free", "condition"), debug = FALSE, outfun = NULL, ...)
obj |
an R vector of class |
param |
numeric, canonical parameter vector. Last component must nonnegative (see Details below). |
nbatch |
the number of batches. |
blen |
the length of batches. |
nspac |
the spacing of iterations that contribute to batches. |
boundary |
type of boundary conditions. The value of this argument can be abbreviated. |
debug |
return additional debugging information. |
outfun |
controls the output. If a function, then the batch means
of |
... |
additional arguments for |
Runs a Swendsen-Wang algorithm producing a Markov chain with equilibrium
distribution having the specified Potts model. The state of a Potts model
is a collection of random variables taking values in a finite set. Here
the finite set is 1, ..., ncolor
and the elements are called
“colors”. The random variables are associated with the nodes of
a rectangular lattice, represented by unpackPotts
as a matrix.
In keeping with calling the values “colors”, the random variables
themselves are often called “pixels”. The probability model is an
exponential family with canonical statistic vector of length ncolor + 1
.
The first ncolor
components are the counts of the number of pixels
of each color. The last component is the number of pairs of neighboring
pixels colored the same. The corresponding canonical parameter, last
component of the canonical parameter vector (argument param
)
must be nonnegative for the Swendsen-Wang algorithm to work (Potts models
are defined for negative dependence parameter, but can't be simulated
using this algorithm).
In the default boundary specification ("torus"
), also called toroidal
or periodic boundary conditions, the vertical edges of the pixel matrix
are considered glued together, as are the horizontal edges.
Thus corresponding pixels in the first and last rows are considered neighbors,
as are corresponding pixels in the first and last columns. In the other
boundary specifications there is no such gluing: pixels in the the relative
interiors of the first and last rows and first and last columns have only
three neighbors, and the four corner pixels have only two neighbors.
In the "torus"
and "free"
boundary specifications, all pixels
are counted in determining the color count canonical statistics, which thus
range from zero to nrow * ncol
, where nrow
and ncol
are
the number of rows and columns of the pixel matrix.
In the "condition"
boundary specification, all pixels in the first
and last rows and first and last columns are fixed (conditioned on), and only
the random pixels are counted in determining the color count canonical
statistics, which thus range from zero to (nrow - 2) * (ncol - 2)
.
In the "torus"
boundary specification, all pixels have four neighbors,
so the neighbor pair canonical statistic ranges from zero
to 2 * nrow * ncol
.
In the "free"
boundary specification, pixels in the interior have four
neighbors, those in the relative interior of edges have three, and those
in the corners have two, so the neighbor pair canonical statistic ranges from
zero to nrow * (ncol - 1) + (nrow - 1) * ncol
.
In the "condition"
boundary specification, only neighbor pairs in which
at least one pixel is random are counted, so the neighbor pair canonical
statistic ranges from zero
to (nrow - 2) * (ncol - 1) + (nrow - 1) * (ncol - 2)
.
an object of class "potts"
,
which is a list containing at least the following components:
initial |
initial state of Markov chain in the format output
by |
final |
final state of Markov chain in the same format. |
initial.seed |
value of |
final.seed |
value of |
time |
running time of Markov chain from |
param |
canonical parameter vector. |
nbatch |
the number of batches. |
blen |
the length of batches. |
nspac |
the spacing of iterations that contribute to batches. |
boundary |
the argument |
batch |
an |
ncolor <- as.integer(4) beta <- log(1 + sqrt(ncolor)) theta <- c(rep(0, ncolor), beta) nrow <- 100 ncol <- 100 x <- matrix(1, nrow = nrow, ncol = ncol) foo <- packPotts(x, ncolor) out <- potts(foo, theta, nbatch = 10) out$batch ## Not run: image(out$final)
ncolor <- as.integer(4) beta <- log(1 + sqrt(ncolor)) theta <- c(rep(0, ncolor), beta) nrow <- 100 ncol <- 100 x <- matrix(1, nrow = nrow, ncol = ncol) foo <- packPotts(x, ncolor) out <- potts(foo, theta, nbatch = 10) out$batch ## Not run: image(out$final)