# R/bc.R In spatstat.core: Core Functionality of the 'spatstat' Family

#### Documented in bcbc.ppmrex

```#'     bc.R
#'
#'  Bias correction techniques
#'
#'  \$Revision: 1.3 \$ \$Date: 2022/01/04 05:30:06 \$

bc <- function(fit, ...) {
UseMethod("bc")
}

bc.ppm <- function(fit, ..., nfine=256) {
stopifnot(is.ppm(fit))
#
theta0 <- coef(fit)
nc <- length(theta0)
#
X <- data.ppm(fit)
# evaluate sufficient statistic at data points
sufX <- model.matrix(fit)[Z, ]
if(ncol(sufX) != nc)
stop("Internal error: model.matrix does not match coef(model)")
# predict on fine grid
if(length(sufF) != nc)
stop("Internal error: model.images does not match coef(model)")
# edge correction
if(fit\$correction == "border" && ((rbord <- fit\$rbord) > 0)) {
bX <- bdist.points(X)
excludeU <- eval.im(b < rbord)
retainX  <- (bX >= rbord)
sufX <- sufX[retainX, , drop=FALSE]
} else {
excludeU <- FALSE
}
# compute fine approximation to score
scoreX <- colSums(sufX)
scoreW <- numeric(nc)
for(k in seq_len(nc)) {
S <- sufF[[k]]
# infinite values of S may occur and correspond to zero cif
Slam <- eval.im(ifelse(is.infinite(S) | excludeU, 0, S * lamF))
scoreW[k] <- integral.im(Slam)
}
score <- scoreX - scoreW
# Newton-Raphson
Iinv <- vcov(fit, hessian=TRUE)
theta <- theta0 + Iinv %*% score
theta <- theta[ , 1L, drop=TRUE]
#
#  return(list(theta0=theta0, theta=theta))
return(theta)
}

# Richardson extrapolation (generic)

rex <- function(x, r=2, k=1, recursive=FALSE) {
# x should be a matrix
# whose columns are successive estimates of a parameter vector
# obtained using "grid step sizes" t, t/r, t/r^2, ...
# Estimate from step size t is assumed to converge at rate t^k
if(!is.matrix(x)) x <- matrix(x, nrow=1)
if(ncol(x) <= 1) return(x)
rk <- r^k
y <- (rk * x[, -1L, drop=FALSE] - x[, -ncol(x), drop=FALSE])/(rk - 1)
if(recursive)
y <- rex(y, r=r, k=k+1, recursive=TRUE)
return(y)
}
```

## Try the spatstat.core package in your browser

Any scripts or data that you put into this service are public.

spatstat.core documentation built on May 18, 2022, 9:05 a.m.