diff options
| -rw-r--r-- | R/mkinfit.R | 50 | ||||
| -rw-r--r-- | tests/testthat/test_error_models.R | 50 |
2 files changed, 86 insertions, 14 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R index 224dc7bd..e0a0e525 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -35,6 +35,7 @@ mkinfit <- function(mkinmod, observed, atol = 1e-8, rtol = 1e-10, n.outtimes = 100,
error_model = c("const", "obs", "tc"),
error_model_algorithm = c("direct", "twostep", "threestep", "fourstep", "IRLS"),
+ reweight.tol = 1e-8, reweight.max.iter = 10,
trace_parms = FALSE,
...)
{
@@ -285,7 +286,7 @@ mkinfit <- function(mkinmod, observed, # Trace parameter values if requested and if we are actually optimising
if(trace_parms & update_data) cat(P, "\n")
- if (fixed_degparms[1] != FALSE) {
+ if (is.numeric(fixed_degparms)) {
degparms <- fixed_degparms
errparms <- P # This version of errparms is local to the function
degparms_fixed = TRUE
@@ -293,7 +294,7 @@ mkinfit <- function(mkinmod, observed, degparms_fixed = FALSE
}
- if (fixed_errparms[1] != FALSE) {
+ if (is.numeric(fixed_errparms)) {
degparms <- P
errparms <- fixed_errparms # Local to the function
errparms_fixed = TRUE
@@ -364,7 +365,8 @@ mkinfit <- function(mkinmod, observed, sum(dnorm(x = value.observed, mean = value.predicted, sd = std, log = TRUE)))
}
- # We update the current likelihood and data during the optimisation, not during hessian calculations
+ # We update the current likelihood and data during the optimisation, not
+ # during hessian calculations
if (update_data) {
assign("out_predicted", out_long, inherits = TRUE)
@@ -434,11 +436,10 @@ mkinfit <- function(mkinmod, observed, # Show parameter names if tracing is requested
if(trace_parms) cat(names_optim, "\n")
+ # browser()
+
# Do the fit and take the time until the hessians are calculated
fit_time <- system.time({
- # In the inital run, we assume a constant standard deviation and do
- # not optimise it, as this is unnecessary and leads to instability of the
- # fit (most obvious when fitting the IORE model)
degparms <- c(state.ini.optim, transparms.optim)
if (err_mod == "const" | error_model_algorithm != "direct") {
@@ -463,7 +464,7 @@ mkinfit <- function(mkinmod, observed, fit <- nlminb(errparms, nlogLik, control = control,
lower = lower[names(errparms)],
upper = upper[names(errparms)],
- degparms_fixed = degparms, ...)
+ fixed_degparms = degparms, ...)
errparms <- fit$par
}
if (error_model_algorithm == "fourstep") {
@@ -471,10 +472,11 @@ mkinfit <- function(mkinmod, observed, fit <- nlminb(degparms, nlogLik, control = control,
lower = lower[names(degparms)],
upper = upper[names(degparms)],
- errparms_fixed = errparms, ...)
+ fixed_errparms = errparms, ...)
degparms <- fit$par
}
- if (error_model_algorithm %in% c("direct", "twostep", "threestep", "fourstep")) {
+ if (error_model_algorithm %in% c("direct", "twostep", "threestep", "fourstep") &
+ err_mod != "const") {
if (!quiet) message("Optimising the complete model")
parms.start <- c(degparms, errparms)
fit <- nlminb(parms.start, nlogLik,
@@ -483,6 +485,36 @@ mkinfit <- function(mkinmod, observed, control = control, ...)
fit$logLik <- - nlogLik.current
}
+ if (err_mod != "const" & error_model_algorithm == "IRLS") {
+ reweight.diff <- 1
+ n.iter <- 0
+ errparms_last <- errparms
+
+ while (reweight.diff > reweight.tol &
+ n.iter < reweight.max.iter) {
+
+ if (!quiet) message("Optimising the error model")
+ fit <- nlminb(errparms, nlogLik, control = control,
+ lower = lower[names(errparms)],
+ upper = upper[names(errparms)],
+ fixed_degparms = degparms, ...)
+ errparms <- fit$par
+
+ if (!quiet) message("Optimising the degradation model")
+ fit <- nlminb(degparms, nlogLik, control = control,
+ lower = lower[names(degparms)],
+ upper = upper[names(degparms)],
+ fixed_errparms = errparms, ...)
+ degparms <- fit$par
+
+ reweight.diff <- dist(rbind(errparms, errparms_last))
+ errparms_last <- errparms
+
+ fit$par <- c(fit$par, errparms)
+ nlogLik.current <- nlogLik(c(degparms, errparms), OLS = FALSE)
+ fit$logLik <- - nlogLik.current
+ }
+ }
# We include the error model in the parameter uncertainty analysis, also
# for constant variance, to get a confidence interval for it
diff --git a/tests/testthat/test_error_models.R b/tests/testthat/test_error_models.R index 94703e7f..c656f7d2 100644 --- a/tests/testthat/test_error_models.R +++ b/tests/testthat/test_error_models.R @@ -180,7 +180,6 @@ test_that("Reweighting method 'tc' produces reasonable variance estimates", { test_that("The two-component error model finds the best known AIC values for parent models", { skip_on_cran() - experimental_data_for_UBA_2019 library(parallel) source("~/git/mkin/R/mkinfit.R") source("~/git/mkin/R/mmkin.R") @@ -195,27 +194,68 @@ test_that("The two-component error model finds the best known AIC values for par error_model = "tc", error_model_algorithm = "fourstep") f_9 <- mkinfit("SFO", experimental_data_for_UBA_2019[[9]]$data, error_model = "tc", error_model_algorithm = "IRLS") + f_9 <- mkinfit("SFO", experimental_data_for_UBA_2019[[9]]$data, + error_model = "tc", error_model_algorithm = "d_3") AIC(f_9) - f_tc_exp <- mmkin(c("SFO"), + f_10 <- mkinfit("DFOP", experimental_data_for_UBA_2019[[10]]$data, + error_model = "tc", error_model_algorithm = "IRLS") + f_tc_exp_direct <- mmkin(c("SFO", "DFOP", "HS"), lapply(experimental_data_for_UBA_2019, function(x) x$data), error_model = "tc", error_model_algorithm = "direct", quiet = TRUE) - f_tc_exp <- mmkin(c("SFO"), + f_tc_exp_twostep <- mmkin(c("SFO", "DFOP", "HS"), lapply(experimental_data_for_UBA_2019, function(x) x$data), error_model = "tc", error_model_algorithm = "twostep", quiet = TRUE) - f_tc_exp <- mmkin(c("SFO"), + f_tc_exp_threestep <- mmkin(c("SFO", "DFOP", "HS"), lapply(experimental_data_for_UBA_2019, function(x) x$data), error_model = "tc", error_model_algorithm = "threestep", quiet = TRUE) + f_tc_exp_fourstep <- mmkin(c("SFO", "DFOP", "HS"), + lapply(experimental_data_for_UBA_2019, function(x) x$data), + error_model = "tc", + error_model_algorithm = "fourstep", + quiet = TRUE) + f_tc_exp_IRLS <- mmkin(c("SFO", "DFOP", "HS"), + lapply(experimental_data_for_UBA_2019, function(x) x$data), + error_model = "tc", + error_model_algorithm = "IRLS", + quiet = TRUE) + + AIC_exp_direct <- lapply(f_tc_exp_direct, AIC) + AIC_exp_direct <- lapply(AIC_exp_direct, round, 1) + dim(AIC_exp_direct) <- dim(f_tc_exp_direct) + dimnames(AIC_exp_direct) <- dimnames(f_tc_exp_direct) + + AIC_exp_twostep <- lapply(f_tc_exp_twostep, AIC) + AIC_exp_twostep <- lapply(AIC_exp_twostep, round, 1) + dim(AIC_exp_twostep) <- dim(f_tc_exp_twostep) + dimnames(AIC_exp_twostep) <- dimnames(f_tc_exp_twostep) + + AIC_exp_threestep <- lapply(f_tc_exp_threestep, AIC) + AIC_exp_threestep <- lapply(AIC_exp_threestep, round, 1) + dim(AIC_exp_threestep) <- dim(f_tc_exp_threestep) + dimnames(AIC_exp_threestep) <- dimnames(f_tc_exp_threestep) + + AIC_exp_fourstep <- lapply(f_tc_exp_fourstep, AIC) + AIC_exp_fourstep <- lapply(AIC_exp_fourstep, round, 1) + dim(AIC_exp_fourstep) <- dim(f_tc_exp_fourstep) + dimnames(AIC_exp_fourstep) <- dimnames(f_tc_exp_fourstep) + + AIC_exp_IRLS <- lapply(f_tc_exp_IRLS, AIC) + AIC_exp_IRLS <- lapply(AIC_exp_IRLS, round, 1) + dim(AIC_exp_IRLS) <- dim(f_tc_exp_IRLS) + dimnames(AIC_exp_IRLS) <- dimnames(f_tc_exp_IRLS) AIC_exp <- lapply(f_tc_exp, AIC) dim(AIC_exp) <- dim(f_tc_exp) dimnames(AIC_exp) <- dimnames(f_tc_exp) - expect_equivalent(round(AIC_exp["SFO", c(9, 11, 12)], 1), c(134.9, 125.5, 82.0)) + unlist(AIC_exp["SFO", c(9, 11, 12)]) + expect_equivalent(round(unlist(AIC_exp["SFO", c(9, 11, 12)]), 1), + c(134.9, 125.5, 82.0)) }) |