|
Readme file
SERIES C
 Applied Statistics
Modelling cell generation times by using the tempered stable distribution, by K. J. Palmer, M. S. Ridout and B. J. T. Morgan, pages 379–397
Programs for the TS distribution and instructions for use are also provided at http://www.kent.ac.uk/ims/personal/msr/TemperedStable.html
#================================================================================================
# The following commands provide examples of how to use the R functions contained in the file
# TSfunctions.txt. The functions provided are:
# dts(x, param)
# pts(x, param)
# qts(x, param, lower = 1/10000, upper = 10000)
# rts(n, param)
# These provide, respectively, the probability density function, the cumulative distribution
# function, the quantile function and a random number generator. For the quantile function,
# lower and upper should be set to values below and above the true quantile; the default values
# will suffice in most examples. The argument param is a vector of length 3, giving the mean,
# coefficient of variation, and further parameter alpha (0 < alpha < 1) of the distribution
#================================================================================================
# first need to source in the file TSfunctions.txt. May need to change the directory here.
# The line below assumes it is in the current directory.
source("TSfunctions.txt")
#-------------------------------------------------------------------------------------
# Examples of how to compute the pdf, cdf, work out quantiles and simulate. Note for
# the quantile function, we need to specify an interval to search within - an error
# message will be given if both these limits give values of the cdf either above or
# below the percentage point we are seeking. We recommend drawing the cdf of the
# relevant TS distribution to determine appropriate limits.
#--------------------------------------------------------------------------------------
par(mfrow=c(2,2))
param <- c(4,0.3,0.8) # set the parameters (mean,CV,alpha)
x <- seq(0,10,by=0.01) # set points at which to evaluate the pdf/cdf
f <- dTS(x,param) # this obtains the pdf ...
plot(x,f,type="l") # ... and plots it
F <- pTS(x,param) # this obtains the cdf ...
plot(x,F,type="l") # ... and plots it
p <- 0.95 # obtain t such that Pr(X<=t) = p,
t <- qTS(p,param) # search for t in (lower,upper); here these
# are left at their default values.
sims <- rTS(1000,param) # simulate some data from TS ...
f <- dTS(x,param) # ... then calculate true density ...
hist(sims,freq=FALSE); lines(x,f) # ... and plot simulations vs true pdf
#-----------------------------------------------------------------------------------
# Example of how to fit the TS distribution to a set of data by maximum likelihood.
# We use simulated data to demonstrate.
#-----------------------------------------------------------------------------------
param <- c(4,0.3,0.8) # set the TS parameters (mean,CV,alpha) which are then
data <- rTS(100,param) # used to simulate a set of data
param0 <- param # param0 is the set of starting values
st <- system.time(v <- optts(param0,data)) # for maximum likelihood, optts carries out ML
cat("\nMaximum likelihood estimation\n\n")
cat("True parameter values\n")
cat(" mu = ",param0[1],"\n")
cat(" CV = ",param0[2],"\n")
cat("alpha = ",param0[3],"\n")
cat("\nParameter estimates\n")
print(v$optparam)
cat("\nElapsed time for estimation = ",st[3],"\n")
f2 <- dTS(x,v$optparam[1,]) # this obtains the pdf ...
plot(x,f,type="l", main="True and estimated densities") # ... and plots it
lines(x,f2,type="l",lty="dashed",col="red") # ... and plots it
legend(6,0.35,lty=c("solid","dashed"), col=c("black","red"),
legend=c("True", "Estimated"), cex=0.75)
# v stores the information from maximum likelihood.
# - v$optparam gives the estimates and standard errors of the MLE's
# - v$est.se.cor.matrix gives a matrix, where the first column gives the MLE's
# [of log(mu), log(CV) and logit(alpha)=log(alpha/(1-alpha))], and the following
# three columns give a 3x3 matrix with the standard errors on the diagonal and
# correlations between tranformed parameters on the lower triangle (se's and cor's
# of the tranformed parameters)
# - v$gradient gives the value of the gradient at the MLE
# - v$hessian.eigenvalues gives the eigenvalues of the hessian at the MLE
# - v$ml gives the normal output obtained from running the optim function in R,
# type ?optim for an explanation (note this output is given in terms of the
# transformed parameters)
Karen Palmer
Institute of Mathematics, Statistics and Actuarial Science
University of Kent
Canterbury
Kent
CT2 7NF
UK
E-mail: kjp24@kent.ac.uk
|
Journals SERIES
A
Statistics
in Society
SERIES B
Statistical
Methodology
SERIES C
Applied
Statistics
SERIES D
The
Statistician
|
