STAR (version 0.3-7)

dinvgauss: The Inverse Gaussian Distribution

Description

Density, distribution function, quantile function, and random generation for the inverse Gaussian.

Usage

dinvgauss(x, mu = 1, sigma2 = 1, boundary = NULL, log = FALSE)
pinvgauss(q, mu = 1, sigma2 = 1, boundary = NULL, lower.tail = TRUE, log.p = FALSE)
qinvgauss(p, mu = 1, sigma2 = 1, boundary = NULL)
rinvgauss(n = 1, mu = 1, sigma2 = 1, boundary = NULL)

Arguments

x, q

vector of quantiles.

p

vector of probabilities.

n

number of observations. If length(n) > 1, the length is taken to be the number required.

mu

mean value of the distribution in the default parameterization, mean value / boundary otherwise. Can also be viewed as the inverse of the drift of the latent Brownian motion.

sigma2

variance of the latent Brownian motion. When this parameterization is used (the default) the distance between the "starting" point and the boundary ("absorbing barrier") is set to 1.

boundary

distance between the starting point and the "absorbing barrier" of the latent Brownian motion. When this parameterization is used the Brownian motion variance is set to 1.

lower.tail

logical; if TRUE (default), probabilities are P[X <= x]<="" code="">, otherwise, P[X > x].

log, log.p

logical; if TRUE, probabilities p are given as log(p).

Value

dinvgauss gives the density, pinvgauss gives the distribution function, qinvgauss gives the quantile function and rinvgauss generates random deviates.

Details

With the default, "sigma2", parameterization (mu = m, sigma2 = s^2) the inverse Gaussian distribution has density: $$% f(x)=\frac{1}{\sqrt{2 \, \pi \, \sigma^2 \, x^3}} \, \exp% (-\frac{1}{2}\frac{(x-\mu)^2}{x \, \sigma^2 \, \mu^2}) $$ with $s^2 > 0$. The theoretical mean is: $m$ and the theoretical variance is: $m^3*s^2$. With the default, "boundary", parameterization (mu = m, boundary = b)the inverse Gaussian distribution has density: $$% f(x)=\frac{b}{\sqrt{2 \, \pi \, x^3}} \, \exp% (-\frac{1}{2}\frac{(x-b \, \mu)^2}{x \, \mu^2}) $$ with $s^2 > 0$. The theoretical mean is: $m * b$ and the theoretical variance is: $m^3*b$. The latent Brownian motion is described in Lindsey (2004) pp 209-213, Whitemore and Seshadri (1987), Aalen and Gjessing (2001) and Gerstein and Mandelbrot (1964).

The expression for the distribution function is given in Eq. 4 of Whitemore and Seshadri (1987).

Initial guesses for the inversion of the distribution function used in qinvgauss are obtained with the transformation of Whitemore and Yalovsky (1978). Random variates are obtained with the method of Michael et al (1976) which is also described by Devroye (1986, p 148) and Gentle (2003, p 193).

References

Gerstein, George L. and Mandelbrot, Benoit (1964) Random Walk Models for the Spike Activity of a Single Neuron. Biophys J. 4: 41--68. http://www.pubmedcentral.nih.gov/articlerender.fcgi?tool=pubmed&pubmedid=14104072.

Whitemore, G. A. and Yalovsky, M. (1978) A normalizing logarithmic transformation for inverse Gaussian random variables. Technometrics 20: 207--208. Whitmore, G. A. and Seshadri, V. (1987) A Heuristic Derivation of the Inverse Gaussian Distribution. The American Statistician 41: 280--281.

Aalen, Odd O. and Gjessing, Hakon K. (2001) Understanding the Shape of the Hazard Rate: A Process Point of View. Statistical Science 16: 1--14.

Lindsey, J.K. (2004) Introduction to Applied Statistics: A Modelling Approach. OUP.

Michael, J. R., Schucany, W. R. and Haas, R. W. (1976) Generating random variates using transformations with multiple roots. The American Statistician 30: 88--90.

Devroye, L. (1986) Non-Uniform Random Variate Generation. Springer-Verlag. http://cg.scs.carleton.ca/~luc/rnbookindex.html.

Gentle, J. E. (2003) Random Number Generation and Monte Carlo Methods. Springer.

See Also

invgaussMLE, Lognormal, hinvgauss

