R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("tweedie-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('tweedie') Loading required package: statmod > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "dtweedie" > > ### * dtweedie > > flush(stderr()); flush(stdout()) > > ### Name: dtweedie > ### Title: Tweedie Distributions > ### Aliases: dtweedie dtweedie.series dtweedie.inversion ptweedie > ### ptweedie.series ptweedie.inversion qtweedie rtweedie > ### Keywords: models > > ### ** Examples > > ### Plot a Tweedie density > power <- 2.5 > mu <- 1 > phi <- 1 > y <- seq(0, 10, length=100) > fy <- dtweedie( y=y, power=power, mu=mu, phi=phi) > plot(y, fy, type="l") > # Compare to the saddlepoint density > f.saddle <- dtweedie.saddle( y=y, power=power, mu=mu, phi=phi) > lines( y, f.saddle, col=2 ) > > ### A histogram of Tweedie random numbers > hist( rtweedie( 1000, power=1.2, mu=1, phi=1) ) > > ### An example of the multimodal feature of the Tweedie > ### family with power near 1 (from the first reference > ### listed above). > y <- seq(0.001,2,len=1000) > mu <- 1 > phi <- 0.1 > p <- 1.02 > f1 <- dtweedie(y,mu=mu,phi=phi,power=p) > plot(y, f1, type="l", xlab="y", ylab="Density") > p <- 1.05 > f2<- dtweedie(y,mu=mu,phi=phi,power=p) > lines(y,f2, col=2) > > > > cleanEx(); ..nameEx <- "dtweedie.saddle" > > ### * dtweedie.saddle > > flush(stderr()); flush(stdout()) > > ### Name: dtweedie.saddle > ### Title: Tweedie Distributions (saddlepoint approximation) > ### Aliases: dtweedie.saddle > ### Keywords: models > > ### ** Examples > > p <- 2.5 > mu <- 1 > phi <- 1 > y <- seq(0, 10, length=100) > fy <- dtweedie( y=y, power=p, mu=mu, phi=phi) > plot(y, fy, type="l") > # Compare to the saddlepoint density > f.saddle <- dtweedie.saddle( y=y, power=p, mu=mu, phi=phi) > lines( y, f.saddle, col=2 ) > > > > cleanEx(); ..nameEx <- "tweedie.profile" > > ### * tweedie.profile > > flush(stderr()); flush(stdout()) > > ### Name: tweedie.profile > ### Title: Tweedie Distributions: mle estimation of p > ### Aliases: tweedie.profile > ### Keywords: models > > ### ** Examples > > library(statmod) # Needed to use tweedie.profile > # Generate some fictitious data > test.data <- rgamma(n=200, scale=1, shape=1) > # The gamma is a Tweedie distribution with power=2; > # let's see if the profile plot shows this > out <- tweedie.profile( test.data ~1, p.vec=seq(1.7, 2.3, length=10), + do.plot=TRUE, method="interpolation", smooth=TRUE, do.ci=TRUE) --- This function may take some time to complete; Please be patient. If it fails, try using method="inversion" rather than the default method="series" Another possible reason for failure is the range of p: Try a different input for p.vec --- p= 1.7 * Fitting initial model: Done * Phi estimation, method: mle (using optimize): Done (phi.vec= 0.8485058 ) * Computing the log-likelihood (method= interpolation ): Done: L= -214.4607 p= 1.766667 * Fitting initial model: Done * Phi estimation, method: mle (using optimize): Done (phi.vec= 0.8898758 ) * Computing the log-likelihood (method= interpolation ): Done: L= -211.0377 p= 1.833333 * Fitting initial model: Done * Phi estimation, method: mle (using optimize): Done (phi.vec= 0.9331202 ) * Computing the log-likelihood (method= interpolation ): Done: L= -208.3912 p= 1.9 * Fitting initial model: Done * Phi estimation, method: mle (using optimize): Done (phi.vec= 0.979525 ) * Computing the log-likelihood (method= interpolation ): Done: L= -206.4052 p= 1.966667 * Fitting initial model: Done * Phi estimation, method: mle (using optimize): Done (phi.vec= 1.030321 ) * Computing the log-likelihood (method= interpolation ): Done: L= -205.0098 p= 2.033333 * Fitting initial model: Done * Phi estimation, method: mle (using optimize): Done (phi.vec= 1.08675 ) * Computing the log-likelihood (method= interpolation ): Done: L= -204.1632 p= 2.1 * Fitting initial model: Done * Phi estimation, method: mle (using optimize): Done (phi.vec= 1.150134 ) * Computing the log-likelihood (method= interpolation ): Done: L= -203.8419 p= 2.166667 * Fitting initial model: Done * Phi estimation, method: mle (using optimize): Done (phi.vec= 1.221961 ) * Computing the log-likelihood (method= interpolation ): Done: L= -204.0342 p= 2.233333 * Fitting initial model: Done * Phi estimation, method: mle (using optimize): Done (phi.vec= 1.303937 ) * Computing the log-likelihood (method= interpolation ): Done: L= -204.7363 p= 2.3 * Fitting initial model: Done * Phi estimation, method: mle (using optimize): Done (phi.vec= 1.398041 ) * Computing the log-likelihood (method= interpolation ): Done: L= -205.9503 --- * Smoothing: Done ML Estimates: p= 2.104082 with phi= 1.154266 giving L= -203.839 --- * Finding confidence interval: Done > # And plot the points to see how the smooth went > points( out$p, out$L) > > > > ### *