# Based on David's "New NN functions.doc" but with the papametres passed explicitly
# I've had to rename "n" as "n.points" for internal consistency
# NN estimation functions for Popas:
# ------------------------------------------
#
# Survey:
#
generate.design.no<-function (reg, n.points=2)
{
if (!is.region(reg))
stop("\n*** The parameter <reg> must be of type 'region'.\n")
if (!is.numeric(n.points))
stop("\n*** The number of points must be numeric.\n")
if (n.points < 1)
stop("\n*** You must sample at least two points (n.points>1).\n")
n <- n.points
nn.x<-runif(n,min=0,max=reg$length)
nn.y<-runif(n,min=0,max=reg$width)
des <- list(region = reg, n.points=n, x=nn.x,y=nn.y)
class(des) <- "design.no"
return(des)
}
plot.design.no<-function (des, ...)
# -------------------------------------------------------------------------------
#
# -------------------------------------------------------------------------------
{
if (!is.design.no(des))
stop("\nThe parameter <des> must be of type 'design.no'.\n")
reg<-des$region
par.was <- par(no.readonly = T)
len <- reg$length
width <- reg$width
margin <- calculate.plot.margin(reg.x = len, reg.y = width,
area.x = par.was$fin[1], area.y = par.was$fin[2])
par(mai = c(margin$y, margin$x, margin$y, margin$x))
par(new = F)
plot(0, 0, type = "n", las = 1, xlab = "", ylab = "", xlim = c(0,
len), ylim = c(0, width), xaxs = "i", yaxs = "i", xpd = T,
...)
points(des$x, des$y, pch = 3)
}
is.design.no<-function (des)
{
inherits(des, "design.no")
}
#
generate.sample.no<-function (pop, des ,with.neighbours=FALSE)
# -------------------------------------------------------------------------------------
# Returns groupIDs of objects nearest design points, and distances to design points.
# If with.neighbours==TRUE, also returns groupIDs of objects nearest these objects (thier
# nearest neighbours), subject to constraint that line from design point to object
# and thence to neighbour is >90 degrees. Also returns distances between the
# objects and nearest neighbours. If a point has no neighbours satisfying the
# constraint, the neighbour groupID is NA and the distance is 0.
# -------------------------------------------------------------------------------------
{
if (!is.population(pop))
stop("\n*** The parameter <pop> must be of type 'population'.\n")
if (!is.design.no(des))
stop("\n*** The parameter <des> must be of type 'design.no'.\n")
if (!equal(pop$region, des$region))
stop(paste("\n*** The given population and design were defined",
"with different regions.\n"))
pos.x <- pop$posx
pos.y <- pop$posy
N <- length(pos.x)
detected <- rep(0, N)
u <- des$n.points
# find distances from points to nearest objects:
point.distance <- rep(NA, u)
point.animal <- rep(NA, u)
for (i in 1:u) {
indx<-which(sqrt((des$x[i] - pos.x)^2 + (des$y[i] -
pos.y)^2) == min(sqrt((des$x[i] - pos.x)^2 + (des$y[i] -
pos.y)^2)))
point.animal[i] <- pop$groupID[indx]
point.distance[i] <- sqrt((des$x[i] - pos.x[indx])^2 + (des$y[i] -
pos.y[indx])^2)
detected[indx] <- 1
# cat("\n Animal: ", animal, " Distance: ", distance[point.animal], "\n")
}
neighbour.distance <- rep(NA, u)
neighbour.animal <- rep(NA, u)
if(with.neighbours) {
# find distances from nearest objects to nearest neighbours at >90 degrees
# from point-to-nearest-object line:
neighbour <- rep(0, N)
skip<-rep(FALSE,u)
for (i in 1:u) {
indx<-which(pop$groupID == point.animal[i])
c.sq<-(pos.x[indx] - pos.x)^2 + (pos.y[indx] - pos.y)^2
b.sq<-(pos.x[indx] - des$x[i])^2 + (pos.y[indx] - des$y[i])^2
a.sq<-(des$x[i] - pos.x)^2 + (des$y[i] - pos.y)^2
cos.A<-(b.sq+c.sq-a.sq)/(2*sqrt(b.sq*c.sq))
cos.A[indx]<-1
# exclude point itself and all points at <90 degrees (i.e. with cos.A>0):
other<-rep(TRUE,N)
other[indx]<-FALSE
other[cos.A>0]<-FALSE
if(sum(other)<1) {
skip[i]<-TRUE
neighbour.distance[i]<-0
neighbour[indx]<-NA
} else {
dist<-sqrt((pos.x[indx] - pos.x)^2 + (pos.y[indx] - pos.y)^2)
dist[!other]<- -99
n.indx<-which(dist>0 & dist == min(sqrt((pos.x[indx] - pos.x[other])^2 + (pos.y[indx] - pos.y[other])^2)))[1]
neighbour.animal[i] <- pop$groupID[n.indx]
neighbour.distance[i] <- sqrt((pos.x[indx] - pos.x[n.indx])^2 + (pos.y[indx] - pos.y[n.indx])^2)
neighbour[indx] <- n.indx
}
}
}
samp <- list(population = pop, design = des, detected = detected,
point.animal=point.animal, point.distance = point.distance, neighbour.animal=neighbour.animal, neighbour.distance = neighbour.distance)
class(samp) <- "sample.no"
return(samp)
}
is.sample.no<-function (samp)
{
inherits(samp, "sample.no")
}
summary.sample.no<-function (samp)
{
if (!is.sample.no(samp))
stop("\nThe parameter <samp> must be of type 'sample.no'.\n")
cat("\n")
cat("SAMPLE SUMMARY (POINT-TO-NEAREST-OBJECT METHOD)\n")
cat("-----------------------------------------\n\n")
cat("Number of points: ", samp$design$n.points, "\n")
cat("Mean of point-to-object distances :",mean(samp$point.distance),"\n")
cat("Mean of squared point-to-object distances :",mean(samp$point.distance^2),"\n")
if(length(!is.na(samp$neighbour.distance))>0) {
cat("Mean of object-to-neighbour distances:",mean(samp$neighbour.distance),"\n")
cat("Mean of squared object-to-neighbour distances:",mean(samp$neighbour.distance^2),"\n")
}
}
plot.sample.no<-function (samp, show.sizes = FALSE, show.exps = FALSE, dsf = 1, whole.population = FALSE,
pnlines = TRUE, nnlines = TRUE, ...)
# -------------------------------------------------------------------------------
# Only does point-to-nearest-object at present.
# -------------------------------------------------------------------------------
{
if (!is.sample.no(samp))
stop("\nThe parameter <samp> must be of type 'sample.no'.\n")
if (!is.numeric(dsf))
stop("\nThe parameter <dsf> must be numeric.\n")
if (dsf <= 0)
stop("\nThe parameter <dsf> must be positive.\n")
if ((show.sizes != TRUE) & (show.sizes != FALSE))
stop("\nThe parameter <show.sizes> must be TRUE or FALSE.\n")
if ((show.exps != TRUE) & (show.exps != FALSE))
stop("\nThe parameter <show.exps> must be TRUE or FALSE.\n")
if ((whole.population != TRUE) & (whole.population != FALSE))
stop("\nThe parameter <whole.population> must be TRUE or FALSE.\n")
pop <- samp$population
des <- samp$design
par.was <- par(no.readonly = TRUE)
len <- pop$reg$length
width <- pop$reg$width
margin <- calculate.plot.margin(reg.x = len, reg.y = width,
area.x = par.was$fin[1], area.y = par.was$fin[2])
par(mai = c(margin$y, margin$x, margin$y, margin$x))
par(new = FALSE)
plot(0, 0, type = "n", las = 1, xlab = "", ylab = "", xlim = c(0,
len), ylim = c(0, width), xaxs = "i", yaxs = "i", xpd = TRUE,
...)
if (whole.population == TRUE)
plot(pop, show.sizes = show.sizes, show.exps = show.exps,
newplot = FALSE, dsf = dsf, ...)
inside <- (samp$detected == 1)
seen <- pop
seen$groupID <- pop$groupID[inside]
seen$posx <- pop$posx[inside]
seen$posy <- pop$posy[inside]
seen$groupsize <- pop$groupsize[inside]
seen$types <- pop$types[inside]
seen$exposure <- pop$exposure[inside]
plot(seen, show.sizes = show.sizes, show.exps = show.exps,
newplot = FALSE, dsf = dsf, group.col = "red", ...)
points(des$x, des$y, pch = 3)
if (pnlines) # point-to-nearest-object lines:
for (i in 1:des$n.points) {
lines(c(des$x[i], pop$posx[which(samp$pop$groupID==samp$point.animal[i])]), c(des$y[i], pop$posy[which(samp$pop$groupID==samp$point.animal[i])]))
}
if (nnlines) # nearest-neighbour lines:
for (i in 1:des$n.points) {
indx <-which(samp$pop$groupID==samp$point.animal[i])
n.indx<-which(samp$pop$groupID==samp$neighbour.animal[i])
lines(c(pop$posx[indx], pop$posx[n.indx]), c(pop$posy[indx], pop$posy[n.indx]),col="red")
}
}
obscure.sample.no<-function(samp)
#----------------------------------------------------------------
# Removes all information about undetected animals from an object
# of class sample.no. Returns object of same class.
#----------------------------------------------------------------
{
if (!is.sample.no(samp))
stop("\n*** <samp> is not an object of type 'sample.no'.\n")
t<-samp
# mark all deteceted animals
detected<-samp$detected
# need to include neighbour animals if they are there:
if(!is.na(samp$neighbour.animal[1])) detected[is.element(samp$population$groupID,samp$neighbour.animal)]<-1
# then filter out all information about others
t$population$groupID<-samp$population$groupID[detected==1]
t$population$posx<-samp$population$posx[detected==1]
t$population$posy<-samp$population$posy[detected==1]
t$population$groupsize<-samp$population$groupsize[detected==1]
t$population$types<-samp$population$types[detected==1]
t$population$exposure<-samp$population$exposure[detected==1]
t$detected<-samp$detected[detected==1]
t$point.animal<-samp$point.animal
t$point.distance<-samp$point.distance
t$neighbour.animal<-samp$neighbour.animal
t$neighbour.distance<-samp$neighbour.distance
t
}
#=================================================
# estimate_nn
# From David's "New NN functions.doc"
#=================================================
#
point.est.no<-function(samp)
{
estimate.no(samp, ci.type = NULL)
}
int.est.no<-function(samp, vlevels = c(0.025, 0.975), ci.type = "boot.nonpar",
nboot = 999, plot = T, ... )
{
estimate.no(samp, vlevels, ci.type, nboot, plot)
}
estimate.no<-function (samp, vlevels = c(0.025, 0.975), ci.type,
nboot = 999, plot = T, ...)
# -------------------------------------------------------------------------------
# Only does point-to-nearest-object at present.
# -------------------------------------------------------------------------------
{
if (!is.sample.no(samp))
stop("\n*** <samp> is not an object of type 'sample.no'.\n")
if (!is.null (ci.type) && !(ci.type %in% c("boot.par", "boot.nonpar")))
stop(paste("\n*** Unrecognised <ci.type>. These are valid: ",
"'boot.par' and 'boot.nonpar'\n"))
if (!is.numeric(vlevels))
stop("\n*** All <vlevels> values must be numeric.\n")
if (any(vlevels < 0) | any(vlevels > 1))
stop("\n*** All <vlevels> values must be between 0 and 1.\n")
if (!is.numeric(nboot))
stop("\n*** <nboot> must be numeric.\n")
if (nboot < 1)
stop("\n*** <nboot> must be at least 1.\n")
if ((plot != T) & (plot != F))
stop("\n*** <plot> must be TRUE or FALSE.\n")
pop <- samp$population
des <- samp$design
reg.length <- des$region$length
reg.width <- des$region$width
n.points <- des$n.points
ra <- reg.length * reg.width
x1 <- samp$point.distance[!is.na(samp$point.distance)]
x2 <- samp$neighbour.distance[!is.na(samp$neighbour.distance)]
if (length(x1) != n.points) {
cat("\n Debug error: length(x)!=n.points.\n")
cat(length(x1))
cat(n.points)
stop("Execution stopped \n")
}
Nmle <- round(n.points * ra/(pi * sum(x1^2)))
if (is.null(ci.type) )
{
pointest <- list(Nhat = Nmle)
class(pointest) <- "point.est.no"
return(pointest)
}
if (ci.type == "boot.par") {
K <- pi * Nmle / ra
nb <- rep(0, nboot)
for (i in 1:nboot)
{
U <- runif(n.points)
xboot <- sqrt( -log(U) / K )
nb[i] <- round(n.points * ra/(pi * sum(xboot^2)))
}
boot.dbn <- list (Nhat.grp = nb)
Nhat.boot <- mean(nb)
boot.mean <- list(Nhat.grp = Nhat.boot)
indices <- round(vlevels * nboot)
indices <- ifelse(indices < 1, 1, indices)
sortedNhats <- sort(nb)
CI.Nhat <- sortedNhats[indices]
ci <- list(Nhat.grp = CI.Nhat)
# x <- seq(sortedNhats[1], sortedNhats[nboot], length = 100)
# y <- density(sortedNhats, kernel = "gaussian", n = 100, from = sortedNhats[1],
# to = sortedNhats[nboot])$y
# pcol <- "blue"
# text2 <- paste("\nParametric Bootstrap with", nboot,
# "replicates")
}
else if (ci.type == "boot.nonpar") {
nb <- numeric(nboot)
for (ib in 1:nboot) {
xb <- sample(x1, n.points, replace = T)
nb[ib] <- round(n.points * ra/(pi * sum(xb^2)))
}
boot.dbn <- list (Nhat.grp = nb)
Nhat.boot <- mean(nb)
boot.mean <- list(Nhat.grp = Nhat.boot)
nb <- round(sort(nb), 0)
cin <- round(nboot * vlevels, 0)
cin <- ifelse(cin < 1, 1, cin)
CI.Nhat <- nb[cin]
ci <- list(Nhat.grp = CI.Nhat)
# x <- seq(nb[1], nb[nboot], length = 100)
# y <- density(nb, kernel = "gaussian", n = 100, from = nb[1],
# to = nb[nboot])$y
# pcol <- "red"
# text2 <- paste("\nNon-parametric Bootstrap with", nboot,
# "replicates")
}
if (plot) {
if (ci.type == "boot.par") {
valid.Nhat <- nb[nb != Inf]
nvalid <- length(valid.Nhat)
main <- paste("Parametric Bootstrap Group abundance distribution:\n",
"(", nboot, "replicates)")
y <- hist(valid.Nhat, prob = T, plot = F, ...)
hist(valid.Nhat, prob = T, xlab = "Estimate", ylab = "Density",
main = main, ...)
yhi <- max(y$intensities * 0.05)
ylo <- 0
CI <- round(CI.Nhat)
segments(CI, ylo, CI, yhi, col = "blue", lwd = 2,
xpd = T)
text(CI, ylo * 1.05, paste(CI), adj = c(0.5, 1),
col = "blue", cex = 0.8, xpd = T)
text(CI, yhi * 1.5, paste(round(100 * vlevels, 1),
"%"), adj = c(0.5, 0), col = "blue", cex = 0.8,
xpd = T)
title(main)
}
if (ci.type == "boot.nonpar") {
valid.Nhat <- nb[nb != Inf]
nvalid <- length(valid.Nhat)
main <- paste("Nonparametric Bootstrap Group abundance distribution:\n",
"(", nboot, "replicates)")
y <- hist(valid.Nhat, prob = T, plot = F, ...)
hist(valid.Nhat, prob = T, xlab = "Estimate", ylab = "Density",
main = main, ...)
yhi <- max(y$intensities * 0.05)
ylo <- 0
CI <- round(CI.Nhat)
segments(CI, ylo, CI, yhi, col = "blue", lwd = 2,
xpd = T)
text(CI, ylo * 1.05, paste(CI), adj = c(0.5, 1),
col = "blue", cex = 0.8, xpd = T)
text(CI, yhi * 1.5, paste(round(100 * vlevels, 1),
"%"), adj = c(0.5, 0), col = "blue", cex = 0.8,
xpd = T)
title(main)
}
# par.was <- par(no.readonly = T)
# plot(x, y, xlab = "population size", ylab = "likelihood",
# type = "n", xaxs = "i", yaxs = "i", ...)
# lines(x, y, type = "l", col = pcol, lwd = 2)
# yhi <- max(y * 0.025)
# ylo <- -yhi
# segments(ci, ylo, ci, yhi, col = pcol, lwd = 2, xpd = T)
# text(ci, ylo * 1.05, paste(ci), adj = c(0.5, 1), col = pcol,
# cex = 0.5, xpd = T)
# text(ci, yhi * 1.5, paste(round(100 * vlevels, 1), "%"),
# adj = c(0.5, 0), col = pcol, cex = 0.5, xpd = T)
# text1 <- paste("Likelihood for N")
# title(paste(text1, text2))
# par(par.was)
}
intest <- list(levels = vlevels, ci = ci, boot.mean = boot.mean, boot.dbn = boot.dbn)
class(intest) <- "int.est.no"
return(intest)
}
is.int.est.no<-function(est)
{
inherits(est, "int.est.no")
}
is.point.est.no<-function(est)
{
inherits(est, "point.est.no")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.