From feef7146a86b0e1b733b6df1df1420dd11101474 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 13 Sep 2018 23:44:42 +0200 Subject: Enable new error model in gmkin Also improve the irls fitting of the error model and add a test for FOCUS_2006_C where the second component of the error model is zero --- GNUmakefile | 1 - NEWS.md | 6 ++ R/mkinfit.R | 27 ++++-- README.html | 207 ++++++++++++++++++++++++++++++++++++++++++--- README.md | 2 +- build.log | 2 +- check.log | 5 +- man/mkinfit.Rd | 5 +- tests/testthat/test_irls.R | 4 + 9 files changed, 234 insertions(+), 25 deletions(-) diff --git a/GNUmakefile b/GNUmakefile index 6584722c..19c74bd8 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -23,7 +23,6 @@ pkgfiles = \ man/* \ NAMESPACE \ NEWS.md \ - README.html \ R/* \ tests/* \ tests/testthat* \ diff --git a/NEWS.md b/NEWS.md index 5af1eadf..55755862 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +# mkin 0.9.47.3 (2018-09-13) + +- 'mkinfit': Improve fitting the error model for reweight.method = 'tc' + +- Test that FOCUS_2006_C can be evaluated with DFOP and reweight.method = 'tc' + # mkin 0.9.47.2 (2018-07-19) - 'sigma_twocomp': Rename 'sigma_rl' to 'sigma_twocomp' as the Rocke and Lorenzato model assumes lognormal distribution for large y. Correct references to the Rocke and Lorenzato model accordingly. diff --git a/R/mkinfit.R b/R/mkinfit.R index 92b455f7..ecbef67e 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -37,7 +37,7 @@ mkinfit <- function(mkinmod, observed, transform_fractions = TRUE, plot = FALSE, quiet = FALSE, err = NULL, - weight = c("none", "std", "mean", "tc"), + weight = c("none", "manual", "std", "mean", "tc"), tc = c(sigma_low = 0.5, rsd_high = 0.07), scaleVar = FALSE, atol = 1e-8, rtol = 1e-10, n.outtimes = 100, @@ -415,10 +415,17 @@ mkinfit <- function(mkinmod, observed, # We need unweighted residuals to update the weighting cost_tmp <- cost(fit$par) - tc_fit <- nls(abs(res.unweighted) ~ sigma_twocomp(obs, sigma_low, rsd_high), - start = list(sigma_low = tc["sigma_low"], rsd_high = tc["rsd_high"]), - data = cost_tmp$residuals, - algorithm = "port") + tc_fit <- try( + nls(abs(res.unweighted) ~ sigma_twocomp(obs, sigma_low, rsd_high), + start = list(sigma_low = tc["sigma_low"], rsd_high = tc["rsd_high"]), + data = cost_tmp$residuals, + lower = 0, + algorithm = "port")) + + if (inherits(tc_fit, "try-error")) { + stop("Estimation of the two error model components failed for the initial fit.\n", + "Try without reweighting or with reweight.method = 'obs'.") + } tc_fitted <- coef(tc_fit) if(!quiet) { @@ -453,10 +460,16 @@ mkinfit <- function(mkinmod, observed, if (reweight.method == "tc") { cost_tmp <- cost(fit$par) - tc_fit <- nls(abs(res.unweighted) ~ sigma_twocomp(obs, sigma_low, rsd_high), + tc_fit <- try(nls(abs(res.unweighted) ~ sigma_twocomp(obs, sigma_low, rsd_high), start = as.list(tc_fitted), data = cost_tmp$residuals, - algorithm = "port") + lower = 0, + algorithm = "port")) + + if (inherits(tc_fit, "try-error")) { + stop("Estimation of the two error model components failed during reweighting.\n", + "Try without reweighting or with reweight.method = 'obs'.") + } tc_fitted <- coef(tc_fit) sr_new <- tc_fitted diff --git a/README.html b/README.html index f71dcfbc..411eea10 100644 --- a/README.html +++ b/README.html @@ -11,17 +11,202 @@ - +README.utf8 - + - - - - - - - + + + + + + +