composite.ll: Composite Log Likelihood for Potts Models

Description

Calculate Composite Log Likelihood (CLL) and the gradient of the CLL for Potts models.

Usage

composite.ll(theta, t_stat, t_cache=NULL, fapply=lapply)
gr.composite.ll(theta, t_stat, t_cache=NULL, fapply=lapply)

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.

Arguments

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. t_cache[[i]][j,] = the value of \(t\) with window \(A_i\) replaced by the \(j^{th}\) element of \(C^A\).
fapply: function. Expected to function as lapply does. Useful for enabling parallel processing. E.g. use the mclapply function from the multicore package.

Details

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.

\(\mathcal{A}\) 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. \(A\) is used to represent a generic window in \(\mathcal{A}\), the code in this package expects that all the windows in \(\mathcal{A}\) have the same size and shape. \(|A|\) is used to denote the size, or number of pixels in a window. Each pixel in a Potts takes on a value in \(C\), the set of possible colors. For simplicity, this implementation takes \(C = \{1,\dots,\texttt{ncolor}\}\). Elements of \(C\) will be referenced using \(c_j\) with \(j \in \{1,\dots,\texttt{ncolor}\}\). \(C^A\) is used to denote all the permutations of \(C\) across the window \(A\), and \(|C|^{|A|}\) is used to denote the size of \(C^A\). In an abuse of notation, we use \(A_a\) to refer to the \(a^{th}\) element of \(\mathcal{A}\). No ordinal or numerical properties of \(\mathcal{A}\), \(C\) or \(C^A\) are used, only that each element in the sets are referenced by one and only one indexing value.

Examples

Run this code


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)

if (FALSE) {
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
}

if (FALSE) {
# 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

}