#' Probability calculator for the normal distribution
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param mean Mean
#' @param stdev Standard deviation
#' @param lb Lower bound (default is -Inf)
#' @param ub Upper bound (default is Inf)
#' @param plb Lower probability bound
#' @param pub Upper probability bound
#' @param dec Number of decimals to show
#'
#' @examples
#' prob_norm(mean = 0, stdev = 1, ub = 0)
#'
#' @seealso \code{\link{summary.prob_norm}} to summarize results
#' @seealso \code{\link{plot.prob_norm}} to plot results
#'
#' @export
prob_norm <- function(mean, stdev, lb = NA, ub = NA,
                      plb = NA, pub = NA, dec = 3) {
  p_ub <- pnorm(ub, mean, stdev)
  p_lb <- pnorm(lb, mean, stdev)
  p_int <- max(p_ub - p_lb, 0) %>% round(dec)

  p_ub %<>% round(dec)
  p_lb %<>% round(dec)

  if (!is.na(pub)) {
    if (pub > 1) pub <- 1
    if (pub < 0) pub <- 0
  }

  if (!is.na(plb)) {
    if (plb > 1) plb <- 1
    if (plb < 0) plb <- 0
  }

  v_ub <- qnorm(pub, mean, stdev) %>% round(dec)
  v_lb <- qnorm(plb, mean, stdev) %>% round(dec)

  if (!is.na(lb) && !is.na(ub)) {
    if (lb > ub) {
      lb <- ub <- NA
      mess_values <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  if (!is.na(plb) && !is.na(pub)) {
    if (plb > pub) {
      plb <- pub <- NA
      mess_probs <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  as.list(environment()) %>% add_class("prob_norm")
}

#' Plot method for the probability calculator (normal)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param x Return value from \code{\link{prob_norm}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @examples
#' result <- prob_norm(mean = 0, stdev = 1, ub = 0)
#' plot(result)
#'
#' @seealso \code{\link{prob_norm}} to calculate results
#' @seealso \code{\link{summary.prob_norm}} to summarize results
#'
#' @importFrom rlang .data
#'
#' @export
plot.prob_norm <- function(x, type = "values", ...) {
  mess <- paste0("mess_", type)
  if (!is.null(x[[mess]])) {
    return(" ")
  }

  if (type == "values") {
    lb <- x$lb
    ub <- x$ub
  } else {
    lb <- x$v_lb
    ub <- x$v_ub
  }

  mean <- x$mean
  stdev <- x$stdev

  limits <- c(mean - 3 * stdev, mean + 3 * stdev)

  dnorm_limit <- function(x) {
    y <- dnorm(x, mean = mean, sd = stdev)
    y[x < lb | x > ub] <- 0
    y
  }

  dnorm_lb <- function(x) {
    if (is.na(lb)) {
      return(0)
    }
    y <- dnorm(x, mean = mean, sd = stdev)
    y[x > lb] <- 0
    y
  }

  dnorm_ub <- function(x) {
    if (is.na(ub)) {
      return(0)
    }
    y <- dnorm(x, mean = mean, sd = stdev)
    y[x < ub] <- 0
    y
  }

  dnorm_lines <- c(ub, lb) %>% na.omit()
  if (length(dnorm_lines) == 0) dnorm_lines <- c(-Inf, Inf)

  ## based on https://rstudio-pubs-static.s3.amazonaws.com/58753_13e35d9c089d4f55b176057235778679.html
  ## and R Graphics Cookbook
  plt <- ggplot(data.frame(x = limits), aes(x = .data$x)) +
    stat_function(fun = stats::dnorm, args = list(mean = mean, sd = stdev)) +
    stat_function(fun = dnorm_limit, geom = "area", fill = "blue", alpha = 0.5, n = 501) +
    stat_function(fun = dnorm_lb, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    stat_function(fun = dnorm_ub, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    geom_vline(xintercept = dnorm_lines, color = "black", linetype = "dashed", linewidth = .5) +
    labs(x = "", y = "")

  sshhr(plt)
}

#' Summary method for the probability calculator (normal)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param object Return value from \code{\link{prob_norm}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @examples
#' result <- prob_norm(mean = 0, stdev = 1, ub = 0)
#' summary(result)
#'
#' @seealso \code{\link{prob_norm}} to calculate results
#' @seealso \code{\link{plot.prob_norm}} to plot results
#'
#' @export
summary.prob_norm <- function(object, type = "values", ...) {
  cat("Probability calculator\n")
  cat("Distribution: Normal\n")

  mess <- object[[paste0("mess_", type)]]
  if (!is.null(mess)) {
    return(mess)
  }
  env <- environment()
  ret <- sapply(names(object), function(x) assign(x, object[[x]], envir = env))

  cat("Mean        :", round(mean, dec), "\n")
  cat("St. dev     :", round(stdev, dec), "\n")

  if (type == "values") {
    cat("Lower bound :", if (is.na(lb)) "-Inf" else lb, "\n")
    cat("Upper bound :", if (is.na(ub)) "Inf" else ub, "\n")

    if (!is.na(ub) || !is.na(lb)) {
      cat("\n")

      if (!is.na(lb)) {
        cat(paste0("P(X < ", lb, ") = ", p_lb, "\n"))
        cat(paste0("P(X > ", lb, ") = ", round(1 - p_lb, dec), "\n"))
      }

      if (!is.na(ub)) {
        cat(paste0("P(X < ", ub, ") = ", p_ub, "\n"))
        cat(paste0("P(X > ", ub, ") = ", round(1 - p_ub, dec), "\n"))
      }

      if (!is.na(lb) && !is.na(ub)) {
        cat(paste0("P(", lb, " < X < ", ub, ")     = ", p_int, "\n"))
        cat(paste0("1 - P(", lb, " < X < ", ub, ") = ", round(1 - p_int, dec), "\n"))
      }
    }
  } else {
    pub <- if (is.na(pub)) 2 else pub
    plb <- if (is.na(plb)) -1 else plb

    cat("Lower bound :", if (plb < 0) "0" else plb, "\n")
    cat("Upper bound :", if (pub > 1) "1" else pub, "\n")

    if (pub <= 1 || plb >= 0) {
      cat("\n")

      if (plb >= 0) {
        cat(paste0("P(X < ", v_lb, ") = ", plb, "\n"))
        cat(paste0("P(X > ", v_lb, ") = ", round(1 - plb, dec), "\n"))
      }

      if (pub <= 1) {
        cat(paste0("P(X < ", v_ub, ") = ", pub, "\n"))
        cat(paste0("P(X > ", v_ub, ") = ", round(1 - pub, dec), "\n"))
      }

      if (pub <= 1 && plb >= 0) {
        cat(paste0("P(", v_lb, " < X < ", v_ub, ")     = ", pub - plb, "\n"))
        cat(paste0("1 - P(", v_lb, " < X < ", v_ub, ") = ", round(1 - (pub - plb), dec), "\n"))
      }
    }
  }
}

#' Probability calculator for the log normal distribution
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param meanlog Mean of the distribution on the log scale
#' @param sdlog Standard deviation of the distribution on the log scale
#' @param lb Lower bound (default is -Inf)
#' @param ub Upper bound (default is Inf)
#' @param plb Lower probability bound
#' @param pub Upper probability bound
#' @param dec Number of decimals to show
#'
#' @seealso \code{\link{summary.prob_lnorm}} to summarize results
#' @seealso \code{\link{plot.prob_lnorm}} to plot results
#'
#' @examples
#' prob_lnorm(meanlog = 0, sdlog = 1, lb = 0, ub = 1)
#'
#' @export
prob_lnorm <- function(meanlog, sdlog, lb = NA, ub = NA,
                       plb = NA, pub = NA, dec = 3) {
  p_ub <- plnorm(ub, meanlog, sdlog)
  p_lb <- plnorm(lb, meanlog, sdlog)
  p_int <- max(p_ub - p_lb, 0) %>% round(dec)

  p_ub %<>% round(dec)
  p_lb %<>% round(dec)

  if (!is.na(pub)) {
    if (pub > 1) pub <- 1
    if (pub < 0) pub <- 0
  }

  if (!is.na(plb)) {
    if (plb > 1) plb <- 1
    if (plb < 0) plb <- 0
  }

  v_ub <- qlnorm(pub, meanlog, sdlog) %>% round(dec)
  v_lb <- qlnorm(plb, meanlog, sdlog) %>% round(dec)

  if (!is.na(lb) && !is.na(ub)) {
    if (lb > ub) {
      lb <- ub <- NA
      mess_values <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  if (!is.na(plb) && !is.na(pub)) {
    if (plb > pub) {
      plb <- pub <- NA
      mess_probs <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  as.list(environment()) %>% add_class("prob_lnorm")
}

#' Plot method for the probability calculator (log normal)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param x Return value from \code{\link{prob_norm}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_lnorm}} to calculate results
#' @seealso \code{\link{plot.prob_lnorm}} to plot results
#'
#' @examples
#' result <- prob_lnorm(meanlog = 0, sdlog = 1, lb = 0, ub = 1)
#' plot(result, type = "values")
#'
#' @importFrom rlang .data
#'
#' @export
plot.prob_lnorm <- function(x, type = "values", ...) {
  mess <- paste0("mess_", type)
  if (!is.null(x[[mess]])) {
    return(" ")
  }

  if (type == "values") {
    lb <- x$lb
    ub <- x$ub
  } else {
    lb <- x$v_lb
    ub <- x$v_ub
  }

  meanlog <- x$meanlog
  sdlog <- x$sdlog

  # limits <- c(meanlog - 3 * sdlog, meanlog + 3 * sdlog)
  limits <- c(0, meanlog + ub * sdlog)

  dlnorm_limit <- function(x) {
    y <- dlnorm(x, meanlog = meanlog, sdlog = sdlog)
    y[x < lb | x > ub] <- 0
    y
  }

  dlnorm_lb <- function(x) {
    if (is.na(lb)) {
      return(0)
    }
    y <- dlnorm(x, meanlog = meanlog, sdlog = sdlog)
    y[x > lb] <- 0
    y
  }

  dlnorm_ub <- function(x) {
    if (is.na(ub)) {
      return(0)
    }
    y <- dlnorm(x, meanlog = meanlog, sdlog = sdlog)
    y[x < ub] <- 0
    y
  }

  dlnorm_lines <- c(ub, lb) %>% na.omit()
  if (length(dlnorm_lines) == 0) dlnorm_lines <- c(-Inf, Inf)

  ## based on https://rstudio-pubs-static.s3.amazonaws.com/58753_13e35d9c089d4f55b176057235778679.html
  ## and R Graphics Cookbook
  plt <- ggplot(data.frame(x = limits), aes(x = .data$x)) +
    stat_function(fun = stats::dlnorm, args = list(meanlog = meanlog, sdlog = sdlog)) +
    stat_function(fun = dlnorm_limit, geom = "area", fill = "blue", alpha = 0.5, n = 501) +
    stat_function(fun = dlnorm_lb, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    stat_function(fun = dlnorm_ub, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    geom_vline(xintercept = dlnorm_lines, color = "black", linetype = "dashed", linewidth = .5) +
    labs(x = "", y = "")

  sshhr(plt)
}

#' Summary method for the probability calculator (log normal)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param object Return value from \code{\link{prob_norm}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_lnorm}} to calculate results
#' @seealso \code{\link{plot.prob_lnorm}} to summarize results
#'
#' @examples
#' result <- prob_lnorm(meanlog = 0, sdlog = 1, lb = 0, ub = 1)
#' summary(result, type = "values")
#'
#' @export
summary.prob_lnorm <- function(object, type = "values", ...) {
  cat("Probability calculator\n")
  cat("Distribution: Log normal\n")

  mess <- object[[paste0("mess_", type)]]
  if (!is.null(mess)) {
    return(mess)
  }
  env <- environment()
  ret <- sapply(names(object), function(x) assign(x, object[[x]], envir = env))

  cat("Mean log    :", round(meanlog, dec), "\n")
  cat("St. dev log :", round(sdlog, dec), "\n")

  if (type == "values") {
    cat("Lower bound :", if (is.na(lb)) "-Inf" else lb, "\n")
    cat("Upper bound :", if (is.na(ub)) "Inf" else ub, "\n")

    if (!is.na(ub) || !is.na(lb)) {
      cat("\n")

      if (!is.na(lb)) {
        cat(paste0("P(X < ", lb, ") = ", p_lb, "\n"))
        cat(paste0("P(X > ", lb, ") = ", round(1 - p_lb, dec), "\n"))
      }

      if (!is.na(ub)) {
        cat(paste0("P(X < ", ub, ") = ", p_ub, "\n"))
        cat(paste0("P(X > ", ub, ") = ", round(1 - p_ub, dec), "\n"))
      }

      if (!is.na(lb) && !is.na(ub)) {
        cat(paste0("P(", lb, " < X < ", ub, ")     = ", p_int, "\n"))
        cat(paste0("1 - P(", lb, " < X < ", ub, ") = ", round(1 - p_int, dec), "\n"))
      }
    }
  } else {
    pub <- if (is.na(pub)) 2 else pub
    plb <- if (is.na(plb)) -1 else plb

    cat(paste0("Lower bound : ", if (plb < 0) "0" else plb, " (", v_lb, ")\n"))
    cat(paste0("Upper bound : ", if (pub > 1) "1" else pub, " (", v_ub, ")\n"))

    if (pub <= 1 || plb >= 0) {
      cat("\n")

      if (plb >= 0) {
        cat(paste0("P(X < ", v_lb, ") = ", plb, "\n"))
        cat(paste0("P(X > ", v_lb, ") = ", round(1 - plb, dec), "\n"))
      }

      if (pub <= 1) {
        cat(paste0("P(X < ", v_ub, ") = ", pub, "\n"))
        cat(paste0("P(X > ", v_ub, ") = ", round(1 - pub, dec), "\n"))
      }

      if (pub <= 1 && plb >= 0) {
        cat(paste0("P(", v_lb, " < X < ", v_ub, ")     = ", pub - plb, "\n"))
        cat(paste0("1 - P(", v_lb, " < X < ", v_ub, ") = ", round(1 - (pub - plb), dec), "\n"))
      }
    }
  }
}

#' Probability calculator for the t-distribution
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param df Degrees of freedom
#' @param lb Lower bound (default is -Inf)
#' @param ub Upper bound (default is Inf)
#' @param plb Lower probability bound
#' @param pub Upper probability bound
#' @param dec Number of decimals to show
#'
#' @seealso \code{\link{summary.prob_tdist}} to summarize results
#' @seealso \code{\link{plot.prob_tdist}} to plot results
#'
#' @examples
#' prob_tdist(df = 10, ub = 2.228)
#'
#' @export
prob_tdist <- function(df, lb = NA, ub = NA,
                       plb = NA, pub = NA, dec = 3) {
  p_ub <- pt(ub, df)
  p_lb <- pt(lb, df)
  p_int <- max(p_ub - p_lb, 0)

  p_ub %<>% round(dec)
  p_lb %<>% round(dec)
  p_int %<>% round(dec)

  if (!is.na(pub)) {
    if (pub > 1) pub <- 1
    if (pub < 0) pub <- 0
  }

  if (!is.na(plb)) {
    if (plb > 1) plb <- 1
    if (plb < 0) plb <- 0
  }

  v_ub <- qt(pub, df)
  v_lb <- qt(plb, df)

  v_ub %<>% round(dec)
  v_lb %<>% round(dec)

  if (!is.na(lb) && !is.na(ub)) {
    if (lb > ub) {
      lb <- ub <- NA
      mess_values <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  if (!is.na(plb) && !is.na(pub)) {
    if (plb > pub) {
      plb <- pub <- NA
      mess_probs <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  as.list(environment()) %>% add_class("prob_tdist")
}

#' Plot method for the probability calculator (t-distribution)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param x Return value from \code{\link{prob_tdist}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_tdist}} to calculate results
#' @seealso \code{\link{summary.prob_tdist}} to summarize results
#'
#' @examples
#' result <- prob_tdist(df = 10, ub = 2.228)
#' plot(result, type = "values")
#'
#' @importFrom rlang .data
#'
#' @export
plot.prob_tdist <- function(x, type = "values", ...) {
  mess <- paste0("mess_", type)
  if (!is.null(x[[mess]])) {
    return(" ")
  }
  if (type == "values") {
    lb <- x$lb
    ub <- x$ub
  } else {
    lb <- x$v_lb
    ub <- x$v_ub
  }

  df <- x$df

  limits <- c(-3, 3)
  dt_limit <- function(x) {
    y <- dt(x, df = df)
    y[x < lb | x > ub] <- 0
    y
  }

  dt_lb <- function(x) {
    if (is.na(lb)) {
      return(0)
    }
    y <- dt(x, df = df)
    y[x > lb] <- 0
    y
  }

  dt_ub <- function(x) {
    if (is.na(ub)) {
      return(0)
    }
    y <- dt(x, df = df)
    y[x < ub] <- 0
    y
  }

  dt_lines <- c(ub, lb) %>% na.omit()
  if (length(dt_lines) == 0) dt_lines <- c(-Inf, Inf)

  ## based on https://rstudio-pubs-static.s3.amazonaws.com/58753_13e35d9c089d4f55b176057235778679.html
  ## and R Graphics Cookbook
  plt <- ggplot(data.frame(x = limits), aes(x = .data$x)) +
    stat_function(fun = stats::dt, args = list(df = df)) +
    stat_function(fun = dt_limit, geom = "area", fill = "blue", alpha = 0.5, n = 501) +
    stat_function(fun = dt_lb, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    stat_function(fun = dt_ub, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    geom_vline(xintercept = dt_lines, color = "black", linetype = "dashed", linewidth = .5) +
    labs(x = "", y = "")

  sshhr(plt)
}

#' Summary method for the probability calculator (t-distribution)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param object Return value from \code{\link{prob_tdist}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_tdist}} to calculate results
#' @seealso \code{\link{plot.prob_tdist}} to plot results
#'
#' @examples
#' result <- prob_tdist(df = 10, ub = 2.228)
#' summary(result, type = "values")
#'
#' @export
summary.prob_tdist <- function(object, type = "values", ...) {
  cat("Probability calculator\n")
  cat("Distribution: t\n")

  mess <- object[[paste0("mess_", type)]]
  if (!is.null(mess)) {
    return(mess)
  }

  env <- environment()
  ret <- sapply(names(object), function(x) assign(x, object[[x]], envir = env))
  n <- df + 1
  cat("Df          :", df, "\n")
  cat("Mean        :", 0, "\n")
  cat("St. dev     :", ifelse(n > 2, round(n / (n - 2), dec), "NA"), "\n")

  if (type == "values") {
    cat("Lower bound :", if (is.na(lb)) "-Inf" else lb, "\n")
    cat("Upper bound :", if (is.na(ub)) "Inf" else ub, "\n")

    if (!is.na(ub) || !is.na(lb)) {
      cat("\n")

      if (!is.na(lb)) {
        cat(paste0("P(X < ", lb, ") = ", p_lb, "\n"))
        cat(paste0("P(X > ", lb, ") = ", round(1 - p_lb, dec), "\n"))
      }

      if (!is.na(ub)) {
        cat(paste0("P(X < ", ub, ") = ", p_ub, "\n"))
        cat(paste0("P(X > ", ub, ") = ", round(1 - p_ub, dec), "\n"))
      }

      if (!is.na(lb) && !is.na(ub)) {
        cat(paste0("P(", lb, " < X < ", ub, ")     = ", p_int, "\n"))
        cat(paste0("1 - P(", lb, " < X < ", ub, ") = ", round(1 - p_int, dec), "\n"))
      }
    }
  } else {
    pub <- if (is.na(pub)) 2 else pub
    plb <- if (is.na(plb)) -1 else plb

    cat("Lower bound :", if (plb < 0) "0" else plb, "\n")
    cat("Upper bound :", if (pub > 1) "1" else pub, "\n")

    if (pub <= 1 || plb >= 0) {
      cat("\n")

      if (plb >= 0) {
        cat(paste0("P(X < ", v_lb, ") = ", plb, "\n"))
        cat(paste0("P(X > ", v_lb, ") = ", round(1 - plb, dec), "\n"))
      }

      if (pub <= 1) {
        cat(paste0("P(X < ", v_ub, ") = ", pub, "\n"))
        cat(paste0("P(X > ", v_ub, ") = ", round(1 - pub, dec), "\n"))
      }

      if (pub <= 1 && plb >= 0) {
        cat(paste0("P(", v_lb, " < X < ", v_ub, ")     = ", pub - plb, "\n"))
        cat(paste0("1 - P(", v_lb, " < X < ", v_ub, ") = ", round(1 - (pub - plb), dec), "\n"))
      }
    }
  }
}

#' Probability calculator for the F-distribution
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param df1 Degrees of freedom
#' @param df2 Degrees of freedom
#' @param lb Lower bound (default is 0)
#' @param ub Upper bound (default is Inf)
#' @param plb Lower probability bound
#' @param pub Upper probability bound
#' @param dec Number of decimals to show
#'
#' @seealso \code{\link{summary.prob_fdist}} to summarize results
#' @seealso \code{\link{plot.prob_fdist}} to plot results
#'
#' @examples
#' prob_fdist(df1 = 10, df2 = 10, ub = 2.978)
#'
#' @export
prob_fdist <- function(df1, df2, lb = NA, ub = NA,
                       plb = NA, pub = NA, dec = 3) {
  if (!is_not(lb) && lb < 0) lb <- 0
  if (!is_not(ub) && ub < 0) ub <- 0

  p_ub <- pf(ub, df1, df2)
  p_lb <- pf(lb, df1, df2)
  p_int <- max(p_ub - p_lb, 0)

  p_ub %<>% round(dec)
  p_lb %<>% round(dec)
  p_int %<>% round(dec)

  if (!is.na(pub)) {
    if (pub > 1) pub <- 1
    if (pub < 0) pub <- 0
  }

  if (!is.na(plb)) {
    if (plb > 1) plb <- 1
    if (plb < 0) plb <- 0
  }

  v_ub <- qf(pub, df1, df2)
  v_lb <- qf(plb, df1, df2)

  v_ub %<>% round(dec)
  v_lb %<>% round(dec)

  if (!is.na(lb) && !is.na(ub)) {
    if (lb > ub) {
      lb <- ub <- NA
      mess_values <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  if (!is.na(plb) && !is.na(pub)) {
    if (plb > pub) {
      plb <- pub <- NA
      mess_probs <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  as.list(environment()) %>% add_class("prob_fdist")
}

#' Plot method for the probability calculator (F-distribution)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param x Return value from \code{\link{prob_fdist}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_fdist}} to calculate results
#' @seealso \code{\link{summary.prob_fdist}} to summarize results
#'
#' @examples
#' result <- prob_fdist(df1 = 10, df2 = 10, ub = 2.978)
#' plot(result, type = "values")
#'
#' @importFrom rlang .data
#'
#' @export
plot.prob_fdist <- function(x, type = "values", ...) {
  mess <- paste0("mess_", type)
  if (!is.null(x[[mess]])) {
    return(" ")
  }
  if (type == "values") {
    lb <- x$lb
    ub <- x$ub
  } else {
    lb <- x$v_lb
    ub <- x$v_ub
  }

  df1 <- x$df1
  df2 <- x$df2

  limits <- c(
    floor(qf(0.01, df1 = df1, df2 = df2)),
    ceiling(qf(1 - 0.01, df1 = df1, df2 = df2))
  )

  dat <- data.frame(
    x = limits,
    Probability = df(limits, df1 = df1, df2 = df2),
    df1 = df1,
    df2 = df2,
    stringsAsFactors = FALSE
  )

  df_line <- function(x) df(x, df1 = df1, df2 = df2)

  df_limit <- function(x) {
    y <- df(x, df1 = df1, df2 = df2)
    y[x < lb | x > ub] <- 0
    y
  }

  df_lb <- function(x) {
    if (is.na(lb)) {
      return(0)
    }
    y <- df(x, df1 = df1, df2 = df2)
    y[x > lb] <- 0
    y
  }

  df_ub <- function(x) {
    if (is.na(ub)) {
      return(0)
    }
    y <- df(x, df1 = df1, df2 = df2)
    y[x < ub] <- 0
    y
  }

  vlines <- c(ub, lb) %>% na.omit()
  if (length(vlines) == 0) vlines <- c(-Inf, Inf)

  ## based on https://rstudio-pubs-static.s3.amazonaws.com/58753_13e35d9c089d4f55b176057235778679.html
  ## and R Graphics Cookbook
  plt <- ggplot(dat, aes(x = .data$x)) +
    stat_function(fun = df_line, geom = "line") +
    stat_function(fun = df_limit, geom = "area", fill = "blue", alpha = 0.5, n = 501) +
    stat_function(fun = df_lb, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    stat_function(fun = df_ub, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    geom_vline(xintercept = vlines, color = "black", linetype = "dashed", linewidth = 0.5) +
    labs(x = "", y = "")

  sshhr(plt)
}

#' Summary method for the probability calculator (F-distribution)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param object Return value from \code{\link{prob_fdist}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_fdist}} to calculate results
#' @seealso \code{\link{plot.prob_fdist}} to plot results
#'
#' @examples
#' result <- prob_fdist(df1 = 10, df2 = 10, ub = 2.978)
#' summary(result, type = "values")
#'
#' @export
summary.prob_fdist <- function(object, type = "values", ...) {
  cat("Probability calculator\n")
  cat("Distribution: F\n")

  mess <- object[[paste0("mess_", type)]]
  if (!is.null(mess)) {
    return(mess)
  }
  env <- environment()
  ret <- sapply(names(object), function(x) assign(x, object[[x]], envir = env))

  cat("Df 1        :", df1, "\n")
  cat("Df 2        :", df2, "\n")
  m <- if (df2 > 2) round(df2 / (df2 - 2), dec) else "NA"
  variance <- if (df2 > 4) {
    round((2 * df2^2 * (df1 + df2 - 2)) / (df1 * (df2 - 2)^2 * (df2 - 4)), dec)
  } else {
    "NA"
  }
  cat("Mean        :", m, "\n")
  cat("Variance    :", variance, "\n")

  if (type == "values") {
    cat("Lower bound :", if (is.na(lb)) "0" else lb, "\n")
    cat("Upper bound :", if (is.na(ub)) "Inf" else ub, "\n")

    if (!is.na(ub) || !is.na(lb)) {
      cat("\n")

      if (!is.na(lb)) {
        cat(paste0("P(X < ", lb, ") = ", p_lb, "\n"))
        cat(paste0("P(X > ", lb, ") = ", round(1 - p_lb, dec), "\n"))
      }

      if (!is.na(ub)) {
        cat(paste0("P(X < ", ub, ") = ", p_ub, "\n"))
        cat(paste0("P(X > ", ub, ") = ", round(1 - p_ub, dec), "\n"))
      }

      if (!is.na(lb) && !is.na(ub)) {
        cat(paste0("P(", lb, " < X < ", ub, ")     = ", p_int, "\n"))
        cat(paste0("1 - P(", lb, " < X < ", ub, ") = ", round(1 - p_int, dec), "\n"))
      }
    }
  } else {
    pub <- if (is.na(pub)) 2 else pub
    plb <- if (is.na(plb)) -1 else plb

    cat("Lower bound :", if (plb < 0) "0" else plb, "\n")
    cat("Upper bound :", if (pub > 1) "1" else pub, "\n")

    if (pub <= 1 || plb >= 0) {
      cat("\n")

      if (plb >= 0) {
        cat(paste0("P(X < ", v_lb, ") = ", plb, "\n"))
        cat(paste0("P(X > ", v_lb, ") = ", round(1 - plb, dec), "\n"))
      }

      if (pub <= 1) {
        cat(paste0("P(X < ", v_ub, ") = ", pub, "\n"))
        cat(paste0("P(X > ", v_ub, ") = ", round(1 - pub, dec), "\n"))
      }

      if (pub <= 1 && plb >= 0) {
        cat(paste0("P(", v_lb, " < X < ", v_ub, ")     = ", pub - plb, "\n"))
        cat(paste0("1 - P(", v_lb, " < X < ", v_ub, ") = ", round(1 - (pub - plb), dec), "\n"))
      }
    }
  }
}

#' Probability calculator for the chi-squared distribution
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param df Degrees of freedom
#' @param lb Lower bound (default is 0)
#' @param ub Upper bound (default is Inf)
#' @param plb Lower probability bound
#' @param pub Upper probability bound
#' @param dec Number of decimals to show
#'
#' @seealso \code{\link{summary.prob_chisq}} to summarize results
#' @seealso \code{\link{plot.prob_chisq}} to plot results
#'
#' @examples
#' prob_chisq(df = 1, ub = 3.841)
#'
#' @export
prob_chisq <- function(df, lb = NA, ub = NA, plb = NA,
                       pub = NA, dec = 3) {
  if (!is_not(lb) && lb < 0) lb <- 0
  if (!is_not(ub) && ub < 0) ub <- 0

  p_ub <- pchisq(ub, df)
  p_lb <- pchisq(lb, df)
  p_int <- max(p_ub - p_lb, 0)

  p_ub %<>% round(dec)
  p_lb %<>% round(dec)
  p_int %<>% round(dec)

  if (!is.na(pub)) {
    if (pub > 1) pub <- 1
    if (pub < 0) pub <- 0
  }

  if (!is.na(plb)) {
    if (plb > 1) plb <- 1
    if (plb < 0) plb <- 0
  }

  v_ub <- qchisq(pub, df)
  v_lb <- qchisq(plb, df)

  v_ub %<>% round(dec)
  v_lb %<>% round(dec)

  if (!is.na(lb) && !is.na(ub)) {
    if (lb > ub) {
      lb <- ub <- NA
      mess_values <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  if (!is.na(plb) && !is.na(pub)) {
    if (plb > pub) {
      plb <- pub <- NA
      mess_probs <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  as.list(environment()) %>% add_class("prob_chisq")
}

#' Plot method for the probability calculator (Chi-squared distribution)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param x Return value from \code{\link{prob_chisq}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_chisq}} to calculate results
#' @seealso \code{\link{summary.prob_chisq}} to summarize results
#'
#' @examples
#' result <- prob_chisq(df = 1, ub = 3.841)
#' plot(result, type = "values")
#'
#' @importFrom rlang .data
#'
#' @export
plot.prob_chisq <- function(x, type = "values", ...) {
  mess <- paste0("mess_", type)
  if (!is.null(x[[mess]])) {
    return(" ")
  }

  if (type == "values") {
    lb <- x$lb
    ub <- x$ub
  } else {
    lb <- x$v_lb
    ub <- x$v_ub
  }

  df <- x$df

  limits <- c(
    floor(qchisq(0.001, df = df)),
    ceiling(qchisq(1 - 0.001, df = df))
  )

  dat <- data.frame(
    x = limits,
    Probability = dchisq(limits, df = df),
    df = df,
    stringsAsFactors = FALSE
  )

  dchisq_limit <- function(x) {
    y <- dchisq(x, df = df)
    y[x < lb | x > ub] <- 0
    y
  }

  dchisq_lb <- function(x) {
    if (is.na(lb)) {
      return(0)
    }
    y <- dchisq(x, df = df)
    y[x > lb] <- 0
    y
  }

  dchisq_ub <- function(x) {
    if (is.na(ub)) {
      return(0)
    }
    y <- dchisq(x, df = df)
    y[x < ub] <- 0
    y
  }

  vlines <- c(ub, lb) %>% na.omit()
  if (length(vlines) == 0) vlines <- c(-Inf, Inf)

  ## based on https://rstudio-pubs-static.s3.amazonaws.com/58753_13e35d9c089d4f55b176057235778679.html
  ## and R Graphics Cookbook
  plt <- ggplot(dat, aes(x = .data$x)) +
    stat_function(fun = stats::dchisq, args = list(df = df)) +
    stat_function(fun = dchisq_limit, geom = "area", fill = "blue", alpha = 0.5, n = 501) +
    stat_function(fun = dchisq_lb, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    stat_function(fun = dchisq_ub, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    geom_vline(xintercept = vlines, color = "black", linetype = "dashed", linewidth = 0.5) +
    labs(x = "", y = "")

  sshhr(plt)
}

#' Summary method for the probability calculator (Chi-squared distribution)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param object Return value from \code{\link{prob_chisq}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_chisq}} to calculate results
#' @seealso \code{\link{plot.prob_chisq}} to plot results
#'
#' @examples
#' result <- prob_chisq(df = 1, ub = 3.841)
#' summary(result, type = "values")
#'
#' @export
summary.prob_chisq <- function(object, type = "values", ...) {
  cat("Probability calculator\n")
  cat("Distribution: Chi-squared\n")

  mess <- object[[paste0("mess_", type)]]
  if (!is.null(mess)) {
    return(mess)
  }
  env <- environment()
  ret <- sapply(names(object), function(x) assign(x, object[[x]], envir = env))

  cat("Df          :", df, "\n")
  cat("Mean        :", df, "\n")
  cat("Variance    :", 2 * df, "\n")

  if (type == "values") {
    cat("Lower bound :", if (is.na(lb)) "0" else lb, "\n")
    cat("Upper bound :", if (is.na(ub)) "Inf" else ub, "\n")

    if (!is.na(ub) || !is.na(lb)) {
      cat("\n")

      if (!is.na(lb)) {
        cat(paste0("P(X < ", lb, ") = ", p_lb, "\n"))
        cat(paste0("P(X > ", lb, ") = ", round(1 - p_lb, dec), "\n"))
      }

      if (!is.na(ub)) {
        cat(paste0("P(X < ", ub, ") = ", p_ub, "\n"))
        cat(paste0("P(X > ", ub, ") = ", round(1 - p_ub, dec), "\n"))
      }

      if (!is.na(lb) && !is.na(ub)) {
        cat(paste0("P(", lb, " < X < ", ub, ")     = ", p_int, "\n"))
        cat(paste0("1 - P(", lb, " < X < ", ub, ") = ", round(1 - p_int, dec), "\n"))
      }
    }
  } else {
    pub <- if (is.na(pub)) 2 else pub
    plb <- if (is.na(plb)) -1 else plb

    cat("Lower bound :", if (plb < 0) "0" else plb, "\n")
    cat("Upper bound :", if (pub > 1) "1" else pub, "\n")

    if (pub <= 1 || plb >= 0) {
      cat("\n")

      if (plb >= 0) {
        cat(paste0("P(X < ", v_lb, ") = ", plb, "\n"))
        cat(paste0("P(X > ", v_lb, ") = ", round(1 - plb, dec), "\n"))
      }

      if (pub <= 1) {
        cat(paste0("P(X < ", v_ub, ") = ", pub, "\n"))
        cat(paste0("P(X > ", v_ub, ") = ", round(1 - pub, dec), "\n"))
      }

      if (pub <= 1 && plb >= 0) {
        cat(paste0("P(", v_lb, " < X < ", v_ub, ")     = ", pub - plb, "\n"))
        cat(paste0("1 - P(", v_lb, " < X < ", v_ub, ") = ", round(1 - (pub - plb), dec), "\n"))
      }
    }
  }
}

#' Probability calculator for the uniform distribution
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param min Minimum value
#' @param max Maximum value
#' @param lb Lower bound (default = 0)
#' @param ub Upper bound (default = 1)
#' @param plb Lower probability bound
#' @param pub Upper probability bound
#' @param dec Number of decimals to show
#'
#' @seealso \code{\link{summary.prob_unif}} to summarize results
#' @seealso \code{\link{plot.prob_unif}} to plot results
#'
#' @examples
#' prob_unif(min = 0, max = 1, ub = 0.3)
#'
#' @export
prob_unif <- function(min, max, lb = NA, ub = NA,
                      plb = NA, pub = NA, dec = 3) {
  if (min > max) {
    mess_values <- "\nThe maximum value must be larger than the minimum value"
    mess_probs <- "\nThe maximum value must be larger than the minimum value"
  }

  if (!is.na(lb) && !is.na(ub)) {
    if (lb > ub) {
      lb <- ub <- NA
      mess_values <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  if (!is.na(plb) && !is.na(pub)) {
    if (plb > pub) {
      plb <- pub <- NA
      mess_probs <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  p_ub <- punif(ub, min, max)
  p_lb <- punif(lb, min, max)
  p_int <- max(p_ub - p_lb, 0) %>% round(dec)

  p_ub %<>% round(dec)
  p_lb %<>% round(dec)

  if (!is.na(pub)) {
    if (pub > 1) pub <- 1
    if (pub < 0) pub <- 0
  }

  if (!is.na(plb)) {
    if (plb > 1) plb <- 1
    if (plb < 0) plb <- 0
  }

  v_ub <- qunif(pub, min, max) %>% round(dec)
  v_lb <- qunif(plb, min, max) %>% round(dec)

  mean <- (max + min) / 2
  stdev <- sqrt((max - min)^2 / 12)

  as.list(environment()) %>% add_class("prob_unif")
}

#' Plot method for the probability calculator (uniform)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param x Return value from \code{\link{prob_unif}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_unif}} to calculate results
#' @seealso \code{\link{summary.prob_unif}} to summarize results
#'
#' @examples
#' result <- prob_unif(min = 0, max = 1, ub = 0.3)
#' plot(result, type = "values")
#'
#' @importFrom rlang .data
#'
#' @export
plot.prob_unif <- function(x, type = "values", ...) {
  mess <- paste0("mess_", type)
  if (!is.null(x[[mess]])) {
    return(" ")
  }

  if (type == "values") {
    lb <- x$lb
    ub <- x$ub
  } else {
    lb <- x$v_lb
    ub <- x$v_ub
  }

  min <- x$min
  max <- x$max

  if (min > max) {
    return(" ")
  }

  limits <- c(min, max)
  dunif_limit <- function(x) {
    y <- dunif(x, min = min, max = max)
    y[x < lb | x > ub] <- 0
    y
  }

  dunif_lb <- function(x) {
    if (is.na(lb)) {
      return(0)
    }
    y <- dunif(x, min = min, max = max)
    y[x > lb] <- 0
    y
  }

  dunif_ub <- function(x) {
    if (is.na(ub)) {
      return(0)
    }
    y <- dunif(x, min = min, max = max)
    y[x < ub] <- 0
    y
  }

  dunif_lines <- c(ub, lb) %>%
    na.omit() %>%
    base::setdiff(c(min, max))
  if (length(dunif_lines) == 0) dunif_lines <- c(-Inf, Inf)

  ## based on https://rstudio-pubs-static.s3.amazonaws.com/58753_13e35d9c089d4f55b176057235778679.html
  ## and R Graphics Cookbook
  plt <- data.frame(x = limits, y = dunif(limits, limits[1], limits[2]), lb = lb, ub = ub) %>%
    ggplot(aes(x = .data$x)) +
    stat_function(fun = dunif_limit, geom = "area", fill = "blue", alpha = 0.5, n = 501) +
    stat_function(fun = dunif_lb, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    stat_function(fun = dunif_ub, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    geom_vline(xintercept = dunif_lines, color = "black", linetype = "dashed", linewidth = 0.5) +
    geom_segment(aes(x = x[1], y = 0, xend = x[1], yend = y[1])) +
    geom_segment(aes(x = x[2], y = 0, xend = x[2], yend = y[2])) +
    geom_segment(aes(x = x[1], y = y[1], xend = x[2], yend = y[2])) +
    labs(x = "", y = "")

  sshhr(plt)
}

#' Summary method for the probability calculator (uniform)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param object Return value from \code{\link{prob_unif}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_unif}} to calculate results
#' @seealso \code{\link{plot.prob_unif}} to plot results
#'
#' @examples
#' result <- prob_unif(min = 0, max = 1, ub = 0.3)
#' summary(result, type = "values")
#'
#' @export
summary.prob_unif <- function(object, type = "values", ...) {
  cat("Probability calculator\n")
  cat("Distribution: Uniform\n")

  mess <- object[[paste0("mess_", type)]]
  if (!is.null(mess)) {
    return(mess)
  }
  env <- environment()
  ret <- sapply(names(object), function(x) assign(x, object[[x]], envir = env))

  cat("Min         :", min, "\n")
  cat("Max         :", max, "\n")
  if (max > min) {
    cat("Mean        :", round(mean, dec), "\n")
    cat("St. dev     :", round(stdev, dec), "\n")
  }

  if (type == "values") {
    cat("Lower bound :", ifelse(is.na(lb), min, lb), "\n")
    cat("Upper bound :", ifelse(is.na(ub), max, ub), "\n")

    if (!is.na(ub) || !is.na(lb)) {
      cat("\n")

      if (!is.na(lb)) {
        cat(paste0("P(X < ", lb, ") = ", p_lb, "\n"))
        cat(paste0("P(X > ", lb, ") = ", round(1 - p_lb, dec), "\n"))
      }

      if (!is.na(ub)) {
        cat(paste0("P(X < ", ub, ") = ", p_ub, "\n"))
        cat(paste0("P(X > ", ub, ") = ", round(1 - p_ub, dec), "\n"))
      }

      if (!is.na(lb) && !is.na(ub)) {
        cat(paste0("P(", lb, " < X < ", ub, ")     = ", p_int, "\n"))
        cat(paste0("1 - P(", lb, " < X < ", ub, ") = ", round(1 - p_int, dec), "\n"))
      }
    }
  } else {
    pub <- if (is.na(pub)) 2 else pub
    plb <- if (is.na(plb)) -1 else plb

    cat("Lower bound :", if (plb < 0) "0" else plb, "\n")
    cat("Upper bound :", if (pub > 1) "1" else pub, "\n")

    if (pub <= 1 || plb >= 0) {
      cat("\n")

      if (plb >= 0) {
        cat(paste0("P(X < ", v_lb, ") = ", plb, "\n"))
        cat(paste0("P(X > ", v_lb, ") = ", round(1 - plb, dec), "\n"))
      }

      if (pub <= 1) {
        cat(paste0("P(X < ", v_ub, ") = ", pub, "\n"))
        cat(paste0("P(X > ", v_ub, ") = ", round(1 - pub, dec), "\n"))
      }

      if (pub <= 1 && plb >= 0) {
        cat(paste0("P(", v_lb, " < X < ", v_ub, ")     = ", pub - plb, "\n"))
        cat(paste0("1 - P(", v_lb, " < X < ", v_ub, ") = ", round(1 - (pub - plb), dec), "\n"))
      }
    }
  }
}

#' Probability calculator for the binomial distribution
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param n Number of trials
#' @param p Probability
#' @param lb Lower bound on the number of successes
#' @param ub Upper bound on the number of successes
#' @param plb Lower probability bound
#' @param pub Upper probability bound
#' @param dec Number of decimals to show
#'
#' @seealso \code{\link{summary.prob_binom}} to summarize results
#' @seealso \code{\link{plot.prob_binom}} to plot results
#'
#' @examples
#' prob_binom(n = 10, p = 0.3, ub = 3)
#'
#' @export
prob_binom <- function(n, p, lb = NA, ub = NA,
                       plb = NA, pub = NA, dec = 3) {

  ## making sure n is integer
  n <- as_integer(n)

  if (!is_not(lb) && lb < 0) lb <- 0
  if (!is_not(ub) && ub < 0) ub <- 0

  if (is.na(lb) || lb < 0) {
    p_elb <- p_lb <- lb <- NA
  } else {
    lb <- as_integer(lb)
    if (lb > n) lb <- n
    p_elb <- dbinom(lb, n, p) %>% round(dec)
    p_lelb <- pbinom(lb, n, p) %>% round(dec)
    if (lb > 0) {
      p_lb <- sum(dbinom(0:max((lb - 1), 0), n, p)) %>% round(dec)
    } else {
      p_lb <- 0
    }
  }

  if (is.na(ub) || ub < 0) {
    p_eub <- p_ub <- ub <- NA
  } else {
    ub <- as_integer(ub)
    if (ub > n) ub <- n
    p_eub <- dbinom(ub, n, p) %>% round(dec)
    p_leub <- pbinom(ub, n, p) %>% round(dec)
    if (ub > 0) {
      p_ub <- sum(dbinom(0:max((ub - 1), 0), n, p)) %>% round(dec)
    } else {
      p_ub <- 0
    }
  }

  if (!is.na(ub) && !is.na(lb)) {
    p_int <- sum(dbinom(lb:ub, n, p)) %>%
      max(0) %>%
      round(dec)
  } else {
    p_int <- NA
  }

  if (is.na(plb)) {
    vlb <- NA
  } else {
    if (plb > 1) plb <- 1
    if (plb < 0) plb <- 0
    vlb <- qbinom(plb, n, p)

    vp_elb <- dbinom(vlb, n, p) %>% round(dec)
    vp_lelb <- pbinom(vlb, n, p) %>% round(dec)
    if (vlb > 0) {
      vp_lb <- sum(dbinom(0:max((vlb - 1), 0), n, p)) %>% round(dec)
    } else {
      vp_lb <- 0
    }
  }

  if (is.na(pub)) {
    vub <- NA
  } else {
    if (pub > 1) pub <- 1
    if (pub < 0) pub <- 0
    vub <- qbinom(pub, n, p)

    vp_eub <- dbinom(vub, n, p) %>% round(dec)
    vp_leub <- pbinom(vub, n, p) %>% round(dec)
    if (vub > 0) {
      vp_ub <- sum(dbinom(0:max((vub - 1), 0), n, p)) %>% round(dec)
    } else {
      vp_ub <- 0
    }
  }

  if (!is.na(pub) && !is.na(plb)) {
    vp_int <- sum(dbinom(vlb:vub, n, p)) %>%
      max(0) %>%
      round(dec)
  } else {
    vp_int <- NA
  }

  if (!is.na(lb) && !is.na(ub)) {
    if (lb > ub) {
      lb <- ub <- NA
      mess_values <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  if (!is.na(vlb) && !is.na(vub)) {
    if (vlb > vub) {
      plb <- pub <- vlb <- vub <- NA
      mess_probs <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  as.list(environment()) %>% add_class("prob_binom")
}


make_colors_discrete <- function(ub, lb, x_range) {
  colors <- factor(rep("blue", length(x_range)), levels = c("red", "green", "blue"))
  if (!is.na(lb) & !is.na(ub)) {
    colors[x_range == lb | x_range == ub] <- "green"
    colors[x_range > lb & x_range < ub] <- "blue"
    colors[x_range > ub | x_range < lb] <- "red"
  } else if (!is.na(lb)) {
    if (lb %in% x_range) colors[x_range == lb] <- "green"
    colors[x_range > lb] <- "blue"
    colors[x_range < lb] <- "red"
  } else if (!is.na(ub)) {
    if (ub %in% x_range) colors[x_range == ub] <- "green"
    colors[x_range > ub] <- "red"
    colors[x_range < ub] <- "blue"
  } else {
    colors[1:length(colors)] <- "blue"
  }
  return(colors)
}

make_bar_plot <- function(ub, lb, x_range, y_range) {
  colors <- make_colors_discrete(ub, lb, x_range)
  dat <- data.frame(x_range = x_range, y_range = y_range, colors = colors)

  if (nrow(dat) < 40) {
    # makes sure each bar has a label
    dat <- dat %>% mutate(x_range = factor(x_range))
  }

  cols <- c(red = "red", green = "green", blue = "blue")
  plt <- ggplot(dat, aes(x = .data$x_range, y = .data$y_range, fill = .data$colors)) +
    geom_bar(stat = "identity", alpha = 0.5) +
    labs(x = "", y = "Probability") +
    scale_fill_manual(values = cols) +
    theme(legend.position = "none")

  sshhr(plt)
}

#' Plot method for the probability calculator (binomial)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param x Return value from \code{\link{prob_binom}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#
#' @seealso \code{\link{prob_binom}} to calculate results
#' @seealso \code{\link{summary.prob_binom}} to summarize results
#'
#' @examples
#' result <- prob_binom(n = 10, p = 0.3, ub = 3)
#' plot(result, type = "values")
#'
#' @importFrom rlang .data
#'
#' @export
plot.prob_binom <- function(x, type = "values", ...) {
  mess <- paste0("mess_", type)
  if (!is.null(x[[mess]])) {
    return(" ")
  }

  if (type == "values") {
    lb <- x$lb
    ub <- x$ub
  } else {
    lb <- x$vlb
    ub <- x$vub
  }

  n <- x$n
  p <- x$p

  limits <- 0:n
  dat <- data.frame(
    x_range = limits,
    y_range = dbinom(limits, size = n, prob = p),
    stringsAsFactors = FALSE
  ) %>%
    filter(., .$y_range > 0.00001)

  make_bar_plot(ub, lb, dat$x_range, dat$y_range)
}

#' Summary method for the probability calculator (binomial)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param object Return value from \code{\link{prob_binom}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_binom}} to calculate results
#' @seealso \code{\link{plot.prob_binom}} to plot results
#'
#' @examples
#' result <- prob_binom(n = 10, p = 0.3, ub = 3)
#' summary(result, type = "values")
#'
#' @export
summary.prob_binom <- function(object, type = "values", ...) {
  cat("Probability calculator\n")
  cat("Distribution: Binomial\n")

  mess <- object[[paste0("mess_", type)]]
  if (!is.null(mess)) {
    return(mess)
  }
  env <- environment()
  ret <- sapply(names(object), function(x) assign(x, object[[x]], envir = env))

  cat("n           :", n, "\n")
  cat("p           :", p, "\n")
  cat("Mean        :", round(n * p, dec), "\n")
  cat("St. dev     :", sqrt(n * p * (1 - p)) %>% round(dec), "\n")

  if (type == "values") {
    cat("Lower bound :", ifelse(is.na(lb), "", lb), "\n")
    cat("Upper bound :", ifelse(is.na(ub), "", ub), "\n")

    if (!is.na(ub) || !is.na(lb)) {
      cat("\n")

      if (!is.na(lb)) {
        cat(paste0("P(X  = ", lb, ") = ", p_elb, "\n"))
        if (lb > 0) {
          cat(paste0("P(X  < ", lb, ") = ", p_lb, "\n"))
          cat(paste0("P(X <= ", lb, ") = ", p_lelb, "\n"))
        }
        if (lb < n) {
          cat(paste0("P(X  > ", lb, ") = ", round(1 - (p_lb + p_elb), dec), "\n"))
          cat(paste0("P(X >= ", lb, ") = ", round(1 - p_lb, dec), "\n"))
        }
      }

      if (!is.na(ub)) {
        cat(paste0("P(X  = ", ub, ") = ", p_eub, "\n"))
        if (ub > 0) {
          cat(paste0("P(X  < ", ub, ") = ", p_ub, "\n"))
          cat(paste0("P(X <= ", ub, ") = ", p_leub, "\n"))
        }
        if (ub < n) {
          cat(paste0("P(X  > ", ub, ") = ", round(1 - (p_ub + p_eub), dec), "\n"))
          cat(paste0("P(X >= ", ub, ") = ", round(1 - p_ub, dec), "\n"))
        }
      }

      if (!is.na(lb) && !is.na(ub)) {
        cat(paste0("P(", lb, " <= X <= ", ub, ")     = ", p_int, "\n"))
        cat(paste0("1 - P(", lb, " <= X <= ", ub, ") = ", round(1 - p_int, dec), "\n"))
      }
    }
  } else {
    cat("Lower bound :", if (is.na(plb)) "\n" else paste0(plb, " (", vlb, ")\n"))
    cat("Upper bound :", if (is.na(pub)) "\n" else paste0(pub, " (", vub, ")\n"))

    if (!is.na(pub) || !is.na(plb)) {
      cat("\n")

      if (!is.na(plb)) {
        cat(paste0("P(X  = ", vlb, ") = ", vp_elb, "\n"))
        if (vlb > 0) {
          cat(paste0("P(X  < ", vlb, ") = ", vp_lb, "\n"))
          cat(paste0("P(X <= ", vlb, ") = ", vp_lelb, "\n"))
        }
        if (vlb < n) {
          cat(paste0("P(X  > ", vlb, ") = ", round(1 - (vp_lb + vp_elb), dec), "\n"))
          cat(paste0("P(X >= ", vlb, ") = ", round(1 - vp_lb, dec), "\n"))
        }
      }

      if (!is.na(pub)) {
        cat(paste0("P(X  = ", vub, ") = ", vp_eub, "\n"))
        if (vub > 0) {
          cat(paste0("P(X  < ", vub, ") = ", vp_ub, "\n"))
          cat(paste0("P(X <= ", vub, ") = ", vp_leub, "\n"))
        }
        if (vub < n) {
          cat(paste0("P(X  > ", vub, ") = ", round(1 - (vp_ub + vp_eub), dec), "\n"))
          cat(paste0("P(X >= ", vub, ") = ", round(1 - vp_ub, dec), "\n"))
        }
      }

      if (!is.na(plb) && !is.na(pub)) {
        cat(paste0("P(", vlb, " <= X <= ", vub, ")     = ", vp_int, "\n"))
        cat(paste0("1 - P(", vlb, " <= X <= ", vub, ") = ", round(1 - vp_int, dec), "\n"))
      }
    }
  }
}

#' Probability calculator for a discrete distribution
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param v Values
#' @param p Probabilities
#' @param lb Lower bound on the number of successes
#' @param ub Upper bound on the number of successes
#' @param plb Lower probability bound
#' @param pub Upper probability bound
#' @param dec Number of decimals to show
#'
#' @seealso \code{\link{summary.prob_disc}} to summarize results
#' @seealso \code{\link{plot.prob_disc}} to plot results
#'
#' @examples
#' prob_disc(v = 1:6, p = 1 / 6, pub = 0.95)
#' prob_disc(v = 1:6, p = c(2 / 6, 2 / 6, 1 / 12, 1 / 12, 1 / 12, 1 / 12), pub = 0.95)
#'
#' @export
prob_disc <- function(v, p, lb = NA, ub = NA,
                      plb = NA, pub = NA, dec = 3) {

  # Think about adding an "expand.grid" setup so you can run this n times. e.g., rolling multiple dice
  # expand.grid(height = 1:6, weight = 1:6)
  rex <- "(\\s*,\\s*|\\s*;\\s*|\\s+)"
  if (is.character(v)) v <- strsplit(v, rex) %>% unlist()
  if (is.character(p)) p <- strsplit(p, rex) %>% unlist()
  rm(rex)

  lp <- length(p)
  lv <- length(v)
  if (lv != lp && lv %% lp == 0) p <- rep(p, lv / lp)

  if (length(v) != length(p)) {
    mess <- "The number of values must be the same or a multiple of the number of probabilities"
    return(list(mess_probs = mess, mess_values = mess) %>% add_class("prob_disc"))
  }

  asNum <- function(x) ifelse(length(x) > 1, as.numeric(x[1]) / as.numeric(x[2]), as.numeric(x[1]))
  if (is.character(v)) v <- sshhr(strsplit(v, "/") %>% sapply(asNum))
  if (is.character(p)) p <- sshhr(strsplit(p, "/") %>% sapply(asNum))

  if (anyNA(p) | anyNA(v)) {
    mess <- "The number of probabilities entered must be a multiple of the number of values"
    mess <- paste0("Invalid inputs:\n\nv: ", paste0(v, collapse = " "), "\np: ", paste0(p, collapse = " "))
    return(list(mess_probs = mess, mess_values = mess) %>% add_class("prob_disc"))
  }

  ## make sure values and probabilities are ordered correctly
  df <- data.frame(v = v, p = p, stringsAsFactors = FALSE) %>%
    arrange(v)
  p <- df$p
  v <- df$v

  if (sum(p) < .99 || sum(p) > 1.01) {
    mess_probs <- mess_values <- paste0("Probabilities for a discrete variable do not sum to 1 (", round(sum(p), 3), ")")
    return(as.list(environment()) %>% add_class("prob_disc"))
  }

  ddisc <- function(b, df) filter(df, v == b)$p
  pdisc <- function(b, df) filter(df, v < b)$p %>% sum()
  ## consistent with http://www.stat.umn.edu/geyer/old/5101/rlook.html#qbinom
  qdisc <- function(prob, df) {
    mutate(df, p = cumsum(df$p)) %>%
      filter(p >= prob) %>%
      .$v %>%
      min()
  }

  if (is.na(lb)) {
    p_elb <- p_lb <- lb <- NA
  } else if (!lb %in% v) {
    p_elb <- 0
    p_lb <- ifelse(lb < min(v), 0, pdisc(lb, df) %>% round(dec))
    p_lelb <- p_elb + p_lb
  } else {
    p_elb <- ddisc(lb, df) %>% round(dec)
    p_lb <- pdisc(lb, df) %>% round(dec)
    p_lelb <- p_elb + p_lb
  }

  if (is.na(ub)) {
    p_eub <- p_ub <- ub <- NA
  } else if (!ub %in% v) {
    p_eub <- 0
    p_ub <- ifelse(ub < min(v), 0, pdisc(ub, df) %>% round(dec))
    p_leub <- p_eub + p_ub
  } else {
    p_eub <- ddisc(ub, df) %>% round(dec)
    p_ub <- pdisc(ub, df) %>% round(dec)
    p_leub <- p_eub + p_ub
  }

  if (!is.na(ub) && !is.na(lb)) {
    p_int <- p_leub - p_lb
  } else {
    p_int <- NA
  }

  if (is.na(plb)) {
    plb <- vlb <- NA
  } else if (length(qdisc(plb, df)) == 0) {
    mess_probs <- "Lower bound is too low"
    return(as.list(environment()) %>% add_class("prob_disc"))
  } else {
    if (plb > 1) plb <- 1
    if (plb < 0) plb <- 0
    vlb <- qdisc(plb, df)
    vp_elb <- ddisc(vlb, df) %>% round(dec)
    vp_lb <- pdisc(vlb, df) %>% round(dec)
    vp_lelb <- vp_elb + vp_lb
  }

  if (is.na(pub)) {
    pub <- vub <- NA
  } else if (length(qdisc(pub, df)) == 0) {
    mess_probs <- "Upper bound is too low"
    return(as.list(environment()) %>% add_class("prob_disc"))
  } else {
    if (pub > 1) pub <- 1
    if (pub < 0) pub <- 0
    vub <- qdisc(pub, df)
    vp_eub <- ddisc(vub, df) %>% round(dec)
    vp_ub <- pdisc(vub, df) %>% round(dec)
    vp_leub <- vp_eub + vp_ub
  }

  if (!is.na(pub) && !is.na(plb)) {
    vp_int <- vp_leub - vp_lb
  } else {
    vp_int <- NA
  }

  if (!is.na(lb) && !is.na(ub)) {
    if (lb > ub) {
      lb <- ub <- NA
      mess_values <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  if (!is.na(vlb) && !is.na(vub)) {
    if (vlb > vub || plb > pub) {
      plb <- pub <- vlb <- vub <- NA
      mess_probs <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  rm(qdisc, pdisc, ddisc, asNum)

  as.list(environment()) %>% add_class("prob_disc")
}

#' Plot method for the probability calculator (discrete)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param x Return value from \code{\link{prob_disc}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_disc}} to calculate results
#' @seealso \code{\link{summary.prob_disc}} to summarize results
#'
#' @examples
#' result <- prob_disc(v = 1:6, p = c(2 / 6, 2 / 6, 1 / 12, 1 / 12, 1 / 12, 1 / 12), pub = 0.95)
#' plot(result, type = "probs")
#'
#' @importFrom rlang .data
#'
#' @export
plot.prob_disc <- function(x, type = "values", ...) {
  mess <- paste0("mess_", type)
  if (!is.null(x[[mess]])) {
    return(" ")
  }

  if (type == "values") {
    lb <- x$lb
    ub <- x$ub
  } else {
    lb <- x$vlb
    ub <- x$vub
  }

  make_bar_plot(ub, lb, x$v, x$p)
}

#' Summary method for the probability calculator (discrete)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param object Return value from \code{\link{prob_disc}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_disc}} to calculate results
#' @seealso \code{\link{plot.prob_disc}} to plot results
#'
#' @examples
#' result <- prob_disc(v = 1:6, p = c(2 / 6, 2 / 6, 1 / 12, 1 / 12, 1 / 12, 1 / 12), pub = 0.95)
#' summary(result, type = "probs")
#'
#' @export
summary.prob_disc <- function(object, type = "values", ...) {
  cat("Probability calculator\n")
  cat("Distribution : Discrete\n")

  mess <- object[[paste0("mess_", type)]]
  if (!is.null(mess)) {
    return(mess)
  }
  env <- environment()
  ret <- sapply(names(object), function(x) assign(x, object[[x]], envir = env))

  cat("Values       :", paste0(v, collapse = " "), "\n")
  cat("Probabilities:", paste0(round(p, dec), collapse = " "), "\n")
  m <- sum(v * p)
  std <- sqrt(sum(p * (v - m)^2))
  cat("Mean         :", round(m, dec), "\n")
  cat("St. dev      :", round(std, dec), "\n")

  if (type == "values") {
    cat("Lower bound  :", ifelse(is.na(lb), "", lb), "\n")
    cat("Upper bound  :", ifelse(is.na(ub), "", ub), "\n")

    if (!is.na(ub) || !is.na(lb)) {
      cat("\n")

      if (!is.na(lb)) {
        cat(paste0("P(X  = ", lb, ") = ", p_elb, "\n"))
        if (lb > min(v)) {
          cat(paste0("P(X  < ", lb, ") = ", p_lb, "\n"))
          cat(paste0("P(X <= ", lb, ") = ", p_lelb, "\n"))
        }
        if (lb < max(v)) {
          cat(paste0("P(X  > ", lb, ") = ", round(1 - (p_lb + p_elb), dec), "\n"))
          cat(paste0("P(X >= ", lb, ") = ", round(1 - p_lb, dec), "\n"))
        }
      }

      if (!is.na(ub)) {
        cat(paste0("P(X  = ", ub, ") = ", p_eub, "\n"))
        if (ub > min(v)) {
          cat(paste0("P(X  < ", ub, ") = ", p_ub, "\n"))
          cat(paste0("P(X <= ", ub, ") = ", p_leub, "\n"))
        }
        if (ub < max(v)) {
          cat(paste0("P(X  > ", ub, ") = ", round(1 - (p_ub + p_eub), dec), "\n"))
          cat(paste0("P(X >= ", ub, ") = ", round(1 - p_ub, dec), "\n"))
        }
      }

      if (!is.na(lb) && !is.na(ub)) {
        cat(paste0("P(", lb, " <= X <= ", ub, ")     = ", p_int, "\n"))
        cat(paste0("1 - P(", lb, " <= X <= ", ub, ") = ", round(1 - p_int, dec), "\n"))
      }
    }
  } else {
    cat("Lower bound  :", if (is.na(plb)) "\n" else paste0(plb, " (", vlb, ")\n"))
    cat("Upper bound  :", if (is.na(pub)) "\n" else paste0(pub, " (", vub, ")\n"))

    if (!is.na(pub) || !is.na(plb)) {
      cat("\n")

      if (!is.na(plb)) {
        cat(paste0("P(X  = ", vlb, ") = ", vp_elb, "\n"))
        if (vlb > min(v)) {
          cat(paste0("P(X  < ", vlb, ") = ", vp_lb, "\n"))
          cat(paste0("P(X <= ", vlb, ") = ", vp_lelb, "\n"))
        }
        if (vlb < max(v)) {
          cat(paste0("P(X  > ", vlb, ") = ", round(1 - (vp_lb + vp_elb), dec), "\n"))
          cat(paste0("P(X >= ", vlb, ") = ", round(1 - vp_lb, dec), "\n"))
        }
      }

      if (!is.na(pub)) {
        cat(paste0("P(X  = ", vub, ") = ", vp_eub, "\n"))
        if (vub > min(v)) {
          cat(paste0("P(X  < ", vub, ") = ", vp_ub, "\n"))
          cat(paste0("P(X <= ", vub, ") = ", vp_leub, "\n"))
        }
        if (vub < max(v)) {
          cat(paste0("P(X  > ", vub, ") = ", round(1 - (vp_ub + vp_eub), dec), "\n"))
          cat(paste0("P(X >= ", vub, ") = ", round(1 - vp_ub, dec), "\n"))
        }
      }

      if (!is.na(plb) && !is.na(pub)) {
        cat(paste0("P(", vlb, " <= X <= ", vub, ")     = ", vp_int, "\n"))
        cat(paste0("1 - P(", vlb, " <= X <= ", vub, ") = ", round(1 - vp_int, dec), "\n"))
      }
    }
  }
}

#' Probability calculator for the exponential distribution
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param rate Rate
#' @param lb Lower bound (default is 0)
#' @param ub Upper bound (default is Inf)
#' @param plb Lower probability bound
#' @param pub Upper probability bound
#' @param dec Number of decimals to show
#'
#' @seealso \code{\link{summary.prob_expo}} to summarize results
#' @seealso \code{\link{plot.prob_expo}} to plot results
#'
#' @examples
#' prob_expo(rate = 1, ub = 2.996)
#'
#' @export
prob_expo <- function(rate, lb = NA, ub = NA,
                      plb = NA, pub = NA, dec = 3) {
  if (!is_not(lb) && lb < 0) lb <- 0
  if (!is_not(ub) && ub < 0) ub <- 0

  p_ub <- pexp(ub, rate)
  p_lb <- pexp(lb, rate)
  p_int <- max(p_ub - p_lb, 0)

  p_ub %<>% round(dec)
  p_lb %<>% round(dec)
  p_int %<>% round(dec)

  if (!is.na(pub)) {
    if (pub > 1) pub <- 1
    if (pub < 0) pub <- 0
  }

  if (!is.na(plb)) {
    if (plb > 1) plb <- 1
    if (plb < 0) plb <- 0
  }

  v_ub <- qexp(pub, rate) %>% round(dec)
  v_lb <- qexp(plb, rate) %>% round(dec)

  if (!is.na(lb) && !is.na(ub)) {
    if (lb > ub) {
      lb <- ub <- NA
      mess_values <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  if (!is.na(plb) && !is.na(pub)) {
    if (plb > pub) {
      plb <- pub <- NA
      mess_probs <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  as.list(environment()) %>% add_class("prob_expo")
}

#' Plot method for the probability calculator (Exponential distribution)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param x Return value from \code{\link{prob_expo}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_expo}} to calculate results
#' @seealso \code{\link{summary.prob_expo}} to summarize results
#'
#' @examples
#' result <- prob_expo(rate = 1, ub = 2.996)
#' plot(result, type = "values")
#'
#' @importFrom rlang .data
#'
#' @export
plot.prob_expo <- function(x, type = "values", ...) {
  mess <- paste0("mess_", type)
  if (!is.null(x[[mess]])) {
    return(" ")
  }

  if (type == "values") {
    lb <- x$lb
    ub <- x$ub
  } else {
    lb <- x$v_lb
    ub <- x$v_ub
  }

  rate <- x$rate

  limits <- c(
    floor(qexp(0.001, rate = rate)),
    ceiling(qexp(1 - 0.001, rate = rate))
  )

  dat <- data.frame(
    x = limits,
    Probability = dexp(limits, rate = rate),
    rate = rate,
    stringsAsFactors = FALSE
  )

  dexp_limit <- function(x) {
    y <- dexp(x, rate = rate)
    y[x < lb | x > ub] <- 0
    y
  }

  dexp_lb <- function(x) {
    if (is.na(lb)) {
      return(0)
    }
    y <- dexp(x, rate = rate)
    y[x > lb] <- 0
    y
  }

  dexp_ub <- function(x) {
    if (is.na(ub)) {
      return(0)
    }
    y <- dexp(x, rate = rate)
    y[x < ub] <- 0
    y
  }

  vlines <- c(ub, lb) %>% na.omit()
  if (length(vlines) == 0) vlines <- c(-Inf, Inf)

  ## based on https://rstudio-pubs-static.s3.amazonaws.com/58753_13e35d9c089d4f55b176057235778679.html
  ## and R Graphics Cookbook
  plt <- ggplot(dat, aes(x = .data$x)) +
    stat_function(fun = stats::dexp, args = list(rate = rate)) +
    stat_function(fun = dexp_limit, geom = "area", fill = "blue", alpha = 0.5, n = 501) +
    stat_function(fun = dexp_lb, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    stat_function(fun = dexp_ub, geom = "area", fill = "red", alpha = 0.5, n = 501) +
    geom_vline(xintercept = vlines, color = "black", linetype = "dashed", linewidth = 0.5) +
    labs(x = "", y = "")

  sshhr(plt)
}


#' Summary method for the probability calculator (exponential)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param object Return value from \code{\link{prob_expo}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_expo}} to calculate results
#' @seealso \code{\link{plot.prob_expo}} to plot results
#'
#' @examples
#' result <- prob_expo(rate = 1, ub = 2.996)
#' summary(result, type = "values")
#'
#' @export
summary.prob_expo <- function(object, type = "values", ...) {
  cat("Probability calculator\n")
  cat("Distribution: Exponential\n")

  mess <- object[[paste0("mess_", type)]]
  if (!is.null(mess)) {
    return(mess)
  }
  env <- environment()
  ret <- sapply(names(object), function(x) assign(x, object[[x]], envir = env))

  cat("Rate        :", rate, "\n")
  cat("Mean        :", round(1 / rate, dec), "\n")
  cat("Variance    :", round(rate^-2, dec), "\n")

  if (type == "values") {
    cat("Lower bound :", if (is.na(lb)) "0" else lb, "\n")
    cat("Upper bound :", if (is.na(ub)) "Inf" else ub, "\n")

    if (!is.na(ub) || !is.na(lb)) {
      cat("\n")

      if (!is.na(lb)) {
        cat(paste0("P(X < ", lb, ") = ", p_lb, "\n"))
        cat(paste0("P(X > ", lb, ") = ", round(1 - p_lb, dec), "\n"))
      }

      if (!is.na(ub)) {
        cat(paste0("P(X < ", ub, ") = ", p_ub, "\n"))
        cat(paste0("P(X > ", ub, ") = ", round(1 - p_ub, dec), "\n"))
      }

      if (!is.na(lb) && !is.na(ub)) {
        cat(paste0("P(", lb, " < X < ", ub, ")     = ", p_int, "\n"))
        cat(paste0("1 - P(", lb, " < X < ", ub, ") = ", round(1 - p_int, dec), "\n"))
      }
    }
  } else {
    pub <- if (is.na(pub)) 2 else pub
    plb <- if (is.na(plb)) -1 else plb

    cat("Lower bound :", if (plb < 0) "0" else plb, "\n")
    cat("Upper bound :", if (pub > 1) "1" else pub, "\n")

    if (pub <= 1 || plb >= 0) {
      cat("\n")

      if (plb >= 0) {
        cat(paste0("P(X < ", v_lb, ") = ", plb, "\n"))
        cat(paste0("P(X > ", v_lb, ") = ", round(1 - plb, dec), "\n"))
      }

      if (pub <= 1) {
        cat(paste0("P(X < ", v_ub, ") = ", pub, "\n"))
        cat(paste0("P(X > ", v_ub, ") = ", round(1 - pub, dec), "\n"))
      }

      if (pub <= 1 && plb >= 0) {
        cat(paste0("P(", v_lb, " < X < ", v_ub, ")     = ", pub - plb, "\n"))
        cat(paste0("1 - P(", v_lb, " < X < ", v_ub, ") = ", round(1 - (pub - plb), dec), "\n"))
      }
    }
  }
}

#' Probability calculator for the poisson distribution
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param lambda Rate
#' @param lb Lower bound (default is 0)
#' @param ub Upper bound (default is Inf)
#' @param plb Lower probability bound
#' @param pub Upper probability bound
#' @param dec Number of decimals to show
#'
#' @seealso \code{\link{summary.prob_pois}} to summarize results
#' @seealso \code{\link{plot.prob_pois}} to plot results
#'
#' @examples
#' prob_pois(lambda = 1, ub = 3)
#'
#' @export
prob_pois <- function(lambda, lb = NA, ub = NA,
                      plb = NA, pub = NA, dec = 3) {
  if (lambda <= 0) mess_values <- "\nLambda must be positive"

  if (!is_not(lb) && lb < 0) lb <- 0
  if (!is_not(ub) && ub < 0) ub <- 0

  if (is.na(lb) || lb < 0) {
    p_elb <- p_lb <- lb <- NA
  } else {
    p_elb <- dpois(lb, lambda) %>% round(dec)
    p_lelb <- ppois(lb, lambda) %>% round(dec)
    if (lb > 0) {
      p_lb <- (ppois(lb, lambda) - dpois(lb, lambda)) %>% round(dec)
    } else {
      p_lb <- 0
    }
  }

  if (is.na(ub) || ub < 0) {
    p_eub <- p_ub <- ub <- NA
  } else {
    p_eub <- dpois(ub, lambda) %>% round(dec)
    p_leub <- ppois(ub, lambda) %>% round(dec)
    if (ub > 0) {
      p_ub <- (ppois(ub, lambda) - dpois(ub, lambda)) %>% round(dec)
    } else {
      p_ub <- 0
    }
  }

  if (!is.na(ub) && !is.na(lb)) {
    p_int <- sum(dpois(lb:ub, lambda)) %>%
      max(0) %>%
      round(dec)
  } else {
    p_int <- NA
  }

  if (is.na(plb)) {
    vlb <- NA
  } else {
    if (plb > 1) plb <- 1
    if (plb < 0) plb <- 0
    vlb <- qpois(plb, lambda)

    vp_elb <- dpois(vlb, lambda) %>% round(dec)
    vp_lelb <- ppois(vlb, lambda) %>% round(dec)
    if (vlb > 0) {
      vp_lb <- (ppois(vlb, lambda) - dpois(vlb, lambda)) %>% round(dec)
    } else {
      vp_lb <- 0
    }
  }

  if (is.na(pub)) {
    vub <- NA
  } else {
    if (pub > 1) pub <- 1
    if (pub < 0) pub <- 0
    vub <- qpois(pub, lambda)

    vp_eub <- dpois(vub, lambda) %>% round(dec)
    vp_leub <- ppois(vub, lambda) %>% round(dec)
    if (vub > 0) {
      vp_ub <- (ppois(vub, lambda) - dpois(vub, lambda)) %>% round(dec)
    } else {
      vp_ub <- 0
    }
  }

  if (!is.na(lb) && !is.na(ub)) {
    if (lb > ub) {
      lb <- ub <- NA
      mess_values <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  if (!is.na(vlb) && !is.na(vub)) {
    if (vlb > vub || plb > pub) {
      plb <- pub <- vlb <- vub <- NA
      mess_probs <- "\nPlease ensure the lower bound is smaller than the upper bound"
    }
  }

  if (!is.na(pub) && !is.na(plb)) {
    vp_int <- sum(dpois(vlb:vub, lambda)) %>%
      max(0) %>%
      round(dec)
  } else {
    vp_int <- NA
  }

  as.list(environment()) %>% add_class("prob_pois")
}

#' Plot method for the probability calculator (poisson)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param x Return value from \code{\link{prob_pois}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_pois}} to calculate results
#' @seealso \code{\link{summary.prob_pois}} to summarize results
#'
#' @examples
#' result <- prob_pois(lambda = 1, ub = 3)
#' plot(result, type = "values")
#'
#' @importFrom rlang .data
#'
#' @export
plot.prob_pois <- function(x, type = "values", ...) {
  mess <- paste0("mess_", type)
  if (!is.null(x[[mess]])) {
    return(" ")
  }

  if (type == "values") {
    lb <- x$lb
    ub <- x$ub
  } else {
    lb <- x$vlb
    ub <- x$vub
  }

  lambda <- x$lambda
  limits <- 0:(ceiling(qpois(1 - 0.00001, lambda)))
  n <- max(limits)

  if (!is.na(lb) && lb > n) {
    limits <- 0:lb
    n <- lb
  }

  if (!is.na(ub) && ub > n) {
    limits <- 0:ub
    n <- ub
  }

  dat <- data.frame(
    x_range = limits,
    y_range = dpois(limits, lambda),
    stringsAsFactors = FALSE
  ) %>% filter(., .$y_range > 0.00001)

  make_bar_plot(ub, lb, dat$x_range, dat$y_range)

}

#' Summary method for the probability calculator (poisson)
#'
#' @details See \url{https://radiant-rstats.github.io/docs/basics/prob_calc.html} for an example in Radiant
#'
#' @param object Return value from \code{\link{prob_pois}}
#' @param type Probabilities ("probs") or values ("values")
#' @param ... further arguments passed to or from other methods
#'
#' @seealso \code{\link{prob_pois}} to calculate results
#' @seealso \code{\link{plot.prob_pois}} to plot results
#'
#' @examples
#' result <- prob_pois(lambda = 1, ub = 3)
#' summary(result, type = "values")
#'
#' @export
summary.prob_pois <- function(object, type = "values", ...) {
  cat("Probability calculator\n")
  cat("Distribution: Poisson\n")

  mess <- object[[paste0("mess_", type)]]
  if (!is.null(mess)) {
    return(mess)
  }
  env <- environment()
  ret <- sapply(names(object), function(x) assign(x, object[[x]], envir = env))

  cat("Lambda      :", lambda, "\n")
  cat("Mean        :", lambda, "\n")
  cat("Variance    :", lambda, "\n")

  if (type == "values") {
    cat("Lower bound :", ifelse(is.na(lb), "", lb), "\n")
    cat("Upper bound :", ifelse(is.na(ub), "", ub), "\n")

    if (!is.na(ub) || !is.na(lb)) {
      cat("\n")

      if (!is.na(lb)) {
        cat(paste0("P(X  = ", lb, ") = ", p_elb, "\n"))
        if (lb > 0) {
          cat(paste0("P(X  < ", lb, ") = ", p_lb, "\n"))
          cat(paste0("P(X <= ", lb, ") = ", p_lelb, "\n"))
        }
        cat(paste0("P(X  > ", lb, ") = ", round(1 - (p_lb + p_elb), dec), "\n"))
        cat(paste0("P(X >= ", lb, ") = ", round(1 - p_lb, dec), "\n"))
      }

      if (!is.na(ub)) {
        cat(paste0("P(X  = ", ub, ") = ", p_eub, "\n"))
        if (ub > 0) {
          cat(paste0("P(X  < ", ub, ") = ", p_ub, "\n"))
          cat(paste0("P(X <= ", ub, ") = ", p_leub, "\n"))
        }
        cat(paste0("P(X  > ", ub, ") = ", round(1 - (p_ub + p_eub), dec), "\n"))
        cat(paste0("P(X >= ", ub, ") = ", round(1 - p_ub, dec), "\n"))
      }

      if (!is.na(lb) && !is.na(ub)) {
        cat(paste0("P(", lb, " <= X <= ", ub, ")     = ", p_int, "\n"))
        cat(paste0("1 - P(", lb, " <= X <= ", ub, ") = ", round(1 - p_int, dec), "\n"))
      }
    }
  } else {
    cat("Lower bound :", if (is.na(plb)) "\n" else paste0(plb, " (", vlb, ")\n"))
    cat("Upper bound :", if (is.na(pub)) "\n" else paste0(pub, " (", vub, ")\n"))

    if (!is.na(pub) || !is.na(plb)) {
      cat("\n")

      if (!is.na(plb)) {
        cat(paste0("P(X  = ", vlb, ") = ", vp_elb, "\n"))
        if (vlb > 0) {
          cat(paste0("P(X  < ", vlb, ") = ", vp_lb, "\n"))
          cat(paste0("P(X <= ", vlb, ") = ", vp_lelb, "\n"))
        }
        cat(paste0("P(X  > ", vlb, ") = ", round(1 - (vp_lb + vp_elb), dec), "\n"))
        cat(paste0("P(X >= ", vlb, ") = ", round(1 - vp_lb, dec), "\n"))
      }

      if (!is.na(pub)) {
        cat(paste0("P(X  = ", vub, ") = ", vp_eub, "\n"))
        if (vub > 0) {
          cat(paste0("P(X  < ", vub, ") = ", vp_ub, "\n"))
          cat(paste0("P(X <= ", vub, ") = ", vp_leub, "\n"))
        }
        cat(paste0("P(X  > ", vub, ") = ", round(1 - (vp_ub + vp_eub), dec), "\n"))
        cat(paste0("P(X >= ", vub, ") = ", round(1 - vp_ub, dec), "\n"))
      }

      if (!is.na(plb) && !is.na(pub)) {
        cat(paste0("P(", vlb, " <= X <= ", vub, ")     = ", vp_int, "\n"))
        cat(paste0("1 - P(", vlb, " <= X <= ", vub, ") = ", round(1 - vp_int, dec), "\n"))
      }
    }
  }
}