R/GeoQQ.R

In GeoModels: Procedures for Gaussian and Non Gaussian Geostatistical (Large) Data Analysis

####################################################
### File name: GeoQQ.r (cleaned; comments in English)
####################################################
GeoQQ <- function(fit, type="Q", add=FALSE, ylim=c(0,1), xlim=NULL, breaks=10, ...)
{
  ###### utilities ##########
  # Log-logistic quantile with shape (>0) and scale (>0)

  # Log-logistic density (fix: density, not quantile)
  ##################################

  if(!inherits(fit,"GeoFit"))  stop("A GeoFit object is needed as input\n")
  if(identical(fit$residual_type, "Pearson")) {
    stop("GeoQQ() is not defined for Pearson-residual objects. Use GeoQQ() or GeoPit() on the original discrete GeoFit object, and use GeoResiduals() output for residual maps or variograms.",
         call. = FALSE)
  }
  if(type!="Q" && type!="D") stop("Type can be Q or D \n")

  model <- fit$model
  discrete_models <- c(
    "Binary", "Bernoulli", "Binomial", "BinomialLogistic",
    "Geom", "Geometric", "BinomialNeg", "Poisson", "PoissonGamma",
    "PoissonZIP", "PoissonZIP1", "PoissonGammaZIP", "PoissonGammaZIP1",
    "BinomialNegZINB", "BinomialNegZINB1"
  )

  pp_tmp <- c(unlist(fit$param), unlist(fit$fixed))
  external_mean <- !is.null(.GeoMean_external(fit$fixed))
  non_intercept_X <- FALSE
  if (!is.null(fit$X)) {
    Xqq <- .GeoMean_design(
      fit$X, n = length(unlist(fit$data, use.names = FALSE)), name = "X"
    )
    non_intercept_X <- ncol(Xqq) != 1L ||
      any(abs(Xqq[, 1L] - 1) > sqrt(.Machine$double.eps))
  }
  if (external_mean || non_intercept_X) {
    if (model %in% discrete_models) {
      stop(
        "GeoQQ() does not define a single theoretical marginal distribution for discrete fits with a nonconstant mean. Use GeoPit(fit) for site-specific marginal diagnostics.",
        call. = FALSE
      )
    }
    stop("GeoQQ() requires an intercept-only marginal fit. ",
         "Please use GeoResiduals(fit) before calling GeoQQ().")
  }

  .GeoQQ_constant_n <- function(n, model_name) {
    nn <- as.numeric(n)
    if (!length(nn) || any(!is.finite(nn)) || any(nn <= 0) || any(nn != round(nn))) {
      stop("For model '", model_name, "', n must contain positive integers.", call. = FALSE)
    }
    if (length(nn) > 1L && any(nn != nn[1L])) {
      stop(
        "GeoQQ() requires a constant n for model '", model_name,
        "'. Use GeoPit(fit) when n is site-specific.",
        call. = FALSE
      )
    }
    nn[1L]
  }

  if (model %in% c("Binomial", "BinomialLogistic", "BinomialNeg")) {
    fit$n <- .GeoQQ_constant_n(fit$n, model)
  }

  copula <- if(!is.null(fit$copula)) fit$copula else NULL

  fit$param <- unlist(fit$param); fit$fixed <- unlist(fit$fixed)
  pp <- c(fit$param, fit$fixed)

  # Safe extraction: single-bracket returns NA when name is absent (no error)
  MM <- as.numeric(pp["mean"])
  VV <- as.numeric(pp["sill"])
  if (is.na(MM)) MM <- 0


  # Which models actually require 'sill' (marginal variance)?
  .requires_VV <- function(m) {
    m %in% c(
      "Gaussian","Gaussian_misp_Binomial","Gaussian_misp_Poisson","Gaussian_misp_BinomialNeg",
      "SkewGaussian","SkewLaplace",
      "StudentT","Gaussian_misp_StudentT",
      "SkewStudentT","Gaussian_misp_SkewStudentT",
      "LogGaussian","Logistic",
      "SinhAsinh","Tukeyh","Tukeyh2",
      "Tukeygh","Gaussian_misp_Tukeygh",
      "TwoPieceGaussian","TwoPieceStudentT","TwoPieceTukeyh"
    )
  }

  # Save and restore par settings
  opar <- par(no.readonly = TRUE); on.exit(par(opar))

  # Extract main from ... and remove it to avoid conflicts
  dots <- list(...)
  main_title <- if("main" %in% names(dots)) dots$main else NULL
  if("main" %in% names(dots)) dots <- dots[names(dots) != "main"]
  
  # Helper to get main title (use default if main_title is NULL)
  .get_main <- function(default) {
    if(is.null(main_title)) default else main_title
  }

  # xlim helpers
  has_xlim <- !is.null(xlim)
  .XL <- function(def) { if (has_xlim) xlim else def }
  .GRID <- function(def_range, n = 512) {
    xr <- if (has_xlim) xlim else def_range
    seq(xr[1], xr[2], length.out = n)
  }

  # Library checks for skew/t families when needed
  need_sn   <- model %in% c("SkewGaussian","SkewStudentT","Gaussian_misp_SkewStudentT")
  need_VGAM <- model %in% c("Tukeyh","Tukeyh2","TwoPieceTukeyh","Tukeygh","Gaussian_misp_Tukeygh")
  if (need_sn && !requireNamespace("sn", quietly = TRUE))
    stop("Package 'sn' is required for skew-* models.")
  if (need_VGAM && !requireNamespace("VGAM", quietly = TRUE))
    stop("Package 'VGAM' is required for Tukey-h/gh models.")

  ########################################
  #### QQ plot
  ########################################
  if(type=="Q") {
    xlab <- "Theoretical Quantiles"
    ylab <- "Sample Quantiles"

    if(!fit$bivariate){
      dd <- if (is.list(fit$coordx_dyn)) unlist(fit$data) else c(t(fit$data))
      dd <- dd[is.finite(dd)]
      if (!length(dd)) stop("empty data")

      N <- length(dd)
      probabilities  <- (1:N)/(N+1)
      probabilities1 <- c(0.25,0.75)
      q_e  <- quantile(dd, probabilities)
      q_e1 <- quantile(dd, probabilities1)

      if (.requires_VV(model) && is.na(VV)) {
        stop("'sill' parameter is required for model '", model,
             "' but is missing in 'fit'. Provide 'sill' in fit$param/fixed.")
      }

      # ---- Models ----
      if(model %in% c("Beta2")){
        mm <- 1/(1+exp(-MM)); sh <- as.numeric(pp["shape"])
        pmin <- as.numeric(pp["min"]); pmax <- as.numeric(pp["max"])
        q_t  <- pmin + (pmax-pmin)*qbeta(probabilities,  shape1=mm*sh, shape2=(1-mm)*sh)
        q_t1 <- pmin + (pmax-pmin)*qbeta(probabilities1, shape1=mm*sh, shape2=(1-mm)*sh)
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Beta qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("Kumaraswamy2")){
        mm <- 1/(1+exp(-MM)); sh <- as.numeric(pp["shape"])
        pmin <- as.numeric(pp["min"]); pmax <- as.numeric(pp["max"])
        shape2 <- log(0.5)/log1p(-(mm^sh))
        q_t  <- pmin + (pmax-pmin)*(1 - (1 - probabilities)^(1/shape2))^(1/sh)
        q_t1 <- pmin + (pmax-pmin)*(1 - (1 - probabilities1)^(1/shape2))^(1/sh)
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Kumaraswamy qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("Gaussian","Gaussian_misp_Binomial","Gaussian_misp_Poisson","Gaussian_misp_BinomialNeg")) {
        q_t  <- qnorm(probabilities,  mean=MM, sd=sqrt(VV))
        q_t1 <- qnorm(probabilities1, mean=MM, sd=sqrt(VV))
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Gaussian qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("Binomial")) {
        q_t  <- qbinom(probabilities,  size=fit$n, prob=pnorm(pp["mean"]))
        q_t1 <- qbinom(probabilities1, size=fit$n, prob=pnorm(pp["mean"]))
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Binomial qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("BinomialLogistic")) {
        q_t  <- qbinom(probabilities,  size=fit$n, prob=plogis(MM))
        q_t1 <- qbinom(probabilities1, size=fit$n, prob=plogis(MM))
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Binomial-Logistic qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("BinomialNeg")) {
        q_t  <- qnbinom(probabilities,  size=fit$n, prob=pnorm(pp["mean"]))
        q_t1 <- qnbinom(probabilities1, size=fit$n, prob=pnorm(pp["mean"]))
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Binomial Neg qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("Poisson")) {
        q_t  <- qpois(probabilities,  lambda=exp(pp["mean"]))
        q_t1 <- qpois(probabilities1, lambda=exp(pp["mean"]))
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Poisson qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("PoissonGamma")) {
        ff <- exp(pp["mean"])
        sh <- as.numeric(pp["shape"])
        q_t  <- qnbinom(probabilities,  size=sh, prob=sh/(ff+sh))
        q_t1 <- qnbinom(probabilities1, size=sh, prob=sh/(ff+sh))
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Poisson qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("SkewGaussian")) {
        omega <- as.numeric(sqrt((pp["skew"]^2 + pp["sill"])/pp["sill"]))
        alpha <- as.numeric(pp["skew"]/sqrt(pp["sill"]))
        q_t  <- MM + sqrt(VV)*sn::qsn(probabilities,  xi=0, omega=omega, alpha=alpha)
        q_t1 <- MM + sqrt(VV)*sn::qsn(probabilities1, xi=0, omega=omega, alpha=alpha)
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Skew Gaussian qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model%in%c("StudentT","Gaussian_misp_StudentT")){
        df <- as.numeric(1/pp["df"])
        q_t  <- MM + sqrt(VV)*qt(probabilities,  df=df)
        q_t1 <- MM + sqrt(VV)*qt(probabilities1, df=df)
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("t qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model%in%c("SkewStudentT","Gaussian_misp_SkewStudentT")){
        alpha <- as.numeric(pp["skew"])
        nu    <- as.numeric(1/pp["df"])
        q_t  <- MM + sqrt(VV)*sn::qst(probabilities,  xi=0, omega=1, alpha=alpha, nu=nu)
        q_t1 <- MM + sqrt(VV)*sn::qst(probabilities1, xi=0, omega=1, alpha=alpha, nu=nu)
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("skewt qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("Weibull")){
        shape <- pp["shape"]
        q_t  <- exp(MM)*qweibull(probabilities,  shape=shape, scale=1/(gamma(1+1/shape)))
        q_t1 <- exp(MM)*qweibull(probabilities1, shape=shape, scale=1/(gamma(1+1/shape)))
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Weibull qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("Gamma")){
        shape <- pp["shape"]
        q_t  <- exp(MM)*qgamma(probabilities,  shape=shape/2, rate=shape/2)
        q_t1 <- exp(MM)*qgamma(probabilities1, shape=shape/2, rate=shape/2)
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Gamma qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("LogGaussian")){
        q_t  <- qlnorm(probabilities,  meanlog = MM - VV/2, sdlog = sqrt(VV))
        q_t1 <- qlnorm(probabilities1, meanlog = MM - VV/2, sdlog = sqrt(VV))
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("LogGaussian qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("LogLogistic")){
        shape <- pp["shape"]
        cc <- gamma(1+1/shape)*gamma(1-1/shape)
        q_t  <- .GeoMarginalQLogLogistic(probabilities,  shape=shape, scale=exp(MM)/cc)
        q_t1 <- .GeoMarginalQLogLogistic(probabilities1, shape=shape, scale=exp(MM)/cc)
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("LogLogistic qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("Logistic")){
        q_t  <- qlogis(probabilities,  location=MM, scale=sqrt(VV))
        q_t1 <- qlogis(probabilities1, location=MM, scale=sqrt(VV))
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Logistic qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("SinhAsinh")){
        tail <- as.numeric(pp["tail"]); skew <- as.numeric(pp["skew"])
        q_t  <- MM + sqrt(VV)*sinh( (asinh(qnorm(probabilities))/tail) + skew/tail )
        q_t1 <- MM + sqrt(VV)*sinh( (asinh(qnorm(probabilities1))/tail) + skew/tail )
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Sas qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("Tukeyh")){
        tail <- as.numeric(pp["tail"])
        uu <- qnorm(probabilities); uu1 <- qnorm(probabilities1)
        q_t  <- MM + sqrt(VV)*uu * exp(0.5*tail*uu^2)
        q_t1 <- MM + sqrt(VV)*uu1* exp(0.5*tail*uu1^2)
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Tukey-h qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("Tukeyh2")){
        tail1 <- as.numeric(pp["tail1"]); tail2 <- as.numeric(pp["tail2"])
        q_t  <- MM + sqrt(VV)*.GeoMarginalQTukeyH2(probabilities, tail1, tail2)
        q_t1 <- MM + sqrt(VV)*.GeoMarginalQTukeyH2(probabilities1, tail1, tail2)
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Tukey-hh qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("SkewLaplace")){
        skew <- as.numeric(pp["skew"])
        q_t  <- MM + sqrt(VV)*.GeoMarginalQSkewLaplace(probabilities, skew)
        q_t1 <- MM + sqrt(VV)*.GeoMarginalQSkewLaplace(probabilities1, skew)
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Skew Laplace qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("Tukeygh","Gaussian_misp_Tukeygh")){
        tail <- as.numeric(pp["tail"]); skew <- as.numeric(pp["skew"])
        uu <- qnorm(probabilities); uu1 <- qnorm(probabilities1)
        if (abs(skew) < 1e-8) {
          q_t  <- MM + sqrt(VV) * (uu  * exp(0.5*tail*uu^2))
          q_t1 <- MM + sqrt(VV) * (uu1 * exp(0.5*tail*uu1^2))
        } else {
          q_t  <- MM + sqrt(VV) * ((exp(skew*uu)  - 1) * exp(0.5*tail*uu^2)  / skew)
          q_t1 <- MM + sqrt(VV) * ((exp(skew*uu1) - 1) * exp(0.5*tail*uu1^2) / skew)
        }
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Tukey-gh qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("TwoPieceGaussian")){
        skew <- as.numeric(pp["skew"])
        q_t  <- MM + sqrt(VV)*.GeoMarginalQTwoPieceGaussian(probabilities, skew)
        q_t1 <- MM + sqrt(VV)*.GeoMarginalQTwoPieceGaussian(probabilities1, skew)
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Two-Piece Gaussian qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("TwoPieceStudentT")){
        skew <- as.numeric(pp["skew"]); df <- 1/as.numeric(pp["df"])
        q_t  <- MM + sqrt(VV)*.GeoMarginalQTwoPieceStudentT(probabilities, skew, df)
        q_t1 <- MM + sqrt(VV)*.GeoMarginalQTwoPieceStudentT(probabilities1, skew, df)
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Two-Piece Student qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      if(model %in% c("TwoPieceTukeyh")){
        skew <- as.numeric(pp["skew"]); tail <- as.numeric(pp["tail"])
        q_t  <- MM + sqrt(VV)*.GeoMarginalQTwoPieceTukeyH(probabilities, skew, tail)
        q_t1 <- MM + sqrt(VV)*.GeoMarginalQTwoPieceTukeyH(probabilities1, skew, tail)
        do.call(plot, c(list(x=q_t, y=q_e, main=.get_main("Two-Piece Tukey-h qq-plot"), 
                             xlab=xlab, ylab=ylab, xlim=.XL(range(q_t, finite=TRUE))), dots))
      }

      x <- q_t1; y <- q_e1
      slope <- diff(y)/diff(x); int <- y[1L] - slope * x[1L]
      abline(int, slope, pch=20)
    }

    if(fit$bivariate){
      par(mfrow=c(1,2))
      if(is.list(fit$coordx_dyn)){ dd1 <- fit$data[[1]]; dd2 <- fit$data[[2]]}
      else  {dd1 <- fit$data[1,]; dd2 <- fit$data[2,];}
      dd1 <- dd1[is.finite(dd1)]; dd2 <- dd2[is.finite(dd2)]

      if(model %in% c("Gaussian")) {
        if (has_xlim) {
          qqnorm(dd1, main=.get_main("First Gaussian qq-plot"), xlim=xlim); abline(0,1)
          qqnorm(dd2, main=.get_main("Second Gaussian qq-plot"), xlim=xlim); abline(0,1)
        } else {
          qqnorm(dd1, main=.get_main("First Gaussian qq-plot")); abline(0,1)
          qqnorm(dd2, main=.get_main("Second Gaussian qq-plot")); abline(0,1)
        }
      }

      if(model %in% c("SkewGaussian")){
        omega1 <- sqrt((pp["skew_1"]^2 + pp["sill_1"])/pp["sill_1"])
        alpha1 <- pp["skew_1"]/sqrt(pp["sill_1"])
        probabilities <- (1:length(dd1))/(length(dd1)+1)
        q_t1 <- sn::qsn(probabilities, xi=0, omega=as.numeric(omega1), alpha=as.numeric(alpha1))
        q_e1 <- quantile(dd1, probabilities)
        do.call(plot, c(list(x=q_t1, y=q_e1, main=.get_main("First Skew Gaussian qq-plot"),
                             xlim=.XL(range(q_t1, finite=TRUE))), dots))
        aa <- lm(q_e1~1+q_t1); abline(as.numeric(aa$coefficients[1]),as.numeric(aa$coefficients[2]))

        omega2 <- sqrt((pp["skew_2"]^2 + pp["sill_2"])/pp["sill_2"])
        alpha2 <- pp["skew_2"]/sqrt(pp["sill_2"])
        probabilities <- (1:length(dd2))/(length(dd2)+1)
        q_t2 <- sn::qsn(probabilities, xi=0, omega=as.numeric(omega2), alpha=as.numeric(alpha2))
        q_e2 <- quantile(dd2, probabilities)
        do.call(plot, c(list(x=q_t2, y=q_e2, main=.get_main("Second Skew Gaussian qq-plot"),
                             xlim=.XL(range(q_t2, finite=TRUE))), dots))
        aa <- lm(q_e2~1+q_t2); abline(as.numeric(aa$coefficients[1]),as.numeric(aa$coefficients[2]))
      }
    }
  }

  ##########################################################
  #### Density / histogram
  ##########################################################
  if(type=="D") {
    if(!fit$bivariate){

      dd <- if (is.list(fit$coordx_dyn)) unlist(fit$data) else c(t(fit$data))
      dd <- dd[is.finite(dd)]
      if (!length(dd)) stop("empty data")

      if (.requires_VV(model) && is.na(VV)) {
        stop("'sill' parameter is required for model '", model,
             "' but is missing in 'fit'. Provide 'sill' in fit$param/fixed.")
      }

      .hist_once <- function(main_title_default, xr_default) {
        if(!add) {
          do.call(hist, c(list(x=dd, freq=FALSE, xlim=.XL(xr_default), ylab="Density", xlab="",
                               main=.get_main(main_title_default), ylim=ylim, breaks=breaks), dots))
        }
      }

      # ------------ Continuous models ------------
      if(model %in% c("Beta2")){
        mm <- 1/(1+exp(-MM)); sh <- as.numeric(pp["shape"])
        pmin <- as.numeric(pp["min"]); pmax <- as.numeric(pp["max"])
        xr <- .XL(c(pmin, pmax))
        ll <- .GRID(xr, n = 800)
        ds <- ifelse(ll>=pmin & ll<=pmax,
                     dbeta((ll-pmin)/(pmax-pmin), shape1=mm*sh, shape2=(1-mm)*sh)/(pmax-pmin), 0)
        .hist_once("Beta Histogram", xr_default = c(pmin,pmax))
        lines(ll, ds, ...)
      }

      if(model %in% c("Kumaraswamy2")){
        sh <- as.numeric(pp["shape"])
        pmin <- as.numeric(pp["min"]); pmax <- as.numeric(pp["max"])
        xr <- .XL(c(pmin,pmax))
        ll <- .GRID(xr, n = 800)
        dkuma <- function(x,MM,sh,pmin,pmax){
          out <- numeric(length(x))
          idx <- which(x>=pmin & x<=pmax)
          if (length(idx)) {
            q <- (x[idx]-pmin)/(pmax-pmin)
            m1 <- 1/(1+exp(-MM))
            shapei <- log(0.5)/log1p(-m1^sh)
            out[idx] <- shapei*sh*q^(sh-1)*(1-q^sh)^(shapei-1)/(pmax-pmin)
          }
          out
        }
        ds <- dkuma(ll, MM, sh, pmin, pmax)
        .hist_once("Kumaraswamy Histogram", xr_default = c(pmin,pmax))
        lines(ll, ds, ...)
      }

      if(model%in%c("Gaussian")){
        xr <- .XL(range(dd, finite=TRUE))
        ll <- .GRID(xr)
        .hist_once("Gaussian Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, dnorm(ll, mean=MM, sd=sqrt(VV)), ...)
      }

      if(model%in%c("StudentT","Gaussian_misp_StudentT")){
        df <- as.numeric(1/pp["df"])
        xr <- .XL(range(dd, finite=TRUE))
        ll <- .GRID(xr)
        .hist_once("Student T Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, dt((ll - MM)/sqrt(VV), df = df) / sqrt(VV), ...)
      }

      if(model%in%c("SkewStudentT","Gaussian_misp_SkewStudentT")){
        alpha <- as.numeric(pp["skew"]); nu <- as.numeric(1/pp["df"])
        xr <- .XL(range(dd, finite=TRUE))
        ll <- .GRID(xr)
        d_st <- sn::dst((ll-MM)/sqrt(VV), xi=0, omega=1, alpha=alpha, nu=nu)/sqrt(VV)
        .hist_once("Skew-T Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, d_st, ...)
      }

      if(model %in% c("SkewGaussian")){
        skew <- as.numeric(pp["skew"]); sill <- as.numeric(pp["sill"])
        omega <- sqrt((skew^2 + sill)/sill); alpha <- skew/sqrt(sill)
        xr <- .XL(range(dd, finite=TRUE))
        ll <- .GRID(xr)
        d_sn <- sn::dsn((ll-MM)/sqrt(VV), xi=0, omega=omega, alpha=alpha)/sqrt(VV)
        .hist_once("Skew Gaussian Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, d_sn, ...)
      }

      if(model %in% c("Weibull")){
        shape <- pp["shape"]
        xr <- .XL(range(dd, finite=TRUE))
        ll <- .GRID(xr)
        d_w <- ifelse(ll>0, dweibull(ll, shape=shape, scale=exp(MM)/(gamma(1+1/shape))), 0)
        .hist_once("Weibull Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, d_w, ...)
      }

     if(model %in% c("Gamma")){
       shape <- as.numeric(pp["shape"])
       xr <- .XL(range(dd, finite = TRUE))
       xr[1] <- max(0, xr[1])
       ll <- seq(from = max(xr[1], .Machine$double.eps), to = xr[2], length.out = 512)
       k    <- shape/2; rate <- shape/(2*exp(MM))  # coerente con qq
       d_g <- dgamma(ll, shape = k, rate = rate)
       .hist_once("Gamma Histogram", xr_default = xr)
       lines(ll, d_g, ...)
}


      if(model %in% c("LogGaussian")){
        xr <- .XL(range(dd, finite=TRUE))
        ll <- .GRID(xr)
        d_l <- ifelse(ll>0, dlnorm(ll, meanlog = MM - VV/2, sdlog = sqrt(VV)), 0)
        .hist_once("LogGaussian Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, d_l, ...)
      }

      if(model %in% c("Logistic")){
        xr <- .XL(range(dd, finite=TRUE))
        ll <- .GRID(xr)
        d_l <- dlogis(ll, location = MM, scale = sqrt(VV))
        .hist_once("Logistic Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, d_l, ...)
      }

      if(model %in% c("LogLogistic")){
        shape <- pp["shape"]
        cc <- gamma(1+1/shape)*gamma(1-1/shape)
        xr <- .XL(range(dd, finite=TRUE))
        ll <- .GRID(xr)
        d_l <- .GeoMarginalDLogLogistic(ll, shape = shape, scale = exp(MM)/cc)
        .hist_once("LogLogistic Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, d_l, ...)
      }

      if(model %in% c("SinhAsinh")){
        tail <- as.numeric(pp["tail"]); skew <- as.numeric(pp["skew"])
        xr <- .XL(range(dd, finite=TRUE))
        ll <- .GRID(xr)
        dsas <- function(z, skew, tail){
          s <- sinh(tail*asinh(z) - skew)
          a <- (2*pi*(1+z^2))^(-0.5)
          c0 <- sqrt(1+s^2)
          d0 <- exp(-0.5*s^2)
          tail*c0*d0*a
        }
        ds <- dsas((ll-MM)/sqrt(VV), skew, tail)/sqrt(VV)
        .hist_once("SAS Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, ds, ...)
      }

      if(model %in% c("SkewLaplace")){
        skew <- as.numeric(pp["skew"])
        xr <- .XL(range(dd, finite=TRUE))
        ll <- .GRID(xr)
        dslaplace <- function(x, sk) {
          if (sk <= 0 || sk >= 1) stop("sk debe estar en (0,1)")
          f <- numeric(length(x))
          sel1 <- x < 0
          f[sel1] <- sk * (1 - sk) * exp((1 - sk) * x[sel1])
          sel2 <- !sel1
          f[sel2] <- sk * (1 - sk) * exp(-sk * x[sel2])
          return(f)
        }
        ds <- dslaplace((ll-MM)/sqrt(VV), skew)/sqrt(VV)
        .hist_once("Skew Laplace Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, ds, ...)
      }

      if(model %in% c("Tukeyh")){
        xr <- .XL(range(dd, finite=TRUE))
        ll <- .GRID(xr)
        ds <- .GeoMarginalDTukeyH((ll-MM)/sqrt(VV), as.numeric(pp["tail"]))/sqrt(VV)
        .hist_once("Tukey-h Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, ds, ...)
      }

      if(model %in% c("Tukeyh2")){
        tail1 <- as.numeric(pp["tail1"]); tail2 <- as.numeric(pp["tail2"])
        xr <- .XL(range(dd, finite=TRUE))
        ll <- .GRID(xr)
        zs <- (ll-MM)/sqrt(VV)
        ds <- numeric(length(zs))
        idx1 <- which(zs>=0); idx2 <- which(zs<0)
        ds[idx1] <- .GeoMarginalDTukeyH(zs[idx1], tail1)
        ds[idx2] <- .GeoMarginalDTukeyH(zs[idx2], tail2)
        ds <- ds / sqrt(VV)
        .hist_once("Tukey-hh Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, ds, ...)
      }

      if(model %in% c("TwoPieceGaussian")){
        skew <- as.numeric(pp["skew"])
        xr <- .XL(range(dd, finite=TRUE)); ll <- .GRID(xr)
        z <- (ll-MM)/sqrt(VV)
        ds <- numeric(length(z))
        idx1 <- which(z>=0); idx2 <- which(z<0)
        ds[idx1] <- dnorm(z[idx1]/(1-skew),0,1) / sqrt(VV)
        ds[idx2] <- dnorm(z[idx2]/(1+skew),0,1) / sqrt(VV)
        .hist_once("Two-Piece Gaussian Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, ds, ...)
      }

      if(model %in% c("TwoPieceStudentT")){
        skew <- as.numeric(pp["skew"]); df <- 1/as.numeric(pp["df"])
        xr <- .XL(range(dd, finite=TRUE)); ll <- .GRID(xr)
        z <- (ll-MM)/sqrt(VV)
        ds <- numeric(length(z))
        idx1 <- which(z>=0); idx2 <- which(z<0)
        ds[idx1] <- dt(z[idx1]/(1-skew), df=df) / sqrt(VV)
        ds[idx2] <- dt(z[idx2]/(1+skew), df=df) / sqrt(VV)
        .hist_once("Two-Piece Student Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, ds, ...)
      }

      if(model %in% c("TwoPieceTukeyh")){
        skew <- as.numeric(pp["skew"]); tail <- as.numeric(pp["tail"])
        xr <- .XL(range(dd, finite=TRUE)); ll <- .GRID(xr)
        z <- (ll-MM)/sqrt(VV)
        ds <- numeric(length(z))
        idx1 <- which(z>=0); idx2 <- which(z<0)
        ds[idx1] <- .GeoMarginalDTukeyH(z[idx1]/(1-skew), tail) / sqrt(VV)
        ds[idx2] <- .GeoMarginalDTukeyH(z[idx2]/(1+skew), tail) / sqrt(VV)
        .hist_once("Two-Piece Tukey-h Histogram", xr_default = range(dd, finite=TRUE))
        lines(ll, ds, ...)
      }

      # ------------ Discrete models ------------
      .disc_support <- function() {
        if (has_xlim) {
          ys <- seq(ceiling(xlim[1]), floor(xlim[2]), by=1)
          if (length(ys) == 0) ys <- sort(unique(as.numeric(dd)))
        } else {
          ys <- sort(unique(as.numeric(dd)))
        }
        ys
      }

      if(model %in% c("BinomialNeg")) {
        count <- as.numeric(table(dd)); ll_emp <- count/sum(count)
        y_emp <- sort(unique(as.numeric(dd)))
        ys <- .disc_support()
        ds <- dnbinom(ys, size=fit$n, prob=pnorm(MM))
        if(!add) {
          do.call(plot, c(list(x=y_emp, y=ll_emp, type="h", col="blue", 
                               main=.get_main("Binomial Negative Histogram"),
                               ylab="Density", xlab="", xlim=.XL(range(c(y_emp, ys))), 
                               ylim=ylim, lwd=2), dots))
        }
        points(ys, ds, type="p", lwd=3)
        lines(ys, ds)
      }

      if(model %in% c("Binomial")) {
        count <- as.numeric(table(dd)); ll_emp <- count/sum(count)
        y_emp <- sort(unique(as.numeric(dd)))
        ys <- .disc_support()
        ds <- dbinom(ys, size=fit$n, prob=pnorm(MM))
        if(!add) {
          do.call(plot, c(list(x=y_emp, y=ll_emp, type="h", col="blue", 
                               main=.get_main("Binomial Histogram"),
                               ylab="Density", xlab="", xlim=.XL(range(c(y_emp, ys))), 
                               ylim=ylim, lwd=2), dots))
        }
        points(ys, ds, type="p", lwd=3)
        lines(ys, ds)
      }

      if(model %in% c("BinomialLogistic")) {
        count <- as.numeric(table(dd)); ll_emp <- count/sum(count)
        y_emp <- sort(unique(as.numeric(dd)))
        ys <- .disc_support()
        ds <- dbinom(ys, size=fit$n, prob=plogis(MM))
        if(!add) {
          do.call(plot, c(list(x=y_emp, y=ll_emp, type="h", col="blue", 
                               main=.get_main("Binomial-Logistic Histogram"),
                               ylab="Density", xlab="", xlim=.XL(range(c(y_emp, ys))), 
                               ylim=ylim, lwd=2), dots))
        }
        points(ys, ds, type="p", lwd=3)
        lines(ys, ds)
      }

      if(model %in% c("Poisson")) {
        count <- as.numeric(table(dd)); ll_emp <- count/sum(count)
        y_emp <- sort(unique(as.numeric(dd)))
        ys <- .disc_support()
        ds <- dpois(ys, lambda=exp(MM))
        if(!add) {
          do.call(plot, c(list(x=y_emp, y=ll_emp, type="h", col="blue", 
                               main=.get_main("Poisson Histogram"),
                               ylab="Density", xlab="", xlim=.XL(range(c(y_emp, ys))), 
                               ylim=ylim, lwd=2), dots))
        }
        points(ys, ds, type="p", lwd=3)
        lines(ys, ds)
      }

      if(model %in% c("PoissonGamma")) {
        count <- as.numeric(table(dd)); ll_emp <- count/sum(count)
        y_emp <- sort(unique(as.numeric(dd)))
        ys <- .disc_support()
        ff <- exp(MM); sh <- as.numeric(pp["shape"])
        ds <- dnbinom(ys, size=sh, prob=sh/(ff+sh))
        if(!add) {
          do.call(plot, c(list(x=y_emp, y=ll_emp, type="h", col="blue", 
                               main=.get_main("Poisson Gamma (NegBin) Histogram"),
                               ylab="Density", xlab="", xlim=.XL(range(c(y_emp, ys))), 
                               ylim=ylim, lwd=2), dots))
        }
        points(ys, ds, type="p", lwd=3)
        lines(ys, ds)
      }
    }
  }

  invisible()
}