Merge pull request #95 from piLaboratory/revert-93-revert-86-distr

Revert "Revert "Fixes redundancy in the distr slot""

.gitattributes

@@ -0,0 +1 @@
+src/gauss* linguist-vendored=true

NAMESPACE

@@ -16,7 +16,7 @@ export(fitbs, fitgamma, fitgeom, fitgs, fitlnorm, fitls, fitmand, fitmzsm,
 ## General fitting and ploting  functions
 export(fitrad, fitsad, octav, octavpred, radpred, pprad, ppsad, qqrad, qqsad, rad)
 ## Accessory functions
-export(plotprofmle, pred.logser, dtrunc, ptrunc, qtrunc, rsad, updatesad, updaterad)
+export(plotprofmle, pred.logser, dtrunc, ptrunc, qtrunc, rsad, updatesad, updaterad, distr)
 ## Explicit classes and methods export
 exportClasses(rad, octav, fitsad, fitrad)
 exportMethods(plot, points, lines, octavpred, radpred, qqsad, qqrad, plotprofmle)

R/distr.R

@@ -0,0 +1,24 @@
+#' Continuous or discrete distributions
+#'
+#' Checks if the distribution is continuous or discrete
+#'
+#' Returns whether a given distribution (used for a sad or rad model) is "discrete"
+#' or "continuous". The name is compared to a list of known distributions, so distributions
+#' not used in the sads package will return "NA" with a warning.
+#' 
+#' In the package sads up to version 0.2.3, the user was required to explicitly set 
+#' a distr argument in some calls to \code{\link{radpred}} and \code{\link{qqsad}}. Now
+#' this is handled automatically, and attempts to set the "distr" argument explicitly are ignored.
+#' 
+#' @param distribution Character. The name of the distribution ("geom" for "fitgeom", "weibull" for "fitweibull", etc.
+distr <- function(distribution) {
+  if (class(distribution)!="character") stop("Distribution must be from class character")
+  if (distribution %in% c("bs", "lnorm", "gamma", "pareto", "weibull"))
+    return("continuous")
+  if (distribution %in% c("gs", "geom", "rbs", "power", "poilog", "nbinom", "mzsm", "mand", 
+                          "ls", "volkov", "zipf"))
+    return("discrete")
+  # if arrived here...
+  warning(paste("Unknown distribution:", distribution))
+  return(NA)
+}

R/fitbs.R

@@ -15,6 +15,6 @@ fitbs <- function(x, trunc, ...){
 	result <- do.call("mle2", c(list(minuslogl=LL, data = list(x = x), fixed=list(N=n, S=s), eval.only=TRUE), dots))
 	#BUGFIX: the show method on mle2 class always expects the "convergence" slot to be set.
 	result@details$convergence = 0 
-    new("fitsad", result, sad = "bs", distr = "C", trunc = ifelse(missing(trunc), NaN, trunc))
+    new("fitsad", result, sad = "bs", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc))
   }
 }

R/fitgamma.R

@@ -42,5 +42,5 @@ fitgamma <- function(x, trunc, start.value, ...){
       LL <- function(shape, rate) -sum(dtrunc("gamma", x = x, coef = list(shape = shape, rate = rate), trunc = trunc, log = TRUE))
   }  
   result <- do.call("mle2", c(list(LL, start = list(shape = ka, rate = 1/theta), data = list(x = x)), dots))
-  new("fitsad", result, sad="gamma", distr = "C", trunc = ifelse(missing(trunc), NaN, trunc)) 
+  new("fitsad", result, sad="gamma", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc)) 
 }

R/fitgeom.R

@@ -15,5 +15,5 @@ fitgeom <- function(x, trunc = 0, start.value, ...){
   }
         ##result <- do.call("mle2", c(list(LL, start = list(prob = phat), data = list(x = x), method = "Brent", lower = 0, upper = 1), dots)) ## Brent method does not converge sometimes.
         result <- do.call("mle2", c(list(LL, start = list(prob = phat), data = list(x = x)), dots))  
-  new("fitsad", result, sad = "geom", distr = "D", trunc = ifelse(is.null(trunc), NaN, trunc))
+  new("fitsad", result, sad = "geom", distr = distr.depr, trunc = ifelse(is.null(trunc), NaN, trunc))
 }

R/fitgs.R

@@ -15,5 +15,5 @@ fitgs <- function(x, trunc, ...){
     LL <- function(S, k) -sum(dgs(y, k, S, log = TRUE))
   }
   result <- do.call("mle2", c(list(LL, start = list(k=0.01), data = list(x = y), fixed=list(S=S), method = "Brent", lower = 1e-16, upper = 1-1e-16), dots))
