R code for chapter 13
Below is the R code used to generate results and figures in chapter 13. The code is in a rather raw format extracted from the book manuscript files – please read the instructions for use if you haven’t done so yet. You can download the script here.
### R code from vignette source '13gibbs'
## Copyright (C) Adrian Baddeley, Ege Rubak and Rolf Turner
###################################################
### code chunk number 1: 13gibbs.Rnw:9-11
###################################################
source("R/startup.R")
requireversion(spatstat, "1.41-1.073")
###################################################
### code chunk number 2: HybridSvensk
###################################################
set.seed(191711)
fitH <- ppm(swedishpines~polynom(x,y,2),
Hybrid(DiggleGatesStibbard(10), AreaInter(8)))
simH <- simulate(fitH, nsim=1)[[1]]
###################################################
### code chunk number 3: Unit2.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0)
###################################################
### code chunk number 4: 13gibbs.Rnw:135-138
###################################################
plot(solist(swedishpines, simH),
main="", main.panel="", nrows=1,
equal.scales=TRUE, mar.panel=0, hsep=1, pch=16)
###################################################
### code chunk number 5: Xhard
###################################################
## generate realisation of hard core process here
## so we can use it for various diagrams
set.seed(424242)
Xhard <- rHardcore(100, 0.1)
###################################################
### code chunk number 6: XhardDemo
###################################################
## Make a diagram showing the simulated pattern (left panel)
## and the corresponding non-overlapping discs (right panel)
XhardX <- layered(Xhard, plotargs=list(pch=16))
XhardDiscs <- layered(Xhard %mark% 0.1,
plotargs=list(markscale=1, legend=FALSE))
XhardDemo <- solist(XhardX, XhardDiscs)
###################################################
### code chunk number 7: BirthDeathDemos
###################################################
## make illustrations for birth-death transitions
## using previous simulated pattern
XX <- Xhard
WW <- Window(XX)
CC <- as.ppp(centroid.owin(WW), WW)
II <- which.min(crossdist(XX, CC))
YY <- XX[-II]
ZZ <- XX[II]
hh <- 0.1
WW <- intersect.owin(Window(XX),
shift(shift(square(0.6), origin="midpoint"), XX[II]))
Window(XX) <- Window(YY) <- Window(ZZ) <- WW
XXr <- intersect.owin(dilation(XX, hh), WW)
YYr <- intersect.owin(dilation(YY, hh), WW)
Less <- layered(YYr, WW, YY, ZZ)
More <- layered(XXr, WW, YY, ZZ)
layerplotargs(Less) <-
list(list(border="grey", hatch=TRUE, hatchargs=list(col=grey(0.3))),
list(),
list(pch=16),
list(pch=1, cex=1.5, col=grey(0.3)))
layerplotargs(More) <-
list(list(border="grey", hatch=TRUE, hatchargs=list(col=grey(0.3))),
list(),
list(pch=16),
list(pch=16))
M <- matrix(c(4,1,1,4)/5, 2, 2)
xx <- M %*% (WW$xrange)
yy <- M %*% (WW$yrange)
Right <- ppp(x=xx, y=rep(yy[2],2), window=WW)
Left <- ppp(x=rev(xx), y=rep(yy[1],2), window=WW)
Arrows <- layered(WW,
onearrow(Right, txt="birth"),
onearrow(Left, txt="death"),
plotargs=list(list(type="n"), list(), list()))
###################################################
### code chunk number 8: 13gibbs.Rnw:211-225
###################################################
## make figure demonstrating 'forbidden zone' for hard core
# First make a yardstick, positioned so that it can be compared with the discs
xr <- WW$xrange
yr <- WW$yrange
xleft <- xr[1] - diff(xr)/6
ytarget <- mean(yr) + hh/2
ybest <- XX$y[which.min(abs(XX$y - ytarget))]
ys.h <- yardstick(xleft, ybest,
xleft, ybest - hh, expression(italic(h)))
# make figure
L3 <- Less[-4]
layerplotargs(L3) <- list(list(col="grey"), list(), list(pch=3))
HardForbid <- layered(L3, ys.h,
plotargs=list(list(), list(angle=90)))
###################################################
### code chunk number 9: StraussDemo
###################################################
Z <- YY[WW]
rr <- 0.11
ZS <- scanmeasure(Z, rr, dimyx=512)[WW]
m <- max(ZS)
ZS <- m - ZS
ZF <- eval.im(factor(as.integer(ZS), levels=0:m))
levels(ZF) <- LZF <- rev(c("beta", paste0("beta/", 2^(1:m))))
EZF <- parse(text=LZF)
xr <- WW$xrange
yr <- WW$yrange
xleft <- xr[1] - diff(xr)/6
ytarget <- mean(yr) + rr/2
ybest <- Z$y[which.min(abs(Z$y - ytarget))]
ys2 <- yardstick(xleft, ybest,
xleft, ybest - rr, expression(italic(R)))
Str <- layered(ZF, Z, ys2,
plotargs=list(
list(col=grey(sqrt((0:m)/m)),
ribargs=list(las=1, labels=EZF),
box=TRUE),
list(pch=3),
list(angle=90, pos=2)))
###################################################
### code chunk number 10: 13gibbs.Rnw:281-288
###################################################
set.seed(424242)
halfnpix <- 10
npix <- 2 * halfnpix + 1
A <- quadrats(owin(), npix, npix)
U <- tiles(A)[[halfnpix * npix + halfnpix + 1]]
Uplus <- grow.rectangle(U, 1/npix)
X <- runifpoint(20, setminus.owin(owin(), Uplus))
###################################################
### code chunk number 11: Unit.Rnw:3-5
###################################################
newplot(6, 0.7)
setmargins(0)
###################################################
### code chunk number 12: 13gibbs.Rnw:293-294
###################################################
zeromargins()
###################################################
### code chunk number 13: 13gibbs.Rnw:298-301
###################################################
plot(Window(A), main="")
plot(X, add=TRUE, pch=16)
plot(U, add=TRUE, lwd=3)
###################################################
### code chunk number 14: Unit2LR.Rnw:3-5
###################################################
newplot(13, 0.9)
setmargins(0,2,0,2)
###################################################
### code chunk number 15: 13gibbs.Rnw:359-362
###################################################
plot(solist(HardForbid,Str),
main="", main.panel="", nrows=1,
equal.scales=TRUE, mar.panel=0, hsep=1)
###################################################
### code chunk number 16: Unit2.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0)
###################################################
### code chunk number 17: 13gibbs.Rnw:510-511
###################################################
zeromargins()
###################################################
### code chunk number 18: 13gibbs.Rnw:515-518
###################################################
plot(XhardDemo,
main="", main.panel="", equal.scales=TRUE,
mar.panel=0, hsep=1.5)
###################################################
### code chunk number 19: Plist
###################################################
set.seed(191710)
Hlist <- Plist <- rpoispp(8, nsim=12)
Rinter <- 0.1
naughty <- (sapply(Hlist, minnndist) <= Rinter)
crossout <- psp(c(0,0), c(0,1), c(1,1), c(1,0), window=square(1), check=FALSE)
Hlist[naughty] <- lapply(Hlist[naughty],
function(z) layered(crossout, z,
plotargs=list(
list(lwd=4, col="grey"),
list(pch=16, cex=0.75))))
Hlist[!naughty] <- lapply(Hlist[!naughty],
function(z) layered(z,
plotargs=list(pch=16, cex=0.75)))
## add yardstick to first panel
ysv <- yardstick(-0.1, 0.5-Rinter/2, -0.1, 0.5+Rinter/2, expression(italic(h)))
Hlist[[1]] <- layered(Hlist[[1]], ysv,
plotargs=list(list(), list(pos=2, angle=90)))
###################################################
### code chunk number 20: Unit3x4y.Rnw:4-5
###################################################
setmargins(0)
###################################################
### code chunk number 21: 13gibbs.Rnw:596-600
###################################################
plot(Hlist, main="", main.panel="",
equal.scales=TRUE, halign=TRUE,
pch=16, cex=0.75,
nrows=3, mar.panel=0.1, hsep=1, vsep=1)
###################################################
### code chunk number 22: 13gibbs.Rnw:626-633
###################################################
intensityStrauss <- function(beta, gamma, R=0.1) {
G <- (1-gamma) * pi * R^2
lam <- LambertW(beta * G)/G
round(lam, 1)
}
lam8 <- intensityStrauss(8, 0)
lam100 <- intensityStrauss(100, 0)
###################################################
### code chunk number 23: birthdeathseq
###################################################
## generate spatial birth-death sequence
Hc <- 0.07
Beta <- 200
mod <- list(cif="hardcore",par=list(beta=Beta,hc=Hc), w=square(1))
X <- rmh(model=mod,
start=list(n.start=0),
control=list(p=0, q=0.35, nrep=1e4, nsave=10, nburn=0, track=TRUE))
Y <- attr(X, "saved")
ntran <- cumsum(attr(X, "history")$accepted)
Z <- list()
Z[[1]] <- ppp(, window=Window(Y[[1]]))
for(i in 1:9)
Z[[i+1]] <- Y[[min(which(ntran == (i * 50))) %/% 10]]
names(Z) <- (0:9) * 50
Z <- as.solist(Z)
###################################################
### code chunk number 24: Unit2x5t.Rnw:4-5
###################################################
zeromargins()
###################################################
### code chunk number 25: 13gibbs.Rnw:872-876
###################################################
plot(Z, main="",
pch=16, cex=0.6, cex.main=0.9, font.main=1,
nrows=2,
equal.scales=TRUE, mar.panel=c(0,0,1,0), hsep=1, vsep=0.75)
###################################################
### code chunk number 26: Unit3.Rnw:3-5
###################################################
newplot(19,0.9)
zeromargins()
###################################################
### code chunk number 27: 13gibbs.Rnw:930-931
###################################################
zeromargins()
###################################################
### code chunk number 28: 13gibbs.Rnw:935-938
###################################################
plot(solist(Less[-1], Arrows, More[-1]),
main="", main.panel="", nrows=1,
equal.scales=TRUE, mar.panel=0, hsep=0.1)
###################################################
### code chunk number 29: 13gibbs.Rnw:1018-1026
###################################################
set.seed(1917)
modIH <- rmhmodel(cif="hardcore",
par=list(beta=1, hc=0.05),
trend = function(x,y) { 500 * exp(-2*(x + y)) },
w=square(1))
XIH <- rmh(modIH, nrep=1e6)
DIH <- density(XIH, bw.scott)
YIH <- layered(XIH, plotargs=list(pch=16))
###################################################
### code chunk number 30: Unit2r.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0,0,0,1)
###################################################
### code chunk number 31: 13gibbs.Rnw:1033-1035
###################################################
plot(solist(YIH, DIH), main="", main.panel="",
equal.scales=TRUE, mar.panel=0, hsep=1)
###################################################
### code chunk number 32: 13gibbs.Rnw:1110-1112
###################################################
set.seed(424242)
SS <- rStrauss(100, 0.5, 0.1)
###################################################
### code chunk number 33: Unit2.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0)
###################################################
### code chunk number 34: 13gibbs.Rnw:1117-1118
###################################################
zeromargins()
###################################################
### code chunk number 35: 13gibbs.Rnw:1122-1127
###################################################
SSS <- layered(SS, plotargs=list(pch=16))
SSM <- layered(SS %mark% 0.1, plotargs=list(markscale=1, legend=FALSE))
plot(solist(SSS, SSM),
main="", main.panel="", equal.scales=TRUE,
mar.panel=0, hsep=1.5)
###################################################
### code chunk number 36: 13gibbs.Rnw:1176-1198
###################################################
sx <- sapply(Plist, function(z, R=Rinter) { d <- pairdist(z)
sum(d[row(d) > col(d)] <= R) })
gamma <- 0.5
accept <- (runif(length(Plist)) < gamma^sx)
Slist <- Plist
Slist[!accept] <- lapply(Slist[!accept],
function(z) layered(owin(),
crossout, z,
plotargs=list(
list(),
list(lwd=4, col="grey"),
list(pch=16, cex=0.75))))
Slist[accept] <- lapply(Slist[accept],
function(z) layered(z,
plotargs=list(pch=16, cex=0.75)))
## add yardstick to first panel
ysvR <- yardstick(-0.1, 0.5-Rinter/2, -0.1, 0.5+Rinter/2, expression(italic(R)))
Slist[[1]] <- layered(Slist[[1]], ysvR,
plotargs=list(list(), list(pos=2, angle=90)))
## main titles
px <- ifelse(sx == 0, "1", paste0("1/", 2^sx))
mains <- parse(text=paste("list(italic(s) == ", sx, ", italic(p) == ", px, ")"))
###################################################
### code chunk number 37: Unit3x4uy.Rnw:5-6
###################################################
setmargins(0, 1, 0, 0)
###################################################
### code chunk number 38: 13gibbs.Rnw:1205-1209
###################################################
plot(Slist, main="", main.panel=mains,
equal.scales=TRUE, halign=TRUE,
pch=16, cex=0.75, panel.args=function(...) list(cex.main=0.9),
nrows=3, mar.panel=0.1 + c(0, 1.5, 0, 0), hsep=1, vsep=1.5)
###################################################
### code chunk number 39: 13gibbs.Rnw:1231-1233
###################################################
lam8s <- intensityStrauss(8, 0.5)
lam100s <- intensityStrauss(100, 0.5)
###################################################
### code chunk number 40: 13gibbs.Rnw:1243-1257
###################################################
set.seed(111222)
rr <- 0.08
lambda <- 50
gamma <- (0:3)/3
beta <- lambda * exp(lambda * pi * rr^2 * (1 - gamma))
HSP <- list()
for(i in seq_along(gamma)) HSP[[i]] <- rStrauss(beta[i], gamma[i], rr)
HSP <- as.solist(HSP)
#ys <- yardstick(0.5-rr/2, -0.1, 0.5+rr/2, -0.1, expression(italic(h)))
HSP[[1]] <- layered(HSP[[1]], #ys,
plotargs=list(list(pch=16, cex=0.5)#,
# list(angle=90, pos=1)))
))
blurb <- parse(text=paste("gamma ==", sapply(gamma, simplenumber)))
###################################################
### code chunk number 41: Unit4u.Rnw:3-5
###################################################
newplot(25, 1)
zeromargins()
###################################################
### code chunk number 42: 13gibbs.Rnw:1269-1272
###################################################
plot(HSP, main="", main.panel=blurb,
equal.scales=TRUE, mar.panel=c(1,0,1,0), hsep=1, nrows=1, valign=TRUE,
pch=16, size=0.5)
###################################################
### code chunk number 43: 13gibbs.Rnw:1541-1595
###################################################
set.seed(67812)
HC <- 0.15
Beta <- 50
mod <- rmhmodel(cif='hardcore', par=list(beta=Beta, hc=HC), w=square(1))
u <- data.frame(x=0.5, y=0.5)
emp <- ppp(window=square(1))
Xu <- rmh(mod, x.cond=u, start=list(x.start=emp), expand=1, nrep=1e6)
CondDemo1 <- layered(owin(),
Xu[1],
Xu[-1],
disc(HC, Xu[1]),
textstring(0.55, 0.55, "u"),
plotargs=list(list(lwd=2),
list(pch=3, cex=1.5),
list(pch=16),
list(lty=2),
list(font=3)))
B <- owin(c(0,1/3), c(0,1))
set.seed(424242)
zz <- rHardcore(Beta, HC, B)
Xz <- rmh(mod, x.cond=zz, start=list(x.start=emp), expand=1, nrep=1e6)
Yz <- Xz[-(1:npoints(zz))]
CondDemo2 <- layered(owin(),
B,
zz,
Yz,
dilation(zz, HC),
textstring(0.1, 0.5, "B"),
textstring(0.75, 0.15, "A"),
plotargs=list(list(lwd=2),
list(lwd=2, lty=3),
list(pch=3),
list(pch=16),
list(lty=2),
list(font=3), list(font=3)))
CW <- owin(c(1/3, 1/3 + HC), c(0,1))
set.seed(999)
cz <- rHardcore(Beta, HC, CW)
Xc <- rmh(mod, x.cond=cz, start=list(x.start=emp), expand=1, nrep=1e6)
Yc <- Xc[-(1:npoints(cz))]
MarkovDemo <- layered(owin(),
CW,
cz,
Yc,
dilation(cz, HC),
textstring(0.1, 0.5, "B"),
textstring(0.4, 0.12, "C"),
textstring(0.75, 0.15, "A"),
plotargs=list(list(lwd=2),
list(lwd=2, lty=3),
list(pch=3),
list(pch=16),
list(lty=2),
list(font=3), list(font=3), list(font=3)))
###################################################
### code chunk number 44: Unit2.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0)
###################################################
### code chunk number 45: 13gibbs.Rnw:1600-1603
###################################################
plot(solist(CondDemo1, CondDemo2),
main="", main.panel="", nrows=1,
equal.scales=TRUE, mar.panel=0, hsep=1, valign=TRUE)
###################################################
### code chunk number 46: Unit.Rnw:3-5
###################################################
newplot(6, 0.7)
setmargins(0)
###################################################
### code chunk number 47: 13gibbs.Rnw:1667-1668
###################################################
plot(MarkovDemo, main="")
###################################################
### code chunk number 48: 13gibbs.Rnw:1729-1734
###################################################
## set digits to 4 for narrower output
Digits <- 4
dopt <- options(digits=Digits)
Terse <- 3
spatstat.options(terse=Terse)
###################################################
### code chunk number 49: fitStrauss-swedishpines5
###################################################
fithom <- ppm(swedishpines~1, Strauss(r=5))
co <- round(exp(coef(fithom)), 3)
###################################################
### code chunk number 50: 13gibbs.Rnw:1767-1768
###################################################
ppm(swedishpines ~ 1, Strauss(r=5))
###################################################
### code chunk number 51: 13gibbs.Rnw:1803-1804
###################################################
ppm(swedishpines ~ polynom(x,y,2), Strauss(5))
###################################################
### code chunk number 52: 13gibbs.Rnw:1833-1835
###################################################
fit0 <- ppm(swedishpines ~ 1, Strauss(5))
coef(fit0)
###################################################
### code chunk number 53: 13gibbs.Rnw:1845-1846
###################################################
unlist(parameters(fit0))
###################################################
### code chunk number 54: fitStrauss2swp
###################################################
fit2 <- ppm(swedishpines ~ polynom(x,y,2), Strauss(5))
###################################################
### code chunk number 55: 13gibbs.Rnw:1880-1882
###################################################
Dactive <- distfun(vesicles.extra$activezone)
ppm(vesicles ~ Dactive, Strauss(40))
###################################################
### code chunk number 56: 13gibbs.Rnw:2048-2049
###################################################
ppm(redwood ~ 1, Strauss(0.1))
###################################################
### code chunk number 57: 13gibbs.Rnw:2055-2056
###################################################
ppm(redwood ~ 1, Strauss(0.1), project=TRUE)
###################################################
### code chunk number 58: 13gibbs.Rnw:2075-2076
###################################################
ppm(cells ~ 1, Strauss(0.00001))
###################################################
### code chunk number 59: 13gibbs.Rnw:2164-2168
###################################################
fit0 <- ppm(swedishpines ~ 1, Strauss(r=5))
fit2 <- update(fit0, . ~ polynom(x,y,2))
fit9 <- update(fit0, Strauss(9))
fit0fine <- update(fit0, nd=256)
###################################################
### code chunk number 60: 13gibbs.Rnw:2184-2186
###################################################
sqrt(diag(vcov(fit0)))
confint(fit0)
###################################################
### code chunk number 61: 13gibbs.Rnw:2189-2190
###################################################
coef(summary(fit0))
###################################################
### code chunk number 62: 13gibbs.Rnw:2198-2199
###################################################
vcov(fit0)
###################################################
### code chunk number 63: 13gibbs.Rnw:2240-2245
###################################################
fit2 <- ppm(swedishpines ~ polynom(x,y,2), Strauss(5))
Z <- solist(predict(fit2, type="trend"),
predict(fit2, type="cif"),
predict(fit2, type="intensity"))
ZZ <- as.anylist(lapply(Z, transect.im))
###################################################
### code chunk number 64: Unit3r.Rnw:3-5
###################################################
newplot(22,0.9)
setmargins(0)
###################################################
### code chunk number 65: 13gibbs.Rnw:2251-2252
###################################################
plot(Z, equal.ribbon=TRUE, main="", main.panel="", mar.panel=0, hsep=1)
###################################################
### code chunk number 66: fv3.Rnw:3-5
###################################################
newplot(12.5, 1.0)
setmargins(0.5+c(3,3,0,1))
###################################################
### code chunk number 67: 13gibbs.Rnw:2283-2284
###################################################
setmargins(0.5+c(2,2,0,0))
###################################################
### code chunk number 68: 13gibbs.Rnw:2288-2289
###################################################
plot(ZZ, equal.scales=TRUE, main="", main.panel="")
###################################################
### code chunk number 69: 13gibbs.Rnw:2301-2303
###################################################
scrambled <- runifpoint(ex=swedishpines)
ll <- predict(fit2, type="cif", X=scrambled)
###################################################
### code chunk number 70: 13gibbs.Rnw:2314-2315
###################################################
inten <- intensity(fit0)
###################################################
### code chunk number 71: 13gibbs.Rnw:2325-2326
###################################################
intensity(fit0)
###################################################
### code chunk number 72: 13gibbs.Rnw:2353-2354
###################################################
fitin(fit2)
###################################################
### code chunk number 73: 13gibbs.Rnw:2368-2369
###################################################
as.interact(fit2)
###################################################
### code chunk number 74: simStrauss2
###################################################
s <- simulate(fit2, nsim=100)
###################################################
### code chunk number 75: Unit4.Rnw:3-5
###################################################
newplot(25, 1)
zeromargins()
###################################################
### code chunk number 76: 13gibbs.Rnw:2392-2394
###################################################
plot(s[1:4], nrows=1, equal.scales=TRUE,
main = "", main.panel="", pch=16, hsep=1, mar.panel=0)
###################################################
### code chunk number 77: 13gibbs.Rnw:2419-2423
###################################################
npoints(swedishpines)
np <- sapply(s, npoints)
mean(np)
sd(np)
###################################################
### code chunk number 78: 13gibbs.Rnw:2442-2443 (eval = FALSE)
###################################################
## E0 <- envelope(fit0, Lest, nsim=39)
###################################################
### code chunk number 79: 13gibbs.Rnw:2485-2486
###################################################
step(fit2, trace=0)
###################################################
### code chunk number 80: 13gibbs.Rnw:2517-2519
###################################################
fit2P <- update(fit2, Poisson())
anova(fit2P, fit2, test="LR")
###################################################
### code chunk number 81: 13gibbs.Rnw:2522-2523
###################################################
anova(fit0, fit2, test="LR")
###################################################
### code chunk number 82: fv4x2.Rnw:3-5
###################################################
newplot(5, 1)
setmargins(0.5+c(3,3,2,1))
###################################################
### code chunk number 83: 13gibbs.Rnw:2637-2638
###################################################
setmargins(0.5+c(2,2,1,0))
###################################################
### code chunk number 84: 13gibbs.Rnw:2642-2714
###################################################
doit <- function(nama, expr, ...) {
height <- 1.3
width <- 1
plot(c(-0.25, width), c(-0.15,height), type="n", main="",
xlab="", ylab="", axes=FALSE)
title(main=nama, cex=1.1, font=2, line=0.5)
AXEX <- 1.2
text(0.9*width, 0.15, expression(bolditalic(d)), cex=AXEX)
## remember to reconcile the next line with \pairpot
text(-0.15, height, expression(bolditalic(c(d))), cex=AXEX)
plot(onearrow(ppp(c(0,0), c(0, height), window=square(height))),
retract=0, headfraction=0.15, col.head=1, add=TRUE)
axis(2, at=c(0,1), tick=TRUE, pos=0, las=1, cex.axis=AXEX)
plot(onearrow(ppp(c(0,width), c(0, 0), window=square(width))),
retract=0, headfraction=0.15, col.head=1, add=TRUE)
axis(1, ..., tick=TRUE, pos=0, cex.axis=AXEX)
segments(0,1,1,1, lty=3, lwd=2)
eval(substitute(curve(E, lwd=3, col=1, add=TRUE, n=513, from=0),
list(E=substitute(expr))))
}
par(mfrow=c(2, 4))
doit("Hard Core",
as.integer(x > 0.5),
at=c(0, 0.5), labels=expression(0, italic(h)))
doit("Strauss",
ifelse(x > 0.5, 1, 0.5),
at=c(0, 0.5), labels=expression(0, italic(R)))
# strauss-hardcore
doit("Strauss-Hard Core",
ifelse(x <= 0.3, 0, ifelse(x <= 0.6, 0.5, 1)),
at=c(0, 0.3, 0.6), labels=expression(0, italic(h), italic(R)))
# pairpiece
doit("Piecewise constant",
c(0.15, 0.4, 0.2, 0.6, 1)[findInterval(x, (0:5)/5, TRUE, TRUE)],
at=(0:4)/5,
labels=expression(0, italic(r[1]), italic(r[2]), italic(r[3]),
italic(r[4])))
# diggle-gratton
delta <- 0.3
rho <- 0.7
kappa <- 1.8
doit("Diggle-Gratton",
ifelse(x <= delta, 0,
ifelse(x >= rho, 1,
((x-delta)/(rho-delta))^kappa)),
at=c(0, delta, rho),
labels=expression(0, delta, rho))
# diggle-gates-stibbard
rho <- 0.7
doit("Diggle-Gates-Stibbard",
ifelse(x >= rho, 1,
sin((pi/2) * x/rho)^2),
at=c(0, rho),
labels=expression(0, rho))
# softcore
sigma <- 0.3
kappa <- 0.5
doit("Soft Core",
exp(-(sigma/x)^(2/kappa)),
at=c(0, sigma),
labels=expression(0, sigma))
# Fiksel
aa <- -3
kappa <- 4
h <- 0.25
r <- 0.75
doit("Fiksel",
ifelse(x <= h, 0,
ifelse(x >= r, 1,
exp( aa * exp(-kappa * x)))),
at=c(0, h, r),
labels=expression(0, h, r))
###################################################
### code chunk number 85: 13gibbs.Rnw:2724-2732
###################################################
set.seed(1000)
Xdgs <- rStraussHard(200, 0.5, 0.08, 0.04)
fakefit <- ppm(Xdgs ~ 1, StraussHard(0.08, 0.04))
cifdgs <- predict(fakefit, type="cif", ngrid=1024,
new.coef=log(c(200, 0.5)))
Xdgs <- layered(Xdgs, plotargs=list(list(pch=16)))
cifXdgs <- layered(cifdgs, plotargs=list(list(col=grey(seq(0.2,1,length=128)))))
Zdgs <- solist(Xdgs, cifXdgs)
###################################################
### code chunk number 86: Unit2r.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0,0,0,1)
###################################################
### code chunk number 87: 13gibbs.Rnw:2738-2739
###################################################
plot(Zdgs, main="", main.panel="", mar.panel=0, hsep=1, equal.scales=TRUE)
###################################################
### code chunk number 88: 13gibbs.Rnw:2983-2985
###################################################
mod <- ppm(swedishpines ~ polynom(x,y,2), PairPiece(1:15))
f <- fitin(mod)
###################################################
### code chunk number 89: 13gibbs.Rnw:3014-3016
###################################################
rr <- data.frame(r=seq(3,14,by=0.05))
p <- profilepl(rr, Strauss, swedishpines~polynom(x,y,2), aic=TRUE)
###################################################
### code chunk number 90: 13gibbs.Rnw:3018-3019
###################################################
p
###################################################
### code chunk number 91: fv2.Rnw:3-5
###################################################
newplot(12, 0.95)
setmargins(0.5+c(3,3,0,1))
###################################################
### code chunk number 92: 13gibbs.Rnw:3032-3034
###################################################
plot(anylist(f, p), main="", main.panel="",
mar.panel=0.2+c(4,4,0,1))
###################################################
### code chunk number 93: 13gibbs.Rnw:3053-3054
###################################################
as.ppm(p)
###################################################
### code chunk number 94: 13gibbs.Rnw:3124-3129
###################################################
set.seed(191919)
ai <- function(beta, eta, r=0.07)
rmhmodel(cif='areaint', par=list(beta=beta, eta=eta, r=r), w=square(1))
XA.02 <- rmh(ai(250, 0.02))
XA50 <- rmh(ai(5, 50))
###################################################
### code chunk number 95: Unit2.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0)
###################################################
### code chunk number 96: 13gibbs.Rnw:3135-3139
###################################################
plot(solist(XA.02, XA50),
main="", main.panel="", nrows=1,
equal.scales=TRUE, mar.panel=0, hsep=1,
pch=16)
###################################################
### code chunk number 97: 13gibbs.Rnw:3184-3207
###################################################
## Make synthetic example for explaining Area Interaction
W <- as.owin(swedishpines)
x <- c(28,29,55,60,66)
y <- c(70,38,32,72,59)
X <- ppp(x=x,y=y, window = W)
u <- list(x=48,y=50)
u <- as.ppp(u, W)
rad <- 14
Xplusr <- dilation(X, rad)
uplusr <- disc(rad, u)
ovlap <- intersect.owin(uplusr, Xplusr)
AIdemo <- layered(W,
ovlap,
Xplusr,
uplusr,
X,
u)
layerplotargs(AIdemo) <- list(list(),
list(col="darkgrey", border=NA),
list(lwd=2),
list(lwd=2, lty=2),
list(pch=16),
list(pch=3))
###################################################
### code chunk number 98: Unit.Rnw:3-5
###################################################
newplot(6, 0.7)
setmargins(0)
###################################################
### code chunk number 99: 13gibbs.Rnw:3214-3215
###################################################
plot(AIdemo, main="")
###################################################
### code chunk number 100: 13gibbs.Rnw:3275-3277 (eval = FALSE)
###################################################
## fitSA <- ppm(swedishpines ~ polynom(x,y,2), AreaInter(4))
## cifSA <- predict(fitSA, type = "cif", ngrid=512)
###################################################
### code chunk number 101: cifSA
###################################################
fitSA <- ppm(swedishpines~polynom(x,y,2), AreaInter(4))
cifSA <- predict(fitSA, type = "cif", ngrid=if(draftversion) 128 else 512)
cifSAdemo <- layered(cifSA, swedishpines,
plotargs=list(list(log=TRUE, ribargs=list(las=1)),
list(col="white", pch=3)))
intSA <- fitin(fitSA)
###################################################
### code chunk number 102: 13gibbs.Rnw:3287-3288
###################################################
fitSA
###################################################
### code chunk number 103: fv.Rnw:3-5
###################################################
newplot(6, 0.5)
setmargins(0.5+c(3,3,1,0))
###################################################
### code chunk number 104: 13gibbs.Rnw:3306-3307
###################################################
plot(intSA, main="", legend=FALSE)
###################################################
### code chunk number 105: UnitR.Rnw:3-6
###################################################
newplot(9, 0.7)
zeromargins() # strip all margins
setmargins(0.1 + c(0,0,0,2)) # back off
###################################################
### code chunk number 106: 13gibbs.Rnw:3311-3312
###################################################
setmargins(0,0,0,3)
###################################################
### code chunk number 107: 13gibbs.Rnw:3314-3315
###################################################
plot(cifSAdemo, main="")
###################################################
### code chunk number 108: 13gibbs.Rnw:3377-3384
###################################################
set.seed(422442)
gs <- function(beta, gamma, r=0.07, sat=2) {
rmhmodel(cif='geyer', par=list(beta=beta, gamma=gamma, sat=sat, r=r),
w=square(1))
}
XG.5 <- rmh(gs(300, 0.5))
XG2 <- rmh(gs(20, 2))
###################################################
### code chunk number 109: Unit2.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0)
###################################################
### code chunk number 110: 13gibbs.Rnw:3391-3395
###################################################
plot(solist(XG.5, XG2),
main="", main.panel="", nrows=1,
equal.scales=TRUE, mar.panel=0, hsep=1,
pch=16)
###################################################
### code chunk number 111: 13gibbs.Rnw:3436-3437
###################################################
ppm(redwood ~ 1, Geyer(r=0.05, sat=2))
###################################################
### code chunk number 112: 13gibbs.Rnw:3441-3444
###################################################
df <- expand.grid(r=seq(0.02, 0.1, by=0.001), sat=c(1,2))
pG <- profilepl(df, Geyer, redwood ~ 1, correction="translate",
aic=TRUE)
###################################################
### code chunk number 113: 13gibbs.Rnw:3446-3447
###################################################
as.ppm(pG)
###################################################
### code chunk number 114: fv.Rnw:3-5
###################################################
newplot(6, 0.5)
setmargins(0.5+c(3,3,1,0))
###################################################
### code chunk number 115: 13gibbs.Rnw:3453-3458
###################################################
plot(pG, main="",
lwd=2,
lty=ifelse(sat == 2, 1, 2), ## undocumented feature of plot.profilepl
pos=2, # text to left of peaks
lwd.opt=2, col.opt=1, tag=FALSE)
###################################################
### code chunk number 116: TripletsCalc
###################################################
set.seed(19501973)
tripmod <- function(gamma) {
rmhmodel(cif="triplets",par=list(beta=200,gamma=gamma,r=0.1), w=square(1))
}
X0 <- rmh(tripmod(0))
X.5 <- rmh(tripmod(0.5))
fitT0 <- ppm(X0 ~ 1, Triplets(0.1))
cifT0 <- predict(fitT0, type="cif", ngrid=256, new.coef=log(c(200, 1e-16)))
facT0 <- eval.im(factor(cifT0 > 100))
faclab <- expression(0, beta)
CifT0 <- layered(facT0, plotargs=list(list(col=grey(c(0.2, 0.9)),
box=TRUE,
ribargs=list(las=1, labels=faclab))))
###################################################
### code chunk number 117: Unit3.Rnw:3-5
###################################################
newplot(19,0.9)
zeromargins()
###################################################
### code chunk number 118: 13gibbs.Rnw:3551-3552
###################################################
setmargins(0,0,0,1.5)
###################################################
### code chunk number 119: 13gibbs.Rnw:3556-3559
###################################################
plot(solist(X0, X.5, CifT0),
main="", main.panel="", nrows=1,
equal.scales=TRUE, mar.panel=0, hsep=1)
###################################################
### code chunk number 120: 13gibbs.Rnw:3588-3590
###################################################
clo <- closepairs(redwoodfull, 0.1)
crx <- with(clo, psp(xi, yi, xj, yj, window=Window(redwoodfull)))
###################################################
### code chunk number 121: 13gibbs.Rnw:3596-3597
###################################################
redcon <- connected(redwoodfull, 0.1)
###################################################
### code chunk number 122: Unit2.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0)
###################################################
### code chunk number 123: 13gibbs.Rnw:3604-3609
###################################################
redcross <- layered(redwoodfull, crx)
layerplotargs(redcon) <- list(legend=FALSE)
plot(solist(redcross, redcon),
main="", main.panel="", nrows=1,
equal.scales=TRUE, mar.panel=0, hsep=1)
###################################################
### code chunk number 124: 13gibbs.Rnw:3646-3647
###################################################
hcg <- hclust(dist(coords(gordon)), "single")
###################################################
### code chunk number 125: Unit2.Rnw:3-5
###################################################
newplot(12.5, 0.9)
setmargins(0)
###################################################
### code chunk number 126: 13gibbs.Rnw:3652-3653
###################################################
setmargins(0.2,3,0,0)
###################################################
### code chunk number 127: 13gibbs.Rnw:3657-3658
###################################################
plot(hcg, main="", xlab="", sub="", labels=FALSE)
###################################################
### code chunk number 128: 13gibbs.Rnw:3667-3672
###################################################
set.seed(191919)
f <- function(X) { coef(ppm(X ~ 1, Concom(2), rbord=0))[["Interaction"]] }
Xsims <- rpoispp(intensity(gordon), win=Window(gordon), nsim=99)
fXsims <- sapply(Xsims, f)
pv <- (1+sum(f(gordon) <= fXsims))/100
###################################################
### code chunk number 129: 13gibbs.Rnw:3685-3686 (eval = FALSE)
###################################################
## plot(hclust(dist(coords(gordon)), "single"))
###################################################
### code chunk number 130: 13gibbs.Rnw:3699-3700
###################################################
ppm(gordon ~ 1, Concom(2), rbord=0)
###################################################
### code chunk number 131: 13gibbs.Rnw:3706-3711
###################################################
etahat <- function(X) {
mod <- ppm(X ~ 1, Concom(2), rbord=0)
return(parameters(mod)$eta)
}
(etaobs <- etahat(gordon))
###################################################
### code chunk number 132: 13gibbs.Rnw:3713-3714
###################################################
set.seed(898292)
###################################################
### code chunk number 133: ConcomMCtest
###################################################
Xsims <- rpoispp(ex=gordon, nsim=99)
etasim <- sapply(Xsims, etahat)
###################################################
### code chunk number 134: 13gibbs.Rnw:3720-3721
###################################################
(pval <- mean(etaobs <= c(etaobs, etasim)))
###################################################
### code chunk number 135: 13gibbs.Rnw:3871-3872 (eval = FALSE)
###################################################
## Hybrid(Strauss(0.1), Geyer(0.2, 3))
###################################################
### code chunk number 136: 13gibbs.Rnw:3883-3886
###################################################
M <- Hybrid(H=Hardcore(), G=Geyer(0.2, 3))
fit <- ppm(redwood ~ 1, M, correction="translate")
fit
###################################################
### code chunk number 137: 13gibbs.Rnw:3927-3929
###################################################
Mplus <- Hybrid(fit, g=Geyer(0.1,1))
update(fit, Mplus)
###################################################
### code chunk number 138: 13gibbs.Rnw:3935-3939 (eval = FALSE)
###################################################
## fit2 <- ppm(swedishpines ~ 1,
## Hybrid(DG=DiggleGratton(2,10), S=Strauss(5)))
## plot(fitin(fit2))
## plot(fitin(fit2), separate=TRUE, mar.panel=rep(4,4))
###################################################
### code chunk number 139: 13gibbs.Rnw:4029-4031
###################################################
fit <- ppm(swedishpines ~ 1, Strauss(8))
Y <- rmh(fit)
###################################################
### code chunk number 140: 13gibbs.Rnw:4072-4076
###################################################
mo <- list(cif="strauss",
par=list(beta=2,gamma=0.2,r=0.7), w=square(10))
X <- rmh(model=mo,
start=list(n.start=42), control=list(nrep=1e6))
###################################################
### code chunk number 141: 13gibbs.Rnw:4096-4098 (eval = FALSE)
###################################################
## mo <- rmhmodel(cif="strauss",
## par=list(beta=2,gamma=0.2,r=0.7), w=square(10))
###################################################
### code chunk number 142: 13gibbs.Rnw:4112-4116 (eval = FALSE)
###################################################
## rmhmodel(cif = c('hardcore', 'strauss'),
## par = list(list(beta = 10, hc = 0.03),
## list(beta = 1, gamma = 0.5, r = 0.07)),
## w = square(1))
###################################################
### code chunk number 143: 13gibbs.Rnw:4164-4165
###################################################
X4 <- rmh(model=mo, nrep=1e4)
###################################################
### code chunk number 144: 13gibbs.Rnw:4269-4272
###################################################
fit <- ppm(swedishpines ~ 1, Strauss(8))
Y <- rmh(fit, track=TRUE)
h <- attr(Y, "history")
###################################################
### code chunk number 145: 13gibbs.Rnw:4274-4276
###################################################
with(h, mean(accepted))
with(h, prop.table(table(proposaltype, accepted), 1))
###################################################
### code chunk number 146: CopperFitS
###################################################
dlin <- distfun(copper$SouthLines)
copfit <- ppm(copper$SouthPoints ~ dlin, Geyer(1, 1))
###################################################
### code chunk number 147: cdftestMCMC
###################################################
coptest <- cdf.test(copfit, dlin, nsim=39)
###################################################
### code chunk number 148: 13gibbs.Rnw:4608-4609
###################################################
coptest
###################################################
### code chunk number 149: SwedFit
###################################################
swedfit <- ppm(swedishpines ~ polynom(x,y,2), Strauss(9))
###################################################
### code chunk number 150: SwedEnv
###################################################
swedenv <- envelope(swedfit, Lest, nsim=39, global=TRUE,
savepatterns=TRUE)
###################################################
### code chunk number 151: fv.Rnw:3-5
###################################################
newplot(6, 0.5)
setmargins(0.5+c(3,3,1,0))
###################################################
### code chunk number 152: 13gibbs.Rnw:4645-4646
###################################################
plot(swedenv, . - r ~ r, main="", legend=FALSE)
###################################################
### code chunk number 153: 13gibbs.Rnw:4763-4764
###################################################
swedR <- residuals(swedfit)
###################################################
### code chunk number 154: Unit2LRboth.Rnw:4-6
###################################################
newplot(13, 1)
setmargins(0,2,0,2)
###################################################
### code chunk number 155: 13gibbs.Rnw:4774-4777
###################################################
plot(solist(swedR, Smooth(swedR)),
main="", main.panel="", nrows=1,
equal.scales=TRUE, mar.panel=c(0,1,0,1), hsep=1)
###################################################
### code chunk number 156: 13gibbs.Rnw:4852-4853
###################################################
swedI <- residuals(swedfit, type="inverse")
###################################################
### code chunk number 157: 13gibbs.Rnw:4863-4864
###################################################
swedE <- eem(swedfit)
###################################################
### code chunk number 158: 13gibbs.Rnw:4872-4873
###################################################
swedP <- residuals(swedfit, type="pearson")
###################################################
### code chunk number 159: Unit2LRboth.Rnw:4-6
###################################################
newplot(13, 1)
setmargins(0,2,0,2)
###################################################
### code chunk number 160: 13gibbs.Rnw:4878-4881
###################################################
plot(solist(swedI, swedP),
main="", main.panel="", nrows=1,
equal.scales=TRUE, mar.panel=c(0,1,0,1), hsep=1)
###################################################
### code chunk number 161: 13gibbs.Rnw:4981-4982 (eval = FALSE)
###################################################
## diagnose.ppm(swedfit, type="Pearson", envelope=TRUE, nsim=39)
###################################################
### code chunk number 162: 13gibbs.Rnw:4985-4986 (eval = FALSE)
###################################################
## diagnose.ppm(swedfit, type="Pearson", envelope=swedenv, nsim=39)
###################################################
### code chunk number 163: 13gibbs.Rnw:4989-4991
###################################################
DS <- diagnose.ppm(swedfit, type="Pearson",
plot.it=FALSE, envelope=swedenv, nsim=39)
###################################################
### code chunk number 164: Unit.Rnw:3-5
###################################################
newplot(6, 0.7)
setmargins(0)
###################################################
### code chunk number 165: 13gibbs.Rnw:4996-4997
###################################################
setmargins(1)
###################################################
### code chunk number 166: 13gibbs.Rnw:5001-5002
###################################################
plot(DS, main="")
###################################################
### code chunk number 167: 13gibbs.Rnw:5037-5039
###################################################
SWP <- Smooth(residuals(swedfit, type="Pearson"))
2/(2 * sqrt(pi) * attr(SWP, "sigma"))
###################################################
### code chunk number 168: 13gibbs.Rnw:5052-5053
###################################################
coppar <- parres(copfit, dlin)
###################################################
### code chunk number 169: fv.Rnw:3-5
###################################################
newplot(6, 0.5)
setmargins(0.5+c(3,3,1,0))
###################################################
### code chunk number 170: 13gibbs.Rnw:5075-5076
###################################################
plot(coppar, main="", legend=FALSE)
###################################################
### code chunk number 171: 13gibbs.Rnw:5131-5133
###################################################
set.seed(191919)
q9 <- qqplot.ppm(swedfit, nsim=39, plot.it=FALSE)
###################################################
### code chunk number 172: fv.Rnw:3-5
###################################################
newplot(6, 0.5)
setmargins(0.5+c(3,3,1,0))
###################################################
### code chunk number 173: 13gibbs.Rnw:5137-5138
###################################################
setmargins(0.5+c(4,3,0,0))
###################################################
### code chunk number 174: 13gibbs.Rnw:5142-5143
###################################################
plot(q9, main="", monochrome=monochrome, pch=3)
###################################################
### code chunk number 175: 13gibbs.Rnw:5188-5189
###################################################
KresSwed <- Kres(swedfit, correction="iso")
###################################################
### code chunk number 176: fv.Rnw:3-5
###################################################
newplot(6, 0.5)
setmargins(0.5+c(3,3,1,0))
###################################################
### code chunk number 177: 13gibbs.Rnw:5196-5197
###################################################
plot(KresSwed, main="", shade=c("ihi", "ilo"), legend=FALSE)
###################################################
### code chunk number 178: chorleyDfun
###################################################
## Code block repeated from Chapter 11
## squared distance to incinerator
d2incin <- function(x, y, xincin=354.5, yincin=413.6) {
(x - xincin)^2 + (y - yincin)^2
}
## Diggle's model of elevated risk as function of distance
fincin <- function(d2, alpha, beta) {
1 + alpha * exp( - beta * d2)
}
## elevated risk as spatial covariate
raisin <- function(x,y, alpha, beta) {
fincin(d2incin(x,y), alpha, beta)
}
###################################################
### code chunk number 179: morescore
###################################################
## Code block repeated from Chapter 11
Zalpha <- function(x,y, alpha, beta) {
expbit <- exp( - beta * d2incin(x,y))
expbit/(1 + alpha * expbit)
}
Zbeta <- function(x,y, alpha, beta) {
d2 <- d2incin(x,y)
topbit <- alpha * exp( - beta * d2)
- d2 * topbit/(1 + topbit)
}
Zscore <- list(alpha=Zalpha, beta=Zbeta)
###################################################
### code chunk number 180: moreHessian
###################################################
## Code block repeated from Chapter 11
Zaa <- function(x,y, alpha, beta) {
expbit <- exp( - beta * d2incin(x,y))
-(expbit/(1 + alpha * expbit))^2
}
Zab <- function(x,y, alpha, beta) {
d2 <- d2incin(x,y)
expbit <- exp( - beta * d2)
- d2 * expbit/(1 + alpha * expbit)^2
}
Zbb <- function(x,y, alpha, beta) {
d2 <- d2incin(x,y)
topbit <- alpha * exp( - beta * d2)
(d2^2) * topbit/(1 + topbit)^2
}
Zhess <- list(aa=Zaa, ab=Zab, ba=Zab, bb=Zbb)
###################################################
### code chunk number 181: chorley
###################################################
lung <- split(chorley)$lung
larynx <- split(chorley)$larynx
smo <- density(lung, sigma=0.15, eps=0.1)
smo <- eval.im(pmax(smo, 1e-10))
###################################################
### code chunk number 182: chorleyDGfit
###################################################
rGey <- 0.25
chorleyDGfit <- ippm(larynx ~ offset(log(smo) + log(raisin)) ,
Geyer(rGey, 1), eps = 0.5, iScore=Zscore,
start=list(alpha=20, beta=1),
nlm.args=list(steptol=0.001))
###################################################
### code chunk number 183: 13gibbs.Rnw:5360-5362
###################################################
unlist(parameters(chorleyDGfit))
chorleyDGopt <- parameters(chorleyDGfit)[c("alpha", "beta")]
###################################################
### code chunk number 184: chorleyDGcalc
###################################################
chorleyDGlev <- leverage(chorleyDGfit,
iScore=Zscore, iHessian=Zhess, iArgs=chorleyDGopt)
chorleyDGinf <- influence(chorleyDGfit,
iScore=Zscore, iHessian=Zhess, iArgs=chorleyDGopt)
chorleyDGdfb <- dfbetas(chorleyDGfit,
iScore=Zscore, iHessian=Zhess, iArgs=chorleyDGopt)
###################################################
### code chunk number 185: 13gibbs.Rnw:5374-5388
###################################################
spatstat.options(image.colfun=greytoblack)
incin <- as.ppp(chorley.extra[["incin"]], W = Window(chorley))
chorleyDGlevDemo <-
layered(chorleyDGlev,
Window(chorley),
incin,
plotargs=list(list(show.all=FALSE, ribbon=TRUE, ribside="bottom"),
list(),
list(pch=10, col="white", etch=TRUE)))
chorleyDGinfDemo <-
layered(chorleyDGinf,
incin,
plotargs=list(list(maxsize=3.0, lwd=3, leg.side="bottom"),
list(pch=10, col="white", etch=TRUE)))
###################################################
### code chunk number 186: Chorley2b.Rnw:3-5
###################################################
newplot(18, 0.8)
setmargins(3,0,0,0)
###################################################
### code chunk number 187: 13gibbs.Rnw:5394-5397
###################################################
plot(solist(chorleyDGlevDemo, chorleyDGinfDemo),
main="", main.panel="", mar.panel=c(1.2,0,0,0), hsep=1,
nrows=1, equal.scales=TRUE)
###################################################
### code chunk number 188: 13gibbs.Rnw:5417-5418
###################################################
spatstat.options(image.colfun=blacktowhite)
###################################################
### code chunk number 189: augmentchorleyDG
###################################################
## This code helps to accelerate the plotting
chorleyDGdfb <- augment.msr(chorleyDGdfb)
###################################################
### code chunk number 190: Chorley2x2LR.Rnw:3-5
###################################################
newplot(20,1)
setmargins(0.5)
###################################################
### code chunk number 191: chorleyDGdfb
###################################################
white134 <- function(i) {
if(i==2) list(cols="black") else
list(cols="white", leg.args=list(cols="black"))
}
plot(chorleyDGdfb,
panel.end=as.owin(chorley),
panel.args=white134, lwd=2, box=FALSE,
mar.panel=0, hsep=2, vsep=1, maxsize=2.5,
main="")
###################################################
### code chunk number 192: chorleyDGsmooth
###################################################
chorleyDGdfbsmo <- Smooth(chorleyDGdfb, sigma=0.8)
###################################################
### code chunk number 193: Chorley2x2R.Rnw:3-5
###################################################
newplot(20,0.95)
setmargins(0, 0.5, 0, 4)
###################################################
### code chunk number 194: chorleyDGdfbsmo
###################################################
plot(chorleyDGdfbsmo,
panel.end=as.owin(chorley),
main="", box=FALSE,
mar.panel=0, hsep=2, vsep=1)
###################################################
### code chunk number 195: 13gibbs.Rnw:5492-5493
###################################################
###################################################
### code chunk number 196: 13gibbs.Rnw:5495-5511
###################################################
## if(!file.exists("data/RedGey.rda")) {
## # get bandwidth value
## if(!file.exists("data/RedBW.rda")) {
## NOT YET RELEASED IN SPATSTAT
## RedBW <- bw.locppm(redwoodfull, ~1, Geyer(0.05, 2),
## correction="isotropic",
## srange=c(0.06, 0.4), ns=32,
## verbose=FALSE)
## save(RedBW, file="data/RedBW.rda", compress=TRUE)
## } else load("data/RedBW.rda")
## # fit model
## NOT YET RELEASED IN SPATSTAT
## RedGey <- locppm(redwoodfull, ~1, Geyer(0.05, 2), sigma=RedBW,
## correction="isotropic",
## vcalc="full", locations="fine",
## verbose=FALSE)
## save(RedGey, file="data/RedGey.rda", compress=TRUE)
## } else load("data/RedGey.rda")
###################################################
### code chunk number 197: 13gibbs.Rnw:5536-5538 (eval = FALSE)
###################################################
## NOT YET RELEASED IN SPATSTAT
## RedGey <- locppm(redwoodfull ~ 1, Geyer(0.05, 2), sigma=bw.locppm,
## correction="isotropic")
###################################################
### code chunk number 198: UnitR.Rnw:3-6
###################################################
newplot(9, 0.7)
zeromargins() # strip all margins
setmargins(0.1 + c(0,0,0,2)) # back off
###################################################
### code chunk number 199: 13gibbs.Rnw:5550-5552
###################################################
## NOT YET RELEASED IN SPATSTAT
## plot(RedGey, main="", which=2,
## style="imagecontour", contourargs=list(col="white"))
###################################################
### code chunk number 200: PromptOff.Rnw:1-2
###################################################
options(prompt=" ")
###################################################
### code chunk number 201: 13gibbs.Rnw:6320-6321 (eval = FALSE)
###################################################
## ppm(X ~ terms, interaction)
###################################################
### code chunk number 202: PromptOn.Rnw:1-2
###################################################
options(prompt="> ")
###################################################
### code chunk number 203: 13gibbs.Rnw:6398-6402
###################################################
## reload.or.compute(datafilepath("ppmStraussSimdatHO.rda"),
## { fit <- ppm(simdat, ~1, Strauss(r=0.275), method="ho") })
set.seed(424242)
fit <- ppm(simdat~1, Strauss(r=0.275), method="ho")
###################################################
### code chunk number 204: 13gibbs.Rnw:6404-6405 (eval = FALSE)
###################################################
## fit <- ppm(simdat ~ 1, Strauss(r=0.275), method="ho")
###################################################
### code chunk number 205: 13gibbs.Rnw:6407-6409
###################################################
fit
vcov(fit)
###################################################
### code chunk number 206: SwedFit
###################################################
p.mean <- rep(0,7)
p.var <- diag(c(rep(10000,6), 0.001))
swedfitVB <- ppm(swedishpines ~ polynom(x,y,2), Strauss(9),
method = "VBlogi",
prior.mean = p.mean, prior.var = p.var)
###################################################
### code chunk number 207: PromptOff.Rnw:1-2
###################################################
options(prompt=" ")
###################################################
### code chunk number 208: 13gibbs.Rnw:6650-6651 (eval = FALSE)
###################################################
## model <- dppGauss(lambda=100, alpha=0.05, d=2)
###################################################
### code chunk number 209: PromptOn.Rnw:1-2
###################################################
options(prompt="> ")
###################################################
### code chunk number 210: PromptOff.Rnw:1-2
###################################################
options(prompt=" ")
###################################################
### code chunk number 211: 13gibbs.Rnw:6663-6664 (eval = FALSE)
###################################################
## simulate(model)
###################################################
### code chunk number 212: PromptOn.Rnw:1-2
###################################################
options(prompt="> ")
###################################################
### code chunk number 213: PromptOff.Rnw:1-2
###################################################
options(prompt=" ")
###################################################
### code chunk number 214: 13gibbs.Rnw:6671-6673 (eval = FALSE)
###################################################
## X <- residualspaper[["Fig1"]]
## jfit <- dppm(X ~ polynom(x,y,3), dppGauss())
###################################################
### code chunk number 215: PromptOn.Rnw:1-2
###################################################
options(prompt="> ")