Add g-and-h distribution

Resolves #84

NAMESPACE

@@ -29,6 +29,7 @@ S3method(cdf,dist_exponential)
 S3method(cdf,dist_f)
 S3method(cdf,dist_gamma)
 S3method(cdf,dist_geometric)
+S3method(cdf,dist_gh)
 S3method(cdf,dist_gk)
 S3method(cdf,dist_gumbel)
 S3method(cdf,dist_hypergeometric)
@@ -109,6 +110,7 @@ S3method(density,dist_exponential)
 S3method(density,dist_f)
 S3method(density,dist_gamma)
 S3method(density,dist_geometric)
+S3method(density,dist_gh)
 S3method(density,dist_gk)
 S3method(density,dist_gumbel)
 S3method(density,dist_hypergeometric)
@@ -156,6 +158,7 @@ S3method(format,dist_exponential)
 S3method(format,dist_f)
 S3method(format,dist_gamma)
 S3method(format,dist_geometric)
+S3method(format,dist_gh)
 S3method(format,dist_gk)
 S3method(format,dist_gumbel)
 S3method(format,dist_hypergeometric)
@@ -202,6 +205,7 @@ S3method(generate,dist_exponential)
 S3method(generate,dist_f)
 S3method(generate,dist_gamma)
 S3method(generate,dist_geometric)
+S3method(generate,dist_gh)
 S3method(generate,dist_gk)
 S3method(generate,dist_gumbel)
 S3method(generate,dist_hypergeometric)
@@ -276,6 +280,7 @@ S3method(log_density,dist_exponential)
 S3method(log_density,dist_f)
 S3method(log_density,dist_gamma)
 S3method(log_density,dist_geometric)
+S3method(log_density,dist_gh)
 S3method(log_density,dist_gk)
 S3method(log_density,dist_gumbel)
 S3method(log_density,dist_hypergeometric)
@@ -365,6 +370,7 @@ S3method(quantile,dist_exponential)
 S3method(quantile,dist_f)
 S3method(quantile,dist_gamma)
 S3method(quantile,dist_geometric)
+S3method(quantile,dist_gh)
 S3method(quantile,dist_gk)
 S3method(quantile,dist_gumbel)
 S3method(quantile,dist_hypergeometric)
@@ -461,6 +467,7 @@ export(dist_exponential)
 export(dist_f)
 export(dist_gamma)
 export(dist_geometric)
+export(dist_gh)
 export(dist_gk)
 export(dist_gumbel)
 export(dist_hypergeometric)

NEWS.md

@@ -16,7 +16,8 @@
 
 ### Probability distributions
 
-* Added `dist_gk()` for the g-and-k distributions.
+* Added `dist_gk()` for g-and-k distributions.
+* Added `dist_gh()` for g-and-h distributions.
 
 ## Improvements
 

R/dist_gh.R

