diff options
| author | Johannes Ranke <jranke@uni-bremen.de> | 2018-09-20 18:19:40 +0200 |
|---|---|---|
| committer | Johannes Ranke <jranke@uni-bremen.de> | 2018-09-20 18:19:40 +0200 |
| commit | 4ba08177631866e8a910f561800f8a45b23e5cb9 (patch) | |
| tree | ceeff2f7600cb9bbc545dc40c6c871e91428a0d7 | |
| parent | 1f100167f4f399062fb077ceb983f50384eccaec (diff) | |
Use the median absolute deviance for fitting
as the absolute value is a biased estimator for the standard deviation
| -rw-r--r-- | R/mkinfit.R | 21 |
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"))
|