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)) 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)) 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 } x * mu }