Skip to content

Commit

Permalink
Improvement to unnormdensity
Browse files Browse the repository at this point in the history
  • Loading branch information
@baddstats
baddstats committed Mar 5, 2023
1 parent 3a0cacf commit 40d6bdd
Show file tree
Hide file tree
Showing 4 changed files with 149 additions and 95 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: spatstat.geom
Version: 3.0-6.013
Date: 2023-03-04
Version: 3.0-6.014
Date: 2023-03-05
Title: Geometrical Functionality of the 'spatstat' Family
Authors@R: c(person("Adrian", "Baddeley",
role = c("aut", "cre", "cph"),
Expand Down
5 changes: 4 additions & 1 deletion NEWS
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
CHANGES IN spatstat.geom VERSION 3.0-6.013
CHANGES IN spatstat.geom VERSION 3.0-6.014

OVERVIEW

Expand Down Expand Up @@ -41,6 +41,9 @@ SIGNIFICANT USER-VISIBLE CHANGES
o unnormdensity
Computation accelerated.

o unnormdensity
Handles datasets containing fewer than 2 values.

BUG FIXES

o unnormdensity
Expand Down
233 changes: 142 additions & 91 deletions R/unnormdensity.R
View file
Original file line number Diff line number Diff line change
@@ -1,109 +1,160 @@
#
# unnormdensity.R
#
# $Revision: 1.17 $ $Date: 2023/03/01 13:33:17 $
# $Revision: 1.18 $ $Date: 2023/03/05 02:04:27 $
#

unnormdensity <- function(x, ..., weights=NULL, defaults=list()) {
if(any(!nzchar(names(list(...)))))
stop("All arguments must be named (tag=value)")
envir.here <- sys.frame(sys.nframe())
if(is.null(weights)) {
## all weights are 1 (not 1/n)
out <- do.call.matched(density.default,
c(list(x=quote(x), ...), defaults),
envir=envir.here)
out$y <- length(x) * out$y
} else if(length(weights) == 1) {
## all weights are equal
out <- do.call.matched(density.default,
c(list(x=quote(x), ...), defaults),
envir=envir.here)
out$y <- weights[1] * length(x) * out$y
} else if(length(weights) != length(x)) {
stop("'x' and 'weights' have unequal length")
} else {
weightrange <- range(weights)
if(all(weightrange == 0)) {
## result is zero
unnormdensity <- local({

unnormdensity <- function(x, ..., weights=NULL, defaults=list()) {
if(any(!nzchar(names(list(...)))))
stop("All arguments must be named (tag=value)")
envir.here <- sys.frame(sys.nframe())
if(length(x) <= 1) {
## density.default does not handle this
out <- do.call(fewdatacase,
resolve.defaults(
list(x=x, ..., weights=weights),
defaults))
} else if(is.null(weights)) {
## all weights are 1 (not 1/n)
out <- do.call.matched(density.default,
c(list(x=quote(x), ...), defaults),
envir=envir.here)
out$y <- 0 * out$y
} else if(all(weightrange >= 0)) {
## all masses are nonnegative, some are positive
out <- do.call.matched(density.default,
c(list(x=quote(x),
weights=quote(weights),
subdensity=TRUE,
...),
defaults),
envir=envir.here)
} else if(all(weightrange <= 0)) {
## all masses are nonpositive, some are negative
w <- (- weights)
out$y <- length(x) * out$y
} else if(length(weights) == 1) {
## all weights are equal
out <- do.call.matched(density.default,
c(list(x=quote(x),
weights=quote(w),
subdensity=TRUE,
...),
defaults),
c(list(x=quote(x), ...), defaults),
envir=envir.here)
out$y <- (- out$y)
out$y <- weights[1] * length(x) * out$y
} else if(length(weights) != length(x)) {
stop("'x' and 'weights' have unequal length")
} else {
## mixture of positive and negative masses
w <- weights
wabs <- abs(w)
wpos <- pmax.int(0, w)
wneg <- - pmin.int(0, w)
## determine bandwidth value
bw <- list(...)$bw # could be NULL
if(is.numeric(bw)) {
## bandwidth is given, as a numeric value
## adjust by factor 'adjust'
adjust <- list(...)$adjust %orifnull% 1
bw <- bw * adjust
weightrange <- range(weights)
if(all(weightrange == 0)) {
## result is zero
out <- do.call.matched(density.default,
c(list(x=quote(x), ...), defaults),
envir=envir.here)
out$y <- 0 * out$y
} else if(all(weightrange >= 0)) {
## all masses are nonnegative, some are positive
out <- do.call.matched(density.default,
c(list(x=quote(x),
weights=quote(weights),
subdensity=TRUE,
...),
defaults),
envir=envir.here)
} else if(all(weightrange <= 0)) {
## all masses are nonpositive, some are negative
w <- (- weights)
out <- do.call.matched(density.default,
c(list(x=quote(x),
weights=quote(w),
subdensity=TRUE,
...),
defaults),
envir=envir.here)
out$y <- (- out$y)
} else {
## compute bandwidth by applying a rule, using absolute masses
dabs <- do.call.matched(density.default,
c(list(x=quote(x),
weights=quote(wabs),
subdensity=TRUE,
...),
defaults),
envir=envir.here)
bw <- dabs$bw
## mixture of positive and negative masses
w <- weights
wabs <- abs(w)
wpos <- pmax.int(0, w)
wneg <- - pmin.int(0, w)
## determine bandwidth value
bw <- list(...)$bw # could be NULL
if(is.numeric(bw)) {
## bandwidth is given, as a numeric value
## adjust by factor 'adjust'
adjust <- list(...)$adjust %orifnull% 1
bw <- bw * adjust
} else {
## compute bandwidth by applying a rule, using absolute masses
dabs <- do.call.matched(density.default,
c(list(x=quote(x),
weights=quote(wabs),
subdensity=TRUE,
...),
defaults),
envir=envir.here)
bw <- dabs$bw
}
## Bandwidth is now determined as a numeric value.
## Compute densities for positive and negative masses separately
outpos <- do.call.matched(density.default,
resolve.defaults(list(x=quote(x),
bw=bw,
adjust=1,
weights=quote(wpos),
subdensity=TRUE,
...),
defaults,
.StripNull=TRUE),
envir=envir.here)
outneg <- do.call.matched(density.default,
resolve.defaults(list(x=quote(x),
bw=bw,
adjust=1,
weights=quote(wneg),
subdensity=TRUE,
...),
defaults,
.StripNull=TRUE),
envir=envir.here)
## combine
out <- outpos
out$y <- outpos$y - outneg$y
}
## Bandwidth is now determined as a numeric value.
## Compute densities for positive and negative masses separately
outpos <- do.call.matched(density.default,
resolve.defaults(list(x=quote(x),
bw=bw,
adjust=1,
weights=quote(wpos),
subdensity=TRUE,
...),
defaults,
.StripNull=TRUE),
envir=envir.here)
outneg <- do.call.matched(density.default,
resolve.defaults(list(x=quote(x),
bw=bw,
adjust=1,
weights=quote(wneg),
subdensity=TRUE,
...),
defaults,
.StripNull=TRUE),
envir=envir.here)
## combine
out <- outpos
out$y <- outpos$y - outneg$y
}
out$call <- match.call()
return(out)
}
out$call <- match.call()
return(out)
}

fewdatacase <- function(x, ...,
weights=NULL, kernel="gaussian",
bw=NULL,
from=NULL, to=NULL, n=512) {
nx <- length(x)
if(nx == 0) {
needed <- list(from=from, to=to, n=n)
} else if(nx == 1) {
needed <- list(bw=bw, from=from, to=to, n=n)
} else stop("Internal function 'fewdatacase' was invoked incorrectly")

if(any(absent <- sapply(needed, is.null))) {
nabsent <- sum(absent)
stop(paste(ngettext(nabsent, "Argument", "Arguments"),
commasep(sQuote(names(needed)[absent])),
ngettext(nabsent, "is", "are"),
"required in density.default when length(x) <= 1"),
call.=FALSE)
}

xx <- seq(from, to, length.out=n)

if(nx == 0) {
yy <- rep(0, n)
} else {
kernel <- match.kernel(kernel)
if(!is.numeric(bw)) {
## Bandwidth selection rules require >= 2 data values
## Use fallback
h <- (to-from)/3
bw <- h / kernel.factor(kernel)
}
yy <- (weights %orifnull% 1) * dkernel(xx, kernel=kernel, mean=x, sd=bw)
}
structure(list(x=xx, y=yy, bw=bw, n=n, call=match.call(),
data.name='x', has.na=FALSE),
class="density")
}

unnormdensity
})


integral.density <- function(f, domain=NULL, weight=NULL, ...) {
x <- f$x
Expand Down
2 changes: 1 addition & 1 deletion inst/doc/packagesizes.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,4 +15,4 @@ date version nhelpfiles nobjects ndatasets Rlines srclines
"2022-10-23" "3.0-3" 442 1182 0 35107 15596
"2023-01-22" "3.0-5" 445 1190 0 35297 15596
"2023-01-30" "3.0-6" 445 1190 0 35307 15596
"2023-03-04" "3.0-6.013" 449 1199 0 35643 15596
"2023-03-05" "3.0-6.014" 449 1199 0 35694 15596

0 comments on commit 40d6bdd

Please sign in to comment.