1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
|
#'
#' bw.pplHeat.R
#'
#' Bandwidth selection for densityHeat.ppp
#' by point process likelihood cross-validation
#'
#' Copyright (c) 2020 Adrian Baddeley, Tilman Davies and Suman Rakshit
#' GNU Public Licence >= 2.0
bw.pplHeat <- function(X, ..., srange=NULL, ns=16, sigma=NULL,
leaveoneout=TRUE, verbose=TRUE) {
#' compute intensity estimates
b <- HeatEstimates.ppp(X, ..., srange=srange, ns=ns, sigma=sigma,
leaveoneout=leaveoneout, verbose=verbose)
lambda <- b$lambda
h <- b$h
hname <- b$hname
#' compute likelihood cross-validation criterion
CV <- rowSums(log(lambda))
iopt <- which.max(CV)
result <- bw.optim(CV, h, iopt,
criterion="Likelihood cross-validation",
hname=hname,
unitname=unitname(X))
return(result)
}
HeatEstimates.ppp <- function(X, ..., srange=NULL, ns=16, sigma=NULL,
leaveoneout=FALSE, verbose=TRUE) {
stopifnot(is.ppp(X))
nX <- npoints(X)
## trap a common error
if(length(argh <- list(...)) &&
(is.null(nama <- names(argh)) || !nzchar(nama[[1L]])) &&
is.numeric(a <- argh[[1L]]) &&
length(a) == 1L)
stop("Use argument 'sigma' to specify the maximum bandwidth!",
call.=FALSE)
## determine candidate bandwidths
if(is.numeric(sigma) && length(sigma)) {
## sigma is a vector of candidate bandwidths, or a maximum bandwidth
sigMax <- max(sigma)
fractions <- if(length(sigma) > 1) sigma/sigMax else
geomseq(from=0.05, to=1, length.out=ns)
} else if(is.im(sigma)) {
## sigma is an image giving the spatially-varying maximum bandwidth
sigMax <- sigma
fractions <- seq_len(ns)/ns
} else if(is.null(sigma)) {
#' make a sequence of candidate bandwidths
if(!is.null(srange)) {
check.range(srange)
} else {
nnd <- nndist(X)
srange <- c(min(nnd[nnd > 0]), diameter(as.owin(X))/2)
}
sigMax <- srange[2]
sigmavalues <- geomseq(from=srange[1L], to=srange[2L], length.out=ns)
fractions <- sigmavalues/sigMax
} else stop("Format of sigma is not understood")
## set up transition matrix and initial state
a <- densityHeat.ppp(X, sigMax, ..., internal=list(setuponly=TRUE))
Y <- a$Y # initial state image
u <- a$u # initial state vector (dropping NA)
Xpos <- a$Xpos # location of each data point, index in 'u'
A <- a$A # transition matrix, operates on 'u;
Nstep <- a$Nstep # total number of iterations
## map desired sigma values to iteration numbers
nits <- pmax(1L, pmin(Nstep, round(Nstep * fractions^2)))
nits <- diff(c(0L,nits))
reciprocalpixelarea <- with(Y, 1/(xstep * ystep))
## compute ....
lambda <- matrix(nrow=ns, ncol=nX)
if(!leaveoneout) {
## usual estimates
for(k in seq_len(ns)) {
for(l in seq_len(nits[k])) u <- u %*% A
lambda[k, ] <- u[Xpos]
}
} else {
## compute leave-one-out estimates
if(verbose) {
cat("Processing", nX, "points ... ")
pstate <- list()
}
for(i in seq_len(nX)) {
## initial state = X[-i]
ui <- u
Xposi <- Xpos[i]
ui[Xposi] <- ui[Xposi] - reciprocalpixelarea
## run iterations, pausing at each sigma value
for(k in seq_len(ns)) {
for(l in seq_len(nits[k])) ui <- ui %*% A
lambda[k, i] <- ui[Xposi]
}
if(verbose) pstate <- progressreport(i, nX, state=pstate)
}
if(verbose) cat("Done.\n")
}
if(!is.im(sigma)) {
h <- sigMax * fractions
hname <- "sigma"
} else {
h <- fractions
hname <- "fract"
}
return(list(lambda=lambda, h=h, hname=hname))
}
|