-  new("fitrad", result, rad = "gs", distr = "D", trunc = ifelse(missing(trunc), NaN, trunc), rad.tab=rad.tab)
+  new("fitrad", result, rad = "gs", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc), rad.tab=rad.tab)
 }

R/fitlnorm.R

@@ -16,5 +16,5 @@ fitlnorm <- function(x, trunc, start.value, ...){
     LL <- function(meanlog, sdlog) -sum(dtrunc("lnorm", x, coef = list(meanlog = meanlog, sdlog = sdlog), trunc = trunc, log = TRUE))
   }  
   result <- do.call("mle2", c(list(LL, start = list(meanlog = meanlog, sdlog = sdlog), data = list(x = x)), dots))
-  new("fitsad", result, sad="lnorm", distr = "C", trunc = ifelse(missing(trunc), NaN, trunc)) 
+  new("fitsad", result, sad="lnorm", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc)) 
 }

R/fitls.R

@@ -24,5 +24,5 @@ fitls <- function(x, trunc, start.value, upper = length(x), ...){
   result <- do.call("mle2", c(list(LL, start = list(alpha = alfa), data = list(x = x), fixed=list(N=N), method = "Brent", lower = 0, upper = upper), dots))
   if(abs(as.numeric(result@coef) - upper) < 0.0000001)
     warning("mle equal to upper bound provided. \n Try new value for the 'upper' argument")
-  new("fitsad", result, sad = "ls", distr = "D", trunc = ifelse(missing(trunc), NaN, trunc))
+  new("fitsad", result, sad = "ls", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc))
 }

R/fitmand.R

@@ -24,5 +24,5 @@ fitmand <- function(x, trunc, start.value, ...){
     LL <- function(N, s, v) -sum(dtrunc("mand", x = y, coef = list(N = N, s = s, v = v), trunc = trunc, log = TRUE))
   }
   result <- do.call("mle2", c(list(LL, start = list(s = shat, v = vhat), fixed=list(N=N), data = list(x = y)), dots))
-  new("fitrad", result, rad="mand", distr = "D", trunc = ifelse(missing(trunc), NaN, trunc), rad.tab=rad.tab)
+  new("fitrad", result, rad="mand", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc), rad.tab=rad.tab)
 }

R/fitmzsm.R

@@ -19,5 +19,5 @@ fitmzsm <- function(x, trunc, start.value, upper = length(x), ...){
     }  
   result <- do.call("mle2", c(list(LL, start = list(theta = thetahat), fixed=list(J=sum(x)), data = list(x = x), method ="Brent", lower=0.001, upper=upper), dots))
   if(abs(as.numeric(result@coef) - upper) < 0.0000001) warning("mle equal to upper bound provided. \n Try value for the 'upper' arguent")
-  new("fitsad", result, sad="mzsm", distr = "D", trunc = ifelse(missing(trunc), NaN, trunc)) 
+  new("fitsad", result, sad="mzsm", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc)) 
 }

R/fitnbinom.R

@@ -17,5 +17,5 @@ fitnbinom <- function(x, trunc=0, start.value, ...){
     LL <- function(size, mu) -sum(dtrunc("nbinom", x = x, coef = list(size=size, mu=mu), trunc = trunc, log = TRUE))
   }
   result <- do.call("mle2", c(list(LL, start = list(size = sizehat, mu = muhat), data = list(x = x)), dots))
-  new("fitsad", result, sad="nbinom", distr = "D", trunc = ifelse(is.null(trunc), NaN, trunc))
+  new("fitsad", result, sad="nbinom", distr = distr.depr, trunc = ifelse(is.null(trunc), NaN, trunc))
 }

R/fitpareto.R

@@ -18,5 +18,5 @@ fitpareto <- function(x, trunc, start.value, upper = 20, ...){
                  data = list(x = x), method = "Brent", lower = 0, upper = upper), dots))
   if(abs(as.numeric(result@coef) - upper) < 0.001) 
     warning("mle equal to upper bound provided. \n Try value for the 'upper' argument")
-  new("fitsad", result, sad="pareto", distr = "C", trunc = ifelse(missing(trunc), NaN, trunc)) 
+  new("fitsad", result, sad="pareto", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc)) 
 }

R/fitpoilog.R

@@ -15,5 +15,5 @@ fitpoilog <- function(x, trunc = 0, ...){
         LL <- function(mu, sig) -sum(dpoilog(x, mu, sig, log = TRUE))
     }
     result <- do.call("mle2", c(list(LL, start = as.list(pl.par), data = list(x = x)), dots))
-    new("fitsad", result, sad="poilog", distr = "D", trunc = ifelse(is.null(trunc), NaN, trunc))
+    new("fitsad", result, sad="poilog", distr = distr.depr, trunc = ifelse(is.null(trunc), NaN, trunc))
 }

R/fitpower.R

