Use relative eqtol

DESCRIPTION

@@ -2,7 +2,7 @@ Package: durmod
 Type: Package
 Title: Mixed Proportional Hazard Competing Risk Model
 Version: 1.1-1
-Date: 2019-07-22
+Date: 2019-07-24
 Authors@R: person("Simen", "Gaure", email="[email protected]", role=c("aut","cre"), 
         comment=c(ORCID="https://orcid.org/0000-0001-7251-8747"))
 URL: https://github.com/sgaure/durmod

R/duration.R

@@ -192,7 +192,7 @@ mphcrm <- function(formula,data,risksets=NULL,
 #' @export
 mphcrm.control <- function(...) {
   ctrl <- list(iters=50,threads=getOption('durmod.threads'),gradient=TRUE, fisher=TRUE,
-               method='BFGS', gdiff=TRUE, minprob=1e-20, eqtol=1e-4, newprob=1e-4, jobname='mphcrm', 
+               method='BFGS', gdiff=TRUE, minprob=1e-20, eqtol=1e-3, newprob=1e-4, jobname='mphcrm', 
                overshoot=0.001,
                startprob=1e-4,
                ll.improve=1e-3, e.improve=1e-3,
@@ -517,7 +517,7 @@ badpoints <- function(pset,control) {
     for(j in seq_len(i-1)) {
       if(!okpt[j]) next
       muj <- sapply(pset$parset, function(pp) pp$mu[j])
-      if(max(abs(exp(muj) - exp(mui))) < control$eqtol) {
+      if(max(abs(exp(muj) - exp(mui))/(1e-9+exp(muj)+exp(mui))) < control$eqtol) {
         mumat <- rbind(muj,mui)
         rownames(mumat) <- c(j,i)
         colnames(mumat) <- names(pset$parset)

R/utils.R

@@ -12,24 +12,21 @@
 #' geninv(x)
 #' @return A matrix of the same dimension as \code{X} is returned, the Moore-Penrose generalized inverse.
 #' @export
-geninv <- function(X, tol=.Machine$double.eps) {
-  if (length(dim(X)) > 2L || !is.numeric(X)) 
-    stop("'X' must be a numeric matrix")
-  if (!is.matrix(X))  X <- as.matrix(X)
+geninv <- function(X, tol=.Machine$double.eps^(2/3)) {
+  stopifnot(is.numeric(X), length(dim(X)) == 2, is.matrix(X))
   nm <- colnames(X)
-  Xsvd <- svd(X)
-  Positive <- Xsvd$d > max(tol * Xsvd$d[1L], 0)
-  inv <- if (all(Positive)) 
-    Xsvd$v %*% (1/Xsvd$d * t(Xsvd$u))
-  else if (!any(Positive)) 
-    array(NA, dim(X)[2L:1L])
+  s <- svd(X)
+  p <- s$d > max(tol * s$d[1L], 0)
+  inv <- if (all(p)) 
+    s$v %*% (1/s$d * t(s$u))
+  else if (!any(p)) 
+    array(0, dim(X)[2L:1L])
   else {
-    im <- Xsvd$v[, Positive, drop = FALSE] %*% ((1/Xsvd$d[Positive]) * 
-                                                 t(Xsvd$u[, Positive, drop = FALSE]))
-    im[!Positive,!Positive] <- NA
-    im
+    s$v[, p, drop = FALSE] %*% ((1/s$d[p]) * 
+                                t(s$u[, p, drop = FALSE]))
   }
-  structure(inv,badvars=nm[!Positive], dimnames=list(nm,nm))
+#  badvars <- if(any(!p)) names(collin(X))
+  structure(inv, dimnames=list(nm,nm))
 }
 
 nazero <- function(x) ifelse(is.na(x),0,x)

vignettes/whatmph.Rnw

@@ -392,10 +392,15 @@ around the program, I have collected them here with their defaults. Some of them
 
 \item{\code{minprob=\Sexpr{def$minprob}}.} A numeric. Masspoints with probability below this are removed. 
 
-\item{\code{eqtol=\Sexpr{def$eqtol}}.} A numeric. It sometimes happens that two location points in the mixed proportional
-  hazard distribution turns out to be equal. One of them can then be removed. \code{eqtol} is the
-  threshold below which the \(\ell^\infty\) distance between two location points is so small that the
-  two points are thought to be equal. Should be lowered if the hazards generally are very small.
+\item{\code{eqtol=\Sexpr{def$eqtol}}.} A numeric. It sometimes happens
+  that two location points in the mixed proportional hazard
+  distribution turns out to be equal. One of them can then be
+  removed. \code{eqtol} is the threshold below which the
+  distance between two location points is so small
+  that the two points are thought to be equal. The ``distance'' is
+  the maximum of pointwise relative differences, i.e. with two
+  vectors \(\{a_i\}_i\) and \(\{b_i\}_i\), the distance is
+  \(\max_i |a_i-b_i|/(|a_i|+|b_i|)\).
 
 \item{\code{newpoint.maxtime=\Sexpr{def$newpoint.maxtime}}.} A numeric. When
   searching for a new location point, a global search algorithm from