-
Notifications
You must be signed in to change notification settings - Fork 5
/
mde.R
143 lines (129 loc) · 4.6 KB
/
mde.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
### actuar: Actuarial Functions and Heavy Tailed Distributions
###
### Minimum distance estimation.
###
### AUTHORS: Mathieu Pigeon, Vincent Goulet <[email protected]>
mde <- function(x, fun, start, measure = c("CvM", "chi-square", "LAS"),
weights = NULL, ...)
{
## General form of the function to minimize.
myfn <- function(parm, x, weights, ...)
{
y <- G(parm, x, ...) - Gn(x)
drop(crossprod(weights * y, y))
}
## Extract call; used to build the call to optim().
Call <- match.call(expand.dots = TRUE)
## Argument checking
if (missing(start) || !is.list(start))
stop(sprintf("%s must be a named list", sQuote("start")))
if (missing(fun) || !(is.function(fun)))
stop(sprintf("%s must be supplied as a function", sQuote("fun")))
grouped <- inherits(x, "grouped.data")
if (!(is.numeric(x) || grouped))
stop(sprintf("%s must be a numeric vector or an object of class %s",
sQuote("x"), dQuote("grouped.data")))
## Make sure that any argument of 'fun' specified in '...' is held
## fixed.
dots <- names(list(...))
dots <- dots[!is.element(dots, c("upper", "lower"))]
start <- start[!is.element(names(start), dots)]
## Adapt 'fun' to our needs; taken from MASS::fitdistr.
nm <- names(start)
f <- formals(fun)
args <- names(f)
m <- match(nm, args)
if (any(is.na(m)))
stop(sprintf("%s specifies names which are not arguments to %s",
sQuote("start"), sQuote("fun")))
formals(fun) <- c(f[c(1, m)], f[-c(1, m)]) # reorder arguments
fn <- function(parm, x, ...) fun(x, parm, ...)
if ((l <- length(nm)) > 1)
body(fn) <- parse(text = paste("fun(x,", paste("parm[", 1:l, "]", collapse = ", "), ")"))
measure <- match.arg(measure)
## Cramer-von Mises. Use the true and empirical cdf for individual
## data, or the true cdf and the ogive for grouped data.
if (measure == "CvM")
{
G <- fn
Gn <- if (grouped) ogive(x) else ecdf(x)
if (is.null(weights)) weights <- 1
Call$x <- knots(Gn)
Call$par <- start
}
## Modified Chi-square.
if (measure == "chi-square")
{
if (!grouped)
stop(sprintf("%s measure requires an object of class %s",
dQuote("chi-square"), dQuote("grouped.data")))
if (any((nj <- x[, 2]) == 0))
stop("frequency must be larger than 0 in all groups")
og <- ogive(x)
x <- knots(og)
n <- sum(nj)
G <- function(...) n * diff(fn(...))
Gn <- function(...) n * diff(og(...))
if (is.null(weights)) weights <- 1/nj
Call$x <- x
Call$par <- start
}
## Layer average severity.
if (measure == "LAS")
{
if (!grouped)
stop(sprintf("%s measure requires an object of class %s",
dQuote("LAS"), dQuote("grouped.data")))
e <- elev(x)
x <- knots(e)
G <- function(...) diff(fn(...))
Gn <- function(...) diff(e(...))
if (is.null(weights)) weights <- 1
Call$x <- x
Call$par <- start
}
## optim() call
Call[[1]] <- as.name("optim")
Call$fun <- Call$start <- Call$measure <- NULL
Call$fn <- myfn
Call$weights <- weights
Call$hessian <- FALSE
if (is.null(Call$method))
{
if (any(c("lower", "upper") %in% names(Call)))
Call$method <- "L-BFGS-B"
else if (length(start) > 1)
Call$method <- "BFGS"
else
Call$method <- "Nelder-Mead"
}
res <- eval(Call)
## Return result
if (res$convergence > 0)
stop("optimization failed")
structure(list(estimate = res$par, distance = res$value),
class = c("mde","list"))
}
print.mde <- function(x, digits = getOption("digits"), ...)
{
ans1 <- format(x$estimate, digits = digits)
ans1 <- sapply(ans1, function(x) paste("", x))
nm1 <- names(ans1)
nm1 <- paste(substring(" ", 1L, (nchar(ans1) - nchar(nm1)) %/% 2),
nm1)
nm1 <- paste(nm1,
substring(" ", 1L, (nchar(ans1) - nchar(nm1)) %/% 2 + 1))
names(ans1) <- nm1
ans2 <- format(x$distance, digits = digits)
ans2 <- sapply(ans2, function(x) paste("", x))
nm2 <- "distance"
nm2 <- paste(substring(" ", 1L, (nchar(ans2) - nchar(nm2)) %/% 2),
nm2)
nm2 <- paste(nm2,
substring(" ", 1L, (nchar(ans2) - nchar(nm2)) %/% 2))
names(ans2) <- nm2
print(ans1, quote = FALSE)
cat("\n")
print(ans2, quote = FALSE)
invisible(x)
}