Merge pull request #55 from poissonconsulting/skew

Skew normal residual functions
poissonconsulting · Jul 27, 2023 · 01c266c · 01c266c
2 parents bcb821f + fbbb6be
commit 01c266c
Show file tree

Hide file tree

Showing 63 changed files with 619 additions and 10 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -43,6 +43,7 @@ export(dev_neg_binom)
 export(dev_norm)
 export(dev_pois)
 export(dev_pois_zi)
+export(dev_skewnorm)
 export(dev_student)
 export(dskewnorm)
 export(exp10)
@@ -69,6 +70,7 @@ export(log_lik_neg_binom)
 export(log_lik_norm)
 export(log_lik_pois)
 export(log_lik_pois_zi)
+export(log_lik_skewnorm)
 export(log_lik_student)
 export(log_odds)
 export(log_odds_ratio)
@@ -105,6 +107,7 @@ export(ran_neg_binom)
 export(ran_norm)
 export(ran_pois)
 export(ran_pois_zi)
+export(ran_skewnorm)
 export(ran_student)
 export(rbern)
 export(res_bern)
@@ -118,6 +121,7 @@ export(res_neg_binom)
 export(res_norm)
 export(res_pois)
 export(res_pois_zi)
+export(res_skewnorm)
 export(res_student)
 export(rskewnorm)
 export(sextreme)

diff --git a/R/dev.R b/R/dev.R
@@ -237,6 +237,37 @@ dev_pois_zi <- function(x, lambda, prob = 0, res = FALSE) {
   dev_res(x, lambda * (1 - prob), dev)
 }
 
+
+#' Skew Normal Deviances
+#'
+#' @inheritParams params
+#' @param x A numeric vector of values.
+#' @param shape A numeric vector of shape.
+#'
+#' @return An numeric vector of the corresponding deviances or deviance residuals.
+#' @family dev_dist
+#' @export
+#'
+#' @examples
+#' dev_skewnorm(c(-2:2))
+#' dev_skewnorm(-2:2, 0, 1, 5)
+#' dev_skewnorm(-2:2, 0, 1, 5, res = TRUE)
+dev_skewnorm <- function(x, mean = 0, sd = 1, shape = 0, res = FALSE) {
+  delta <- shape / sqrt(1 + shape^2)
+  mu_z <- sqrt(2 / pi) * delta
+  sig_z <- sqrt(1 - mu_z^2)
+  gam_1 <- ((4 - pi) / 2) * ((delta * sqrt(2 / pi))^3 / (1 - (2 * delta^2) / pi)^(3/2))
+  m_o <- mu_z - (gam_1 * sig_z / 2) - (sign(shape) / 2) * exp(-2 * pi / abs(shape))
+  mode_sat <- mean + sd * m_o
+  dev <- log_lik_skewnorm(mode_sat, mean = mean, sd = sd, shape = shape) -
+    log_lik_skewnorm(x, mean = mean, sd = sd, shape = shape)
+  neg <- dev < 0
+  dev[neg] <- 0
+  dev <- dev * 2
+  if(vld_false(res)) return(dev)
+  dev_res(x, mean + sd * (shape / sqrt(1 + shape^2)) * sqrt(2 / pi), dev)
+}
+
 #' Student's t Deviances
 #'
 #' @inheritParams params

