# from http://www.statsci.org/s/invgauss.html
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)
}

# generate data
set.seed(5732)
test <- function(x)
        {(1e6*(x^11*(1-x)^6)+1e4*(x^3*(1-x)^10))+.1}
x <- (0:100)/100
mu <- test(x)
y <- rinvgauss(x,mu,1)

# fit the model and evaluate
ig.fit <- gssanova(y~x,family="inverse.gaussian")
est <- predict(ig.fit,data.frame(x=x),se=TRUE)

# plot the fit
plot(x,y,log="y",ylab=expression(mu(x)))
lines(x,mu,col=6)
lines(x,exp(est$fit),col=4)
lines(x,exp(est$fit+1.96*est$se),col=5)
lines(x,exp(est$fit-1.96*est$se),col=5)