/
np.den.R
executable file
·228 lines (207 loc) · 11.5 KB
/
np.den.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
#····································································
# np.den.R (npsp package)
#····································································
# bin.den S3 class and methods
# bin.den(x, nbin)
# as.bin.den()
# as.bin.den.bin.data(object, ...)
# np.den S3 class and methods
# np.den() S3 generic
# np.den.default(x, nbin, h, degree, drv, ncv, ...)
# np.den.bin.den(x, h, degree, drv, ncv, ...)
# np.den.bin.data(x, h, degree, drv, ncv, ...)
# np.den.svar.bin(x, h, degree, drv, ncv, ...)
#
# (c) Ruben Fernandez-Casal
# Created: Oct 2013
#
# NOTE: Press Ctrl + Shift + O to show document outline in RStudio
#····································································
#····································································
# bin.den(x, nbin = NULL) {
#····································································
#' Linear binning for density estimation
#'
#' Creates a \code{bin.den}-\code{\link{class}} (gridded binned density) object
#' with linear binning counts.
#' @aliases bin.den-class
#' @param x vector or matrix of covariates (e.g. spatial coordinates).
#' Columns correspond with dimensions and rows with observations.
#' @param nbin vector with the number of bins on each dimension.
#' @details If parameter \code{nbin} is not specified is set to \code{pmax(25, rule.binning(x))}.
#' @return Returns an S3 object of \code{\link{class}} \code{bin.den} (extends \code{\link{data.grid}}).
#' A list with the following 3 components:
#' \item{binw}{vector or array (dimension \code{nbin}) with the bin counts (weights).}
#' \item{grid}{a \code{\link{grid.par}}-\code{\link{class}} object with the grid parameters.}
#' \item{data}{a list with a component \code{$x} with argument \code{x}.}
#' @seealso \code{\link{np.den}}, \code{\link{h.cv}}, \code{\link{bin.data}},
#' \code{\link{locpol}}, \code{\link{rule.binning}}.
#' @examples
#' binden <- bin.den(earthquakes[, c("lon", "lat")], nbin = c(30,30))
#' bindat <- binning(earthquakes[, c("lon", "lat")], earthquakes$mag, nbin = c(30,30))
#' all.equal(binden, as.bin.den(bindat))
#' @export
# Interface to the Fortran routine "bin_den"
bin.den <- function(x, nbin = NULL) {
#····································································
x <- as.matrix(x)
ny <- nrow(x) # number of data
# Remove missing values
ok <- complete.cases(x) # observations having no missing values across x
if (any(!ok)) {
warning("missing values removed")
x <- x[ok,]
ny <- nrow(x)
}
nd <- ncol(x) # number of dimensions
if(is.null(nbin)) nbin <- pmax(25L, rule.binning(x)) else # dimensions of the binning grid
if (nd != length(nbin)) stop("arguments 'x' and 'nbin' have incompatible dimensions.")
nt <- prod(nbin)
# Let's go FORTRAN!
# subroutine bin_den(nd, nbin, x, ny, bin_min, bin_max, bin_w)
ret <-.Fortran( "bin_den", nd = as.integer(nd), nbin = as.integer(nbin),
xt = as.double(t(x)), ny = as.integer(ny),
min = double(nd), max = double(nd),
binw = double(nt), PACKAGE = "npsp")
# Construir o resultado
result <- with( ret, data.grid( binw = binw,
grid = grid.par(n = nbin, min = min, max = max,
dimnames = dimnames(x)[[2]])) )
result$data <- list(x = x)
oldClass(result) <- c("bin.den", "data.grid")
return(result)
#····································································
} # bin.den
#····································································
#' @rdname bin.den
#' @param object (gridded data) used to select a method.
#' @param ... further arguments passed to or from other methods.
#' @export
#····································································
as.bin.den <- function(object, ...) {
UseMethod("as.bin.den")
} # S3 generic function as.bin.den
#····································································
#' @rdname bin.den
#' @method as.bin.den data.grid
#' @param weights.ind integer or character with the index or name of the component
#' containing the bin counts/weights.
#' @export
as.bin.den.data.grid <- function(object, weights.ind = 1, ...) {
#····································································
if (!inherits(object, "data.grid"))
stop("function only works for objects of class (or extending) 'data.grid'")
x <- coords(object$grid)
binw <- object[[weights.ind]]
index <- is.finite(binw)
binw[!index] <- 0
if (!all(index <- binw > 0)) x <- x[index, ]
result <- list(binw = binw, grid = object$grid, data = list(x = x))
oldClass(result) <- c("bin.den", "data.grid")
return(result)
}
#····································································
#' @rdname bin.den
#' @method as.bin.den bin.den
#' @export
as.bin.den.bin.den <- function(object, ...) {
#····································································
if (inherits(object, "svar.bin")) {
warning("Conversion not yet implemented; using 'as.bin.den.data.grid()'...")
return(as.bin.den.data.grid(object, weights.ind = 'binw'))
}
if (!is.null(object$mask)) object$binw[!object$mask] <- 0
result <- with( object, list(binw = binw, grid = grid, data = list(x = data$x)))
oldClass(result) <- c("bin.den", "data.grid")
return(result)
}
#····································································
# np.den(x, ...)
#····································································
#' Local polynomial density estimation
#'
#' Estimates a multidimensional probability density function (and its first derivatives)
#' using local polynomial kernel smoothing of linearly binned data.
#' @aliases np.den-class
#' @param x a (data) object used to select a method.
#' @param ... further arguments passed to or from other methods.
#' @return Returns an S3 object of class \code{np.den} (locpol den + bin den + grid par.).
#' A \code{\link{bin.den}} object with the additional (some optional) 3 components:
#' \item{est}{vector or array (dimension \code{nbin}) with the local polynomial density estimates. }
#' \item{locpol}{a list with 6 components:
#' \itemize{
#' \item{\code{degree} degree of the polinomial.}
#' \item{\code{h} bandwidth matrix.}
#' \item{\code{rm} residual mean (of the escaled bin counts).}
#' \item{\code{rss} sum of squared residuals (of the escaled bin counts).}
#' \item{\code{ncv} number of cells ignored (in each dimension).}
#' }}
#' \item{deriv}{(if requested) matrix of first derivatives.}
#' @seealso \code{\link{bin.den}}, \code{\link{binning}}, \code{\link{h.cv}},
#' \code{\link{data.grid}}.
#' @examples
#' bin.den <- binning(earthquakes[, c("lon", "lat")], nbin = c(30,30))
#' h.den <- h.cv(bin.den)
#' den <- np.den(bin.den, h = h.den$h)
#' plot(den, main = 'Estimated log(density)')
#' @export
#····································································
np.den <- function(x, ...) {
UseMethod("np.den")
}
#····································································
# np.den.default(x, nbin = NULL, h = NULL, degree = 1 + as.numeric(drv), drv = FALSE, ncv = 0, ...) {
#····································································
#' @rdname np.den
#' @method np.den default
#' @inheritParams locpol.default
#' @details Standard generic function with a default method (interface to the
#' fortran routine \code{lp_data_grid}), in which argument \code{x}
#' is a vector or matrix of covariates (e.g. spatial coordinates).
#' In this case, the data are binned (calls \code{\link{bin.den}}) and the local fitting
#' procedure is applied to the scaled bin counts (calls \code{\link{np.den.bin.den}}).
#'
#' If parameter \code{nbim} is not specified is set to \code{rep(25, ncol(x))}.
#'
#' A multiplicative triweight kernel is used to compute the weights.
#'
#' If \code{ncv > 1}, estimates are computed by leaving out cells with indexes within
#' the intervals \eqn{[x_i - ncv + 1, x_i + ncv - 1]}, at each dimension i, where \eqn{x}
#' denotes the index of the estimation position.
#' @references
#' Wand, M.P. and Jones, M.C. (1995) \emph{Kernel Smoothing}. Chapman and Hall, London.
#' @export
np.den.default <- function(x, nbin = NULL, h = NULL, degree = 1 + as.numeric(drv),
drv = FALSE, ncv = 0, ...) {
xbin <- bin.den(x, nbin)
return(locpol.bin.den(xbin, h = h, degree = degree, drv = drv, ncv = ncv, ...))
}
#····································································
# np.den.bin.den(x, h, degree, drv, ncv, ...) ----
#····································································
#' @rdname np.den
#' @method np.den bin.den
#' @export
np.den.bin.den <- locpol.bin.den
#····································································
#' @rdname np.den
#' @method np.den bin.data
#' @export
np.den.bin.data <- function(x, h = NULL, degree = 1 + as.numeric(drv), drv = FALSE,
ncv = 0, ...) {
#····································································
result <- np.den.bin.den(x, h = h, degree = degree, drv = drv, ncv = ncv, ...)
oldClass(result) <- c("np.den", "bin.data", "bin.den", "data.grid")
return(result)
}
#····································································
#' @rdname np.den
#' @method np.den svar.bin
#' @export
np.den.svar.bin <- function(x, h = NULL, degree = 1 + as.numeric(drv), drv = FALSE,
ncv = 0, ...) {
#····································································
result <- np.den.bin.den(x, h = h, degree = degree, drv = drv, ncv = ncv, ...)
oldClass(result) <- c("np.den", "svar.bin", "bin.data", "bin.den", "data.grid")
return(result)
}