Revision 69b0f9dca8eb051f132725ecc679fe1997246e50 authored by Adrian Baddeley on 18 January 2006, 21:47:25 UTC, committed by cran-robot on 18 January 2006, 21:47:25 UTC
1 parent cb2215f
multistrhard.S
#
#
# multistrhard.S
#
# $Revision: 2.8 $ $Date: 2005/07/27 06:26:01 $
#
# The multitype Strauss/hardcore process
#
# MultiStraussHard()
# create an instance of the multitype Strauss/ harcore
# point process
# [an object of class 'interact']
#
# -------------------------------------------------------------------
#
MultiStraussHard <- function(types, iradii, hradii) {
if(length(types) == 1)
stop("The \`types\' argument should be a vector of all possible types")
if(is.factor(types)) {
types <- levels(types)
} else {
types <- levels(factor(types, levels=types))
}
dimnames(iradii) <- dimnames(hradii) <- list(types, types)
out <-
list(
name = "Multitype Strauss Hardcore process",
family = pairwise.family,
pot = function(d, tx, tu, par) {
# arguments:
# d[i,j] distance between points X[i] and U[j]
# tx[i] type (mark) of point X[i]
# tu[i] type (mark) of point U[j]
#
# get matrices of interaction radii
r <- par$iradii
h <- par$hradii
# get possible marks and validate
if(!is.factor(tx) || !is.factor(tu))
stop("marks of data and dummy points must be factor variables")
lx <- levels(tx)
lu <- levels(tu)
if(length(lx) != length(lu) || any(lx != lu))
stop("marks of data and dummy points do not have same possible levels")
if(!identical(lx, par$types))
stop("data and model do not have the same possible levels of marks")
if(!identical(lu, par$types))
stop("dummy points and model do not have the same possible levels of marks")
# list all UNORDERED pairs of types to be checked
# (the interaction must be symmetric in type, and scored as such)
uptri <- (row(r) <= col(r)) & (!is.na(r) | !is.na(h))
mark1 <- (lx[row(r)])[uptri]
mark2 <- (lx[col(r)])[uptri]
vname <- apply(cbind(mark1,mark2), 1, paste, collapse="x")
vname <- paste("mark", vname, sep="")
npairs <- length(vname)
# list all ORDERED pairs of types to be checked
# (to save writing the same code twice)
different <- mark1 != mark2
mark1o <- c(mark1, mark2[different])
mark2o <- c(mark2, mark1[different])
nordpairs <- length(mark1o)
# unordered pair corresponding to each ordered pair
ucode <- c(1:npairs, (1:npairs)[different])
#
# go....
# apply the relevant interaction distance to each pair of points
rxu <- r[ tx, tu ]
str <- (d < rxu)
str[is.na(str)] <- FALSE
# and the relevant hard core distance
hxu <- h[ tx, tu ]
forbid <- (d < hxu)
forbid[is.na(forbid)] <- FALSE
# form the potential
value <- ifelse(forbid, -Inf, str)
# create numeric array for result
z <- array(0, dim=c(dim(d), npairs),
dimnames=list(character(0), character(0), vname))
# assign value[i,j] -> z[i,j,k] where k is relevant interaction code
for(i in 1:nordpairs) {
# data points with mark m1
Xsub <- (tx == mark1o[i])
# quadrature points with mark m2
Qsub <- (tu == mark2o[i])
# assign
z[Xsub, Qsub, ucode[i]] <- value[Xsub, Qsub]
}
return(z)
},
#### end of 'pot' function ####
#
par = list(types=types, iradii = iradii, hradii = hradii),
parnames = c("possible types", "interaction distances", "hardcore distances"),
init = function(self) {
r <- self$par$iradii
h <- self$par$hradii
nt <- length(self$par$types)
MultiPair.checkmatrix(r, nt, "\`iradii\'")
MultiPair.checkmatrix(h, nt, "\`hradii\'")
ina <- is.na(iradii)
hna <- is.na(hradii)
if(all(ina))
stop("All entries of \`iradii\' are NA")
both <- !ina & !hna
if(any(iradii[both] <= hradii[both]))
stop("iradii must be larger than hradii")
},
update = NULL, # default OK
print = function(self) {
print.isf(self$family)
cat(paste("Interaction:\t", self$name, "\n"))
cat(paste(length(self$par$types), "types of points\n"))
cat("Possible types: \n")
print(self$par$types)
cat("Interaction radii:\n")
print(self$par$iradii)
cat("Hardcore radii:\n")
print(self$par$hradii)
invisible()
},
interpret = function(coeffs, self) {
# get possible types
typ <- self$par$types
ntypes <- length(typ)
# get matrices of interaction radii
r <- self$par$iradii
h <- self$par$hradii
# list all unordered pairs of types
uptri <- (row(r) <= col(r)) & (!is.na(r) | !is.na(h))
mark1 <- (typ[row(r)])[uptri]
mark2 <- (typ[col(r)])[uptri]
# names of coefficients
basename <- apply(cbind(mark1,mark2), 1, paste, collapse="x")
thetaname <- paste("mark", basename, sep="")
npairs <- length(basename)
# extract canonical parameters; shape them into a matrix
gammas <- matrix(, ntypes, ntypes)
dimnames(gammas) <- list(typ, typ)
for(i in 1:npairs) {
theta <- coeffs[thetaname[i]]
gamma <- exp(theta)
gammas[mark1[i], mark2[i]] <- gamma
gammas[mark2[i], mark1[i]] <- gamma
}
#
return(list(param=list(gammas=gammas),
inames="interaction parameters gamma_ij",
printable=round(gammas,4)))
},
valid = function(coeffs, self) {
# interaction radii r[i,j]
iradii <- self$par$iradii
# hard core radii r[i,j]
hradii <- self$par$hradii
# interaction parameters gamma[i,j]
gamma <- (self$interpret)(coeffs, self)$param$gammas
# Check that we managed to estimate all required parameters
required <- !is.na(iradii)
if(!all(is.finite(gamma[required])))
return(FALSE)
# Check that the model is integrable
# inactive hard cores ...
ihc <- (is.na(hradii) | hradii == 0)
# .. must have gamma <= 1
return(all(gamma[required & ihc] <= 1))
},
project = function(coeffs, self) {
# interaction radii r[i,j]
iradii <- self$par$iradii
# hard core radii r[i,j]
hradii <- self$par$hradii
# interaction parameters gamma[i,j]
gamma <- (self$interpret)(coeffs, self)$param$gammas
# remove NA's
gamma[is.na(gamma)] <- 1
# inactive hard cores?
ihc <- (is.na(hradii) | hradii == 0)
if(any(ihc & (gamma > 1)))
gamma[ihc] <- pmin(gamma[ihc], 1)
# now put them back... %^[
# get possible types
typ <- self$par$types
ntypes <- length(typ)
# get matrices of interaction radii
r <- self$par$iradii
h <- self$par$hradii
# list all unordered pairs of types
uptri <- (row(r) <= col(r)) & (!is.na(r) | !is.na(h))
mark1 <- (typ[row(r)])[uptri]
mark2 <- (typ[col(r)])[uptri]
# names of coefficients
basename <- apply(cbind(mark1,mark2), 1, paste, collapse="x")
thetaname <- paste("mark", basename, sep="")
npairs <- length(basename)
# reassign
coeffs[thetaname] <- log(gamma[uptri])
return(coeffs)
},
irange = function(self, coeffs=NA, epsilon=0, ...) {
r <- self$par$iradii
h <- self$par$hradii
ractive <- !is.na(r)
hactive <- !is.na(h)
if(any(!is.na(coeffs))) {
gamma <- (self$interpret)(coeffs, self)$param$gammas
gamma[is.na(gamma)] <- 1
ractive <- ractive & (abs(log(gamma)) > epsilon)
}
if(!any(c(ractive,hactive)))
return(0)
else
return(max(c(r[ractive],h[hactive])))
}
)
class(out) <- "interact"
out$init(out)
return(out)
}
Computing file changes ...