-
Notifications
You must be signed in to change notification settings - Fork 0
/
test-logLik-mev.R
128 lines (115 loc) · 4.84 KB
/
test-logLik-mev.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
# Check that logLik(object) and logLik(logLikVec(object)) agree
if (requireNamespace("mev", quietly = TRUE)) {
library(mev)
# mev::fit.gev
# The example from the mev::gev.fit documentation
gev_fit <- mev::fit.gev(revdbayes::portpirie, show = FALSE)
temp <- gev_fit
adj_gev_fit <- alogLik(gev_fit)
class(temp) <- "laxmev_gev"
test_that("mev::fit.gev: logLik() vs. logLik(logLikVec)", {
testthat::expect_equal(logLik(gev_fit), logLik(logLikVec(temp)))
})
# Check that alogLik also returned the correct maximised log-likelihood
test_that("mev::fit.gev: logLik() vs. logLik(logLikVec)", {
testthat::expect_equal(logLik(gev_fit), logLik(adj_gev_fit))
})
# Check logLik.gev.fit: trivially correct
test_that("mev::fit.gev: logLik() vs. logLik(logLikVec)", {
testthat::expect_equal(logLik(temp), logLik(logLikVec(temp)))
})
# mev::fit.gpd
# An example from the mev::fit.gpd documentation
data(eskrain)
gpd_mev <- fit.gpd(eskrain, threshold = 35, method = 'Grimshaw')
temp <- gpd_mev
adj_gpd_mev <- alogLik(gpd_mev)
class(temp) <- "laxmev_gpd"
test_that("mev::fit.gpd: logLik() vs. logLik(logLikVec)", {
testthat::expect_equal(logLik(gpd_mev), logLik(logLikVec(temp)))
})
# Check that alogLik also returned the correct maximised log-likelihood
test_that("mev::fit.gpd: logLik() vs. logLik(logLikVec)", {
testthat::expect_equal(logLik(gpd_mev), logLik(adj_gpd_mev))
})
# Check logLik.evd_fgev: trivially correct
test_that("mev::fit.gpd: logLik() vs. logLik(logLikVec)", {
testthat::expect_equal(logLik(temp), logLik(logLikVec(temp)))
})
# mev::fit.pp
# Use simulated data
set.seed(1112019)
x <- revdbayes::rgp(365 * 10, loc = 0, scale = 1, shape = 0.1)
pfit <- mev::fit.pp(x, threshold = 1, npp = 365)
temp <- pfit
adj_pfit <- alogLik(pfit)
class(temp) <- "laxmev_pp"
test_that("mev::fit.pp: logLik() vs. logLik(logLikVec)", {
testthat::expect_equal(logLik(pfit), logLik(logLikVec(temp)))
})
# Check that alogLik also returned the correct maximised log-likelihood
test_that("mev::fit.pp: logLik() vs. logLik(logLikVec)", {
testthat::expect_equal(logLik(pfit), logLik(adj_pfit))
})
# Check logLik.gev.fit: trivially correct
test_that("mev::fit.pp: logLik() vs. logLik(logLikVec)", {
testthat::expect_equal(logLik(temp), logLik(logLikVec(temp)))
})
# Check that if only threshold exceedances are supplied to mev::fit.pp()
# then an error is thrown
data(eskrain)
pp_mle <- fit.pp(eskrain, threshold = 30, np = 6201)
adj_pp_mle <- try(alogLik(pp_mle), silent = TRUE)
test_that("mev::fit.pp: only exceedances throws an error", {
testthat::expect_identical(class(adj_pp_mle), "try-error")
})
# ismev::fit.rlarg
# An example from the mev::fit.rlarg documentation
set.seed(31102019)
xdat <- rrlarg(n = 10, loc = 0, scale = 1, shape = 0.1, r = 4)
rfit <- fit.rlarg(xdat)
temp <- rfit
adj_rfit <- alogLik(rfit)
class(temp) <- "laxmev_rlarg"
# logLik.mev_rlarg() produces an object with an extra attribute names
# equal to "scale". logLik.mev_rlarg() and logLik.laxmev_rlarg() also
# disagree on nobs. Therefore, use equivalent in the first two tests
test_that("mev::fit.rlarg: logLik() vs. logLik(logLikVec)", {
testthat::expect_equivalent(logLik(rfit), logLik(logLikVec(temp)))
})
# Check that alogLik also returned the correct maximised log-likelihood
test_that("mev::fit.rlarg: logLik() vs. logLik(logLikVec)", {
testthat::expect_equivalent(logLik(rfit), logLik(adj_rfit))
})
# Check logLik.rlarg.fit: trivially correct
test_that("mev::fit.rlarg: logLik() vs. logLik(logLikVec)", {
testthat::expect_equal(logLik(temp), logLik(logLikVec(temp)))
})
# mev::fit.egp
# An example from the mev::fit.egp documentation
if (requireNamespace("evd", quietly = TRUE)) {
models <- c("egp1", "egp2", "egp3")
set.seed(7102019)
xdat <- revdbayes::rgp(n = 100, loc = 0, scale = 1, shape = 0.5)
for (i in 1:3) {
fitted <- fit.egp(xdat = xdat, thresh = 1, model = models[i],
show = FALSE)
temp <- fitted
adj_fitted <- alogLik(fitted)
class(temp) <- "laxmev_egp"
# logLik.mev_egp() produces an object with an extra attribute names
# equal to "scale". Therefore, use equivalent in the first two tests
test_that("mev::fit.egp: logLik() vs. logLik(logLikVec)", {
testthat::expect_equivalent(logLik(fitted), logLik(logLikVec(temp)))
})
# Check that alogLik also returned the correct maximised log-likelihood
test_that("mev::fit.egp: logLik() vs. logLik(logLikVec)", {
testthat::expect_equivalent(logLik(fitted), logLik(adj_fitted))
})
# Check logLik.evd_fgev: trivially correct (up to naming)
test_that("mev::fit.egp: logLik() vs. logLik(logLikVec)", {
testthat::expect_equal(logLik(temp), logLik(logLikVec(temp)))
})
}
}
}