Merge pull request #51 from tripartio/master

Add support for discrete Poisson distribution + plotting of CDFs and quantile functions.

DESCRIPTION

@@ -29,7 +29,7 @@ Suggests:
     copula,
     dplyr,
     covr
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 VignetteBuilder: knitr
 URL: https://github.com/JonasMoss/univariateML, https://jonasmoss.github.io/univariateML/
 BugReports: https://github.com/JonasMoss/univariateML/issues

NAMESPACE

@@ -31,6 +31,7 @@ export(mllomax)
 export(mlnaka)
 export(mlnorm)
 export(mlpareto)
+export(mlpois)
 export(mlpower)
 export(mlrayleigh)
 export(mlsged)

R/generics.R

@@ -5,35 +5,60 @@
 #' @param points Boolean; should points be plotted by default?
 #' @keywords internal
 
-plot_wrangler <- function(x, range, points = FALSE, ...) {
+plot_wrangler <- function(x, range, points = FALSE, kind, ...) {
+  continuous <- if(is.null(attr(x,"continuous"))) TRUE else FALSE
   support <- attr(x, "support")
   if (is.null(range)) {
     if (abs(support[1]) + abs(support[2]) < Inf) {
       limits <- support
     } else if (abs(support[1]) == 0 & abs(support[2]) == Inf) {
-      limits <- c(0, qml(0.99, x))
+      limits <- c(0, qml(0.995, x))
     } else {
-      limits <- qml(c(0.01, 0.99), x)
+      limits <- qml(c(0.005, 0.995), x)
     }
 
-    range <- seq(limits[1], limits[2], length.out = 1000)
+    limits_untransformed <- limits
+    if(kind == "q") limits <- pml(limits, x)
+
+    range <- if (continuous) {
+      seq(limits[1], limits[2], length.out = 1000)
+    } else {
+      if(kind == "q") {
+        pml(seq(limits_untransformed[1], limits_untransformed[2]), x)
+      } else {
+        seq(limits[1], limits[2])
+      }
+    }
   }
 
+  ylab <- list(d="Density",p="Cumulative probability",q="Quantile")
+  xlab <- list(d="x",p="Quantile",q="Cumulative probability")
   defaults <- list(
     type = if (points) "p" else "l",
     main = paste0(attr(x, "model"), " model"),
-    ylab = "Density",
-    xlab = "x",
+    ylab = ylab[[kind]],
+    xlab = xlab[[kind]],
     lwd = 1
   )
 
+  if(!continuous) {
+    defaults$pch = 20
+    defaults$type = if (points) "p" else "b"
+  }
+
   args <- listmerge(
     x = defaults,
     y = list(...)
   )
 
   args$x <- range
-  args$y <- dml(args$x, x)
+  args$y <- if(kind == "d") {
+    dml(args$x, x)
+  } else if (kind == "p") {
+    pml(args$x, x)
+  } else {
+    qml(args$x, x)
+  }
   args
 }
 
@@ -45,6 +70,7 @@ plot_wrangler <- function(x, range, points = FALSE, ...) {
 #' @export
 #' @param x a `univariateML` object.
 #' @param range range of `x` values to plot, i.e. `c(lower, upper)`.
+#' @param kind can be `density`, `probability`, or `quantile`.
 #' @param ... parameters passed to `plot`, `lines`, or `points`.
 #' @return An invisible copy of `x`.
 #' @examples
@@ -53,24 +79,27 @@ plot_wrangler <- function(x, range, points = FALSE, ...) {
 #' rug(datasets::precip)
 #' @export
 #'
-plot.univariateML <- function(x, range = NULL, ...) {
-  args <- plot_wrangler(x, range, points = FALSE, ...)
+plot.univariateML <- function(x, range = NULL, kind = c("d", "p", "q"), ...) {
+  kind <- match.arg(kind)
+  args <- plot_wrangler(x, range, points = FALSE, kind = kind, ...)
   do.call(graphics::plot, args)
   invisible(x)
 }
 
 #' @export
 #' @rdname plot.univariateML
-lines.univariateML <- function(x, range = NULL, ...) {
-  args <- plot_wrangler(x, range, points = FALSE, ...)
+lines.univariateML <- function(x, range = NULL, kind = c("d", "p", "q"), ...) {
+  kind <- match.arg(kind)
+  args <- plot_wrangler(x, range, points = FALSE, kind = kind, ...)
   do.call(graphics::lines, args)
   invisible(x)
 }
 
 #' @export
 #' @rdname plot.univariateML
-points.univariateML <- function(x, range = NULL, ...) {
-  args <- plot_wrangler(x, range, points = TRUE, ...)
+points.univariateML <- function(x, range = NULL, kind = c("d", "p", "q"), ...) {
+  kind <- match.arg(kind)
+  args <- plot_wrangler(x, range, points = TRUE, kind = kind,...)
   do.call(graphics::points, args)
   invisible(x)
 }

R/helpers.R

@@ -64,12 +64,16 @@ to_univariateML <- function(y, obj) {
 #' Wrangles arguments for use in ppml and qqml functions.
 #'
 #' @param y The input data.
-#' @param obj Function or `"univariateML"` object.
+#' @param obj Function or continuous `"univariateML"` object.
 #' @param datax logical; if true, plots the data on the x axis.
 #' @param ... Arguments passed to `plot` or `points` down the line.
 #' @keywords internal
 
 ppqq_wrangler <- function(y, obj, datax, pp, ...) {
+
+  if (!is.null(attr(obj, "continuous"))) {
+    stop("QQ and PP plots are only supported for continuous distributions.")
+  }
   ## Nas are removed by default in this function, following qqplot.
 
   y <- y[!is.na(y)]

R/mllnorm.R

@@ -11,7 +11,7 @@
 #' @param x a (non-empty) numeric vector of data values.
 #' @param na.rm logical. Should missing values be removed?
 #' @param ... currently affects nothing.
-#' @return `mllonorm` returns an object of [class][base::class] `univariateML`.
+#' @return `mllnorm` returns an object of [class][base::class] `univariateML`.
 #'    This is a named numeric vector with maximum likelihood estimates for
 #'    `meanlog` and `sdlog` and the following attributes:
 #'     \item{`model`}{The name of the model.}

R/mllogitnorm.R

@@ -22,7 +22,7 @@
 #'     \item{`call`}{The call as captured my `match.call`}
 #' @examples
 #' AIC(mllogitnorm(USArrests$Rape / 100))
-#' @seealso link[dlogitnorm]{dlogitnorm}for the normal density.
+#' @seealso [Normal][stats::dnorm] for the normal density.
 #' @references Atchison, J., & Shen, S. M. (1980). Logistic-normal
 #' distributions: Some properties and uses. Biometrika, 67(2), 261-272.
 #' @export

R/mlpois.R

@@ -0,0 +1,42 @@
+#' Poisson distribution maximum likelihood estimation
+#'
+#' The maximum likelihood estimate of `lambda` is the empirical mean.
+#'
+#' For the density function of the Poisson distribution see
+#' [Poisson][stats::Poisson].
+#'
+#' @param x a (non-empty) numeric vector of data values.
+#' @param na.rm logical. Should missing values be removed?
+#' @param ... currently affects nothing.
+#' @return `mlpois` returns an object of [class][base::class] `univariateML`.
+#'    This is a named numeric vector with maximum likelihood estimates for
+#'    `mean` and `sd` and the following attributes:
+#'     \item{`model`}{The name of the model.}
+#'     \item{`density`}{The density associated with the estimates.}
+#'     \item{`logLik`}{The loglikelihood at the maximum.}
+#'     \item{`support`}{The support of the density.}
+#'     \item{`n`}{The number of observations.}
+#'     \item{`call`}{The call as captured my `match.call`}
+#' @examples
+#' mlpois(precip)
+#' @seealso [Poisson][stats::Poisson] for the Poisson density.
+#' @export
+
+mlpois <- function(x, na.rm = FALSE, ...) {
+  if (na.rm) x <- x[!is.na(x)] else assertthat::assert_that(!anyNA(x))
+  ml_input_checker(x)
+  n <- length(x)
+
+  lambda <- mean(x)
+
+  object <- c(lambda = lambda)
+  class(object) <- "univariateML"
+  attr(object, "model") <- "Poisson"
+  attr(object, "continuous") = FALSE
+  attr(object, "density") <- "stats::dpois"
+  attr(object, "logLik") <- -n * lambda + sum(x) * log(lambda) - sum(lgamma(x + 1))
+  attr(object, "support") <- c(0, Inf)
+  attr(object, "n") <- length(x)
+  attr(object, "call") <- match.call()
+  object
+}

R/probability_plots.R

@@ -18,6 +18,9 @@
 #'
 #' This function is modeled after [qqnorm][stats::qqnorm].
 #'
+#' Quantile-quantile plots and probability-probability plots are only supported
+#' for continuous distributions.
+#'
 #' Graphical parameters may be given as arguments to all the functions below.
 #'
 #' @param y Numeric vector; The data to plot on the `y` axis when

README.Rmd

@@ -30,7 +30,7 @@ set.seed(313)
 ## Overview
 [`univariateML`](https://jonasmoss.github.io/univariateML/index.html) is an `R`-package for
 user-friendly maximum likelihood estimation of a
-[selection](https://jonasmoss.github.io/univariateML/articles/distributions.html) of parametric univariate densities. In addition to basic estimation capabilities,
+[selection](https://jonasmoss.github.io/univariateML/articles/distributions.html) of parametric univariate densities and probability mass functions. In addition to basic estimation capabilities,
 this package support visualization through `plot` and `qqmlplot`, model selection
 by `AIC` and `BIC`, confidence sets through the parametric bootstrap with
 `bootstrapml`, and convenience functions such as the density, distribution 
@@ -99,6 +99,7 @@ lines(mlweibull(egypt$age))
 | Logit-normal                        | `mllogitnorm`          | logitnorm  |
 | Uniform distribution                | `mlunif`               | stats      |
 | Power distribution                  | `mlpower`              | extraDistr |
+| Poisson distribution                | `mlpois`               | stats      |
 
 ## Implementations
 

README.md

@@ -17,12 +17,13 @@ developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.re
 [`univariateML`](https://jonasmoss.github.io/univariateML/index.html) is
 an `R`-package for user-friendly maximum likelihood estimation of a
 [selection](https://jonasmoss.github.io/univariateML/articles/distributions.html)
-of parametric univariate densities. In addition to basic estimation
-capabilities, this package support visualization through `plot` and
-`qqmlplot`, model selection by `AIC` and `BIC`, confidence sets through
-the parametric bootstrap with `bootstrapml`, and convenience functions
-such as the density, distribution function, quantile function, and
-random sampling at the estimated distribution parameters.
+of parametric univariate densities and probability mass functions. In
+addition to basic estimation capabilities, this package support
+visualization through `plot` and `qqmlplot`, model selection by `AIC`
+and `BIC`, confidence sets through the parametric bootstrap with
+`bootstrapml`, and convenience functions such as the density,
+distribution function, quantile function, and random sampling at the
+estimated distribution parameters.
 
 ## Installation
 
@@ -94,6 +95,7 @@ lines(mlweibull(egypt$age))
 | Logit-normal                        | `mllogitnorm`         | logitnorm  |
 | Uniform distribution                | `mlunif`              | stats      |
 | Power distribution                  | `mlpower`             | extraDistr |
+| Poisson distribution                | `mlpois`              | stats      |
 
 ## Implementations
 
@@ -110,11 +112,12 @@ x <- rbeta(500, 2, 7)
 
 microbenchmark::microbenchmark(
   univariateML = univariateML::mlbeta(x),
-  naive = nlm(function(p) -sum(dbeta(x, p[1], p[2], log = TRUE)), p = c(1, 1)))
+  naive = nlm(function(p) -sum(dbeta(x, p[1], p[2], log = TRUE)), p = c(1, 1))
+)
 #> Unit: microseconds
-#>          expr     min       lq      mean   median       uq     max neval
-#>  univariateML   259.2   348.75   557.959   447.05   536.40  5103.5   100
-#>         naive 15349.1 15978.35 16955.165 16365.45 17082.25 48941.4   100
+#>          expr     min       lq      mean   median       uq      max neval
+#>  univariateML   205.1   329.85   800.077   464.35   730.35  13938.4   100
+#>         naive 10665.9 11773.30 21428.104 15233.70 28628.65 103372.1   100
 ```
 
 The maximum likelihood estimators in this package have all been subject