# R/BMT.Psi.R In BMT: The BMT Distribution

#' @title The BMT-Psi Distribution.
#' @description Density, distribution function, quantile function, random number
#'   generation for the BMT-Psi distribution with mean equal to \code{mean} and
#'   standard deviation equal to \code{sd}.
#' @rdname BMT.Psi
#' @name BMT.Psi
#' @aliases dBMT.Psi
#' @aliases pBMT.Psi
#' @aliases qBMT.Psi
#' @aliases rBMT.Psi
#'
#' @details If \code{mean} or \code{sd} are not specified they assume the
#'   default values of 0 and 1, respectively.
#'
#'   The BMT-Psi distribution is the BMT distribution with \eqn{\kappa_l =
#'   \kappa_r = 0.63355781127887611515}. The BMT-Psi cumulative distribution
#'   function (cdf) is the closest BMT cdf to the logistic cdf with scale =
#'   1 / d and d = 1.70174439 (Camilli, 1994, p. 295).
#'
#' @param x,q vector of quantiles.
#' @param p vector of probabilities.
#' @param n number of observations. If \code{length(n) > 1}, the lenght is taken
#'   to be the number required
#' @param mean vector of means.
#' @param sd vector of standard deviations.
#' @param log,log.p logical; if TRUE, probabilities p are given as log(p).
#' @param lower.tail logical; if TRUE (default), probabilities are \eqn{P[X \le
#'   x]}, otherwise, \eqn{P[X > x]}.
#'
#' @return \code{dBMT.Psi} gives the density, \code{pBMT.Psi} the distribution
#'   function, \code{qBMT.Psi} the quantile function, and \code{rBMT.Psi}
#'   generates random deviates.
#'
#'   The length of the result is determined by \code{n} for \code{rBMT.Psi}, and
#'   is the maximum of the lengths of the numerical arguments for the other
#'   functions.
#'
#'   The numerical arguments other than \code{n} are recycled to the length of
#'   the result. Only the first elements of the logical arguments are used.
#'
#'   \code{sd <= 0} is an error and returns \code{NaN}.
#'
#' @references Torres-Jimenez, C. J. (2018), \emph{The BMT Item Response Theory
#'   model: A new skewed distribution family with bounded domain and an IRT
#'   model based on it}, PhD thesis, Doctorado en ciencias - Estadistica,
#'   Universidad Nacional de Colombia, Sede Bogota.
#'
#'   Camilli, G. (1994). Teacher's corner: origin of the scaling constant d= 1.7
#'   in item response theory. Journal of Educational Statistics, 19(3), 293-295.
#'
#' @seealso \link{Distributions} for other standard distributions.
#'   the BMT-Phi distribution.
#'
#' @author Camilo Jose Torres-Jimenez [aut,cre] \email{cjtorresj@unal.edu.co}
#'
#' @examples
#'
#' layout(matrix(1:4, 2, 2))
#'
#' curve(plogis(x, scale = 1 / 1.70174439), -4, 4, col = "red", lty = 2, ylab = "cdf")
#' curve(pBMT.Psi(x), add = TRUE, col = "blue", lty = 3)
#' legend("topleft", legend = c("logis(0, 1 / 1.70174439)","BMT-Psi(0,1)"),
#'        bty = "n", col = c("red","blue"), lty = 2:3)
#'
#' curve(plogis(x, scale = 1 / 1.70174439)-pBMT.Psi(x), -4, 4)
#'
#' curve(qlogis(x, scale = 1 / 1.70174439), col = "red", lty = 2, xlab = "p", ylab = "qf")
#' curve(qBMT.Psi(x), add = TRUE, col = "blue", lty = 3)
#'
#' hist(rBMT.Psi(10000), freq = FALSE, breaks = seq(-4, 4, 0.25), border = "blue")
#' curve(dlogis(x, scale = 1 / 1.70174439), add = TRUE, col = "red", lty = 2)
#' curve(dBMT.Psi(x), add = TRUE, col = "blue", lty = 3)

