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
|
#' @title
#' The Grenander mode estimator
#'
#' @description
#' This function computes the Grenander mode estimator.
#'
#' @details
#' The Grenander estimate is defined by
#' \deqn{ \frac{ \sum_{j=1}^{n-k} \frac{(x_{j+k} + x_{j})}{2(x_{j+k} - x_{j})^p} }
#' { \sum_{j=1}^{n-k} \frac{1}{(x_{j+k} - x_{j})^p} } }{ ( sum_{j=1}^{n-k} (x_{j+k} + x_{j})/(2(x_{j+k} - x_{j})^p) ) / ( sum_{j=1}^{n-k} 1/((x_{j+k} - x_{j})^p) ) }
#'
#' If \eqn{p}{p} tends to infinity, this estimate tends to the Venter mode estimate;
#' this justifies to call \code{\link[modeest]{venter}} if \code{p = Inf}.
#'
#' The user should either give the bandwidth \code{bw} or the argument \code{k},
#' \code{k} being taken equal to \code{ceiling(bw*n) - 1} if missing.
#'
#' @references
#' \itemize{
#' \item Grenander U. (1965).
#' Some direct estimates of the mode.
#' \emph{Ann. Math. Statist.}, \bold{36}:131-138.
#'
#' \item Dalenius T. (1965).
#' The Mode - A Negleted Statistical Parameter.
#' \emph{J. Royal Statist. Soc. A}, \emph{128}:110-117.
#'
#' \item Adriano K.N., Gentle J.E. and Sposito V.A. (1977).
#' On the asymptotic bias of Grenander's mode estimator.
#' \emph{Commun. Statist.-Theor. Meth. A}, \bold{6}:773-776.
#'
#' \item Hall P. (1982).
#' Asymptotic Theory of Grenander's Mode Estimator.
#' \emph{Z. Wahrsch. Verw. Gebiete}, \bold{60}:315-334.
#' }
#'
#' @note
#' The user may call \code{grenander} through
#' \code{mlv(x, method = "grenander", bw, k, p, ...)}.
#'
#' @param x
#' numeric. Vector of observations.
#'
#' @param bw
#' numeric. The bandwidth to be used. Should belong to (0, 1].
#'
#' @param k
#' numeric. Paramater 'k' in Grenander's mode estimate, see below.
#'
#' @param p
#' numeric. Paramater 'p' in Grenander's mode estimate, see below.
#' If \code{p = Inf}, the function \code{\link[modeest]{venter}} is used.
#'
#' @param ...
#' Additional arguments to be passed to \code{\link[modeest]{venter}}.
#'
#' @return
#' A numeric value is returned, the mode estimate.
#' If \code{p = Inf}, the \code{\link[modeest]{venter}} mode estimator is returned.
#'
#' @author D.R. Bickel for the original code,
#' P. Poncet for the slight modifications introduced.
#'
#' @seealso
#' \code{\link[modeest]{mlv}} for general mode estimation;
#' \code{\link[modeest]{venter}} for the Venter mode estimate.
#'
#' @export
#' @aliases Grenander
#'
#' @examples
#' # Unimodal distribution
#' x <- rnorm(1000, mean = 23, sd = 0.5)
#'
#' ## True mode
#' normMode(mean = 23, sd = 0.5) # (!)
#'
#' ## Parameter 'k'
#' k <- 5
#'
#' ## Many values of parameter 'p'
#' ps <- seq(0.1, 4, 0.01)
#'
#' ## Estimate of the mode with these parameters
#' M <- sapply(ps, function(p) grenander(x, p = p, k = k))
#'
#' ## Distribution obtained
#' plot(density(M), xlim = c(22.5, 23.5))
#'
grenander <-
function(x,
bw = NULL,
k,
p,
...)
{
if (p == Inf) {
message("argument 'p' is infinite. Venter's mode estimator is used")
return(venter(x = x, bw = bw, k = k, ...))
}
ny <- length(x)
if (missing(k) && !is.null(bw)) {
if (bw <= 0 || bw > 1) {
stop("argument 'bw' must belong to (0, 1]",
call. = FALSE)
}
k <- ceiling(bw*ny) - 1
} else if (missing(k) && is.null(bw)) {
k <- ceiling(ny/2) - 1
}
if (k < 0 || k >= ny) {
stop("argument 'k' must belong to [0, length('x'))",
call. = FALSE)
}
y <- sort(x)
inf <- y[1:(ny-k)]
sup <- y[(k+1):ny]
diff <- sup - inf
tot <- inf + sup
if (any(diff==0)) {
warning("limiting value of Grenander mode used",
call. = FALSE) # TODO: be more specific
M <- mean(ifelse(diff==0, tot, NA), na.rm = TRUE)/2
} else {
b <- sum(tot/diff^p)/2
a <- sum(1/diff^p)
if (is.finite(b/a)) {
M <- b/a
} else {
stop("function 'grenander' failed. Argument 'p' may be too large",
call. = FALSE)
}
}
M
}
|
