From 70591022c07f0e8fb4dd67789b7c8d78af8ebc18 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Thu, 2 May 2019 13:17:05 +0200
Subject: Better initials for error model parameters

- Also make it possible to specify initial values for error model
parameters.
- Run tests
- Rebuild docs
---
 R/mkinfit.R | 29 +++++++++++++++++++++--------
 1 file changed, 21 insertions(+), 8 deletions(-)

(limited to 'R/mkinfit.R')

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 754e72b8..dca71ecf 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -21,6 +21,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value"))
 mkinfit <- function(mkinmod, observed,
   parms.ini = "auto",
   state.ini = "auto",
+  err.ini = "auto",
   fixed_parms = NULL,
   fixed_initials = names(mkinmod$diffs)[-1],
   from_max_mean = FALSE,
@@ -247,6 +248,26 @@ mkinfit <- function(mkinmod, observed,
     "obs" = paste0("sigma_", obs_vars),
     "tc" = c("sigma_low", "rsd_high"))
 
+  # Define starting values for the error model
+  if (err.ini[1] != "auto") {
+    if (!identical(names(err.ini), errparm_names)) {
+      stop("Please supply initial values for error model components ", paste(errparm_names, collapse = ", "))
+    } else {
+      errparms = err.ini
+    }
+  } else {
+    if (err_mod == "const") {
+      errparms = 3
+    }
+    if (err_mod == "obs") {
+      errparms = rep(3, length(obs_vars))
+    }
+    if (err_mod == "tc") {
+      errparms <- c(sigma_low = 3, rsd_high = 0.01)
+    }
+    names(errparms) <- errparm_names
+  }
+
   # Define outtimes for model solution.
   # Include time points at which observed data are available
   outtimes = sort(unique(c(observed$time, seq(min(observed$time),
@@ -407,14 +428,6 @@ mkinfit <- function(mkinmod, observed,
       fit <- fit.ols
       fit$logLik <- - nlogLik(c(fit$par, sigma = sigma_mle), OLS = FALSE)
     } else {
-      if (err_mod == "obs") {
-        errparms = rep(3, length(obs_vars))
-      }
-      if (err_mod == "tc") {
-        errparms <- c(sigma_low = 0.5, rsd_high = 0.07)
-      }
-      names(errparms) <- errparm_names
-
       fit <- nlminb(c(state.ini.optim, transparms.optim, errparms),
                     nlogLik, control = control,
                     lower = lower, upper = upper, ...)
-- 
cgit v1.2.1