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
|
#' @title
#' The Venter / Dalenius / LMS mode estimator
#'
#' @description
#' This function computes the Venter mode estimator, also called the Dalenius,
#' or LMS (Least Median Square) mode estimator.
#'
#' @details
#' The modal interval, i.e. the shortest interval among intervals containing
#' \code{k+1} observations, is first computed. (In dimension > 1, this question
#' is known as a 'k-enclosing problem'.)
#' 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, so
#' \code{bw} can be seen as the fraction of the observations to be considered
#' for the shortest interval.
#'
#' If \code{type = 1}, the midpoint of the modal interval is returned.
#' If \code{type = 2}, the \code{floor((k+1)/2)}th element of the modal
#' interval is returned.
#' If \code{type = 3} or \code{type = "dalenius"}, the median of the modal
#' interval is returned.
#' If \code{type = 4} or \code{type = "shorth"}, the mean of the modal interval
#' is returned.
#' If \code{type = 5} or \code{type = "ekblom"}, Ekblom's
#' \eqn{L_{-\infty}}{L_{-infinity}} estimate is returned, see Ekblom (1972).
#' If \code{type = 6} or \code{type = "hsm"}, the half sample mode (hsm) is
#' computed, see \code{\link{hsm}}.
#'
#' @note
#' The user may call \code{venter} through
#' \code{mlv(x, method = "venter", ...)}.
#'
#' @references
#' \itemize{
#' \item Dalenius T. (1965).
#' The Mode - A Negleted Statistical Parameter.
#' \emph{J. Royal Statist. Soc. A}, \emph{128}:110-117.
#'
#' \item Venter J.H. (1967).
#' On estimation of the mode.
#' \emph{Ann. Math. Statist.}, \bold{38}(5):1446-1455.
#'
#' \item Ekblom H. (1972).
#' A Monte Carlo investigation of mode estimators in small samples.
#' \emph{Applied Statistics}, \bold{21}:177-184.
#
#\item Rousseeuw and Leroy, 1987 #(ou bien Andrews ?)
#'
#' \item Leclerc J. (1997).
#' Comportement limite fort de deux estimateurs du mode : le shorth et l'estimateur naif.
#' \emph{C. R. Acad. Sci. Paris, Serie I}, \bold{325}(11):1207-1210.
#' }
#'
#' @param x
#' numeric. Vector of observations.
#'
#' @param bw
#' numeric. The bandwidth to be used. Should belong to (0, 1]. See 'Details'.
#'
#' @param k
#' numeric. See 'Details'.
#'
#' @param iter
#' numeric. Number of iterations.
#'
#' @param type
#' numeric or character. The type of Venter estimate to be computed. See 'Details'.
#'
#' @param tie.action
#' character. The action to take if a tie is encountered.
#'
#' @param tie.limit
#' numeric. A limit deciding whether or not a warning is given when a tie is
#' encountered.
#'
#' @param warn
#' logical. If \code{TRUE}, a warning is thrown when a tie is encountered.
#'
#' @param ...
#' Further arguments.
#'
#' @return
#' A numeric value is returned, the mode estimate.
#'
#' @seealso
#' \code{\link[modeest]{mlv}} for general mode estimation,
#' \code{\link[modeest]{hsm}} for the half sample mode.
#'
#' @importFrom stats median
#' @export
#' @aliases Venter
#'
#' @examples
#' library(evd)
#'
#' # Unimodal distribution
#' x <- rgev(1000, loc = 23, scale = 1.5, shape = 0)
#'
#' ## True mode
#' gevMode(loc = 23, scale = 1.5, shape = 0)
#'
#' ## Estimate of the mode
#' venter(x, bw = 1/3)
#' mlv(x, method = "venter", bw = 1/3)
#'
venter <-
function(x,
bw = NULL,
k,
iter = 1,
type = 1,
tie.action = "mean",
tie.limit = 0.05,
warn = FALSE)
{
ny <- length(x)
## Initialization
type <- match.arg(tolower(as.character(type)),
c("-inf", "1", "2", "3", "dalenius", "4", "shorth", "5", "ekblom", "6", "hsm"))
if (type == "3") type <- "dalenius"
if (type == "4") type <- "shorth"
if (type == "-Inf" || type == "5") type <- "ekblom"
if (type == "6") type <- "hsm"
if (type == "hsm") {
return(hsm(x = x, bw = bw, k = k, tie.action = tie.action,
tie.limit = tie.limit))
}
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)) {
if (type == "ekblom") {
k <- 1
} else {
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]
diffs <- sup - inf
i <- which(diffs == min(diffs))
## Ties?
if (length(i) > 1) i <- .deal.ties(ny, i, tie.action, tie.limit, warn = warn)
## Output
M <- switch(type,
"1" = (y[i] + y[i+k])/2,
"2" = y[i+floor((k+1)/2)],
"dalenius" = stats::median(y[i:(i+k)]),
"shorth" = mean(y[i:(i+k)]),
"ekblom" = ifelse(y[i+2]-y[i+1]>y[i]-y[i-1], y[i], y[i+1]))
if (iter > 1) {
M <- Recall(x = y[i:(i+k)], bw = (k + 1)/ny, iter = iter-1, type = type,
tie.action = tie.action, tie.limit = tie.limit)
}
#attr(M, "inf") <- y[i]
#attr(M, "sup") <- y[i+k]
M
}
#' @export
#' @rdname venter
#'
shorth <- function(x, ...) venter(x, type = "shorth", ...)
|