@@ -16,5 +16,5 @@ fitpower <- function(x, trunc, start.value, upper = 20, ...){
   }
   result <- do.call("mle2", c(list(LL, start = list(s = shat), data = list(x = x), method = "Brent", lower = 1, upper = upper), dots))
   if(abs(as.numeric(result@coef) - upper) < 0.0000001) warning("mle equal to upper bound provided. \n Try value for the 'upper' arguent")
-  new("fitsad", result, sad = "power", distr = "D", trunc = ifelse(missing(trunc), NaN, trunc))
+  new("fitsad", result, sad = "power", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc))
 }

R/fitrbs.R

@@ -20,5 +20,5 @@ fitrbs <- function(x, trunc, ...){
   result <- do.call("mle2", c(list(LL, start = list(N=N, S = S), data = list(x = y), fixed=list(N=N, S=S), eval.only=T), dots))
   #BUGFIX: the show method on mle2 class always expects the "convergence" slot to be set.
   result@details$convergence = 0 
-  new("fitrad", result, rad="rbs", distr = "D", trunc = ifelse(missing(trunc), NaN, trunc), rad.tab=rad.tab)
+  new("fitrad", result, rad="rbs", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc), rad.tab=rad.tab)
 }

R/fitvolkov.R

@@ -26,5 +26,5 @@ fitvolkov <- function(x, trunc, start.value, ...){
     }
   }
   result <- do.call("mle2", c(list(minuslogl=LL, start = list(theta = thetahat, m = mhat), fixed=list(J=sum(x)), data = list(x = x)), dots))  
-  new("fitsad", result, sad="volkov", distr = "D", trunc = ifelse(missing(trunc), NaN, trunc))
+  new("fitsad", result, sad="volkov", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc))
 }

R/fitweibull.R

@@ -20,5 +20,5 @@ fitweibull <- function(x, trunc, start.value, ...){
     LL <- function(shape, scale) -sum(dtrunc("weibull", x = x, coef = list(shape, scale), trunc = trunc, log = TRUE))
   }  
   result <- do.call("mle2", c(list(LL, start = list(shape = ka, scale = theta), data = list(x = x)), dots))
-  new("fitsad", result, sad="weibull", distr = "C", trunc = ifelse(missing(trunc), NaN, trunc)) 
+  new("fitsad", result, sad="weibull", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc)) 
 }

R/fitzipf.R

@@ -28,5 +28,5 @@ fitzipf <- function(x, N, trunc, start.value, upper = 20, ...){
   result <- do.call("mle2", c(list(LL, start = list(s = sss), data = list(x = y), fixed=list(N=N), method = "Brent", lower = 0, upper = upper), dots))
   if(abs(as.numeric(result@coef) - upper) < 0.001)
     warning("mle equal to upper bound provided. \n Try increase value for the 'upper' argument")
-  new("fitrad", result, rad="zipf", distr = "D", trunc = ifelse(missing(trunc), NaN, trunc), rad.tab=rad.tab)
+  new("fitrad", result, rad="zipf", distr = distr.depr, trunc = ifelse(missing(trunc), NaN, trunc), rad.tab=rad.tab)
 }

R/octav.R

@@ -20,12 +20,18 @@ octav <- function(x, oct, preston=FALSE){
   N <- 2^(oct)
   oct.hist <- hist(y, breaks=c(0,N), plot=FALSE)
   res <- data.frame(octave = oct, upper = oct.hist$breaks[-1], Freq = oct.hist$counts)
-  if(preston){
-    j <- N[-length(N)]
-    w <- y[y%in%j]
-    ties <- table(factor(w, levels=j))
-    res[-1, 3] <- res[-1, 3]+ties/2
-    res[-length(N), 3] <- res[-length(N), 3]-ties/2
-  }
+  if(preston) res <- prestonfy(res, y)
   new("octav", res)
 }
+
+# Helper function to create Preston octaves. Is also used by octavpred
+prestonfy <- function(res, y) {
+  N <- 2^(res$octave)
+  j <- N[-length(N)]
+  w <- y[y%in%j]
+  ties <- table(factor(w, levels=j))
+  res[-1, 3] <- res[-1, 3]+ties/2
+  res[-length(N), 3] <- res[-length(N), 3]-ties/2
+  return(res)
+}
+

R/sads-classes.R

@@ -22,3 +22,5 @@ setClass("rad", representation("data.frame"), validity = function(object) {
 setClass("fitsad", representation("mle2", sad="character", distr="character", trunc="numeric"))
 setClass("fitrad", representation("mle2", rad="character", distr="character", trunc="numeric", rad.tab="rad"))
 #setClass("fitsadlist", representation("list"))
+
+distr.depr <- "The 'distr' slot of fitrad and fitsad objects have been deprecated. Please see ?distr"