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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
#
#
# satpiece.S
#
# $Revision: 1.17 $ $Date: 2018/03/15 07:37:41 $
#
# Saturated pairwise interaction process with piecewise constant potential
#
# SatPiece() create an instance of the process
# [an object of class 'interact']
#
#
# -------------------------------------------------------------------
#
SatPiece <- local({
# ..... auxiliary functions ......
delSP <- function(i, r, sat) {
r <- r[-i]
sat <- sat[-i]
nr <- length(r)
if(nr == 0) return(Poisson())
if(nr == 1) return(Geyer(r, sat))
return(SatPiece(r, sat))
}
# ....... template object ..........
BlankSatPiece <-
list(
name = "piecewise constant Saturated pairwise interaction process",
creator = "SatPiece",
family = "pairsat.family", # evaluated later
pot = function(d, par) {
r <- par$r
nr <- length(r)
out <- array(FALSE, dim=c(dim(d), nr))
out[,,1] <- (d < r[1])
if(nr > 1) {
for(i in 2:nr)
out[,,i] <- (d >= r[i-1]) & (d < r[i])
}
out
},
par = list(r = NULL, sat=NULL), # filled in later
parnames = c("interaction thresholds", "saturation parameters"),
hasInf = FALSE,
init = function(self) {
r <- self$par$r
sat <- self$par$sat
if(!is.numeric(r) || !all(r > 0))
stop("interaction thresholds r must be positive numbers")
if(length(r) > 1 && !all(diff(r) > 0))
stop("interaction thresholds r must be strictly increasing")
if(!is.numeric(sat) || any(sat < 0))
stop("saturation parameters must be nonnegative numbers")
if(any(ceiling(sat) != floor(sat)))
warning("saturation parameter has a non-integer value")
if(length(sat) != length(r) && length(sat) != 1)
stop("vectors r and sat must have equal length")
},
update = NULL, # default OK
print = NULL, # default OK
interpret = function(coeffs, self) {
r <- self$par$r
npiece <- length(r)
# extract coefficients
gammas <- exp(as.numeric(coeffs))
# name them
gn <- gammas
names(gn) <- paste("[", c(0,r[-npiece]),",", r, ")", sep="")
#
return(list(param=list(gammas=gammas),
inames="interaction parameters gamma_i",
printable=dround(gn)))
},
valid = function(coeffs, self) {
# interaction parameters gamma must be
# non-NA
# finite, if sat > 0
# less than 1, if sat = Inf
gamma <- (self$interpret)(coeffs, self)$param$gammas
sat <- self$par$sat
if(anyNA(gamma))
return(FALSE)
return(all((is.finite(gamma) | sat == 0)
& (gamma <= 1 | sat != Inf)))
},
project = function(coeffs, self){
loggammas <- as.numeric(coeffs)
sat <- self$par$sat
r <- self$par$r
ok <- is.finite(loggammas) & (is.finite(sat) | loggammas <= 0)
if(all(ok))
return(NULL)
if(!any(ok))
return(Poisson())
bad <- !ok
if(spatstat.options("project.fast") || sum(bad) == 1) {
# remove smallest threshold with an unidentifiable parameter
firstbad <- min(which(bad))
return(delSP(firstbad, r, sat))
} else {
# consider all candidate submodels
subs <- lapply(which(bad), delSP, r=r, sat=sat)
return(subs)
}
},
irange = function(self, coeffs=NA, epsilon=0, ...) {
r <- self$par$r
sat <- self$par$sat
if(all(is.na(coeffs)))
return(2 * max(r))
gamma <- (self$interpret)(coeffs, self)$param$gammas
gamma[is.na(gamma)] <- 1
active <- (abs(log(gamma)) > epsilon) & (sat > 0)
if(!any(active))
return(0)
else return(2 * max(r[active]))
},
version=NULL # added later
)
class(BlankSatPiece) <- "interact"
SatPiece <- function(r, sat) {
instantiate.interact(BlankSatPiece, list(r=r, sat=sat))
}
SatPiece <- intermaker(SatPiece, BlankSatPiece)
SatPiece
})
|
