# R/qqplot.R In distr6: The Complete R6 Probability Distributions Interface

#### Documented in qqplot

```#' @title Quantile-Quantile Plots for distr6 Objects
#' @author Chijing Zeng
#'
#' @description Quantile-quantile plots are used to compare a "theoretical" or
#' empirical distribution to a reference distribution. They can also compare the quantiles of two
#' reference distributions.
#'
#' @param x \code{distr6} object or numeric vector.
#' @param y \code{distr6} object or numeric vector.
#' @param npoints number of evaluation points.
#' @param plot logical; if TRUE (default), figures are displayed in the plot window; otherwise a
#' [data.table::data.table] of points and calculated values is returned.
#' @param idline logical; if TRUE (default), the line \eqn{y = x} is plotted
#' @param ... graphical parameters.
#'
#' @details If \code{x} or \code{y} are given as numeric vectors then they are first passed to the
#' [Empirical] distribution. The [Empirical] distribution is a discrete distribution so quantiles
#' are equivalent to the the Type 1 method in \code{\link[stats]{quantile}}.
#'
#'
#' @seealso \code{\link{plot.Distribution}} for plotting a \code{distr6} object.
#'
#' @examples
#' qqplot(Normal\$new(mean = 15, sd = sqrt(30)), ChiSquared\$new(df = 15))
#' qqplot(rt(200, df = 5), rt(300, df = 5),
#'   main = "QQ-Plot", xlab = "t-200",
#'   ylab = "t-300"
#' )
#' qqplot(Normal\$new(mean = 2), rnorm(100, mean = 3))
#' @export
qqplot <- function(x, y, npoints = 3000, idline = TRUE, plot = TRUE, ...) {

# random number generator: sample quantiles from 0 to 1 (cdf)
points <- seq(0, 1, length.out = npoints)

if (is.numeric(x)) {
x <- Empirical\$new(x)
}
if (is.numeric(y)) {
y <- Empirical\$new(y)
}

assertDistributionList(list(x, y), "x and y should both be distributions")

quantile.X <- x\$quantile(points)
quantile.Y <- y\$quantile(points)
#
#   nx <- length(quantile.X)
#   ny <- length(quantile.Y)
#
#   if(nx == ny)
#     print(TRUE)
#   else
#     print(FALSE)
#
#   if (nx != ny) {
#     ppoints.x <- ppoints(nx, a = pos)
#     ppoints.y <- ppoints(ny, a = pos)
#     if (nx > ny)
#       quantile.X <- approx(ppoints.x, quantile.X, xout = ppoints.y, rule = 2)\$y
#     else
#       quantile.Y <- approx(ppoints.y, quantile.Y, xout = ppoints.x, rule = 2)\$y
#   }

# Now plot!

## Quantile vs Quantile##
if (plot) {
graphics::plot(x = quantile.X, y = quantile.Y, type = "p", ...)

## Identity Line##
if (idline) {
graphics::abline(a = 0, b = 1, col = "blue", lty = 2, lwd = 2)
}
}

#
#     ##Confidence Lines##
#     x0 <- quantile.X
#     y0 <- quantile.Y
#
#     x0[is.infinite(x0)] <- NA
#     y0[is.infinite(y0)] <- NA
#
#     quantile.df <- data.frame(x0, y0)
#     quantile.df_no_na <- na.omit(quantile.df)
#
#     coef <- coef(MASS::rlm(y0 ~ x0, data = quantile.df_no_na))
#     a <- coef[1]
#     b <- coef[2]
#
#     conf <- if (conf == FALSE) 0.95 else conf
#     xx <- Normal\$new()\$quantile(1 - (1 - conf)/2)
#
#     SE <- (b / Normal\$new()\$pdf(Normal\$new()\$quantile(ppoints(npoints)))) *
#      sqrt(ppoints(npoints)*(1 - ppoints(npoints))/npoints)
#
#     fit.value <- a + b*x0
#     upper <- fit.value + (xx*SE)
#     lower <- fit.value - (xx*SE)

# if (withConf != FALSE) {
#   lines(x0, upper, lty=3, col = "orange", lwd = 2)
#   lines(x0, lower, lty=3, col = "orange", lwd = 2)
# }

invisible(list(x = quantile.X, y = quantile.Y))

}
```

## Try the distr6 package in your browser

Any scripts or data that you put into this service are public.

distr6 documentation built on March 28, 2022, 1:05 a.m.