diff --git a/R/log-lik.R b/R/log-lik.R
@@ -155,8 +155,8 @@ log_lik_neg_binom <- function(x, lambda = 1, theta = 0) {
 #' @export
 #'
 #' @examples
-#' dev_norm(c(-2:2))
-log_lik_norm <- function(x,  mean = 0, sd = 1) {
+#' log_lik_norm(c(-2:2))
+log_lik_norm <- function(x, mean = 0, sd = 1) {
   dnorm(x, mean = mean, sd = sd, log = TRUE)
 }
 
@@ -194,6 +194,29 @@ log_lik_pois_zi <- function(x, lambda = 1, prob = 0) {
   log(lpois)
 }
 
+#' Skew Normal Log-Likelihood
+#'
+#' @inheritParams params
+#' @param x A numeric vector of values.
+#' @param shape A numeric vector of shape.
+#'
+#' @return An numeric vector of the corresponding log-likelihoods.
+#' @family log_lik_dist
+#' @export
+#'
+#' @examples
+#' log_lik_skewnorm(c(-2:2))
+#' log_lik_skewnorm(c(-2:2), shape = -2)
+#' log_lik_skewnorm(c(-2:2), shape = 2)
+log_lik_skewnorm <- function(x, mean = 0, sd = 1, shape = 0) {
+  log_lik <- dskewnorm(x = x, mean = mean, sd = sd, shape = shape, log = TRUE)
+  use_norm <- !is.na(shape) & shape == 0
+  lnorm <- log_lik_norm(x = x, mean = mean, sd = sd)
+  lengths <- as.logical(length(x)) + as.logical(length(mean)) + as.logical(length(sd)) + as.logical(length(shape))
+  if (lengths >= 4)  log_lik[use_norm] <- lnorm[use_norm]
+  log_lik
+}
+
 #' Student's t Log-Likelihood
 #'
 #' @inheritParams params

diff --git a/R/ran.R b/R/ran.R
@@ -164,6 +164,24 @@ ran_pois_zi <- function(n = 1, lambda = 1, prob = 0) {
   stats::rpois(n, lambda = lambda) * ran_bern(n, prob = 1 - prob)
 }
 
+#' Skew Normal Random Samples
+#'
+#' @inheritParams params
+#' @param shape A numeric vector of shape.
+#' @return A numeric vector of the random samples.
+#' @family ran_dist
+#' @export
+#'
+#' @examples
+#' ran_skewnorm(10, shape = -1)
+#' ran_skewnorm(10, shape = 0)
+#' ran_skewnorm(10, shape = 1)
+ran_skewnorm <- function(n = 1, mean = 0, sd = 1, shape = 0) {
+  chk_whole_number(n)
+  chk_gte(n)
+  rskewnorm(n = n, mean = mean, sd = sd, shape = shape)
+}
+
 #' Student's t Random Samples
 #'
 #' @inheritParams params

diff --git a/R/res.R b/R/res.R
@@ -215,7 +215,7 @@ res_neg_binom <- function(x, lambda = 1, theta = 0, type = "dev", simulate = FAL
 #' @export
 #'
 #' @examples
-#' dev_norm(c(-2:2))
+#' res_norm(c(-2:2))
 res_norm <- function(x,  mean = 0, sd = 1, type = "dev", simulate = FALSE) {
   chk_string(type)
   if(!vld_false(simulate)) {
@@ -296,6 +296,32 @@ res_student_standardized <- function(x, mean, sd, theta) {
   return(res)
 }
 
+#' Skew Normal Residuals
+#'
+#' @inheritParams params
+#' @param x A numeric vector of values.
+#' @param shape A numeric vector of shape.
+#'
+#' @return An numeric vector of the corresponding residuals.
+#' @family res_dist
+#' @export
+#'
+#' @examples
+#' res_skewnorm(c(-2:2))
+res_skewnorm <- function(x, mean = 0, sd = 1, shape = 0, type = "dev", simulate = FALSE) {
+  chk_string(type)
+  if(!vld_false(simulate)) {
+    x <- ran_skewnorm(length(x), mean = mean, sd = sd, shape = shape)
+  }
+  switch(type,
+         data = x,
+         raw = x - mean + sd * (shape / sqrt(1 + shape^2)) * sqrt(2 / pi),
+         standardized = (x - (mean + sd * (shape / sqrt(1 + shape^2)) * sqrt(2 / pi))) /
+           (sd^2 * (1 - ((2 * (shape / sqrt(1 + shape^2))^2) / pi))),
+         dev = dev_skewnorm(x, mean = mean, sd = sd, shape = shape, res = TRUE),
+         chk_subset(x, c("data", "raw", "dev", "standardized")))
+}
+
 #' Student's t Residuals
 #'
 #' @inheritParams params

diff --git a/R/skewnorm.R b/R/skewnorm.R
@@ -27,6 +27,8 @@ dskewnorm <- function(x, mean = 0, sd = 1, shape = 0, log = FALSE) {
     )
   }
   chk_gte(sd)
+  nulls <- any(is.null(x), is.null(mean), is.null(sd), is.null(shape))
+  if (nulls) stop("invalid arguments")
   lengths <- as.logical(length(x)) + as.logical(length(mean)) + as.logical(length(sd)) + as.logical(length(shape))
   if (lengths >= 4) {
     nas <- any(is.na(x), is.na(mean), is.na(sd), is.na(shape))
@@ -54,6 +56,8 @@ pskewnorm <- function(q, mean = 0, sd = 1, shape = 0) {
     )
   }
   chk_gte(sd)
+  nulls <- any(is.null(q), is.null(mean), is.null(sd), is.null(shape))
+  if (nulls) stop("invalid arguments")
   lengths <- as.logical(length(q)) + as.logical(length(mean)) + as.logical(length(sd)) + as.logical(length(shape))
   if (lengths >= 4) {
     nas <- any(is.na(q), is.na(mean), is.na(sd), is.na(shape))
@@ -83,6 +87,8 @@ qskewnorm <- function(p, mean = 0, sd = 1, shape = 0) {
   chk_gte(sd)
   chk_gte(p)
   chk_lte(p, 1)
+  nulls <- any(is.null(p), is.null(mean), is.null(sd), is.null(shape))
+  if (nulls) stop("invalid arguments")
   lengths <- as.logical(length(p)) + as.logical(length(mean)) + as.logical(length(sd)) + as.logical(length(shape))
   if (lengths >= 4) {
     nas <- any(is.na(p), is.na(mean), is.na(sd), is.na(shape))
@@ -116,13 +122,15 @@ rskewnorm <- function(n = 1, mean = 0, sd = 1, shape = 0) {
   chk_lt(n, Inf)
   chk_not_any_na(n)
   chk_gte(sd)
+  nulls <- any(is.null(n), is.null(mean), is.null(sd), is.null(shape))
+  if (nulls) stop("invalid arguments")
   lengths <- as.logical(length(n)) + as.logical(length(mean)) + as.logical(length(sd)) + as.logical(length(shape))
   character <- any(is.character(n), is.character(mean), is.character(sd), is.character(shape))
   if (lengths < 4 & !character) {
     return(vector(mode = "numeric"))
   }
   chk_whole_number(n)
-  if (lengths >= 4) {
+  if (lengths >= 4 & n != 0L) {
     nas <- any(is.na(n), is.na(mean), is.na(sd), is.na(shape))
     if (!nas) {
       chk_compatible_lengths(rep(1, n), mean, sd, shape)
@@ -131,5 +139,6 @@ rskewnorm <- function(n = 1, mean = 0, sd = 1, shape = 0) {
   chk_false(character)
   ran <- sn::rsn(n, xi = mean, omega = sd, alpha = shape)
   attributes(ran) <- NULL
+  if (n == 0L) return(ran)
   ran[1:n]
 }
diff --git a/_pkgdown.yml b/_pkgdown.yml
@@ -97,6 +97,7 @@ reference:
   - '`dev_norm`'
   - '`dev_pois`'
   - '`dev_pois_zi`'
+  - '`dev_skewnorm`'
   - '`dev_student`'
 - title: Residual
   desc: Raw and Deviance Residuals
@@ -112,6 +113,7 @@ reference:
   - '`res_norm`'
   - '`res_pois`'
   - '`res_pois_zi`'
+  - '`res_skewnorm`'
   - '`res_student`'
 - title: Log-Likelihood
   desc: Log-likelihood functions
@@ -127,6 +129,7 @@ reference:
   - '`log_lik_norm`'
   - '`log_lik_pois`'
   - '`log_lik_pois_zi`'
+  - '`log_lik_skewnorm`'
   - '`log_lik_student`'
 - title: Random
   desc: Random sample functions
@@ -142,6 +145,7 @@ reference:
   - '`ran_norm`'
   - '`ran_pois`'
   - '`ran_pois_zi`'
+  - '`ran_skewnorm`'
   - '`ran_student`'
 - title: Bernoulli
   desc: Bernoulli distribution functions

diff --git a/man/dev_bern.Rd b/man/dev_bern.Rd
diff --git a/man/dev_beta_binom.Rd b/man/dev_beta_binom.Rd
diff --git a/man/dev_binom.Rd b/man/dev_binom.Rd
diff --git a/man/dev_gamma.Rd b/man/dev_gamma.Rd
diff --git a/man/dev_gamma_pois.Rd b/man/dev_gamma_pois.Rd
diff --git a/man/dev_lnorm.Rd b/man/dev_lnorm.Rd
diff --git a/man/dev_neg_binom.Rd b/man/dev_neg_binom.Rd
diff --git a/man/dev_norm.Rd b/man/dev_norm.Rd
diff --git a/man/dev_pois.Rd b/man/dev_pois.Rd
diff --git a/man/dev_pois_zi.Rd b/man/dev_pois_zi.Rd
diff --git a/man/dev_skewnorm.Rd b/man/dev_skewnorm.Rd
diff --git a/man/dev_student.Rd b/man/dev_student.Rd