## 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)
Run the code above in your browser using DataLab