aboutsummaryrefslogtreecommitdiff
path: root/R/mkinfit.R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2018-09-20 18:19:40 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2018-09-20 18:19:40 +0200
commit4ba08177631866e8a910f561800f8a45b23e5cb9 (patch)
treeceeff2f7600cb9bbc545dc40c6c871e91428a0d7 /R/mkinfit.R
parent1f100167f4f399062fb077ceb983f50384eccaec (diff)
Use the median absolute deviance for fitting
as the absolute value is a biased estimator for the standard deviation
Diffstat (limited to 'R/mkinfit.R')
-rw-r--r--R/mkinfit.R21
1 files changed, 16 insertions, 5 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R
index 2f05e81a..d56c663a 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -415,10 +415,16 @@ mkinfit <- function(mkinmod, observed,
# We need unweighted residuals to update the weighting
tmp_res <- cost(fit$par)$residuals
+ mad_agg <- aggregate(tmp_res$res.unweighted,
+ by = list(name = tmp_res$name, x_res = tmp_res$x),
+ FUN = function(x) mad(x, center = 0))
+ names(mad_agg) <- c("name", "x", "mad")
+ tmp_res_mad <- merge(tmp_res, mad_agg)
+
tc_fit <- try(
- nls(abs(res.unweighted) ~ sigma_twocomp(obs, sigma_low, rsd_high),
+ nls(mad ~ sigma_twocomp(mod, sigma_low, rsd_high),
start = list(sigma_low = tc["sigma_low"], rsd_high = tc["rsd_high"]),
- data = tmp_res,
+ data = tmp_res_mad,
lower = 0,
algorithm = "port"))
@@ -429,7 +435,7 @@ mkinfit <- function(mkinmod, observed,
tc_fitted <- coef(tc_fit)
if(!quiet) {
- cat("IRLS based on variance estimates according to the Rocke and Lorenzato model\n")
+ cat("IRLS based on variance estimates according to the two component error model\n")
cat("Initial variance components are:\n")
print(signif(tc_fitted))
}
@@ -459,11 +465,16 @@ mkinfit <- function(mkinmod, observed,
}
if (reweight.method == "tc") {
tmp_res <- cost(fit$par)$residuals
+ mad_agg <- aggregate(tmp_res$res.unweighted,
+ by = list(name = tmp_res$name, x_res = tmp_res$x),
+ FUN = function(x) mad(x, center = 0))
+ names(mad_agg) <- c("name", "x", "mad")
+ tmp_res_mad <- merge(tmp_res, mad_agg)
tc_fit <- try(
- nls(abs(res.unweighted) ~ sigma_twocomp(obs, sigma_low, rsd_high),
+ nls(mad ~ sigma_twocomp(mod, sigma_low, rsd_high),
start = list(sigma_low = tc["sigma_low"], rsd_high = tc["rsd_high"]),
- data = tmp_res,
+ data = tmp_res_mad,
lower = 0,
algorithm = "port"))

Contact - Imprint