-
Notifications
You must be signed in to change notification settings - Fork 1
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
10 changed files
with
511 additions
and
6 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,173 @@ | ||
| #' | ||
| #' Spatially weighted quantile | ||
| #' | ||
|
|
||
| SpatialMedian <- function(X, ...) { | ||
| UseMethod("SpatialMedian") | ||
| } | ||
|
|
||
| SpatialQuantile <- function(X, prob=0.5, ...) { | ||
| UseMethod("SpatialQuantile") | ||
| } | ||
|
|
||
| #' methods for 'ppp' class | ||
|
|
||
| SpatialMedian.ppp <- function(X, sigma=NULL, ..., | ||
| type=1, | ||
| at=c("pixels", "points"), leaveoneout=TRUE, | ||
| weights=NULL, | ||
| edge=TRUE, diggle=FALSE, verbose=FALSE) { | ||
| SpatialQuantile.ppp(X, sigma=sigma, prob=0.5, | ||
| ..., | ||
| type=type, | ||
| at=at, leaveoneout=leaveoneout, | ||
| weights=weights, | ||
| edge=edge, diggle=diggle, verbose=verbose) | ||
| } | ||
|
|
||
| SpatialQuantile.ppp <- function(X, prob=0.5, sigma=NULL, ..., | ||
| type=1, | ||
| at=c("pixels", "points"), leaveoneout=TRUE, | ||
| weights=NULL, | ||
| edge=TRUE, diggle=FALSE, verbose=FALSE) { | ||
| if(!is.ppp(X)) stop("X should be a point pattern") | ||
| if(!is.marked(X)) stop("The point pattern X should have marks") | ||
| check.1.real(prob) | ||
| stopifnot(prob >= 0) | ||
| stopifnot(prob <= 1) | ||
| at <- match.arg(at) | ||
| atName <- switch(at, pixels="pixels", points="data points") | ||
| check.1.integer(type) | ||
| type <- as.integer(type) | ||
| if(!any(type == c(1L,4L))) | ||
| stop(paste("Quantiles of type", type, "are not supported"), call.=FALSE) | ||
| ## extract marks | ||
| X <- coerce.marks.numeric(X) | ||
| m <- marks(X) | ||
| ## multiple columns of marks? | ||
| if(!is.null(dim(m)) && ncol(m) > 1) { | ||
| ## compute separately for each column | ||
| Xlist <- unstack(X) | ||
| Zlist <- lapply(Xlist, SpatialQuantile, prob=prob, sigma=sigma, ..., | ||
| type=type, | ||
| at=at, leaveoneout=leaveoneout, weights=weights, | ||
| edge=edge, diggle=diggle, verbose=verbose) | ||
| ZZ <- switch(at, | ||
| pixels = as.imlist(Zlist), | ||
| points = do.call(data.frame, Zlist)) | ||
| return(ZZ) | ||
| } | ||
| ## single column of marks | ||
| m <- as.numeric(m) | ||
| nX <- npoints(X) | ||
| #' unique mark values | ||
| um <- sort(unique(m)) | ||
| Num <- length(um) | ||
| #' trivial cases | ||
| if(nX == 0 || ((Num == 1) && leaveoneout)) { | ||
| Z <- switch(at, | ||
| pixels = as.im(NA_real_, W=Window(X), ...), | ||
| points = rep(NA_real_, nX)) | ||
| attr(Z, "sigma") <- sigma | ||
| return(Z) | ||
| } | ||
| if(Num == 1) { | ||
| Z <- switch(at, | ||
| pixels = as.im(um[1], W=Window(X), ...), | ||
| points = rep(um[1], nX)) | ||
| attr(Z, "sigma") <- sigma | ||
| return(Z) | ||
| } | ||
| #' numerical weights | ||
| if(!is.null(weights)) { | ||
| check.nvector(weights, nX) | ||
| if(any(weights < 0)) stop("Negative weights are not permitted") | ||
| if(sum(weights) < .Machine$double.eps) | ||
| stop("Weights are numerically zero; quantiles are undefined", call.=FALSE) | ||
| } | ||
| #' start main calculation | ||
| ## bandwidth selector | ||
| if(is.function(sigma)) | ||
| sigma <- sigma(X, ...) | ||
| #' edge correction has no effect if diggle=FALSE | ||
| #' (because uniform edge correction cancels) | ||
| edge <- edge && diggle | ||
| #' smoothed intensity of entire pattern | ||
| UX <- unmark(X) | ||
| LX <- density(UX, ..., | ||
| at=at, leaveoneout=leaveoneout, | ||
| weights=weights, | ||
| edge=edge, diggle=diggle, positive=TRUE) | ||
| #' extract smoothing bandwidth actually used | ||
| sigma <- attr(LX, "sigma") | ||
| varcov <- attr(LX, "varcov") | ||
| #' initialise result | ||
| Z <- LX | ||
| Z[] <- NA | ||
| #' guard against underflow | ||
| tinythresh <- 8 * .Machine$double.eps | ||
| if(underflow <- (min(LX) < tinythresh)) { | ||
| Bad <- (LX < tinythresh) | ||
| warning(paste("Numerical underflow detected at", | ||
| percentage(Bad, 1), "of", paste0(atName, ";"), | ||
| "sigma is probably too small"), | ||
| call.=FALSE) | ||
| #' apply l'Hopital's Rule at the problem locations | ||
| Z[Bad] <- nnmark(X, at=at, xy=LX)[Bad] | ||
| Good <- !Bad | ||
| } | ||
| #' compute | ||
| for(k in 1:Num) { | ||
| #' cumulative spatial weight of points with marks <= m_[k] | ||
| if(k == Num) { | ||
| Acum.k <- 1 | ||
| } else { | ||
| w.k <- (m <= um[k]) * (weights %orifnull% 1) | ||
| Lcum.k <- density(UX, | ||
| weights=w.k, | ||
| sigma=sigma, varcov=varcov, | ||
| xy=LX, | ||
| at=at, leaveoneout=leaveoneout, | ||
| edge=edge, diggle=diggle, positive=TRUE) | ||
| Acum.k <- Lcum.k/LX | ||
| } | ||
| if(k == 1) { | ||
| #' region where quantile is um[1] | ||
| relevant <- (Acum.k >= prob) | ||
| if(underflow) relevant <- relevant & Good | ||
| if(any(relevant)) { | ||
| Z[relevant] <- um[1] | ||
| if(verbose) | ||
| splat("value um[1] =", um[1], "assigned to", sum(relevant), atName) | ||
| } | ||
| } else { | ||
| #' region where quantile is between um[k-1] and um[k] | ||
| unassigned <- (Acum.kprev < prob) | ||
| if(underflow) unassigned <- unassigned & Good | ||
| if(!any(unassigned)) break | ||
| relevant <- unassigned & (Acum.k >= prob) | ||
| if(any(relevant)) { | ||
| if(type == 1) { | ||
| ## left-continuous inverse | ||
| left <- (Acum.k > prob) | ||
| Z[relevant & left] <- um[k-1] | ||
| Z[relevant & !left] <- um[k] | ||
| } else if(type == 4) { | ||
| ## linear interpolation | ||
| Z[relevant] <- um[k-1] + | ||
| (um[k] - um[k-1]) * ((prob - Acum.kprev)/(Acum.k - Acum.kprev))[relevant] | ||
| } | ||
| if(verbose) | ||
| splat("values between", | ||
| paste0("um", paren(k-1, "[")), "=", um[k-1], | ||
| "and", | ||
| paste0("um", paren(k, "[")), "=", um[k], | ||
| "assigned to", sum(relevant), atName) | ||
| } | ||
| } | ||
| Acum.kprev <- Acum.k | ||
| } | ||
| attr(Z, "sigma") <- sigma | ||
| attr(Z, "varcov") <- varcov | ||
| return(Z) | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,130 @@ | ||
| \name{SpatialMedian.ppp} | ||
| \alias{SpatialMedian.ppp} | ||
| \title{ | ||
| Spatially Weighted Median of Values at Points | ||
| } | ||
| \description{ | ||
| Given a spatial point pattern with numeric marks, | ||
| compute a weighted median of the mark values, | ||
| with spatially-varying weights that depend on distance to the data points. | ||
| } | ||
| \usage{ | ||
| \method{SpatialMedian}{ppp}(X, sigma = NULL, \dots, | ||
| type = 1, at = c("pixels", "points"), leaveoneout = TRUE, | ||
| weights = NULL, edge = TRUE, diggle = FALSE, verbose = FALSE) | ||
| } | ||
| \arguments{ | ||
| \item{X}{ | ||
| A spatial point pattern (object of class \code{"ppp"}) | ||
| with numeric marks. | ||
| } | ||
| \item{sigma}{ | ||
| Smoothing bandwidth, passed to \code{\link{density.ppp}}. | ||
| } | ||
| \item{\dots}{ | ||
| Further arguments passed to \code{\link{density.ppp}} controlling the | ||
| spatial smoothing. | ||
| } | ||
| \item{type}{ | ||
| Type of median (see \code{\link[stats]{quantile.default}}. | ||
| Only types 1 and 4 are currently implemented. | ||
| } | ||
| \item{at}{ | ||
| Character string indicating whether to compute the median | ||
| at every pixel of a pixel image (\code{at="pixels"}, the default) | ||
| or at every data point of \code{X} (\code{at="points"}). | ||
| } | ||
| \item{leaveoneout}{ | ||
| Logical value indicating whether to compute a leave-one-out | ||
| estimator. Applicable only when \code{at="points"}. | ||
| } | ||
| \item{weights}{ | ||
| Optional vector of numeric weights attached to the points of \code{X}. | ||
| } | ||
| \item{edge,diggle}{ | ||
| Arguments passed to \code{\link{density.ppp}} to | ||
| determine the edge correction. | ||
| } | ||
| \item{verbose}{ | ||
| Logical value specifying whether to print progress reports | ||
| during the calculation. | ||
| } | ||
| } | ||
| \details{ | ||
| The argument \code{X} should be a spatial point pattern | ||
| (object of class \code{"ppp"}) with numeric marks. | ||
|
|
||
| The algorithm computes the weighted median of the mark values | ||
| at each desired spatial location, using spatially-varying weights | ||
| which depend on distance to the data points. | ||
|
|
||
| Suppose the data points are at spatial locations | ||
| \eqn{x_1,\ldots,x_n}{x[1], ..., x[n]} | ||
| and have mark values | ||
| \eqn{y_1,\ldots,y_n}{y[1], ..., y[n]}. | ||
| For a query location \eqn{u}, the smoothed median is defined | ||
| as the weighted median of the mark values | ||
| \eqn{y_1,\ldots,y_n}{y[1], ..., y[n]} with weights | ||
| \eqn{w_1,\ldots,w_n}{w[1], ..., w[n]}, | ||
| where | ||
| \deqn{ | ||
| w_i = \frac{k(u,x_i)}{\sum_{j=1}^n k(u,x_j)} | ||
| }{ | ||
| w[i] = k(u,x[i])/(k(u, x[1]) + ... + k(u, x[n])) | ||
| } | ||
| where \eqn{k(u,v)} is the smoothing kernel with bandwidth \code{sigma} | ||
|
|
||
| If \code{at="points"} and \code{leaveoneout=TRUE}, then | ||
| a leave-one-out calculation is performed, which means that | ||
| when the query location is a data point \eqn{x_i}{x[i]}, | ||
| the value at the data point is ignored, and | ||
| the weighted median is computed from the values \eqn{y_j}{y[j]} | ||
| for all \eqn{j} not equal to \eqn{i}. | ||
| } | ||
| \value{ | ||
| \emph{If \code{X} has a single column of marks:} | ||
| \itemize{ | ||
| \item | ||
| If \code{at="pixels"} (the default), the result is | ||
| a pixel image (object of class \code{"im"}). | ||
| \item | ||
| If \code{at="points"}, the result is a numeric vector | ||
| of length equal to the number of points in \code{X}. | ||
| } | ||
| \emph{If \code{X} has a data frame of marks:} | ||
| \itemize{ | ||
| \item | ||
| If \code{at="pixels"} (the default), the result is a named list of | ||
| pixel images (object of class \code{"im"}). There is one | ||
| image for each column of marks. This list also belongs to | ||
| the class \code{"solist"}, for which there is a plot method. | ||
| \item | ||
| If \code{at="points"}, the result is a data frame | ||
| with one row for each point of \code{X}, | ||
| and one column for each column of marks. | ||
| Entries are values of the interpolated function at the points of \code{X}. | ||
| } | ||
| The return value has attributes | ||
| \code{"sigma"} and \code{"varcov"} which report the smoothing | ||
| bandwidth that was used. | ||
| } | ||
| \author{ | ||
| \adrian. | ||
| } | ||
| \seealso{ | ||
| \code{\link{SpatialQuantile.ppp}} | ||
|
|
||
| \code{\link{Smooth.ppp}} | ||
| } | ||
| \examples{ | ||
| X <- longleaf | ||
| if(!interactive()) { | ||
| ## mark values rounded to nearest multiple of 10 to reduce check time | ||
| marks(X) <- round(marks(X), -1) | ||
| } | ||
| Z <- SpatialMedian(X, sigma=30) | ||
| ZX <- SpatialMedian(X, sigma=30, at="points") | ||
| } | ||
| \keyword{spatial} | ||
| \keyword{methods} | ||
| \keyword{smooth} |
Oops, something went wrong.