#' @rdname BMT.Psi
#' @export dBMT.Psi
dBMT.Psi <- function(x, mean = 0, sd = 1, log = FALSE) {
# The length of the result is determined by the maximum of the lengths of the
# numerical arguments. The numerical arguments are recycled to the length of
# the result.
len <- max(length(x), length(mean), length(sd))
x <- rep(x, len = len)
mean <- rep(mean, len = len)
sd <- rep(sd, len = len)
# Control location-scale parameters
sd <- replace(sd, sd <= 0, NaN)
# Transform x to 0,1 given location-scale parameters
range <- sd / .sd.BMT.Psi.DD01
x <- (x - mean) / range + .mean.BMT.Psi.DD01
# NaNs
y <- x + mean + sd
if (any(is.nan(y))) {
warning("NaNs founded or produced")
}
ind <- !is.na(y)
# For transformed x outside 0,1
y[ind] <- 0
# For transformed x inside 0,1
ind[ind] <- x[ind] > 0 & x[ind] < 1
if (any(ind)) {
# inv.x.t.Psi
t <- .inv.x.t.Psi(x[ind])
# yf.t.Psi
y[ind] <- .yf.t.Psi(t) / range[ind]
}
# density values y are given as log(y)
if (log)
y <- log(y)
return(y)
}

#' @rdname BMT.Psi
#' @export pBMT.Psi
pBMT.Psi <- function(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE){
# The length of the result is determined by the maximum of the lengths of the
# numerical arguments. The numerical arguments are recycled to the length of
# the result.
len <- max(length(q), length(mean), length(sd))
q <- rep(q, len=len)
mean <- rep(mean, len=len)
sd <- rep(sd, len=len)
# Control location-scale parameters
sd <- replace(sd, sd <= 0, NaN)
# Transform q to 0,1 given location-scale parameters
range <- sd / .sd.BMT.Psi.DD01
q <- (q - mean) / range + .mean.BMT.Psi.DD01
# NaNs
p <- q + mean + sd
if(any(is.nan(p)))
warning("NaNs founded or produced")
ind <- !is.na(p)
# For transformed q outside 0,1
p[ind] <- 0
ind2 <- ind
ind2[ind] <- q[ind] >= 1
p[ind2] <- 1
# For transformed q inside 0,1
ind[ind] <- q[ind] > 0 & q[ind] < 1
if(any(ind)){
# inv.x.t.Psi
t <- .inv.x.t.Psi(q[ind])
# yF.t.Psi
p[ind] <- .yF.t.Psi(t)
}
# probabilities are \eqn{P[X > x]}
if(!lower.tail)
p <- 1 - p
# probabilities p are given as log(p)
if(log.p)
p <- log(p)
return(p)
}

#' @rdname BMT.Psi
#' @export qBMT.Psi
qBMT.Psi <- function(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE){
# probabilities p are given as log(p)
if(log.p)
p <- exp(p)
# probabilities are \eqn{P[X > x]}
if(!lower.tail)
p <- 1 - p
# The length of the result is determined by the maximum of the lengths of the
# numerical arguments. The numerical arguments are recycled to the length of
# the result.
len <- max(length(p),length(mean),length(sd))
p <- rep(p, len=len)
mean <- rep(mean, len=len)
sd <- rep(sd, len=len)
# Control location-scale parameters
sd <- replace(sd, sd <= 0, NaN)
# NaNs
q <- p + mean + sd
if(any(is.nan(q)))
warning("NaNs founded or produced")
ind <- !is.na(q)
# For p outside (0,1)
q <- rep(NaN,len)
ind2 <- ind
ind2[ind] <- p[ind] == 0
q[ind2] <- 0
ind2 <- ind
ind2[ind] <- p[ind] == 1
q[ind2] <- 1
# For p inside (0,1)
ind[ind] <- p[ind] > 0 & p[ind] < 1
if(any(ind)){
# inv.yF.t
t <- .inv.yF.t.Psi(p[ind])
# x.t
q[ind] <- .x.t.Psi(t)
}
# Transform q
range <- sd / .sd.BMT.Psi.DD01
q <- (q - .mean.BMT.Psi.DD01) * range + mean
return(q)
}

#' @rdname BMT.Psi
#' @export rBMT.Psi
rBMT.Psi <- function(n, mean = 0, sd = 1){
#
len <- length(n)
if(len > 1)
n <- len
else
n <- trunc(n)
if(n < 1)
stop("invalid arguments")
# Method of inversion
p <- runif(n)
x <- qBMT.Psi(p, mean, sd)
return(x)
}

