R/invgauss.R

In cwhmisc: Miscellaneous Functions for Math, Plotting, Printing, Statistics, Strings, and Tools

dinvgauss <- function(x, mu = stop("no shape arg"), lambda = 1) {
#  Density of inverse Gaussian distribution  GKS  15 Jan 98
  if(any(mu <= 0)) stop("mu must be positive")
  if(any(lambda <= 0)) stop("lambda must be positive")
  d <- ifelse(x > 0, sqrt(lambda/(2 * pi * x ^ 3))*exp(-lambda * (x - mu) ^ 2 /
             (2 * mu ^ 2 * x)), 0)
  if(!is.null(Names <- names(x)))
         names(d) <- rep(Names, length = length(d))
  return(d)
}

pinvgauss <- function(q, mu = stop("no shape arg"), lambda = 1) {
#  Inverse Gaussian distribution function  GKS  15 Jan 98
  if(any(mu <= 0)) stop("mu must be positive")
  if(any(lambda <= 0)) stop("lambda must be positive")
  n <- length(q)
  if(length(mu) > 1 && length(mu) != n) mu <- rep(mu, length = n)
  if(length(lambda) > 1 && length(lambda) != n) lambda <- rep(lambda, length = n)
  lq <- sqrt(lambda / q)
  qm <- q / mu
  p <- ifelse(q > 0, pnorm(lq * (qm - 1)) + exp(2 * lambda / mu) *
            pnorm(-lq * (qm + 1)), 0)
  if(!is.null(Names <- names(q)))
            names(p) <- rep(Names, length = length(p))
   return(p)
}

rinvgauss <- function(n, mu = stop("no shape arg"), lambda = 1) {
#  Random variates from inverse Gaussian distribution
#  Reference:
#      Chhikara and Folks, The Inverse Gaussian Distribution,
#      Marcel Dekker, 1989, page 53.
#  GKS  15 Jan 98
#
  if(any(mu <= 0)) stop("mu must be positive")
  if(any(lambda <= 0)) stop("lambda must be positive")
  if(length(n) > 1) n <- length(n)
  if(length(mu) > 1 && length(mu) != n) mu <- rep(mu, length=n)
  if(length(lambda) > 1 && length(lambda) != n)
       lambda <- rep(lambda,length = n)
  y2 <- rchisq(n, 1)
  u <- runif(n)
  r1 <- mu / (2 * lambda) * (2 * lambda + mu * y2 - sqrt(4 * lambda *
         mu * y2 + mu ^ 2 * y2 ^ 2))
  r2 <- mu ^ 2 / r1
  ifelse(u < mu/(mu + r1), r1, r2)
}

qinvgauss  <- function(p, mu = stop("no mean arg"), lambda = 1) {
#  Quantiles of the inverse Gaussian distribution
#  Dr Paul Bagshaw
#  Centre National d'Etudes des Telecommunications (DIH/DIPS)
#  Technopole Anticipa, France
#  paul.bagshaw@cnet.francetelecom.fr 23 Dec 98
  if(any(mu <= 0.)) stop("mu must be positive")
  if(any(lambda <= 0.)) stop("lambda must be positive")
  n <- length(p)
  if(length(mu) > 1 && length(mu) != n)
    mu <- rep(mu, length = n)
  if(length(lambda) > 1 && length(lambda) != n)
    lambda <- rep(lambda, length = n)
  thi <- lambda / mu
  U <- qnorm (p)
  r1 <- 1 + U / sqrt (thi) + U ^ 2 / (2 * thi) + U ^ 3 / (8 * thi * sqrt(thi))
  x <- r1
  for (i in 1:10) {
    cum <- pinvgauss (x, 1., thi)
    dx <- (cum - p) / dinvgauss (x, 1., thi)
    dx <- ifelse (is.finite(dx), dx, ifelse (p > cum, -1, 1))
    dx[dx < -1] <- -1
    if (all(dx == 0.)) break
    x <- x - dx
  }
  return(x * mu)
}