Update pb-infinite-NaN-gof-metrics.R

lbbe-software · Oct 22, 2024 · 438dae6 · 438dae6
1 parent ebbab5a
commit 438dae6
Showing 1 changed file with 46 additions and 18 deletions.
diff --git a/share/pb-infinite-NaN-gof-metrics.R b/share/pb-infinite-NaN-gof-metrics.R
@@ -7,7 +7,7 @@ mygof <- function(gof, echo=FALSE)
     { 
       n <- length(obs)
       s <- sort(obs)
-      theop <- do.call(pdistnam,c(list(s),as.list(par),as.list(fix.arg)))
+      theop <- do.call(pdistnam, c(list(s), as.list(par), as.list(fix.arg)))
       1/(12*n) + sum( ( theop - (2 * 1:n - 1)/(2 * n) )^2 )
     }
   }else if (gof == "KS")
@@ -18,70 +18,90 @@ mygof <- function(gof, echo=FALSE)
       s <- sort(obs)
       obspu <- seq(1,n)/n
       obspl <- seq(0,n-1)/n
-      theop <- do.call(pdistnam,c(list(s),as.list(par),as.list(fix.arg)))
-      max(pmax(abs(theop-obspu),abs(theop-obspl)))
+      theop <- do.call(pdistnam, c(list(s), as.list(par), as.list(fix.arg)))
+      max(pmax(abs(theop-obspu), abs(theop-obspl)))
     }
   }else if (gof == "AD")
   {
     fnobj <- function(par, fix.arg, obs, pdistnam)
     { 
       n <- length(obs)
       s <- sort(obs)
-      theop <- do.call(pdistnam,c(list(s),as.list(par),as.list(fix.arg)))
+      theop <- do.call(pdistnam, c(list(s), as.list(par), as.list(fix.arg)))
       if(echo)
         print(cbind(2 * 1:n - 1, s, theop, 1-rev(theop), log(theop) + log(1 - rev(theop)), (2 * 1:n - 1) * (log(theop) + log(1 - rev(theop))) ),
               digit=22)
-      - n - mean( (2 * 1:n - 1) * (log(theop) + log(1 - rev(theop))) ) 
+      ilogpi <- log(theop * (1 - rev(theop))) * (2 * 1:n - 1)
+      idx <- is.finite(ilogpi)
+      - sum(idx) - mean( ilogpi[idx] ) 
     }
   }else if (gof == "ADR")
   {
     fnobj <- function(par, fix.arg, obs, pdistnam)
     { 
       n <- length(obs)
       s <- sort(obs)
-      theop <- do.call(pdistnam,c(list(s),as.list(par),as.list(fix.arg)))
+      theop <- do.call(pdistnam, c(list(s), as.list(par), as.list(fix.arg)))
       if(echo)
         print(cbind(2 * 1:n - 1, s, 1-rev(theop), log(1 - rev(theop)), (2 * 1:n - 1) * log(1 - rev(theop)) ))
-      n/2 - 2 * sum(theop) - mean ( (2 * 1:n - 1) * log(1 - rev(theop)) )
+
+      ilogpi <- log(1 - rev(theop)) * (2 * 1:n - 1)
+      idx <- is.finite(ilogpi)
+      sum(idx)/2 - 2 * sum(theop[idx]) - mean ( ilogpi[idx] )
     }
   }else if (gof == "ADL")
   {
     fnobj <- function(par, fix.arg, obs, pdistnam)
     { 
       n <- length(obs)
       s <- sort(obs)
-      theop <- do.call(pdistnam,c(list(s),as.list(par),as.list(fix.arg)))
-      -3*n/2 + 2 * sum(theop) - mean ( (2 * 1:n - 1) * log(theop) )
+      theop <- do.call(pdistnam, c(list(s), as.list(par), as.list(fix.arg)))
+      ilogpi <- (2 * 1:n - 1) * log(theop)
+      idx <- is.finite(ilogpi)
+      -3*sum(idx)/2 + 2 * sum(theop[idx]) - mean ( ilogpi[idx] )
     }
   }else if (gof == "AD2R")
   {
     fnobj <- function(par, fix.arg, obs, pdistnam)
     { 
       n <- length(obs)
       s <- sort(obs)
-      theop <- do.call(pdistnam,c(list(s),as.list(par),as.list(fix.arg)))
-      2 * sum(log(1 - theop)) + mean ( (2 * 1:n - 1) / (1 - rev(theop)) )
+      theop <- do.call(pdistnam, c(list(s), as.list(par), as.list(fix.arg)))
+
+      logpi <- log(1 - theop)
+      i1pi2 <- (2 * 1:n - 1) / (1 - rev(theop))
+      idx <- is.finite(logpi) & is.finite(i1pi2)
+
+      2 * sum(logpi[idx]) + mean( i1pi2[idx] )
     }
   }else if (gof == "AD2L")
   {
     fnobj <- function(par, fix.arg, obs, pdistnam)
     { 
       n <- length(obs)
       s <- sort(obs)
-      theop <- do.call(pdistnam,c(list(s),as.list(par),as.list(fix.arg)))
-      2 * sum(log(theop)) + mean ( (2 * 1:n - 1) / theop )
+      theop <- do.call(pdistnam, c(list(s), as.list(par), as.list(fix.arg)))
+      logpi <- log(theop)
+      i1pi <- (2 * 1:n - 1) / theop 
+      idx <- is.finite(logpi) & is.finite(i1pi)
+      2 * sum( logpi[idx] ) + mean ( i1pi[idx] )
     }
   }else if (gof == "AD2")
   {
     fnobj <- function(par, fix.arg, obs, pdistnam)
     { 
       n <- length(obs)
       s <- sort(obs)
-      theop <- do.call(pdistnam,c(list(s),as.list(par),as.list(fix.arg)))
+      theop <- do.call(pdistnam, c(list(s), as.list(par), as.list(fix.arg)))
       if(echo)
         print(cbind(theop, log(theop), log(1-theop), log(theop) + log(1 - theop), ((2 * 1:n - 1) / (1 - rev(theop)))))
-      2 * sum(log(theop) + log(1 - theop) ) + 
-        mean ( ((2 * 1:n - 1) / theop) + ((2 * 1:n - 1) / (1 - rev(theop))) )
+
+      logpi <- log(theop * (1 - theop)) 
+      i1pi <- (2 * 1:n - 1) / theop
+      i1pi2 <- (2 * 1:n - 1) / (1 - rev(theop))
+      idx <- is.finite(logpi) & is.finite(i1pi) & is.finite(i1pi2)
+
+      2 * sum( logpi[idx] ) + mean ( i1pi[idx]  + i1pi2[idx] )
     }
   }else
     fnobj <- NULL
