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.

> ### * <HEADER>
> ###
> 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("Davies-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('Davies')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "Davies"
> 
> ### * Davies
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Davies
> ### Title: The Davies distribution
> ### Aliases: Davies ddavies pdavies qdavies rdavies ddavies.p
> ### Keywords: distribution
> 
> ### ** Examples
> 
> params <- c(10,0.1,-0.1)
> x <- seq(from=4,to=20,by=0.2)
> p <- seq(from=1e-3,to=1-1e-3,len=50)
> 
> rdavies(n=5,params)
[1]  9.032526  9.490331 10.297852 12.575867  8.714633
> least.squares(rdavies(100,params))
[1] 10.03596537  0.08224260 -0.09182977
> plot(pdavies(x,params))
> 
> plot(p,qdavies(p,params))
> plot(x,ddavies(x,params),type="b")
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "Gld"
> 
> ### * Gld
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Gld
> ### Title: The Generalized Lambda Distribution
> ### Aliases: Gld dgld dgld.p pgld qgld rgld
> ### Keywords: distribution
> 
> ### ** Examples
> 
> params <- c(4.114,0.1333,0.0193,0.1588)  #taken straight from some paper
> 
> gld.rv <- rgld(100,params)
> 
> hist(gld.rv)
> fit.davies.q(gld.rv)  #remember the Davies distn has 3 DF and the GLD 4...
> 
> 
> 
> cleanEx(); ..nameEx <- "davies.moment"
> 
> ### * davies.moment
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: davies.moment
> ### Title: Moments of the Davies distribution
> ### Aliases: davies.moment kurtosis skewness expected.value
> ###   expected.value.approx mu variance M
> ### Keywords: distribution
> 
> ### ** Examples
> 
> params <- c(10,0.1,-0.1)
> davies.moment(n=100,i=99,2,params) # ie the second moment of the 99th smallest
[1] 233.2944
>                             # observation of 100 drawn from a Davies
>                             # distribution with parameters p
> 
> mean(rdavies(1e6,params))-mu(params)
[1] -0.001530719
> 
> #now reproduce the S-K graph:
> 
> f <- function(x,y){c(skewness(c(1,x,y)),kurtosis(c(1,x,y)))}
> g <- function(j,vector,pp,qq=1){points(t(sapply(vector,f,y=j)),type="l",col="black",lty=qq)}
> 
> vector <- c((0:300)/100 , (0:300)/10000 , seq(from=3,to=10,len=100))
> vector <- sort(unique(vector))
> 
> plot(t(sapply((0:10)/10,f,y=0)),xlim=c(-3,3),ylim=c(0,10),type="n",xlab="skewness",ylab="kurtosis")
> g(-0.001,vector,"red",qq=1)
> g(-0.01,vector,"yellow",qq=2)
> g(-0.02,vector,"green",qq=3)
> g(-0.05,vector,"blue",qq=4)
> g(-0.1 ,vector,"purple",qq=5)
> g(-0.14,vector,"black",qq=6)
> 
> x <- seq(from=-3,to=3,len=30)
> points(x,x^2+1,type="l",lwd=2)
> 
> leg.txt <- expression(lambda[2]==-0.001,lambda[2]==-0.01,lambda[2]==-0.02,lambda[2]==-0.05,lambda[2]==-0.1,lambda[2]==-0.14)
> legend(-1.1,10,leg.txt,col="black",lty=1:6)
> 
> 
> 
> cleanEx(); ..nameEx <- "davies.start"
> 
> ### * davies.start
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: davies.start
> ### Title: start value for Davies minimization routines
> ### Aliases: davies.start
> ### Keywords: distribution
> 
> ### ** Examples
> 
> data <- rnorm(100)^2
> davies.start(data)
[1]  3.014035  3.094997 -0.010000
> least.squares(data)
[1]  2.1191420  3.3448700 -0.2026601
> 
> params <- c(10 , 0.1 , -0.1)
> x <- rdavies(100 , params)
> davies.start(x)
[1]  9.51454821  0.06579446 -0.12307374
> 
> f <- function(threeps){objective(davies.start(x,threeps),x)}
> 
> (jj<-optim(c(0.1,0.5,0.9),f))
$par
[1] 0.05050504 0.46694944 0.90704276

$value
[1] 7.883017

$counts
function gradient 
     451       NA 

$convergence
[1] 0

$message
NULL

> davies.start(x,jj$par)
[1]  9.41223146  0.06930978 -0.10272929
> least.squares(x)
[1]  9.43420781  0.07206733 -0.10207140
> 
> #not bad at all.
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "expected.gld"
> 
> ### * expected.gld
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: expected.gld
> ### Title: expected value of the Generalized Lambda Distribution
> ### Aliases: expected.gld expected.gld.approx
> ### Keywords: distribution
> 
> ### ** Examples
> 
> params <- c(4.114,0.1333,0.0193,0.1588)
> mean(rgld(1000,params))
[1] 5.008641
> expected.gld(n=1,i=1,params)
[1] 5
> expected.gld.approx(n=1,i=1,params)
[1] 4.796232
> 
> f <- function(n){apply(matrix(rgld(n+n,params),2,n),2,min)}
> #ie f(n) gives the smaller of 2 rgld RVs, n times.
> 
> mean(f(1000))
[1] 4.397669
> expected.gld(n=2,i=1,params)
[1] 4.453445
> expected.gld.approx(n=2,i=1,params)
[1] 4.424418
> 
> plot(1:100,expected.gld.approx(n=100,i=1:100,params)-expected.gld(n=100,i=1:100,params))
> # not bad, eh? ....yyyeeeeesss, but the parameters given by Ozturk give
> #an almost zero second derivative for d(qgld)/dp, so the good agreement
> #isn't surprising really.  Observe that the error is minimized at about
> #p=0.2, where the point of inflection is.
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "fit.davies.p"
> 
> ### * fit.davies.p
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fit.davies.p
> ### Title: Fits and plots Davies distributions to datasets
> ### Aliases: fit.davies.p fit.davies.q
> ### Keywords: distribution
> 
> ### ** Examples
> 
>   fit.davies.q(rnorm(100)^2)
>   fit.davies.p(exp(rnorm(100))) 
> 
>   data(x00m700p4)
>   fit.davies.q(x00m700p4)
> 
> 
> 
> cleanEx(); ..nameEx <- "hypergeo"
> 
> ### * hypergeo
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: hypergeo
> ### Title: The hypergeometric function
> ### Aliases: hypergeo
> ### Keywords: distribution
> 
> ### ** Examples
> 
> hypergeo(1,1,1,0.234)
[1] 1.305483
> hypergeo(1,1,1,(1:10)/10)
 [1]  1.111111  1.250000  1.428571  1.666667  2.000000  2.500000  3.333333
 [8]  5.000000 10.000000        NA
> 
> 
> 
> cleanEx(); ..nameEx <- "least.squares"
> 
> ### * least.squares
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: least.squares
> ### Title: Finds the optimal Davies distribution for a dataset
> ### Aliases: least.squares maximum.likelihood
> ### Keywords: distribution
> 
> ### ** Examples
> 
> p <- c(10 , 0.1 , -0.1)
> data <-rdavies(50,p)
> system.time(print(maximum.likelihood(data)))
[1] 10.40230021  0.10087827 -0.07581067
[1] 13.17  0.01 13.23  0.00  0.00
>                            #observe how long this takes.
>                            #The time is taken in repeated calls
>                            #to pdavies(), which uses uniroot().
> 
> system.time(print(least.squares(data)))
[1] 10.22134109  0.09564707 -0.08767713
[1] 0.07 0.00 0.07 0.00 0.00
>                            #Much faster.
> 
> 
> 
> cleanEx(); ..nameEx <- "likelihood"
> 
> ### * likelihood
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: likelihood
> ### Title: likelihood for the Davies distribution
> ### Aliases: likelihood neg.log.likelihood
> ### Keywords: distribution
> 
> ### ** Examples
> 
> p1 <- c(10,0.1,-0.1)
> p2 <- c(10,0.4,-0.1)
> data <- rdavies(100,p1)
> likelihood(p1,data)
[1] 3.193886e-81
> likelihood(p2,data)                 #should be smaller.
[1] 8.383199e-105
> neg.log.likelihood(p1,rstupid(100)) #should be large negative.
[1] Inf
> 
> 
> 
> cleanEx(); ..nameEx <- "normsquare"
> 
> ### * normsquare
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: normsquare
> ### Title: Square normal distribution
> ### Aliases: normsquare do.normsquare do.many.normsquare lognorm do.lognorm
> ###   do.many.lognorm gld.ozturk do.gld.ozturk do.many.gld.ozturk
> ### Keywords: distribution
> 
> ### ** Examples
> 
>  do.many.normsquare(c(20,30),4)
     vector                    
[1,]     20 0.41 0.76 0.56 0.39
[2,]     30 0.06 0.73 0.80 0.68
> 
> #carries out resampling on a
> #square-normal random dataset of size 20 and 30, four times.
> #Takes a long time.
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "objective"
> 
> ### * objective
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: objective
> ### Title: The objective function for fitting the Davies distribution
> ### Aliases: objective objective.approx
> ### Keywords: distribution
> 
> ### ** Examples
> 
> params <- c(10,0.1,-0.1)
> x <- rdavies(100,params)
> objective(params,x)
[1] 11.75608
> objective.approx(params,x)
[1] 7.672145
> 
> objective(least.squares(x),x)
[1] 1.864535
> objective(davies.start(x),x)
[1] 2.912496
> 
>   
> 
> 
> cleanEx(); ..nameEx <- "plotcf"
> 
> ### * plotcf
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plotcf
> ### Title: p-value investigation
> ### Aliases: plotcf
> ### Keywords: distribution
> 
> ### ** Examples
> 
> f.H0.T <- function(n,free=5){t.test(rt(n,df=free))$p.value}
> f.H0.F <- function(n,free=5){t.test(rf(n,df1=free,df2=free))$p.value}
> 
> plotcf(sapply(rep(10,100),f.H0.T))  # should reject about 5: thus
>                                      # probability of a type I error is
>                                      # about 0.05 (as it should be; this
>                                      # is an exact test
> plotcf(sapply(rep(10,100),f.H0.F))  # should reject about 80: thus
>                                      # probability of a type II error is
>                                      # about 0.2 for this H_A.
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "resample.ls"
> 
> ### * resample.ls
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: resample.ls
> ### Title: Parametric resampling from a Davies distribution
> ### Aliases: resample.ls resample.ml
> ### Keywords: distribution
> 
> ### ** Examples
> 
> hist(resample.ls(rnorm(50)^2)$errors)
> fit.davies.q(resample.ls(rnorm(80)^2)$errors)
> 
> 
> 
> cleanEx(); ..nameEx <- "rstupid"
> 
> ### * rstupid
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: rstupid
> ### Title: A stupid PDF
> ### Aliases: rstupid
> ### Keywords: distribution
> 
> ### ** Examples
> 
> stupid <-rstupid(500)
> fit.davies.q(stupid)
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "twolines.vert"
> 
> ### * twolines.vert
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: twolines.vert
> ### Title: Order statistic comparison
> ### Aliases: twolines.vert
> ### Keywords: distribution
> 
> ### ** Examples
> 
> twolines.vert(1:100,sort(rnorm(100)),sort(rnorm(100)))
> params <- c(10 , 0.1 , -0.1)
> twolines.vert(1:100 , sort(rdavies(100,params)) , sort(rdavies(100,params)))
> 
> 
> 
> cleanEx(); ..nameEx <- "x00m700p4"
> 
> ### * x00m700p4
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: x00m700p4
> ### Title: Peak concentration for 100 instantaneous releases
> ### Aliases: x00m700p4
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(x00m700p4)
> fit.davies.q(x00m700p4)
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  119.08 0.38 120.42 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')