dinvgauss: The Inverse Gaussian Distribution In STAR: Spike Train Analysis with R

Description

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

Usage

 ```1 2 3 4 5 6 7 8``` ```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]`, otherwise, `P[X > x]`. `log, log.p` logical; if `TRUE`, probabilities p are given as log(p).

Details

With the default, `"sigma2"`, parameterization (```mu = m, sigma2 = s^2```) the inverse Gaussian distribution has density:

f(x) = 1/sqrt(2*pi*x^3*s^2) * exp(-0.5*(x-mu)^2/(x*s^2*m^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) = (b/sqrt(2*pi*x^3)) * exp(-0.5*(x-b*mu)^2/(x*m^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).

Value

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

Author(s)

Christophe Pouzat [email protected]

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.

`invgaussMLE`, `Lognormal`, `hinvgauss`

Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71``` ```## 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-3*sampIGmleIG\$se, sampIGmleIG\$estimate+3*sampIGmleIG\$se, sampIGmleIG\$se/10) Sigma2 <- seq(sampIGmleIG\$estimate-7*sampIGmleIG\$se, sampIGmleIG\$estimate+7*sampIGmleIG\$se, sampIGmleIG\$se/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,sampIG2mleIG\$estimate,lower.tail=FALSE),col=2) ## End(Not run) ```

Example output      ```Loading required package: survival
This is mgcv 1.8-20. For overview type 'help("mgcv-package")'.
mu    sigma2
est  0.082531749 3.0257187
se   0.004124258 0.4279012
true 0.075000000 3.0000000
 0.8697612
Warning message:
In rep(yi[si > 0], each = si[si > 0]) :
first element used of 'each' argument
mu    sigma2
est  0.084037995 2.8688045
se   0.005374312 0.4903304
true 0.075000000 3.0000000
 0.8507114
```

STAR documentation built on May 2, 2019, 4:53 p.m.