@@ -90,7 +110,7 @@ mygof <- function(gof, echo=FALSE)
 }
 
 pdistnam <- pweibull
-x <- c(rep(1, 13), 10)
+x <- c(rep(1, 14), 10)
 par <- fitdistrplus:::startarg_transgamma_family(x, "weibull")
 fix.arg <- NULL
 
@@ -100,7 +120,15 @@ mygof("AD", TRUE)(par, fix.arg, x, pdistnam)
 mygof("ADR", TRUE)(par, fix.arg, x, pdistnam)
 mygof("AD2", TRUE)(par, fix.arg, x, pdistnam)
 
+optim(par, function(theta) mygof("AD", FALSE)(theta, fix.arg, x, pdistnam), control=list(trace=1, REPORT=1))
 
 optim(par, function(theta) mygof("AD2", FALSE)(theta, fix.arg, x, pdistnam), control=list(trace=1, REPORT=1))
 
-abs(mygof("AD")(as.list(unlist(par)+eps), fix.arg, x, pdistnam) - mygof("AD")(par, fix.arg, x, pdistnam))/eps
+optim(par, function(theta) mygof("ADR", FALSE)(theta, fix.arg, x, pdistnam), control=list(trace=1, REPORT=1))
+
+optim(par, function(theta) mygof("AD2R", FALSE)(theta, fix.arg, x, pdistnam), control=list(trace=1, REPORT=1))
+
+optim(par, function(theta) mygof("AD2L", FALSE)(theta, fix.arg, x, pdistnam), control=list(trace=1, REPORT=1))
+
+
+