From 59e2c3b5ce470ace53b93ad5e6e7d9acbfe86f21 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Fri, 25 Apr 2014 18:43:51 +0200
Subject: McCall 2,4,5-T dataset, lower bound for transform_rates=FALSE

---
 R/mkinfit.R            | 13 +++++++++++--
 R/transform_odeparms.R |  4 ++--
 2 files changed, 13 insertions(+), 4 deletions(-)

(limited to 'R')

diff --git a/R/mkinfit.R b/R/mkinfit.R
index a623146..16b4ff3 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -208,8 +208,17 @@ mkinfit <- function(mkinmod, observed,
     return(mC)
   }
 
+  lower <- rep(-Inf, length(c(state.ini.optim, parms.optim)))
+  names(lower) <- c(names(state.ini.optim), names(parms.optim))
+  if (!transform_rates) {
+    index_k <- grep("^k_", names(lower))
+    lower[index_k] <- 0
+    other_rate_parms <- intersect(c("alpha", "beta", "k1", "k2"), names(lower))
+    lower[other_rate_parms] <- 0
+  }
+
   fit <- modFit(cost, c(state.ini.optim, parms.optim), 
-                method = method.modFit, control = control.modFit, ...)
+                method = method.modFit, control = control.modFit, lower = lower, ...)
 
   # Reiterate the fit until convergence of the variance components (IRLS)
   # if requested by the user
@@ -234,7 +243,7 @@ mkinfit <- function(mkinmod, observed,
       sigma.old <- sqrt(fit$var_ms_unweighted)
       observed[err] <- sqrt(fit$var_ms_unweighted)[as.character(observed$name)]
       fit <- modFit(cost, fit$par, method = method.modFit,
-                    control = control.modFit, ...)
+                    control = control.modFit, lower = lower, ...)
       reweight.diff = sum((sqrt(fit$var_ms_unweighted) - sigma.old)^2)
       if (!quiet) {
         cat("Iteration", n.iter, "yields variance estimates:\n")
diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
index 31200c7..24de074 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -44,8 +44,8 @@ transform_odeparms <- function(parms, mod_vars,
     }
   }
 
-  # Transform parameters also for FOMC, DFOP and HS models if requested
-  for (pname in c("alpha", "beta", "k1", "k2", "tb")) {
+  # Transform rate and fraction parameters also for FOMC, DFOP and HS models if requested
+  for (pname in c("alpha", "beta", "k1", "k2")) {
     if (!is.na(parms[pname])) {
       transparms[pname] <- ifelse(transform_rates, log(parms[pname]), parms[pname])
       transparms[pname] <- ifelse(transform_rates, log(parms[pname]), parms[pname])
-- 
cgit v1.2.1