## This file is part of the 'agop' library.
##
## Copyleft (c) 2013-2023, Marek Gagolewski <https://www.gagolewski.com/>
##
##
## 'agop' is free software: you can redistribute it and/or modify it under
## the terms of the GNU Lesser General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## 'agop' is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU Lesser General Public License for more details.
##
## A copy of the GNU Lesser General Public License can be downloaded
## from <http://www.gnu.org/licenses/>.
#' @title
#' Pareto Type-II (Lomax) Distribution
#'
#' @description
#' Density, cumulative distribution function,
#' quantile function, and random generation for the
#' Pareto Type-II (Lomax)
#' distribution with shape
#' parameter \eqn{k>0} and scale parameter \eqn{s>0}.
#'
#' [TO DO: rewrite in C, add NA handling]
#'
#' @details
#' If \eqn{X\sim\mathrm{P2}(k,s)}{X~P2(k,s)},
#' then \eqn{\mathrm{supp}\,X=[0,\infty)}{supp X=[0,\infty)}.
#' The c.d.f. for \eqn{x\ge 0} is given by \deqn{F(x)=1-s^k/(s+x)^k}
#' and the density by \deqn{f(x)=k s^k/(s+x)^{k+1}.}
#'
#' @param x,q vector of quantiles
#' @param p vector of probabilities
#' @param n integer; number of observations
#' @param k vector of shape parameters, \eqn{k>0}
#' @param s vector of scale parameters, \eqn{s>0}
#' @param lower.tail logical; if \code{TRUE} (default),
#' probabilities are \eqn{P(X \le x)}, and \eqn{P(X > x)} otherwise
#' @return
#' numeric vector;
#' \code{dpareto2} gives the density,
#' \code{ppareto2} gives the cumulative distribution function,
#' \code{qpareto2} calculates the quantile function,
#' and \code{rpareto2} generates random deviates.
#'
#' @export
#' @rdname Pareto2
#' @family distributions
#' @family Pareto2
rpareto2 <- function(n, k=1, s=1)
{
stopifnot(is.numeric(k), k > 0)
stopifnot(is.numeric(s), s > 0)
# n checked by runif
s*((runif(n)^(-1.0/k)) - 1.0)
}
ppareto2_old <- function(q, k=1, s=1, lower.tail=TRUE)
{
stopifnot(is.numeric(k), k > 0)
stopifnot(is.numeric(s), s > 0)
stopifnot(is.numeric(q))
ret <- ifelse(q<0, 0, (1-(s/(s+q))^k))
if (identical(lower.tail[1], FALSE))
ret <- 1-ret
else
ret
}
#' @export
#' @rdname Pareto2
ppareto2 <- function(q, k=1, s=1, lower.tail=TRUE)
{
.Call("ppareto2", q, k, s, lower.tail, PACKAGE="agop") # args checked internally
}
#' @export
#' @rdname Pareto2
qpareto2 <- function(p, k=1, s=1, lower.tail=TRUE)
{
stopifnot(is.numeric(k), k > 0)
stopifnot(is.numeric(s), s > 0)
stopifnot(is.numeric(p))
if (identical(lower.tail[1], FALSE))
p <- 1-p
ifelse(p<0 | p>1, NaN, s*((1-p)^(-1/k)-1))
}
#' @export
#' @rdname Pareto2
dpareto2 <- function(x, k=1, s=1)
{
stopifnot(is.numeric(k), k > 0)
stopifnot(is.numeric(s), s > 0)
stopifnot(is.numeric(x))
ifelse(x<=0, 0, k/(s+x)*(s/(s+x))^k)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.