Examples

## Not run: 
# ## Start with the inverse Gauss
# ## Define standard mu and sigma
# mu.true <- 0.075 ## a mean ISI of 75 ms
# sigma2.true <- 3
# ## Define a sequence of points on the time axis
# X <- seq(0.001,0.3,0.001)
# ## look at the density
# plot(X,dinvgauss(X,mu.true,sigma2.true),type="l",xlab="ISI (s)",ylab="Density")
# 
# ## Generate a sample of 100 ISI from this distribution
# sampleSize <- 100
# sampIG <- rinvgauss(sampleSize,mu=mu.true,sigma2=sigma2.true)
# ## check out the empirical survival function (obtained with the Kaplan-Meyer
# ## estimator) against the true one
# library(survival)
# sampIG.KMfit <- survfit(Surv(sampIG,1+numeric(length(sampIG))) ~1)
# plot(sampIG.KMfit,log=TRUE)
# lines(X,pinvgauss(X,mu.true,sigma2.true,lower.tail=FALSE),col=2)
# 
# ## Get a ML fit
# sampIGmleIG <- invgaussMLE(sampIG)
# ## compare true and estimated parameters
# rbind(est = sampIGmleIG$estimate,se = sampIGmleIG$se,true = c(mu.true,sigma2.true))
# ## plot contours of the log relative likelihood function
# Mu <- seq(sampIGmleIG$estimate[1]-3*sampIGmleIG$se[1],
#           sampIGmleIG$estimate[1]+3*sampIGmleIG$se[1],
#           sampIGmleIG$se[1]/10)
# Sigma2 <- seq(sampIGmleIG$estimate[2]-7*sampIGmleIG$se[2],
#               sampIGmleIG$estimate[2]+7*sampIGmleIG$se[2],
#               sampIGmleIG$se[2]/10)
# sampIGmleIGcontour <- sapply(Mu, function(mu) sapply(Sigma2, function(s2) sampIGmleIG$r(mu,s2))) 
# contour(Mu,Sigma2,t(sampIGmleIGcontour),
#         levels=c(log(c(0.5,0.1)),-0.5*qchisq(c(0.95,0.99),df=2)),
#         labels=c("log(0.5)",
#           "log(0.1)",
#           "-1/2*P(Chi2=0.95)",
#           "-1/2*P(Chi2=0.99)"),
#         xlab=expression(mu),ylab=expression(sigma^2))
# points(mu.true,sigma2.true,pch=16,col=2)
# ## We can see that the contours are more parabola like on a log scale
# contour(log(Mu),log(Sigma2),t(sampIGmleIGcontour),
#         levels=c(log(c(0.5,0.1)),-0.5*qchisq(c(0.95,0.99),df=2)),
#         labels=c("log(0.5)",
#           "log(0.1)",
#           "-1/2*P(Chi2=0.95)",
#           "-1/2*P(Chi2=0.99)"),
#         xlab=expression(log(mu)),ylab=expression(log(sigma^2)))
# points(log(mu.true),log(sigma2.true),pch=16,col=2)
# ## make a deviance test for the true parameters
# pchisq(-2*sampIGmleIG$r(mu.true,sigma2.true),df=2)
# ## check fit with a QQ plot
# qqDuration(sampIGmleIG,log="xy")
# 
# ## Generate a censored sample using an exponential distribution
# sampEXP <- rexp(sampleSize,1/(2*mu.true))
# sampIGtime <- pmin(sampIG,sampEXP)
# sampIGstatus <- as.numeric(sampIG <= sampEXP)
# ## fit the censored sample
# sampIG2mleIG <- invgaussMLE(sampIGtime,,sampIGstatus)
# ## look at the results
# rbind(est = sampIG2mleIG$estimate,
#       se = sampIG2mleIG$se,
#       true = c(mu.true,sigma2.true))
# pchisq(-2*sampIG2mleIG$r(mu.true,sigma2.true),df=2)
# ## repeat the survival function estimation
# sampIG2.KMfit <- survfit(Surv(sampIGtime,sampIGstatus) ~1)
# plot(sampIG2.KMfit,log=TRUE)
# lines(X,pinvgauss(X,sampIG2mleIG$estimate[1],sampIG2mleIG$estimate[2],lower.tail=FALSE),col=2)
# ## End(Not run)