@@ -0,0 +1,117 @@
+#' The generalised g-and-h Distribution
+#'
+#' @description
+#' `r lifecycle::badge('stable')`
+#'
+#' The generalised g-and-h distribution is a flexible distribution used to model univariate data, similar to the g-k distribution.
+#' It is known for its ability to handle skewness and heavy-tailed behavior.
+#'
+#' @inheritParams gk::dgh
+#'
+#' @details
+#'
+#'   We recommend reading this documentation on
+#'   <https://pkg.mitchelloharawild.com/distributional/>, where the math
+#'   will render nicely.
+#'
+#'   In the following, let \eqn{X} be a g-and-h random variable with parameters
+#'   `A`, `B`, `g`, `h`, and `c`.
+#'
+#'   **Support**: \eqn{(-\infty, \infty)}
+#'
+#'   **Mean**: Not available in closed form.
+#'
+#'   **Variance**: Not available in closed form.
+#'
+#'   **Probability density function (p.d.f)**:
+#'
+#'   The g-and-h distribution does not have a closed-form expression for its density. Instead,
+#'   it is defined through its quantile function:
+#'
+#'   \deqn{
+#'     Q(u) = A + B \left( 1 + c \frac{1 - \exp(-gz(u))}{1 + \exp(-gz(u))} \right) \exp(h z(u)^2/2) z(u)
+#'   }{
+#'     Q(u) = A + B * (1 + c * ((1 - exp(-g * z(u))) / (1 + exp(-g * z(u))))) * exp(h * z(u)^2/2) * z(u)
+#'   }
+#'
+#'   where \eqn{z(u) = \Phi^{-1}(u)}
+#'
+#'   **Cumulative distribution function (c.d.f)**:
+#'
+#'   The cumulative distribution function is typically evaluated numerically due to the lack
+#'   of a closed-form expression.
+#'
+#' @seealso [gk::dgh], [distributional::dist_gk]
+#'
+#' @examples
+#' dist <- dist_gh(A = 0, B = 1, g = 0, h = 0.5)
+#' dist
+#'
+#' @examplesIf requireNamespace("gh", quietly = TRUE)
+#' mean(dist)
+#' variance(dist)
+#' support(dist)
+#' generate(dist, 10)
+#'
+#' density(dist, 2)
+#' density(dist, 2, log = TRUE)
+#'
+#' cdf(dist, 4)
+#'
+#' quantile(dist, 0.7)
+#'
+#' @name dist_gh
+#' @export
+dist_gh <- function(A, B, g, h, c = 0.8){
+  A <- vec_cast(A, double())
+  B <- vec_cast(B, double())
+  g <- vec_cast(g, double())
+  h <- vec_cast(h, double())
+  c <- vec_cast(c, double())
+  if(any(B <= 0)){
+    abort("The B parameter (scale) of the g-and-h distribution must be strictly positive.")
+  }
+  new_dist(A = A, B = B, g = g, h = h, c = c, class = "dist_gh")
+}
+
+#' @export
+format.dist_gh <- function(x, digits = 2, ...){
+  sprintf(
+    "gh(A = %s, B = %s, g = %s, h = %s%s)",
+    format(x[["A"]], digits = digits, ...),
+    format(x[["B"]], digits = digits, ...),
+    format(x[["g"]], digits = digits, ...),
+    format(x[["h"]], digits = digits, ...),
+    if (x[["c"]]==0.8) "" else paste0(", c = ", format(x[["c"]], digits = digits, ...))
+  )
+}
+
+#' @export
+density.dist_gh <- function(x, at, ...){
+  require_package("gh")
+  gk::dgh(at, x[["A"]], x[["B"]], x[["g"]], x[["h"]], x[["c"]])
+}
+
+#' @export
+log_density.dist_gh <- function(x, at, ...){
+  require_package("gh")
+  gk::dgh(at, x[["A"]], x[["B"]], x[["g"]], x[["h"]], x[["c"]], log = TRUE)
+}
+
+#' @export
+quantile.dist_gh <- function(x, p, ...){
+  require_package("gh")
+  gk::qgh(p, x[["A"]], x[["B"]], x[["g"]], x[["h"]], x[["c"]])
+}
+
+#' @export
+cdf.dist_gh <- function(x, q, ...){
+  require_package("gh")
+  gk::pgh(q, x[["A"]], x[["B"]], x[["g"]], x[["h"]], x[["c"]])
+}
+
+#' @export
+generate.dist_gh <- function(x, times, ...){
+  require_package("gh")
+  gk::rgh(times, x[["A"]], x[["B"]], x[["g"]], x[["h"]], x[["c"]])
+}

R/dist_gk.R

@@ -41,7 +41,7 @@
 #'   The cumulative distribution function is typically evaluated numerically due to the lack
 #'   of a closed-form expression.
 #'
-#' @seealso [gk::gk]
+#' @seealso [gk::dgk], [distributional::dist_gh]
 #'
 #' @examples
 #' dist <- dist_gk(A = 0, B = 1, g = 0, k = 0.5)

tests/testthat/test-dist-gh.R

@@ -0,0 +1,36 @@
+test_that("g-and-h distribution", {
+  # Define g-and-h distribution parameters
+  A <- 0
+  B <- 1
+  g <- 0
+  h <- 0.5
+  c <- 0.8
+  dist <- dist_gh(A, B, g, h, c)
+
+  # Check formatting
+  expect_equal(format(dist), "gh(A = 0, B = 1, g = 0, h = 0.5)")
+
+  # Require package installed
+  skip_if_not_installed("gh", "0.1.0")
+
+  # quantiles
+  expect_equal(quantile(dist, 0.1), gk::qgh(0.1, A, B, g, h, c))
+  expect_equal(quantile(dist, 0.5), gk::qgh(0.5, A, B, g, h, c))
+
+  # pdf
+  expect_equal(density(dist, 0), gk::dgh(0, A, B, g, h, c))
+  expect_equal(density(dist, 3), gk::dgh(3, A, B, g, h, c))
+
+  # cdf
+  expect_equal(cdf(dist, 0), gk::pgh(0, A, B, g, h, c))
+  expect_equal(cdf(dist, 3), gk::pgh(3, A, B, g, h, c))
+
+  # F(Finv(a)) ~= a
+  expect_equal(cdf(dist, quantile(dist, 0.4)), 0.4, tolerance = 1e-3)
+
+  # Generate random samples
+  set.seed(123)
+  samples <- generate(dist, 10)
+  set.seed(123)
+  expect_equal(samples[[1L]], gk::rgh(10, A, B, g, h, c))
+})