## Global constants
.epsilon <- 1e-10
.zero <- 0 - .epsilon
.one <- 1 + .epsilon
# X \sim BMT.Psi(mean, sd) is the same X \sim BMT(.D.Psi, .D.Psi, "t w", mean, sd, "l-s")
# .D.Psi <- BMTfit.mge(qlogis(1:1e6/(1e6+1), scale=1/1.70174439),
#                  "KS", start=list(p4=0.5), fix.arg=list(p1=0, p2=1, p3=0),
#                  type.p.3.4="a-s", type.p.1.2="l-s", custom.optim="nlminb")$estimate .D.Psi <- 0.63355781127887611515 # Mean of X \sim BMT(.D.Psi, .D.Psi, "t w", 0, 1, "c-d") # .mean.BMT.Psi.DD01 <- BMTmean(.D.Psi, .D.Psi) .mean.BMT.Psi.DD01 <- 0.5 # Standard deviation of X \sim BMT(.D.Psi, .D.Psi, "t w", 0, 1, "c-d") # .sd.BMT.Psi.DD01 <- BMTsd(.D.Psi, .D.Psi) .sd.BMT.Psi.DD01 <- 0.16771818811837588270 # Maximum of X \sim BMT.Psi(0, 1) # .d.BMT.Psi <- BMTchangepars(.D.Psi, .D.Psi, "t w", 0, 1, "l-s")$p2
.d.BMT.Psi <-         2.98119128050142556674
# Coefficients of polynomial x.t = a_3.Psi * t^3 + a_2.Psi * t^2 + a_1.Psi * t
# .a_3.Psi <- 6*.D.Psi - 2
.a_3.Psi <-           1.80134686767325646883
# .a_2.Psi <- -9*.D.Psi + 3
.a_2.Psi <-           -2.70202030150988470325
# .a_1.Psi <- 3*.D.Psi
.a_1.Psi <-           1.90067343383662823442
# For the real root of x.t = x
# .a.Psi <- .a_2.Psi/.a_3.Psi = -1.5
.a.Psi <-             -1.5
# .b.Psi <- .a_1.Psi/.a_3.Psi
.b.Psi <-             1.05514016647536013060
# .Q.Psi <- (.a.Psi*.a.Psi - 3*.b.Psi)/9
.Q.Psi <-             -0.10171338882512001578
# .auxR.Psi <- (2*.a.Psi*.a.Psi - 9*.b.Psi)*.a.Psi / 54
.auxR.Psi <-          0.13878504161884000490

# x.t.Psi
.x.t.Psi <- function(t){
x <- ((.a_3.Psi*t + .a_2.Psi)*t + .a_1.Psi)*t
return(x)
}

# Inverse function for x.t.Psi
.inv.x.t.Psi <- function(x){
# Press W.H., Teukolsky S.A., Vetterling W.T. & Flannery B.P. 2007.
# Numerical recipes: The art of scientific computing
# Section 5.6: Quadratic and Cubic Equations. Page 228.
len <- length(x)
c <- -x/.a_3.Psi
R <- .auxR.Psi + 0.5*c
# One real root
A <- - sign(R) * (abs(R) + sqrt(R*R-.Q.Psi*.Q.Psi*.Q.Psi))^(1/3)
B <- rep(0, len)
B[A!=0] <- .Q.Psi/A[A!=0]
r <- (A + B) - .a.Psi/3
r[r < 0] <- 0
r[r > 1] <- 1
return(r)
}

# yF.t.Psi (same yF.t)
.yF.t.Psi <- function(t){
yF <- (-2*t + 3)*t*t
return(yF)
}

# inv.yF.t.Psi (same inv.yF.t)
.inv.yF.t.Psi <- function(yF){
t <- 0.5-cos((acos(2*yF-1)-2*pi)/3)
return(t)
}

# yf.t.Psi
.yf.t.Psi <- function(t){
yf <- ((-6*t + 6)*t) / ((3*.a_3.Psi*t + 2*.a_2.Psi)*t + .a_1.Psi)
return(yf)
}


## Try the BMT package in your browser

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

BMT documentation built on May 2, 2019, 5:41 a.m.