Nothing
R/refclass.r
In palm: Fitting Point Process Models via the Palm Likelihood
Defines functions
create.obj
set.giveparent.class
set.totaldeletion.class
set.void.class
set.matern.class
set.thomas.class
set.twocamerachild.class
set.binomchild.class
set.poischild.class
set.sibling.class
set.ns.class
set.buffer.class
set.pbc.class
set.nlminb.class
set.bobyqa.class
set.fit.class
set.CLASSNAME.class
######
## General base class.
######
base.class.R6 <- R6Class("palm",
public = list(
## Setting fields.
n.patterns = NULL,
lims = NULL,
lims.list = NULL,
vol = NULL,
vol.list = NULL,
dim = NULL,
dim.list = NULL,
classes = NULL,
## Initialisation method.
initialize = function(lims.list, classes, ...){
self$lims.list <- lims.list
self$vol.list <- lapply(self$lims.list,
function(x) prod(apply(x, 1, diff)))
self$dim.list <- lapply(self$lims.list,
function(x) nrow(x))
self$n.patterns <- length(lims.list)
self$classes <- classes
},
## A method to set up a new pattern.
setup.pattern = function(pattern){
if (pattern > self$n.patterns){
stop("Pattern index exceeds the number of patterns")
}
self$lims <- self$lims.list[[pattern]]
self$vol <- self$vol.list[[pattern]]
self$dim <- self$dim.list[[pattern]]
},
## A method to simulate multiple patterns.
simulate = function(pars = self$par.fitted){
out <- vector(mode = "list", length = self$n.patterns)
for (i in 1:self$n.patterns){
self$setup.pattern(i)
out[[i]] <- self$simulate.pattern(pars)
}
if (self$n.patterns == 1){
out <- out[[1]]
}
out
},
## A method to simulate a single pattern.
simulate.pattern = function(pars){},
## A method to trim points to the observation window.
trim.points = function(points, output.indices = FALSE){
in.window <- rep(TRUE, nrow(points))
for (i in 1:self$dim){
in.window <- in.window & (self$lims[i, 1] <= points[, i] & self$lims[i, 2] >= points[, i])
}
if (output.indices){
out <- which(in.window)
} else {
out <- points[in.window, , drop = FALSE]
}
out
}
))
######
## Template for new classes.
######
## Replace CLASSNAME with class name, then add fields and methods.
set.CLASSNAME.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("CLASSNAME.inherit", class, envir = class.env)
R6Class("palm_CLASSNAME",
inherit = class.env$CLASSNAME.inherit,
public = list(
))
}
######
## Class for models fitted to data.
######
set.fit.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("fit.inherit", class, envir = class.env)
R6Class("palm_fit",
inherit = class.env$fit.inherit,
public = list(
points.list = NULL,
points = NULL,
n.points = NULL,
n.points.list = NULL,
contrasts = NULL,
contrasts.list = NULL,
n.contrasts = NULL,
n.contrasts.list = NULL,
contrast.pairs = NULL,
contrast.pairs.list = NULL,
R = NULL,
boots = NULL,
trace = NULL,
conv.code = NULL,
converged = NULL,
n.par = NULL,
set.start = NULL,
set.bounds = NULL,
par.names = NULL,
fixed.names = NULL,
par.names.link = NULL,
par.start = NULL,
par.start.link = NULL,
par.fixed = NULL,
par.fixed.link = NULL,
par.links = NULL,
par.invlinks = NULL,
par.fitted = NULL,
par.fitted.link = NULL,
par.lower = NULL,
par.lower.link = NULL,
par.upper = NULL,
par.upper.link = NULL,
pi.multiplier = NULL,
## Initialisation method.
initialize = function(points.list, R, trace, start, bounds, ...){
super$initialize(...)
self$points.list <- points.list
## Checking list of points and list of lims have the same length.
if (length(self$points.list) != length(self$lims.list)){
stop("The list 'points' must have the same number of components as 'lims'.")
}
self$R <- R
## Creating numbers of points and contrasts for each pattern.
self$n.points.list <- lapply(points.list, nrow)
self$contrasts.list <- list()
self$contrast.pairs.list <- list()
self$n.contrasts.list <- list()
for (i in 1:self$n.patterns){
self$get.contrasts(i)
}
## Starting out with the first pattern for start values and such.
self$setup.pattern(1)
self$trace <- trace
self$set.start <- start
self$set.bounds <- bounds
self$get.pars()
self$n.par <- length(self$par.start)
self$par.start.link <- self$link.pars(self$par.start)
self$get.link.bounds()
},
## Overwriting the method to set up a new pattern.
setup.pattern = function(pattern, do.contrasts = TRUE){
super$setup.pattern(pattern)
self$points <- self$points.list[[pattern]]
self$n.points <- self$n.points.list[[pattern]]
if (do.contrasts){
self$contrast.pairs <- self$contrast.pairs.list[[pattern]]
self$contrasts <- self$contrasts.list[[pattern]]
self$n.contrasts <- self$n.contrasts.list[[pattern]]
}
},
## An empty method for getting contrasts.
get.contrasts = function(pattern){
self$n.contrasts.list[[pattern]] <- length(self$contrasts.list[[pattern]])
},
## A method for getting the parameters across all classes.
get.pars = function(){
self$fetch.pars()
},
## An empty method for fetching parameters from a particular class.
fetch.pars = function(){
},
## A method for adding new parameters.
add.pars = function(name, link.name = NULL, link, invlink = NULL,
start, lower, upper){
self$par.names <- c(self$par.names, name)
## Sorting out inverse link.
if (identical(link, log)){
full.link <- function(pars) log(pars[name])
## Default name and inverse for log link.
if (is.null(link.name)){
link.name <- paste("log", name, sep = ".")
}
if (is.null(invlink)){
full.invlink <- function(pars) exp(pars[link.name])
}
} else if (identical(link, logit)){
full.link <- function(pars) link(pars[name])
## Default name and inverse for logit link.
if (is.null(link.name)){
link.name <- paste("logit", name, sep = ".")
}
if (is.null(invlink)){
full.invlink <- function(pars) invlogit(pars[link.name])
}
} else {
if (is.null(link.name)){
stop("Please provide link name.")
}
if (is.null(invlink)){
stop("Please provide inverse link function.")
}
full.link <- link
full.invlink <- invlink
}
self$par.links <- c(self$par.links, full.link)
names(self$par.links) <- self$par.names
self$par.invlinks <- c(self$par.invlinks, full.invlink)
names(self$par.invlinks) <- self$par.names
self$par.names.link <- c(self$par.names.link, link.name)
## Overwriting start parameter, if one provided.
if (any(name == names(self$set.start))){
start <- self$set.start[name]
}
if (any(name == names(self$set.bounds))){
lower <- self$set.bounds[[name]][1]
upper <- self$set.bounds[[name]][2]
}
self$par.upper <- c(self$par.upper, upper)
names(self$par.upper) <- self$par.names
self$par.lower <- c(self$par.lower, lower)
## Adjustment for start values on the bounds.
if (start >= upper){
start <- lower + 0.9*(upper - lower)
} else if (start <= lower){
start <- lower + 0.9*(upper - lower)
}
self$par.start <- c(self$par.start, start)
names(self$par.start) <- self$par.names
names(self$par.lower) <- self$par.names
},
## A method to set the upper and lower parameter bounds on the link scale.
get.link.bounds = function(){
self$par.lower.link <- self$link.pars(self$par.lower)
self$par.upper.link <- self$link.pars(self$par.upper)
},
## A method for converting parameters to their link scales.
link.pars = function(pars){
out <- numeric(self$n.par)
for (i in 1:self$n.par){
out[i] <- self$par.links[[i]](pars)
}
names(out) <- self$par.names.link
out
},
## A method for converting parameters to their real scale.
invlink.pars = function(pars){
out <- numeric(self$n.par)
for (i in 1:self$n.par){
out[i] <- self$par.invlinks[[i]](pars)
}
names(out) <- self$par.names
out
},
## A method for the objective function.
obj.fun = function(link.pars, fixed.link.pars = NULL,
est.names = NULL, fixed.names = NULL){
combined.pars <- c(link.pars, fixed.link.pars)
names(combined.pars) <- c(est.names, fixed.names)
link.pars <- combined.pars[self$par.names.link]
pars <- self$invlink.pars(link.pars)
obj.fun.components <- numeric(self$n.patterns)
## Summing over all the patterns.
for (i in 1:self$n.patterns){
self$setup.pattern(i)
obj.fun.components[i] <- self$neg.log.palm.likelihood(pars)
}
out <- sum(obj.fun.components)
ll <- -out
if (self$trace){
for (i in 1:self$n.par){
cat(self$par.names[i], ": ", sep = "")
cat(pars[i], ", ", sep = "")
}
cat("Log-lik: ", ll, "\n", sep = "")
}
out
},
## An empty method for the Palm intensity.
palm.intensity = function(r, pars){},
## An empty method for the sum of the log Palm intensities.
sum.log.intensities = function(pars){
sum(log(self$pi.multiplier*self$palm.intensity(self$contrasts, pars)))
},
## A default method for the integral in the Palm likelihood.
palm.likelihood.integral = function(pars){
f <- function(r, pars){
self$palm.intensity(r, pars)*Sd(r, self$dim)
}
-self$pi.multiplier*integrate(f, lower = 0, upper = self$R, pars = pars)$value
},
## A default method for the log of the Palm likelihood function.
log.palm.likelihood = function(pars){
self$sum.log.intensities(pars) + self$palm.likelihood.integral(pars)
},
## A method for the negative Palm likelihood function.
neg.log.palm.likelihood = function(pars){
-self$log.palm.likelihood(pars)
},
## An empty method for optimisation.
palm.optimise = function(){},
## A method for non-optimisation.
palm.not.optimise = function(){
self$converged <- FALSE
self$par.fitted.link <- rep(NA, self$n.par)
names(self$par.fitted.link) <- self$par.names.link
self$par.fitted <- rep(NA, self$n.par)
self$conv.code <- 5
},
## A method for model fitting.
fit = function(){
## If there are contrasts, fit the model. Otherwise, return something else.
if (self$n.contrasts <= 1){
warning("Not enough observed points to compute parameter estimates.")
self$palm.not.optimise()
} else {
self$palm.optimise()
}
},
## A method for bootstrapping.
boot = function(N, prog = TRUE){
boots <- matrix(0, nrow = N, ncol = self$n.par)
## Setting up progress bar.
if (prog){
pb <- txtProgressBar(min = 0, max = N, style = 3)
}
for (i in 1:N){
zero.points <- TRUE
while (zero.points){
sim.obj <- self$simulate()
if (self$n.patterns > 1){
sim.obj$points <- lapply(sim.obj, function(x) x$points)
sim.obj$sibling.list <- lapply(sim.obj, function(x) x$sibling.list)
zero.points <- sum(sapply(sim.obj$points, nrow)) == 0
} else {
zero.points <- nrow(sim.obj$points) == 0
}
}
## Doesn't matter that sibling.list is non-null below, as sibling class is not passed.
obj.boot <- create.obj(classes = self$classes, points = sim.obj$points, lims = self$lims,
R = self$R, child.list = self$child.list,
sibling.list = sim.obj$sibling.list, trace = FALSE,
start = self$par.fitted, bounds = self$bounds)
obj.boot$fit()
if (obj.boot$converged){
boots[i, ] <- obj.boot$par.fitted
} else {
boots[i, ] <- rep(NA, self$n.par)
}
## Updating progress bar.
if (prog){
setTxtProgressBar(pb, i)
}
}
cat("\n")
self$boots <- boots
colnames(self$boots) <- self$par.names
},
## A method for plotting.
plot = function(xlim = NULL, ylim = NULL, show.empirical = TRUE, breaks = 50, ...){
if (show.empirical){
emp <- self$empirical(xlim, breaks)
}
if (is.null(xlim)){
xlim <- self$get.xlim()
}
xx <- seq(xlim[1], xlim[2], length.out = 1000)
yy <- self$palm.intensity(xx, self$par.fitted)
if (is.null(ylim)){
if (show.empirical){
ylim <- c(0, max(c(emp$y, yy)))
} else {
ylim <- c(0, max(yy))
}
}
par(xaxs = "i")
plot(xx, yy, xlab = "", ylab = "", xlim = range(xx), ylim = ylim, type = "n", ...)
title(xlab = "r", ylab = "Palm intensity")
lines(xx, yy)
if (show.empirical){
lines(emp$x, emp$y, lty = "dashed")
}
},
## A method for default xlim.
get.xlim = function(){
c(0, self$R)
},
## A method for empirical Palm intensity.
empirical = function(xlim = NULL, breaks = 50){
if (is.null(xlim)){
xlim <- self$get.xlim()
}
dists <- pbc_distances(self$points, self$lims)
midpoints <- seq(0, xlim[2], length.out = breaks)
midpoints <- midpoints[-length(midpoints)]
h <- diff(midpoints[c(1, 2)])
midpoints[1] <- midpoints[1] + h/2
intensities <- numeric(length(midpoints))
for (i in 1:length(midpoints)){
halfwidth <- ifelse(i == 1, 0.25*h, 0.5*h)
n.interval <- sum(dists <= (midpoints[i] + halfwidth)) -
sum(dists <= (midpoints[i] - halfwidth))
area <- Vd(midpoints[i] + halfwidth, self$dim) - Vd(midpoints[i] - halfwidth, self$dim)
intensities[i] <- n.interval/(self$pi.multiplier*area)
}
list(x = midpoints, y = intensities)
}
))
}
######
## Class for bobyqa() optimisation.
######
set.bobyqa.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("bobyqa.inherit", class, envir = class.env)
R6Class("palm_bobyqa",
inherit = class.env$bobyqa.inherit,
public = list(
## A method to minimise the negative Palm likelihood.
palm.optimise = function(){
out <- try(bobyqa(self$par.start.link, self$obj.fun,
lower = self$par.lower.link, upper = self$par.upper.link,
fixed.link.pars = self$par.fixed.link,
est.names = self$par.names.link,
fixed.names = self$fixed.names.link,
control = list(maxfun = 2000, npt = 2*self$n.par + 1)), silent = TRUE)
if (class(out)[1] == "try-error"){
self$palm.not.optimise()
} else {
if (out$ierr > 0){
warning("Failed convergence.")
self$converged <- FALSE
} else {
self$converged <- TRUE
}
self$par.fitted.link <- out$par
names(self$par.fitted.link) <- self$par.names.link
self$par.fitted <- self$invlink.pars(self$par.fitted.link)
self$conv.code <- out$ierr
}
}
))
}
######
## Class for nlminb() optimisation.
######
set.nlminb.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("nlminb.inherit", class, envir = class.env)
R6Class("palm_nlminb",
inherit = class.env$nlminb.inherit,
public = list(
## A method to minimise the negative Palm likelihood.
palm.optimise = function(){
out <- nlminb(self$par.start.link, self$obj.fun,
fixed.link.pars = self$par.fixed.link,
est.names = self$par.names.link, fixed.names = self$fixed.names.link,
control = list(eval.max = 2000, iter.max = 1500),
lower = self$par.lower.link, upper = self$par.upper.link)
if (out$convergence > 1){
warning("Failed convergence.")
self$converged <- FALSE
} else {
self$converged <- TRUE
}
self$par.fitted.link <- out$par
names(self$par.fitted.link) <- self$par.names.link
self$par.fitted <- self$invlink.pars(self$par.fitted.link)
self$conv.code <- out$convergence
}
))
}
######
## Class for periodic boundary conditions.
######
set.pbc.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("pbc.inherit", class, envir = class.env)
R6Class("palm_pbc",
inherit = class.env$pbc.inherit,
public = list(
## Overwriting the method to set up a new pattern.
setup.pattern = function(pattern, do.contrasts = TRUE){
super$setup.pattern(pattern, do.contrasts)
self$pi.multiplier <- self$n.points/2
},
## A method to generate contrasts.
get.contrasts = function(pattern){
self$setup.pattern(pattern, do.contrasts = FALSE)
## Saving which contrast applies to which pair of observations.
contrast.pairs <- matrix(0, nrow = (self$n.points^2 - self$n.points)/2, ncol = 2)
k <- 1
if (self$n.points > 1){
for (i in 1:(self$n.points - 1)){
for (j in (i + 1):self$n.points){
contrast.pairs[k, ] <- c(i, j)
k <- k + 1
}
}
contrasts <- pbc_distances(points = self$points, lims = self$lims)
self$contrast.pairs.list[[pattern]] <- contrast.pairs[contrasts <= self$R, ]
self$contrasts.list[[pattern]] <- contrasts[contrasts <= self$R]
} else {
self$contrast.pairs.list[[pattern]] <- contrast.pairs
self$contrasts[[pattern]] <- numeric(0)
}
super$get.contrasts(pattern)
}
))
}
######
## Class for buffer-zone boundary conditions.
######
set.buffer.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("buffer.inherit", class, envir = class.env)
R6Class("palm_buffer",
inherit = class.env$buffer.inherit,
public = list(
## Overwriting the method to set up a new pattern.
setup.pattern = function(pattern, do.contrasts = TRUE){
super$setup.pattern(pattern, do.contrasts)
self$pi.multiplier <- self$n.points
},
## A method to generate contrasts.
get.contrasts = function(pattern){
self$setup.pattern(pattern, do.contrasts = FALSE)
## Getting rid of contrasts between two external points.
buffer.keep <- buffer_keep(points = self$points, lims = self$lims,
R = self$R)
contrasts <- as.vector(as.matrix(dist(self$points))[buffer.keep])
contrast.pairs.1 <- matrix(rep(1:self$n.points, self$n.points), nrow = self$n.points)
contrast.pairs.2 <- matrix(rep(1:self$n.points, self$n.points), nrow = self$n.points, byrow = TRUE)
contrast.pairs.1 <- as.vector(contrast.pairs.1[buffer.keep])
contrast.pairs.2 <- as.vector(contrast.pairs.2[buffer.keep])
contrast.pairs <- cbind(contrast.pairs.1, contrast.pairs.2)
## Now truncating to contrasts less than R.
self$contrasts.list[[pattern]] <- contrasts[contrasts <= self$R]
self$contrast.pairs.list[[pattern]] <- contrast.pairs[contrasts <= self$R, ]
super$get.contrasts(pattern)
}
))
}
######
## Class for Neyman-Scott processes.
######
set.ns.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("ns.inherit", class, envir = class.env)
R6Class("palm_ns",
inherit = class.env$ns.inherit,
public = list(
child.list = NULL,
## Adding density parameter.
fetch.pars = function(){
super$fetch.pars()
link <- function(pars){
log(pars["D"]*self$child.expectation(pars))
}
## An inverse link that uses the other unlinked paramters.
invlink <- function(pars, out){
exp(pars["log.Dc"])/self$child.expectation(out)
}
self$add.pars(name = "D", link.name = "log.Dc", link = link, invlink = invlink,
start = sqrt(self$n.points/self$vol), lower = 0, upper = Inf)
},
## Overwriting the method for converting parameters to their real scale.
invlink.pars = function(pars){
out <- numeric(self$n.par)
names(out) <- self$par.names
which.D <- which(self$par.names == "D")
for (i in (1:self$n.par)[-which.D]){
out[i] <- self$par.invlinks[[i]](pars)
}
out[which.D] <- self$par.invlinks[[which.D]](pars, out)
out
},
## Overwriting simulation methods.
simulate.pattern = function(pars = self$par.fitted){
parent.locs <- self$get.parents(pars)
n.parents <- nrow(parent.locs)
sim.n.children <- self$simulate.n.children(n.parents, pars)
n.children <- sim.n.children$n.children
sibling.list <- sim.n.children$sibling.list
child.locs <- matrix(0, nrow = sum(n.children), ncol = self$dim)
j <- 0
for (i in seq_along(n.children)){
if (n.children[i] > 0){
child.locs[(j + 1):(j + n.children[i]), ] <-
self$simulate.children(n.children[i], parent.locs[i, ], pars)
j <- j + n.children[i]
}
}
parent.ids <- rep(seq_along(n.children), times = n.children)
trimmed <- self$trim.points(child.locs, output.indices = TRUE)
list(points = child.locs[trimmed, , drop = FALSE],
parents = parent.locs,
parent.ids = parent.ids[trimmed],
child.ys = sim.n.children$child.ys[trimmed],
sibling.list = self$trim.siblings(sibling.list, trimmed))
},
## A method to trim the sibling list.
trim.siblings = function(sibling.list, trimmed){
sibling.list
},
## A method to get parents by simulation.
get.parents = function(pars){
expected.parents <- pars["D"]*self$vol
n.parents <- rpois(1, expected.parents)
parent.locs <- matrix(0, nrow = n.parents, ncol = self$dim)
for (i in 1:self$dim){
parent.locs[, i] <- runif(n.parents, self$lims[i, 1], self$lims[i, 2])
}
parent.locs
},
## Overwriting method for the Palm intensity.
palm.intensity = function(r, pars){
self$sibling.pi(r, pars) + self$nonsibling.pi(pars)
},
## An empty method for the expectation of the child distribution.
child.expectation = function(pars){},
## An empty method for the variance of the child distribution.
child.variance = function(pars){},
## A method for the expected number of siblings of a randomly chosen point.
sibling.expectation = function(pars){
(self$child.variance(pars) + self$child.expectation(pars)^2)/self$child.expectation(pars) - 1
},
## An empty method for the PDF of Q, the between-sibling distances.
fq = function(r, pars){},
## An empty method for the CDF of Q.
Fq = function(r, pars){},
## A default method for the quotient of the PDF of Q and the surface volume.
## Numerically unstable; better to replace in child classes.
fq.over.s = function(r, pars){
self$fq(r, pars)/Sd(r, dim)
},
## A method for the Palm intensity of nonsibling points.
nonsibling.pi = function(pars){
pars["D"]*self$child.expectation(pars)
},
## A method for the PI of sibling points.
sibling.pi = function(r, pars){
self$sibling.expectation(pars)*self$fq.over.s(r, pars)
}
))
}
######
## Class for known sibling relationships.
######
set.sibling.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("sibling.inherit", class, envir = class.env)
R6Class("palm_sibling",
inherit = class.env$sibling.inherit,
public = list(
siblings = NULL,
## Lol one day you'll get a bug because you're an
## idiot and you have objects called 'sibling.list'
## and 'siblings.list' and you'll be so mad.
siblings.list = NULL,
sibling.list = NULL,
sibling.mat = NULL,
sibling.alpha = NULL,
sibling.beta = NULL,
## Initialisation method.
initialize = function(sibling.list, ...){
self$sibling.list <- sibling.list
super$initialize(...)
self$siblings.list <- list()
for (i in 1:self$n.patterns){
self$get.siblings(i)
}
},
## Overwriting the method to set up a new pattern.
setup.pattern = function(pattern, do.contrasts = TRUE, do.siblings = TRUE){
super$setup.pattern(pattern, do.contrasts)
self$sibling.mat <- self$sibling.list[[pattern]]$sibling.mat
self$sibling.alpha <- self$sibling.list[[pattern]]$alpha
self$sibling.beta <- self$sibling.list[[pattern]]$beta
if (do.siblings){
self$siblings <- self$siblings.list[[pattern]]
}
},
## A method to get vector of sibling relationships
## that matches with the contrasts.
get.siblings = function(pattern){
self$setup.pattern(pattern, do.siblings = FALSE)
siblings <- numeric(self$n.contrasts)
for (i in 1:self$n.contrasts){
pair <- self$contrast.pairs[i, ]
siblings[i] <- self$sibling.mat[pair[1], pair[2]]
}
## 0 for known nonsiblings.
## 1 for known siblings.
## 2 for unknown sibling status.
siblings <- as.numeric(siblings)
siblings[is.na(siblings)] <- 2
self$siblings.list[[pattern]] <- siblings
},
## Overwriting the sum of the log intensities.
sum.log.intensities = function(pars){
all.palm.intensities <- numeric(self$n.contrasts)
all.palm.intensities[self$siblings == 0] <-
self$sibling.beta*self$nonsibling.pi(pars)
all.palm.intensities[self$siblings == 1] <-
self$sibling.alpha*
self$sibling.pi(self$contrasts[self$siblings == 1], pars)
all.palm.intensities[self$siblings == 2] <-
(1 - self$sibling.beta)*self$nonsibling.pi(pars) +
(1 - self$sibling.alpha)*
self$sibling.pi(self$contrasts[self$siblings == 2], pars)
sum(log(self$n.points*all.palm.intensities))
}
))
}
######
## Class for Poisson number of children.
######
set.poischild.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("poischild.inherit", class, envir = class.env)
R6Class("palm_poischild",
inherit = class.env$poischild.inherit,
public = list(
## Adding lambda parameter.
fetch.pars = function(){
super$fetch.pars()
self$add.pars(name = "lambda", link = log, start = sqrt(self$n.points/self$vol),
lower = 0, upper = self$n.points)
},
## Simulation method for the number of children per parent.
simulate.n.children = function(n, pars){
list(n.children = rpois(n, pars["lambda"]))
},
## A method for the expectation of the child distribution.
child.expectation = function(pars){
pars["lambda"]
},
## A method for the variance of the child distribution.
child.variance = function(pars){
pars["lambda"]
}
))
}
######
## Class for binomial number of children.
######
set.binomchild.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("binomchild.inherit", class, envir = class.env)
R6Class("palm_binomchild",
inherit = class.env$binomchild.inherit,
public = list(
binom.size = NULL,
initialize = function(child.list, ...){
self$child.list <- child.list
self$binom.size <- child.list$size
super$initialize(...)
},
## Adding p parameter.
fetch.pars = function(){
super$fetch.pars()
p.start <- sqrt(self$n.points/self$vol)/self$binom.size
p.start <- ifelse(p.start > 1, 0.9, p.start)
p.start <- ifelse(p.start < 0, 0.1, p.start)
self$add.pars(name = "p", link = logit, start = p.start, lower = 0, upper = 1)
},
## Simulation method for the number of children per parent.
simulate.n.children = function(n, pars){
list(n.children = rbinom(n, self$binom.size, pars["p"]))
},
## A method for the expectation of the child distribution.
child.expectation = function(pars){
self$binom.size*pars["p"]
},
## A method for the variance of the child distribution.
child.variance = function(pars){
self$binom.size*pars["p"]*(1 - pars["p"])
}
))
}
######
## Class for two-camera children distribution.
######
set.twocamerachild.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("twocamerachild.inherit", class, envir = class.env)
R6Class("palm_twocamerachild",
inherit = class.env$twocamerachild.inherit,
public = list(
## Detection zone halfwidth.
twocamera.w = NULL,
## Survey area halfwidth.
twocamera.b = NULL,
## Lag between cameras.
twocamera.l = NULL,
## Mean dive-cycle duration.
twocamera.tau = NULL,
initialize = function(child.list, ...){
self$child.list <- child.list
self$twocamera.w <- child.list$twocamera.w
self$twocamera.b <- child.list$twocamera.b
self$twocamera.l <- child.list$twocamera.l
self$twocamera.tau <- child.list$twocamera.tau
super$initialize(...)
self$setup.pattern(1)
if (self$dim != 1){
stop("Analysis of two-camera surveys is only implemented for one-dimensional processes.")
}
},
## Adding kappa parameter.
fetch.pars = function(){
super$fetch.pars()
self$add.pars(name = "kappa", link = log, start = 0.1*self$twocamera.tau, lower = 0,
upper = self$twocamera.tau)
},
## Overwriting base simulation method in case D.2D is provided.
simulate.pattern = function(pars = self$par.fitted){
if (any(names(pars) == "D.2D")){
which.D <- which(names(pars) == "D.2D")
pars[which.D] <- 2*pars[which.D]*self$twocamera.b
names(pars)[which.D] <- "D"
}
super$simulate.pattern(pars)
},
## Simulation method for the number of children per parent.
simulate.n.children = function(n, pars){
probs <- twocamera.probs(self$twocamera.l, self$twocamera.tau, self$twocamera.w,
self$twocamera.b, pars["kappa"], pars["sigma"])
## Generating some y values.
parent.ys <- runif(n, -self$twocamera.b, self$twocamera.b)
child.ys <- matrix(nrow = n, ncol = 2)
for (i in 1:n){
child.ys[i, ] <- rnorm(2, parent.ys[i], pars["sigma"])
}
camera1.in <- ifelse(child.ys[, 1] < self$twocamera.w & child.ys[, 1] > -self$twocamera.w,
TRUE, FALSE)
camera2.in <- ifelse(child.ys[, 2] < self$twocamera.w & child.ys[, 2] > -self$twocamera.w,
TRUE, FALSE)
camera1.up <- sample(c(TRUE, FALSE), n, replace = TRUE, prob = c(probs$p.up, probs$p.down))
camera2.up <- logical(n)
for (i in 1:n){
p.up <- ifelse(camera1.up[i], 1 - probs$p.down.up, probs$p.up.down)
camera2.up[i] <- sample(c(TRUE, FALSE), 1, prob = c(p.up, 1 - p.up))
}
camera1.det <- camera1.in & camera1.up
camera2.det <- camera2.in & camera2.up
child.ys <- t(child.ys)[t(cbind(camera1.det, camera2.det))]
n.children <- camera1.det + camera2.det
cameras <- numeric(sum(n.children))
j <- 1
for (i in 1:n){
if (n.children[i] > 0){
cameras[j:(j + n.children[i] - 1)] <- c(1[camera1.det[i]], 2[camera2.det[i]])
}
j <- j + n.children[i]
}
sibling.list <- siblings.twocamera(cameras)
sibling.list$cameras <- cameras
list(n.children = n.children, sibling.list = sibling.list, child.ys = child.ys)
},
## Overwriting the method to trim the sibling list for children outside the window.
trim.siblings = function(sibling.list, trimmed){
sibling.list$sibling.mat <- sibling.list$sibling.mat[trimmed, trimmed, drop = FALSE]
sibling.list$cameras <- sibling.list$cameras[trimmed]
super$trim.siblings(sibling.list, trimmed)
},
## A method for the expectation of the child distribution.
child.expectation = function(pars){
probs <- twocamera.probs(self$twocamera.l, self$twocamera.tau, self$twocamera.w,
self$twocamera.b, pars["kappa"], pars["sigma"])
2*probs$p.10/(probs$p.10 + probs$p.01)
},
## A method for the variance of the child distribution.
child.variance = function(pars){
probs <- twocamera.probs(self$twocamera.l, self$twocamera.tau, self$twocamera.w,
self$twocamera.b, pars["kappa"], pars["sigma"])
2*probs$p.10*probs$p.01*(2 - probs$p.10 - probs$p.01)/(probs$p.10 + probs$p.01)^2
}
))
}
######
## Class for Thomas processes.
######
set.thomas.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("thomas.inherit", class, envir = class.env)
R6Class("palm_thomas",
inherit = class.env$thomas.inherit,
public = list(
## Adding sigma paremter.
fetch.pars = function(){
super$fetch.pars()
self$add.pars(name = "sigma", link = log, start = 0.1*self$R, lower = 0,
upper = self$R)
},
## Simulation method for children locations.
simulate.children = function(n, parent.loc, pars){
rmvnorm(n, parent.loc, pars["sigma"]^2*diag(self$dim))
},
## Overwriting method for the integral in the Palm likelihood.
palm.likelihood.integral = function(pars){
-self$pi.multiplier*(pars["D"]*self$child.expectation(pars)*Vd(self$R, self$dim) +
self$sibling.expectation(pars)*self$Fq(self$R, pars))
},
## Overwriting method for the PDF of Q.
fq = function(r, pars){
2^(1 - self$dim/2)*r^(self$dim - 1)*exp(-r^2/(4*pars["sigma"]^2))/
((pars["sigma"]*sqrt(2))^self$dim*gamma(self$dim/2))
},
## Overwriting method for the CDF of Q.
Fq = function(r, pars){
pgamma(r^2/(4*pars["sigma"]^2), self$dim/2)
},
## Overwriting method for the quotient of the PDF of Q and the surface volume.
fq.over.s = function(r, pars){
exp(-r^2/(4*pars["sigma"]^2))/((2*pars["sigma"])^self$dim*pi^(self$dim/2))
},
## A method for xlim.
get.xlim = function(){
c(0, min(self$R, 10*self$par.fitted["sigma"]))
}
))
}
######
## Class for Matern processes.
######
set.matern.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("matern.inherit", class, envir = class.env)
R6Class("palm_matern",
inherit = class.env$matern.inherit,
public = list(
## Adding tau parameter.
fetch.pars = function(){
super$fetch.pars()
self$add.pars(name = "tau", link = log, start = 0.1*self$R, lower = 0, upper = self$R)
},
## Simulation method for children locations via rejection.
simulate.children = function(n, parent.loc, pars){
out <- matrix(0, nrow = n, ncol = self$dim)
for (i in 1:n){
keep <- FALSE
while (!keep){
proposal <- runif(self$dim, -pars["tau"], pars["tau"])
if (sqrt(sum(proposal^2)) <= pars["tau"]){
proposal <- proposal + parent.loc
out[i, ] <- proposal
keep <- TRUE
}
}
}
out
},
## Overwriting method for the PDF of Q.
fq = function(r, pars){
ifelse(r > 2*pars["tau"], 0,
2*self$dim*r^(self$dim - 1)*(pars["tau"]*hyperg_2F1(0.5, 0.5 - self$dim/2, 1.5, 1) -
r/2*hyperg_2F1(0.5, 0.5 - self$dim/2, 1.5, r^2/(4*pars["tau"]^2)))/
(beta(self$dim/2 + 0.5, 0.5)*pars["tau"]^(self$dim + 1)))
},
## Overwriting method for the CDF of Q.
Fq = function(r, pars){
alpha <- r^2/(4*pars["tau"]^2)
r^2/pars["tau"]^2*(1 - pbeta(alpha, 0.5, self$dim/2 + 0.5)) +
2^self$dim*incomplete.beta(alpha, self$dim/2 + 0.5, self$dim/2 + 0.5)/
beta(0.5, self$dim/2 + 0.5)
},
## Overwriting method for the quotient of the PDF of Q and the surface volume.
fq.over.s = function(r, pars){
ifelse(r > 2*pars["tau"], 0,
2*(pars["tau"]*hyperg_2F1(0.5, 0.5 - self$dim/2, 1.5, 1) -
r/2*hyperg_2F1(0.5, 0.5 - self$dim/2, 1.5, r^2/(4*pars["tau"]^2)))*gamma(self$dim/2 + 1)/
(beta(self$dim/2 + 0.5, 0.5)*pars["tau"]^(self$dim + 1)*pi^(self$dim/2)))
},
## A method for xlim.
get.xlim = function(){
c(0, min(self$R, 10*self$par.fitted["tau"]))
}
))
}
######
## Class for void processes.
######
set.void.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("void.inherit", class, envir = class.env)
R6Class("palm_void",
inherit = class.env$void.inherit,
public = list(
## Adding children and parent density parameters.
fetch.pars = function(){
super$fetch.pars()
self$add.pars(name = "Dp", link = log, start = 10/self$vol, lower = 0, upper = Inf)
self$add.pars(name = "Dc", link = log, start = self$n.points/self$vol, lower = 0, upper = Inf)
},
## Overwriting simulation method.
simulate.pattern = function(pars = self$par.fitted){
## Generating children.
expected.children <- pars["Dc"]*self$vol
n.children <- rpois(1, expected.children)
child.locs <- matrix(0, nrow = n.children, ncol = self$dim)
for (i in 1:self$dim){
child.locs[, i] <- runif(n.children, self$lims[i, 1], self$lims[i, 2])
}
## Generating parents.
parent.locs <- self$get.parents(pars)
list(points = self$delete.points(child.locs, parent.locs, pars), parents = parent.locs)
},
## A method to get parents by simulation.
get.parents = function(pars){
expected.parents <- pars["Dp"]*self$vol
n.parents <- rpois(1, expected.parents)
parent.locs <- matrix(0, nrow = n.parents, ncol = self$dim)
for (i in 1:self$dim){
parent.locs[, i] <- runif(n.parents, self$lims[i, 1], self$lims[i, 2])
}
parent.locs
},
## Overwriting method for the Palm intensity.
palm.intensity = function(r, pars){
pars["Dc"]*self$prob.safe(r, pars)
}
))
}
set.totaldeletion.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("totaldeletion.inherit", class, envir = class.env)
R6Class("palm_totaldeletion",
inherit = class.env$totaldeletion.inherit,
public = list(
fetch.pars = function(){
super$fetch.pars()
self$add.pars(name = "tau", link = log, start = 0.1*self$R, lower = 0, upper = self$R)
},
## Probability of being safe, given distance r from a known point.
prob.safe = function(r, pars){
exp(-pars["Dp"]*Vd(pars["tau"], self$dim)*(1 - pbeta(1 - (r/(2*pars["tau"])), (self$dim + 1)/2, 0.5)))
},
## Function to delete children, given children and parent locations.
delete.points = function(child.locs, parent.locs, pars){
dists <- euc_distances(child.locs[, 1], child.locs[, 2], parent.locs[, 1], parent.locs[, 2])
child.locs[apply(dists, 1, min) > pars["tau"], ]
},
## A method for xlim.
get.xlim = function(){
c(0, min(self$R, 10*self$par.fitted["tau"]))
}
))
}
######
## Class for simulating given known parent locations.
######
set.giveparent.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("giveparent.inherit", class, envir = class.env)
R6Class("palm_giveparent",
inherit = class.env$giveparent.inherit,
public = list(
parent.locs = NULL,
initialize = function(parent.locs, ...){
self$parent.locs <- parent.locs
super$initialize(...)
},
## Overwriting method for getting parents.
get.parents = function(pars){
if (!(ncol(self$parent.locs) == self$dim)){
stop("Incorrection dimensions of parent locations.")
}
self$parent.locs
}
))
}
## Function to create R6 object with correct class hierarchy.
create.obj <- function(classes, points, lims, R, child.list, parent.locs, sibling.list, trace, start, bounds){
class <- base.class.R6
n.classes <- length(classes)
class.env <- new.env()
for (i in 1:n.classes){
set.class <- get(paste("set", classes[i], "class", sep = "."))
class <- set.class(class, class.env)
}
if (any(classes == "twocamerachild") & !any(classes == "thomas")){
stop("Analysis of two-camera surveys is only implemented for Thomas processes.")
}
if (!is.null(points) & !(is.matrix(points) | is.list(points))){
stop("The argument 'points' must be a matrix, or a list of matrices.")
}
if (is.matrix(points)){
points <- list(points)
sibling.list <- list(sibling.list)
}
if (is.matrix(lims)){
lims <- list(lims)
}
class$new(points.list = points, lims.list = lims, R = R,
child.list = child.list, parent.locs = parent.locs,
sibling.list = sibling.list, trace = trace,
classes = classes, start = start, bounds = bounds)
}
## Some objects to get around R CMD check.
super <- NULL
self <- NULL
Try the palm package in your browser
Any scripts or data that you put into this service are public.
palm documentation built on Sept. 22, 2023, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions create.obj set.giveparent.class set.totaldeletion.class set.void.class set.matern.class set.thomas.class set.twocamerachild.class set.binomchild.class set.poischild.class set.sibling.class set.ns.class set.buffer.class set.pbc.class set.nlminb.class set.bobyqa.class set.fit.class set.CLASSNAME.class
######
## General base class.
######
base.class.R6 <- R6Class("palm",
public = list(
## Setting fields.
n.patterns = NULL,
lims = NULL,
lims.list = NULL,
vol = NULL,
vol.list = NULL,
dim = NULL,
dim.list = NULL,
classes = NULL,
## Initialisation method.
initialize = function(lims.list, classes, ...){
self$lims.list <- lims.list
self$vol.list <- lapply(self$lims.list,
function(x) prod(apply(x, 1, diff)))
self$dim.list <- lapply(self$lims.list,
function(x) nrow(x))
self$n.patterns <- length(lims.list)
self$classes <- classes
},
## A method to set up a new pattern.
setup.pattern = function(pattern){
if (pattern > self$n.patterns){
stop("Pattern index exceeds the number of patterns")
}
self$lims <- self$lims.list[[pattern]]
self$vol <- self$vol.list[[pattern]]
self$dim <- self$dim.list[[pattern]]
},
## A method to simulate multiple patterns.
simulate = function(pars = self$par.fitted){
out <- vector(mode = "list", length = self$n.patterns)
for (i in 1:self$n.patterns){
self$setup.pattern(i)
out[[i]] <- self$simulate.pattern(pars)
}
if (self$n.patterns == 1){
out <- out[[1]]
}
out
},
## A method to simulate a single pattern.
simulate.pattern = function(pars){},
## A method to trim points to the observation window.
trim.points = function(points, output.indices = FALSE){
in.window <- rep(TRUE, nrow(points))
for (i in 1:self$dim){
in.window <- in.window & (self$lims[i, 1] <= points[, i] & self$lims[i, 2] >= points[, i])
}
if (output.indices){
out <- which(in.window)
} else {
out <- points[in.window, , drop = FALSE]
}
out
}
))
######
## Template for new classes.
######
## Replace CLASSNAME with class name, then add fields and methods.
set.CLASSNAME.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("CLASSNAME.inherit", class, envir = class.env)
R6Class("palm_CLASSNAME",
inherit = class.env$CLASSNAME.inherit,
public = list(
))
}
######
## Class for models fitted to data.
######
set.fit.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("fit.inherit", class, envir = class.env)
R6Class("palm_fit",
inherit = class.env$fit.inherit,
public = list(
points.list = NULL,
points = NULL,
n.points = NULL,
n.points.list = NULL,
contrasts = NULL,
contrasts.list = NULL,
n.contrasts = NULL,
n.contrasts.list = NULL,
contrast.pairs = NULL,
contrast.pairs.list = NULL,
R = NULL,
boots = NULL,
trace = NULL,
conv.code = NULL,
converged = NULL,
n.par = NULL,
set.start = NULL,
set.bounds = NULL,
par.names = NULL,
fixed.names = NULL,
par.names.link = NULL,
par.start = NULL,
par.start.link = NULL,
par.fixed = NULL,
par.fixed.link = NULL,
par.links = NULL,
par.invlinks = NULL,
par.fitted = NULL,
par.fitted.link = NULL,
par.lower = NULL,
par.lower.link = NULL,
par.upper = NULL,
par.upper.link = NULL,
pi.multiplier = NULL,
## Initialisation method.
initialize = function(points.list, R, trace, start, bounds, ...){
super$initialize(...)
self$points.list <- points.list
## Checking list of points and list of lims have the same length.
if (length(self$points.list) != length(self$lims.list)){
stop("The list 'points' must have the same number of components as 'lims'.")
}
self$R <- R
## Creating numbers of points and contrasts for each pattern.
self$n.points.list <- lapply(points.list, nrow)
self$contrasts.list <- list()
self$contrast.pairs.list <- list()
self$n.contrasts.list <- list()
for (i in 1:self$n.patterns){
self$get.contrasts(i)
}
## Starting out with the first pattern for start values and such.
self$setup.pattern(1)
self$trace <- trace
self$set.start <- start
self$set.bounds <- bounds
self$get.pars()
self$n.par <- length(self$par.start)
self$par.start.link <- self$link.pars(self$par.start)
self$get.link.bounds()
},
## Overwriting the method to set up a new pattern.
setup.pattern = function(pattern, do.contrasts = TRUE){
super$setup.pattern(pattern)
self$points <- self$points.list[[pattern]]
self$n.points <- self$n.points.list[[pattern]]
if (do.contrasts){
self$contrast.pairs <- self$contrast.pairs.list[[pattern]]
self$contrasts <- self$contrasts.list[[pattern]]
self$n.contrasts <- self$n.contrasts.list[[pattern]]
}
},
## An empty method for getting contrasts.
get.contrasts = function(pattern){
self$n.contrasts.list[[pattern]] <- length(self$contrasts.list[[pattern]])
},
## A method for getting the parameters across all classes.
get.pars = function(){
self$fetch.pars()
},
## An empty method for fetching parameters from a particular class.
fetch.pars = function(){
},
## A method for adding new parameters.
add.pars = function(name, link.name = NULL, link, invlink = NULL,
start, lower, upper){
self$par.names <- c(self$par.names, name)
## Sorting out inverse link.
if (identical(link, log)){
full.link <- function(pars) log(pars[name])
## Default name and inverse for log link.
if (is.null(link.name)){
link.name <- paste("log", name, sep = ".")
}
if (is.null(invlink)){
full.invlink <- function(pars) exp(pars[link.name])
}
} else if (identical(link, logit)){
full.link <- function(pars) link(pars[name])
## Default name and inverse for logit link.
if (is.null(link.name)){
link.name <- paste("logit", name, sep = ".")
}
if (is.null(invlink)){
full.invlink <- function(pars) invlogit(pars[link.name])
}
} else {
if (is.null(link.name)){
stop("Please provide link name.")
}
if (is.null(invlink)){
stop("Please provide inverse link function.")
}
full.link <- link
full.invlink <- invlink
}
self$par.links <- c(self$par.links, full.link)
names(self$par.links) <- self$par.names
self$par.invlinks <- c(self$par.invlinks, full.invlink)
names(self$par.invlinks) <- self$par.names
self$par.names.link <- c(self$par.names.link, link.name)
## Overwriting start parameter, if one provided.
if (any(name == names(self$set.start))){
start <- self$set.start[name]
}
if (any(name == names(self$set.bounds))){
lower <- self$set.bounds[[name]][1]
upper <- self$set.bounds[[name]][2]
}
self$par.upper <- c(self$par.upper, upper)
names(self$par.upper) <- self$par.names
self$par.lower <- c(self$par.lower, lower)
## Adjustment for start values on the bounds.
if (start >= upper){
start <- lower + 0.9*(upper - lower)
} else if (start <= lower){
start <- lower + 0.9*(upper - lower)
}
self$par.start <- c(self$par.start, start)
names(self$par.start) <- self$par.names
names(self$par.lower) <- self$par.names
},
## A method to set the upper and lower parameter bounds on the link scale.
get.link.bounds = function(){
self$par.lower.link <- self$link.pars(self$par.lower)
self$par.upper.link <- self$link.pars(self$par.upper)
},
## A method for converting parameters to their link scales.
link.pars = function(pars){
out <- numeric(self$n.par)
for (i in 1:self$n.par){
out[i] <- self$par.links[[i]](pars)
}
names(out) <- self$par.names.link
out
},
## A method for converting parameters to their real scale.
invlink.pars = function(pars){
out <- numeric(self$n.par)
for (i in 1:self$n.par){
out[i] <- self$par.invlinks[[i]](pars)
}
names(out) <- self$par.names
out
},
## A method for the objective function.
obj.fun = function(link.pars, fixed.link.pars = NULL,
est.names = NULL, fixed.names = NULL){
combined.pars <- c(link.pars, fixed.link.pars)
names(combined.pars) <- c(est.names, fixed.names)
link.pars <- combined.pars[self$par.names.link]
pars <- self$invlink.pars(link.pars)
obj.fun.components <- numeric(self$n.patterns)
## Summing over all the patterns.
for (i in 1:self$n.patterns){
self$setup.pattern(i)
obj.fun.components[i] <- self$neg.log.palm.likelihood(pars)
}
out <- sum(obj.fun.components)
ll <- -out
if (self$trace){
for (i in 1:self$n.par){
cat(self$par.names[i], ": ", sep = "")
cat(pars[i], ", ", sep = "")
}
cat("Log-lik: ", ll, "\n", sep = "")
}
out
},
## An empty method for the Palm intensity.
palm.intensity = function(r, pars){},
## An empty method for the sum of the log Palm intensities.
sum.log.intensities = function(pars){
sum(log(self$pi.multiplier*self$palm.intensity(self$contrasts, pars)))
},
## A default method for the integral in the Palm likelihood.
palm.likelihood.integral = function(pars){
f <- function(r, pars){
self$palm.intensity(r, pars)*Sd(r, self$dim)
}
-self$pi.multiplier*integrate(f, lower = 0, upper = self$R, pars = pars)$value
},
## A default method for the log of the Palm likelihood function.
log.palm.likelihood = function(pars){
self$sum.log.intensities(pars) + self$palm.likelihood.integral(pars)
},
## A method for the negative Palm likelihood function.
neg.log.palm.likelihood = function(pars){
-self$log.palm.likelihood(pars)
},
## An empty method for optimisation.
palm.optimise = function(){},
## A method for non-optimisation.
palm.not.optimise = function(){
self$converged <- FALSE
self$par.fitted.link <- rep(NA, self$n.par)
names(self$par.fitted.link) <- self$par.names.link
self$par.fitted <- rep(NA, self$n.par)
self$conv.code <- 5
},
## A method for model fitting.
fit = function(){
## If there are contrasts, fit the model. Otherwise, return something else.
if (self$n.contrasts <= 1){
warning("Not enough observed points to compute parameter estimates.")
self$palm.not.optimise()
} else {
self$palm.optimise()
}
},
## A method for bootstrapping.
boot = function(N, prog = TRUE){
boots <- matrix(0, nrow = N, ncol = self$n.par)
## Setting up progress bar.
if (prog){
pb <- txtProgressBar(min = 0, max = N, style = 3)
}
for (i in 1:N){
zero.points <- TRUE
while (zero.points){
sim.obj <- self$simulate()
if (self$n.patterns > 1){
sim.obj$points <- lapply(sim.obj, function(x) x$points)
sim.obj$sibling.list <- lapply(sim.obj, function(x) x$sibling.list)
zero.points <- sum(sapply(sim.obj$points, nrow)) == 0
} else {
zero.points <- nrow(sim.obj$points) == 0
}
}
## Doesn't matter that sibling.list is non-null below, as sibling class is not passed.
obj.boot <- create.obj(classes = self$classes, points = sim.obj$points, lims = self$lims,
R = self$R, child.list = self$child.list,
sibling.list = sim.obj$sibling.list, trace = FALSE,
start = self$par.fitted, bounds = self$bounds)
obj.boot$fit()
if (obj.boot$converged){
boots[i, ] <- obj.boot$par.fitted
} else {
boots[i, ] <- rep(NA, self$n.par)
}
## Updating progress bar.
if (prog){
setTxtProgressBar(pb, i)
}
}
cat("\n")
self$boots <- boots
colnames(self$boots) <- self$par.names
},
## A method for plotting.
plot = function(xlim = NULL, ylim = NULL, show.empirical = TRUE, breaks = 50, ...){
if (show.empirical){
emp <- self$empirical(xlim, breaks)
}
if (is.null(xlim)){
xlim <- self$get.xlim()
}
xx <- seq(xlim[1], xlim[2], length.out = 1000)
yy <- self$palm.intensity(xx, self$par.fitted)
if (is.null(ylim)){
if (show.empirical){
ylim <- c(0, max(c(emp$y, yy)))
} else {
ylim <- c(0, max(yy))
}
}
par(xaxs = "i")
plot(xx, yy, xlab = "", ylab = "", xlim = range(xx), ylim = ylim, type = "n", ...)
title(xlab = "r", ylab = "Palm intensity")
lines(xx, yy)
if (show.empirical){
lines(emp$x, emp$y, lty = "dashed")
}
},
## A method for default xlim.
get.xlim = function(){
c(0, self$R)
},
## A method for empirical Palm intensity.
empirical = function(xlim = NULL, breaks = 50){
if (is.null(xlim)){
xlim <- self$get.xlim()
}
dists <- pbc_distances(self$points, self$lims)
midpoints <- seq(0, xlim[2], length.out = breaks)
midpoints <- midpoints[-length(midpoints)]
h <- diff(midpoints[c(1, 2)])
midpoints[1] <- midpoints[1] + h/2
intensities <- numeric(length(midpoints))
for (i in 1:length(midpoints)){
halfwidth <- ifelse(i == 1, 0.25*h, 0.5*h)
n.interval <- sum(dists <= (midpoints[i] + halfwidth)) -
sum(dists <= (midpoints[i] - halfwidth))
area <- Vd(midpoints[i] + halfwidth, self$dim) - Vd(midpoints[i] - halfwidth, self$dim)
intensities[i] <- n.interval/(self$pi.multiplier*area)
}
list(x = midpoints, y = intensities)
}
))
}
######
## Class for bobyqa() optimisation.
######
set.bobyqa.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("bobyqa.inherit", class, envir = class.env)
R6Class("palm_bobyqa",
inherit = class.env$bobyqa.inherit,
public = list(
## A method to minimise the negative Palm likelihood.
palm.optimise = function(){
out <- try(bobyqa(self$par.start.link, self$obj.fun,
lower = self$par.lower.link, upper = self$par.upper.link,
fixed.link.pars = self$par.fixed.link,
est.names = self$par.names.link,
fixed.names = self$fixed.names.link,
control = list(maxfun = 2000, npt = 2*self$n.par + 1)), silent = TRUE)
if (class(out)[1] == "try-error"){
self$palm.not.optimise()
} else {
if (out$ierr > 0){
warning("Failed convergence.")
self$converged <- FALSE
} else {
self$converged <- TRUE
}
self$par.fitted.link <- out$par
names(self$par.fitted.link) <- self$par.names.link
self$par.fitted <- self$invlink.pars(self$par.fitted.link)
self$conv.code <- out$ierr
}
}
))
}
######
## Class for nlminb() optimisation.
######
set.nlminb.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("nlminb.inherit", class, envir = class.env)
R6Class("palm_nlminb",
inherit = class.env$nlminb.inherit,
public = list(
## A method to minimise the negative Palm likelihood.
palm.optimise = function(){
out <- nlminb(self$par.start.link, self$obj.fun,
fixed.link.pars = self$par.fixed.link,
est.names = self$par.names.link, fixed.names = self$fixed.names.link,
control = list(eval.max = 2000, iter.max = 1500),
lower = self$par.lower.link, upper = self$par.upper.link)
if (out$convergence > 1){
warning("Failed convergence.")
self$converged <- FALSE
} else {
self$converged <- TRUE
}
self$par.fitted.link <- out$par
names(self$par.fitted.link) <- self$par.names.link
self$par.fitted <- self$invlink.pars(self$par.fitted.link)
self$conv.code <- out$convergence
}
))
}
######
## Class for periodic boundary conditions.
######
set.pbc.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("pbc.inherit", class, envir = class.env)
R6Class("palm_pbc",
inherit = class.env$pbc.inherit,
public = list(
## Overwriting the method to set up a new pattern.
setup.pattern = function(pattern, do.contrasts = TRUE){
super$setup.pattern(pattern, do.contrasts)
self$pi.multiplier <- self$n.points/2
},
## A method to generate contrasts.
get.contrasts = function(pattern){
self$setup.pattern(pattern, do.contrasts = FALSE)
## Saving which contrast applies to which pair of observations.
contrast.pairs <- matrix(0, nrow = (self$n.points^2 - self$n.points)/2, ncol = 2)
k <- 1
if (self$n.points > 1){
for (i in 1:(self$n.points - 1)){
for (j in (i + 1):self$n.points){
contrast.pairs[k, ] <- c(i, j)
k <- k + 1
}
}
contrasts <- pbc_distances(points = self$points, lims = self$lims)
self$contrast.pairs.list[[pattern]] <- contrast.pairs[contrasts <= self$R, ]
self$contrasts.list[[pattern]] <- contrasts[contrasts <= self$R]
} else {
self$contrast.pairs.list[[pattern]] <- contrast.pairs
self$contrasts[[pattern]] <- numeric(0)
}
super$get.contrasts(pattern)
}
))
}
######
## Class for buffer-zone boundary conditions.
######
set.buffer.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("buffer.inherit", class, envir = class.env)
R6Class("palm_buffer",
inherit = class.env$buffer.inherit,
public = list(
## Overwriting the method to set up a new pattern.
setup.pattern = function(pattern, do.contrasts = TRUE){
super$setup.pattern(pattern, do.contrasts)
self$pi.multiplier <- self$n.points
},
## A method to generate contrasts.
get.contrasts = function(pattern){
self$setup.pattern(pattern, do.contrasts = FALSE)
## Getting rid of contrasts between two external points.
buffer.keep <- buffer_keep(points = self$points, lims = self$lims,
R = self$R)
contrasts <- as.vector(as.matrix(dist(self$points))[buffer.keep])
contrast.pairs.1 <- matrix(rep(1:self$n.points, self$n.points), nrow = self$n.points)
contrast.pairs.2 <- matrix(rep(1:self$n.points, self$n.points), nrow = self$n.points, byrow = TRUE)
contrast.pairs.1 <- as.vector(contrast.pairs.1[buffer.keep])
contrast.pairs.2 <- as.vector(contrast.pairs.2[buffer.keep])
contrast.pairs <- cbind(contrast.pairs.1, contrast.pairs.2)
## Now truncating to contrasts less than R.
self$contrasts.list[[pattern]] <- contrasts[contrasts <= self$R]
self$contrast.pairs.list[[pattern]] <- contrast.pairs[contrasts <= self$R, ]
super$get.contrasts(pattern)
}
))
}
######
## Class for Neyman-Scott processes.
######
set.ns.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("ns.inherit", class, envir = class.env)
R6Class("palm_ns",
inherit = class.env$ns.inherit,
public = list(
child.list = NULL,
## Adding density parameter.
fetch.pars = function(){
super$fetch.pars()
link <- function(pars){
log(pars["D"]*self$child.expectation(pars))
}
## An inverse link that uses the other unlinked paramters.
invlink <- function(pars, out){
exp(pars["log.Dc"])/self$child.expectation(out)
}
self$add.pars(name = "D", link.name = "log.Dc", link = link, invlink = invlink,
start = sqrt(self$n.points/self$vol), lower = 0, upper = Inf)
},
## Overwriting the method for converting parameters to their real scale.
invlink.pars = function(pars){
out <- numeric(self$n.par)
names(out) <- self$par.names
which.D <- which(self$par.names == "D")
for (i in (1:self$n.par)[-which.D]){
out[i] <- self$par.invlinks[[i]](pars)
}
out[which.D] <- self$par.invlinks[[which.D]](pars, out)
out
},
## Overwriting simulation methods.
simulate.pattern = function(pars = self$par.fitted){
parent.locs <- self$get.parents(pars)
n.parents <- nrow(parent.locs)
sim.n.children <- self$simulate.n.children(n.parents, pars)
n.children <- sim.n.children$n.children
sibling.list <- sim.n.children$sibling.list
child.locs <- matrix(0, nrow = sum(n.children), ncol = self$dim)
j <- 0
for (i in seq_along(n.children)){
if (n.children[i] > 0){
child.locs[(j + 1):(j + n.children[i]), ] <-
self$simulate.children(n.children[i], parent.locs[i, ], pars)
j <- j + n.children[i]
}
}
parent.ids <- rep(seq_along(n.children), times = n.children)
trimmed <- self$trim.points(child.locs, output.indices = TRUE)
list(points = child.locs[trimmed, , drop = FALSE],
parents = parent.locs,
parent.ids = parent.ids[trimmed],
child.ys = sim.n.children$child.ys[trimmed],
sibling.list = self$trim.siblings(sibling.list, trimmed))
},
## A method to trim the sibling list.
trim.siblings = function(sibling.list, trimmed){
sibling.list
},
## A method to get parents by simulation.
get.parents = function(pars){
expected.parents <- pars["D"]*self$vol
n.parents <- rpois(1, expected.parents)
parent.locs <- matrix(0, nrow = n.parents, ncol = self$dim)
for (i in 1:self$dim){
parent.locs[, i] <- runif(n.parents, self$lims[i, 1], self$lims[i, 2])
}
parent.locs
},
## Overwriting method for the Palm intensity.
palm.intensity = function(r, pars){
self$sibling.pi(r, pars) + self$nonsibling.pi(pars)
},
## An empty method for the expectation of the child distribution.
child.expectation = function(pars){},
## An empty method for the variance of the child distribution.
child.variance = function(pars){},
## A method for the expected number of siblings of a randomly chosen point.
sibling.expectation = function(pars){
(self$child.variance(pars) + self$child.expectation(pars)^2)/self$child.expectation(pars) - 1
},
## An empty method for the PDF of Q, the between-sibling distances.
fq = function(r, pars){},
## An empty method for the CDF of Q.
Fq = function(r, pars){},
## A default method for the quotient of the PDF of Q and the surface volume.
## Numerically unstable; better to replace in child classes.
fq.over.s = function(r, pars){
self$fq(r, pars)/Sd(r, dim)
},
## A method for the Palm intensity of nonsibling points.
nonsibling.pi = function(pars){
pars["D"]*self$child.expectation(pars)
},
## A method for the PI of sibling points.
sibling.pi = function(r, pars){
self$sibling.expectation(pars)*self$fq.over.s(r, pars)
}
))
}
######
## Class for known sibling relationships.
######
set.sibling.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("sibling.inherit", class, envir = class.env)
R6Class("palm_sibling",
inherit = class.env$sibling.inherit,
public = list(
siblings = NULL,
## Lol one day you'll get a bug because you're an
## idiot and you have objects called 'sibling.list'
## and 'siblings.list' and you'll be so mad.
siblings.list = NULL,
sibling.list = NULL,
sibling.mat = NULL,
sibling.alpha = NULL,
sibling.beta = NULL,
## Initialisation method.
initialize = function(sibling.list, ...){
self$sibling.list <- sibling.list
super$initialize(...)
self$siblings.list <- list()
for (i in 1:self$n.patterns){
self$get.siblings(i)
}
},
## Overwriting the method to set up a new pattern.
setup.pattern = function(pattern, do.contrasts = TRUE, do.siblings = TRUE){
super$setup.pattern(pattern, do.contrasts)
self$sibling.mat <- self$sibling.list[[pattern]]$sibling.mat
self$sibling.alpha <- self$sibling.list[[pattern]]$alpha
self$sibling.beta <- self$sibling.list[[pattern]]$beta
if (do.siblings){
self$siblings <- self$siblings.list[[pattern]]
}
},
## A method to get vector of sibling relationships
## that matches with the contrasts.
get.siblings = function(pattern){
self$setup.pattern(pattern, do.siblings = FALSE)
siblings <- numeric(self$n.contrasts)
for (i in 1:self$n.contrasts){
pair <- self$contrast.pairs[i, ]
siblings[i] <- self$sibling.mat[pair[1], pair[2]]
}
## 0 for known nonsiblings.
## 1 for known siblings.
## 2 for unknown sibling status.
siblings <- as.numeric(siblings)
siblings[is.na(siblings)] <- 2
self$siblings.list[[pattern]] <- siblings
},
## Overwriting the sum of the log intensities.
sum.log.intensities = function(pars){
all.palm.intensities <- numeric(self$n.contrasts)
all.palm.intensities[self$siblings == 0] <-
self$sibling.beta*self$nonsibling.pi(pars)
all.palm.intensities[self$siblings == 1] <-
self$sibling.alpha*
self$sibling.pi(self$contrasts[self$siblings == 1], pars)
all.palm.intensities[self$siblings == 2] <-
(1 - self$sibling.beta)*self$nonsibling.pi(pars) +
(1 - self$sibling.alpha)*
self$sibling.pi(self$contrasts[self$siblings == 2], pars)
sum(log(self$n.points*all.palm.intensities))
}
))
}
######
## Class for Poisson number of children.
######
set.poischild.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("poischild.inherit", class, envir = class.env)
R6Class("palm_poischild",
inherit = class.env$poischild.inherit,
public = list(
## Adding lambda parameter.
fetch.pars = function(){
super$fetch.pars()
self$add.pars(name = "lambda", link = log, start = sqrt(self$n.points/self$vol),
lower = 0, upper = self$n.points)
},
## Simulation method for the number of children per parent.
simulate.n.children = function(n, pars){
list(n.children = rpois(n, pars["lambda"]))
},
## A method for the expectation of the child distribution.
child.expectation = function(pars){
pars["lambda"]
},
## A method for the variance of the child distribution.
child.variance = function(pars){
pars["lambda"]
}
))
}
######
## Class for binomial number of children.
######
set.binomchild.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("binomchild.inherit", class, envir = class.env)
R6Class("palm_binomchild",
inherit = class.env$binomchild.inherit,
public = list(
binom.size = NULL,
initialize = function(child.list, ...){
self$child.list <- child.list
self$binom.size <- child.list$size
super$initialize(...)
},
## Adding p parameter.
fetch.pars = function(){
super$fetch.pars()
p.start <- sqrt(self$n.points/self$vol)/self$binom.size
p.start <- ifelse(p.start > 1, 0.9, p.start)
p.start <- ifelse(p.start < 0, 0.1, p.start)
self$add.pars(name = "p", link = logit, start = p.start, lower = 0, upper = 1)
},
## Simulation method for the number of children per parent.
simulate.n.children = function(n, pars){
list(n.children = rbinom(n, self$binom.size, pars["p"]))
},
## A method for the expectation of the child distribution.
child.expectation = function(pars){
self$binom.size*pars["p"]
},
## A method for the variance of the child distribution.
child.variance = function(pars){
self$binom.size*pars["p"]*(1 - pars["p"])
}
))
}
######
## Class for two-camera children distribution.
######
set.twocamerachild.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("twocamerachild.inherit", class, envir = class.env)
R6Class("palm_twocamerachild",
inherit = class.env$twocamerachild.inherit,
public = list(
## Detection zone halfwidth.
twocamera.w = NULL,
## Survey area halfwidth.
twocamera.b = NULL,
## Lag between cameras.
twocamera.l = NULL,
## Mean dive-cycle duration.
twocamera.tau = NULL,
initialize = function(child.list, ...){
self$child.list <- child.list
self$twocamera.w <- child.list$twocamera.w
self$twocamera.b <- child.list$twocamera.b
self$twocamera.l <- child.list$twocamera.l
self$twocamera.tau <- child.list$twocamera.tau
super$initialize(...)
self$setup.pattern(1)
if (self$dim != 1){
stop("Analysis of two-camera surveys is only implemented for one-dimensional processes.")
}
},
## Adding kappa parameter.
fetch.pars = function(){
super$fetch.pars()
self$add.pars(name = "kappa", link = log, start = 0.1*self$twocamera.tau, lower = 0,
upper = self$twocamera.tau)
},
## Overwriting base simulation method in case D.2D is provided.
simulate.pattern = function(pars = self$par.fitted){
if (any(names(pars) == "D.2D")){
which.D <- which(names(pars) == "D.2D")
pars[which.D] <- 2*pars[which.D]*self$twocamera.b
names(pars)[which.D] <- "D"
}
super$simulate.pattern(pars)
},
## Simulation method for the number of children per parent.
simulate.n.children = function(n, pars){
probs <- twocamera.probs(self$twocamera.l, self$twocamera.tau, self$twocamera.w,
self$twocamera.b, pars["kappa"], pars["sigma"])
## Generating some y values.
parent.ys <- runif(n, -self$twocamera.b, self$twocamera.b)
child.ys <- matrix(nrow = n, ncol = 2)
for (i in 1:n){
child.ys[i, ] <- rnorm(2, parent.ys[i], pars["sigma"])
}
camera1.in <- ifelse(child.ys[, 1] < self$twocamera.w & child.ys[, 1] > -self$twocamera.w,
TRUE, FALSE)
camera2.in <- ifelse(child.ys[, 2] < self$twocamera.w & child.ys[, 2] > -self$twocamera.w,
TRUE, FALSE)
camera1.up <- sample(c(TRUE, FALSE), n, replace = TRUE, prob = c(probs$p.up, probs$p.down))
camera2.up <- logical(n)
for (i in 1:n){
p.up <- ifelse(camera1.up[i], 1 - probs$p.down.up, probs$p.up.down)
camera2.up[i] <- sample(c(TRUE, FALSE), 1, prob = c(p.up, 1 - p.up))
}
camera1.det <- camera1.in & camera1.up
camera2.det <- camera2.in & camera2.up
child.ys <- t(child.ys)[t(cbind(camera1.det, camera2.det))]
n.children <- camera1.det + camera2.det
cameras <- numeric(sum(n.children))
j <- 1
for (i in 1:n){
if (n.children[i] > 0){
cameras[j:(j + n.children[i] - 1)] <- c(1[camera1.det[i]], 2[camera2.det[i]])
}
j <- j + n.children[i]
}
sibling.list <- siblings.twocamera(cameras)
sibling.list$cameras <- cameras
list(n.children = n.children, sibling.list = sibling.list, child.ys = child.ys)
},
## Overwriting the method to trim the sibling list for children outside the window.
trim.siblings = function(sibling.list, trimmed){
sibling.list$sibling.mat <- sibling.list$sibling.mat[trimmed, trimmed, drop = FALSE]
sibling.list$cameras <- sibling.list$cameras[trimmed]
super$trim.siblings(sibling.list, trimmed)
},
## A method for the expectation of the child distribution.
child.expectation = function(pars){
probs <- twocamera.probs(self$twocamera.l, self$twocamera.tau, self$twocamera.w,
self$twocamera.b, pars["kappa"], pars["sigma"])
2*probs$p.10/(probs$p.10 + probs$p.01)
},
## A method for the variance of the child distribution.
child.variance = function(pars){
probs <- twocamera.probs(self$twocamera.l, self$twocamera.tau, self$twocamera.w,
self$twocamera.b, pars["kappa"], pars["sigma"])
2*probs$p.10*probs$p.01*(2 - probs$p.10 - probs$p.01)/(probs$p.10 + probs$p.01)^2
}
))
}
######
## Class for Thomas processes.
######
set.thomas.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("thomas.inherit", class, envir = class.env)
R6Class("palm_thomas",
inherit = class.env$thomas.inherit,
public = list(
## Adding sigma paremter.
fetch.pars = function(){
super$fetch.pars()
self$add.pars(name = "sigma", link = log, start = 0.1*self$R, lower = 0,
upper = self$R)
},
## Simulation method for children locations.
simulate.children = function(n, parent.loc, pars){
rmvnorm(n, parent.loc, pars["sigma"]^2*diag(self$dim))
},
## Overwriting method for the integral in the Palm likelihood.
palm.likelihood.integral = function(pars){
-self$pi.multiplier*(pars["D"]*self$child.expectation(pars)*Vd(self$R, self$dim) +
self$sibling.expectation(pars)*self$Fq(self$R, pars))
},
## Overwriting method for the PDF of Q.
fq = function(r, pars){
2^(1 - self$dim/2)*r^(self$dim - 1)*exp(-r^2/(4*pars["sigma"]^2))/
((pars["sigma"]*sqrt(2))^self$dim*gamma(self$dim/2))
},
## Overwriting method for the CDF of Q.
Fq = function(r, pars){
pgamma(r^2/(4*pars["sigma"]^2), self$dim/2)
},
## Overwriting method for the quotient of the PDF of Q and the surface volume.
fq.over.s = function(r, pars){
exp(-r^2/(4*pars["sigma"]^2))/((2*pars["sigma"])^self$dim*pi^(self$dim/2))
},
## A method for xlim.
get.xlim = function(){
c(0, min(self$R, 10*self$par.fitted["sigma"]))
}
))
}
######
## Class for Matern processes.
######
set.matern.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("matern.inherit", class, envir = class.env)
R6Class("palm_matern",
inherit = class.env$matern.inherit,
public = list(
## Adding tau parameter.
fetch.pars = function(){
super$fetch.pars()
self$add.pars(name = "tau", link = log, start = 0.1*self$R, lower = 0, upper = self$R)
},
## Simulation method for children locations via rejection.
simulate.children = function(n, parent.loc, pars){
out <- matrix(0, nrow = n, ncol = self$dim)
for (i in 1:n){
keep <- FALSE
while (!keep){
proposal <- runif(self$dim, -pars["tau"], pars["tau"])
if (sqrt(sum(proposal^2)) <= pars["tau"]){
proposal <- proposal + parent.loc
out[i, ] <- proposal
keep <- TRUE
}
}
}
out
},
## Overwriting method for the PDF of Q.
fq = function(r, pars){
ifelse(r > 2*pars["tau"], 0,
2*self$dim*r^(self$dim - 1)*(pars["tau"]*hyperg_2F1(0.5, 0.5 - self$dim/2, 1.5, 1) -
r/2*hyperg_2F1(0.5, 0.5 - self$dim/2, 1.5, r^2/(4*pars["tau"]^2)))/
(beta(self$dim/2 + 0.5, 0.5)*pars["tau"]^(self$dim + 1)))
},
## Overwriting method for the CDF of Q.
Fq = function(r, pars){
alpha <- r^2/(4*pars["tau"]^2)
r^2/pars["tau"]^2*(1 - pbeta(alpha, 0.5, self$dim/2 + 0.5)) +
2^self$dim*incomplete.beta(alpha, self$dim/2 + 0.5, self$dim/2 + 0.5)/
beta(0.5, self$dim/2 + 0.5)
},
## Overwriting method for the quotient of the PDF of Q and the surface volume.
fq.over.s = function(r, pars){
ifelse(r > 2*pars["tau"], 0,
2*(pars["tau"]*hyperg_2F1(0.5, 0.5 - self$dim/2, 1.5, 1) -
r/2*hyperg_2F1(0.5, 0.5 - self$dim/2, 1.5, r^2/(4*pars["tau"]^2)))*gamma(self$dim/2 + 1)/
(beta(self$dim/2 + 0.5, 0.5)*pars["tau"]^(self$dim + 1)*pi^(self$dim/2)))
},
## A method for xlim.
get.xlim = function(){
c(0, min(self$R, 10*self$par.fitted["tau"]))
}
))
}
######
## Class for void processes.
######
set.void.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("void.inherit", class, envir = class.env)
R6Class("palm_void",
inherit = class.env$void.inherit,
public = list(
## Adding children and parent density parameters.
fetch.pars = function(){
super$fetch.pars()
self$add.pars(name = "Dp", link = log, start = 10/self$vol, lower = 0, upper = Inf)
self$add.pars(name = "Dc", link = log, start = self$n.points/self$vol, lower = 0, upper = Inf)
},
## Overwriting simulation method.
simulate.pattern = function(pars = self$par.fitted){
## Generating children.
expected.children <- pars["Dc"]*self$vol
n.children <- rpois(1, expected.children)
child.locs <- matrix(0, nrow = n.children, ncol = self$dim)
for (i in 1:self$dim){
child.locs[, i] <- runif(n.children, self$lims[i, 1], self$lims[i, 2])
}
## Generating parents.
parent.locs <- self$get.parents(pars)
list(points = self$delete.points(child.locs, parent.locs, pars), parents = parent.locs)
},
## A method to get parents by simulation.
get.parents = function(pars){
expected.parents <- pars["Dp"]*self$vol
n.parents <- rpois(1, expected.parents)
parent.locs <- matrix(0, nrow = n.parents, ncol = self$dim)
for (i in 1:self$dim){
parent.locs[, i] <- runif(n.parents, self$lims[i, 1], self$lims[i, 2])
}
parent.locs
},
## Overwriting method for the Palm intensity.
palm.intensity = function(r, pars){
pars["Dc"]*self$prob.safe(r, pars)
}
))
}
set.totaldeletion.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("totaldeletion.inherit", class, envir = class.env)
R6Class("palm_totaldeletion",
inherit = class.env$totaldeletion.inherit,
public = list(
fetch.pars = function(){
super$fetch.pars()
self$add.pars(name = "tau", link = log, start = 0.1*self$R, lower = 0, upper = self$R)
},
## Probability of being safe, given distance r from a known point.
prob.safe = function(r, pars){
exp(-pars["Dp"]*Vd(pars["tau"], self$dim)*(1 - pbeta(1 - (r/(2*pars["tau"])), (self$dim + 1)/2, 0.5)))
},
## Function to delete children, given children and parent locations.
delete.points = function(child.locs, parent.locs, pars){
dists <- euc_distances(child.locs[, 1], child.locs[, 2], parent.locs[, 1], parent.locs[, 2])
child.locs[apply(dists, 1, min) > pars["tau"], ]
},
## A method for xlim.
get.xlim = function(){
c(0, min(self$R, 10*self$par.fitted["tau"]))
}
))
}
######
## Class for simulating given known parent locations.
######
set.giveparent.class <- function(class, class.env){
## Saving inherited class to class.env.
assign("giveparent.inherit", class, envir = class.env)
R6Class("palm_giveparent",
inherit = class.env$giveparent.inherit,
public = list(
parent.locs = NULL,
initialize = function(parent.locs, ...){
self$parent.locs <- parent.locs
super$initialize(...)
},
## Overwriting method for getting parents.
get.parents = function(pars){
if (!(ncol(self$parent.locs) == self$dim)){
stop("Incorrection dimensions of parent locations.")
}
self$parent.locs
}
))
}
## Function to create R6 object with correct class hierarchy.
create.obj <- function(classes, points, lims, R, child.list, parent.locs, sibling.list, trace, start, bounds){
class <- base.class.R6
n.classes <- length(classes)
class.env <- new.env()
for (i in 1:n.classes){
set.class <- get(paste("set", classes[i], "class", sep = "."))
class <- set.class(class, class.env)
}
if (any(classes == "twocamerachild") & !any(classes == "thomas")){
stop("Analysis of two-camera surveys is only implemented for Thomas processes.")
}
if (!is.null(points) & !(is.matrix(points) | is.list(points))){
stop("The argument 'points' must be a matrix, or a list of matrices.")
}
if (is.matrix(points)){
points <- list(points)
sibling.list <- list(sibling.list)
}
if (is.matrix(lims)){
lims <- list(lims)
}
class$new(points.list = points, lims.list = lims, R = R,
child.list = child.list, parent.locs = parent.locs,
sibling.list = sibling.list, trace = trace,
classes = classes, start = start, bounds = bounds)
}
## Some objects to get around R CMD check.
super <- NULL
self <- NULL
Try the palm package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.