Nothing
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General
# Public License along with this library; if not, write to the
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
# MA 02111-1307 USA
################################################################################
# FUNCTION:
# .pgld
# .qgld
# .qgl.fmkl
# .qgl.fm5
# .qgl.rs
# .qdgld
# .qdgl.rs
# .qdgl.fmkl
# .qdgl.fm5
# .rgld
# .gl.parameter.tidy
# .gl.check.lambda
################################################################################
# Code borrowed from
# R's contributed package "gld" written by Robert King
# Rmetrics:
# Note that gmm is not available on Debian as of 2009-11-11.
# To run these functions under Debian/Rmetrics we have them
# implemented here as a builtin.
# Package: gld
# Version: 1.8.4
# Date: 2008/10/01
# Title: Estimation and use of the generalised (Tukey) lambda distribution
# Author: Robert King <Robert.King@newcastle.edu.au>
# Maintainer: Robert King <Robert.King@newcastle.edu.au>
# Description: The generalised lambda distribution, or Tukey lambda
# distribution, provides a wide variety of shapes with one
# functional form. This package provides random numbers,
# quantiles, probabilities, densities and plots. It also
# includes an implementation of the starship estimation method
# for the distribution.
# License: GPL (>= 2)
# URL: http://tolstoy.newcastle.edu.au/~rking/publ/rprogs/information.html
# Packaged: 2009-10-14
################################################################################
.dgld <-
function(x, lambda1=0, lambda2=NULL, lambda3=NULL, lambda4=NULL, param="fmkl",
lambda5=NULL, inverse.eps=1e-8, max.iterations=500)
{
# Tidy the parameters so gl.check.lambda will work
lambdas <-
.gl.parameter.tidy(lambda1, lambda2, lambda3, lambda4, param, lambda5)
# Check the parameters
if(!.gl.check.lambda(lambdas, param=param, vect=TRUE)) {
stop(paste("The parameter values", lambdas,
"\ndo not produce a proper distribution with the", param,
"parameterisation - see \ndocumentation for .gl.check.lambda"))
}
# calculate u=F(x) numerically, then use qdgl
# Unless x is outside the range, then density should be zero
extreme <- .qgld(c(0, 1), lambda1=lambdas, param=param)
# It may be better to change this to simply
# (x <= extreme[2])*(x >= extreme[1])
outside.range <- !as.logical((x<=extreme[2])*(x>=extreme[1]))
u <- .pgld(x, lambdas, param=param, inverse.eps=inverse.eps,
max.iterations=max.iterations)
dens <- .qdgld(u, lambda1=lambdas, param=param)
dens[outside.range] <- 0
dens
}
################################################################################
.pgld <-
function(q, lambda1=0, lambda2=NULL, lambda3=NULL, lambda4=NULL,
param="fmkl", lambda5=NULL, inverse.eps=1e-8, max.iterations=500)
{
# Thanks to Steve Su, <s.su@qut.edu.au>, for improvements to this code
# If lambda1 is a vector, the default value for lambda2 will cause a
# problem.
# I did have a warning about this, but it will occur too often to make
# up for the benefit, so I've deleted it.
# Tidy the parameters so gl.check.lambda will work
lambdas <-
.gl.parameter.tidy(lambda1, lambda2, lambda3, lambda4, param, lambda5)
# Check the parameters
if(!.gl.check.lambda(lambdas, param=param, vect=TRUE)) {
stop(paste("The parameter values", lambda1, lambda2, lambda3, lambda4,
"\ndo not produce a proper distribution with the", param,
"parameterisation - see \ndocumentation for gl.check.lambda"))
}
jr <- q; jr[sort.list(q)] <- seq(along=q)
order.x<-order(q)
xx<-sort(q)
# xx has the sorted data, and jr & order.x the information to get back
# to theoriginal order.
extreme <- .qgld(c(inverse.eps, 1-inverse.eps), lambda1=lambdas,
param=param)
max.e<-extreme[2]
min.e<-extreme[1]
ind.min<-xx<=min.e
ind.max<-xx>=max.e
# This simpler comparison works here because we are using inverse.eps as our
# tolerance
q<-xx[as.logical((xx<max.e)*(xx>min.e))]
# We only want to calculate the probabilities for q values inside
# the support
length.of.vector <- length(q)
# Need a blank u to send to C
u <- 0*q
result <- switch(param,
freimer=, # allows for alternate expressions
frm=, # allows for alternate expressions
FMKL=, # Notes on .C call - the "numerics", lambdas and inverse.eps
# don't need the as.??? call as they are implicitly double
fmkl=.C(
C_gl_fmkl_distfunc, lambdas[1], lambdas[2], lambdas[3], lambdas[4],
as.double(0), as.double(1), inverse.eps,
as.integer(max.iterations), as.double(q), as.double(u),
as.integer(length.of.vector)),
ramberg=, # Ramberg & Schmeiser
ram=,
RS=,
rs=.C(
C_gl_rs_distfunc, lambdas[1], lambdas[2], lambdas[3], lambdas[4],
as.double(0), as.double(1), inverse.eps,
as.integer(max.iterations), as.double(q), as.double(u),
as.integer(length.of.vector)),
fm5=.C(C_gl_fm5_distfunc, lambdas[1], lambdas[2], lambdas[3],
lambdas[4], lambdas[5],
as.double(0), as.double(1), inverse.eps,
as.integer(max.iterations), as.double(q), as.double(u),
as.integer(length.of.vector)),
stop("Error: Parameterisation must be one of fmkl, fm5 or rs")
) # closes "switch"
if (!(is.numeric(result[[1]]))){
stop("Values for quantiles outside range. This shouldn't happen")
}
if (param=="fm5") {
u <- result[[11]]
} else {
u <- result[[10]]
}
xx[as.logical((xx<max.e)*(xx>min.e))]<-u
xx[ind.min]<-0
xx[ind.max]<-1
# Revert to the original order of the dataset:
xx[jr]
}
################################################################################
.qgld <-
function(p, lambda1, lambda2=NULL, lambda3=NULL, lambda4=NULL,
param="fmkl", lambda5=NULL)
{
lambdas <-
.gl.parameter.tidy(lambda1, lambda2, lambda3, lambda4, param, lambda5)
# Check the values are OK
if(!.gl.check.lambda(lambdas, param=param, vect=TRUE)) {
stop(paste("The parameter values", paste(lambdas, collapse=" "),
"\ndo not produce a proper distribution for the", param,
"parameterisation \n - see documentation for gl.check.lambda"))
}
result <- switch(param,
freimer=, # allows for alternate expressions
frm=, # allows for alternate expressions
FMKL=,
fmkl=.qgl.fmkl(p, lambdas),
ramberg=, # Ramberg & Schmeiser
ram=,
RS=,
rs=.qgl.rs(p, lambdas),
fm5 = .qgl.fm5(p, lambdas),
stop("Error: Parameterisation must be fmkl, fm5 or rs")
) # closes "switch"
result
}
# ------------------------------------------------------------------------------
.qgl.fmkl <-
function(p, lambdas)
{
# No checking - use qgl if you want that
lambda4 = lambdas[4]
lambda3 = lambdas[3]
lambda2 = lambdas[2]
lambda1 = lambdas[1]
p <- as.double(p)
# abandoned this for the simpler one below
# outside.range <- !as.logical(((p<1)*(p>0))|
# (sapply(p, all.equal, 1)=="TRUE")| (sapply(p, all.equal, 0)=="TRUE"))
outside.range <- !as.logical((p<=1)*(p>=0))
# u gets only the probabilities in [0, 1]
u <- p[!outside.range]
# If OK, determine special cases
if (lambda3 == 0) {
if (lambda4 == 0) { # both log
quants <- lambda1 + (log(u) - log(1 - u))/lambda2
}
else { # l3 zero, l4 non-zero
quants <- lambda1 +
(log(u) - ((1 - u)^lambda4 - 1)/lambda4)/lambda2
}
} else { # lambda3 non-zero
if (lambda4 == 0) { # non-zero, l4 zero
quants <- lambda1 +
((u^lambda3 - 1)/lambda3 - log(1 - u))/lambda2
} else { # both non-zero - use usual formula
quants <- lambda1 + ( ( u ^ lambda3 - 1) / lambda3 -
( (1-u)^lambda4 - 1) /lambda4 ) / lambda2
}
}
# Now we have the quantiles for p values inside [0, 1], put them in
# the right place in the result vector
result <- p*NaN
result[!outside.range] <- quants
# The remaining "quantiles" are NaN
result
}
# ------------------------------------------------------------------------------
.qgl.fm5 <-
function(p, lambdas)
{
# No parameter value checking. If you want that, use qgl!
lambda5 = as.double(lambdas[5])
lambda4 = as.double(lambdas[4])
lambda3 = as.double(lambdas[3])
lambda2 = as.double(lambdas[2])
lambda1 = as.double(lambdas[1])
p <- as.double(p)
# abandoned this for the simpler
# outside.range <- !as.logical(((p<1)*(p>0))|
# (sapply(p, all.equal, 1)=="TRUE")| (sapply(p, all.equal, 0)=="TRUE"))
outside.range <- !as.logical((p<=1)*(p>=0))
# u gets only the probabilities in [0, 1]
u <- p[!outside.range]
# If OK, determine special cases
if (lambda3 == 0) {
if (lambda4 == 0) { # both log
quants <- lambda1 + ((1-lambda5)*log(u) -
(1+lambda5)*log(1 - u))/lambda2
} else { # l3 zero, l4 non-zero
quants <- lambda1 +
((1-lambda5)*log(u) -
(1+lambda5)*((1 - u)^lambda4 - 1)/lambda4)/lambda2
}
} else { # lambda3 non-zero
if (lambda4 == 0) { # non-zero, l4 zero
quants <- lambda1 +
((1-lambda5)*(u^lambda3 - 1)/lambda3 -
(1+lambda5)*log(1 - u))/lambda2
} else { # both non-zero - use usual formula
quants <- lambda1 + ((1-lambda5)* ( u ^ lambda3 - 1) / lambda3
- (1+lambda5)*( (1-u)^lambda4 - 1) /lambda4 ) / lambda2
}
}
# Now we have the quantiles for p values inside [0, 1], put them in
# the right place in the result vector
result <- p*NaN
result[!outside.range] <- quants
# The remaining "quantiles" are NaN
result
}
# ------------------------------------------------------------------------------
.qgl.rs <-
function(p, lambdas)
{
u <- p
# No parameter value checking - use qgl!
lambda4 = lambdas[4]
lambda3 = lambdas[3]
lambda2 = lambdas[2]
lambda1 = lambdas[1]
# At present, I'm rejecting zero values for l3 and l4, though I think
# there are limit results, so one functional form.
quants <- lambda1 + ( u ^ lambda3 - (1-u)^lambda4 ) / lambda2
quants
}
# ------------------------------------------------------------------------------
.qdgld <-
function(p, lambda1, lambda2=NULL, lambda3=NULL, lambda4=NULL,
param="fmkl", lambda5=NULL)
{
# Don't allow characters in lambda5 -
# common error with parameterisation stuff
if(is.character(lambda5)) {
stop(paste("lambda5=", lambda5,
"It should be a number between -1 and 1"))
}
lambdas <-
.gl.parameter.tidy(lambda1, lambda2, lambda3, lambda4, param, lambda5)
# Check the values are OK
if(!.gl.check.lambda(lambdas, param=param, vect=TRUE)) {
stop(paste("The parameter values", paste(lambdas, collapse=" "),
"\ndo not produce a proper distribution for the", param,
"parameterisation \n - see documentation for gl.check.lambda"))
}
result <- switch(param,
# Different tests apply for each parameterisation
freimer=, # allows for alternate expressions
frm=, # allows for alternate expressions
FMKL=,
fmkl=.qdgl.fmkl(p, lambdas),
ramberg=, # Ramberg & Schmeiser
ram=,
RS=,
rs=.qdgl.rs(p, lambdas),
fm5 = .qdgl.fm5(p, lambdas),
stop("Error: Parameterisation must be fmkl, fm5 or rs")
) # closes "switch"
result
}
# ------------------------------------------------------------------------------
.qdgl.rs <-
function(p, lambdas)
{
# Check the values are OK)
if(!.gl.check.lambda(lambdas, param="rs", vect=TRUE)) {
stop(paste("The parameter values", paste(lambdas, collapse=" "),
"\ndo not produce a proper distribution with the RS parameterisation",
" - see \ndocumentation for gl.check.lambda"))
}
outside.range <- !as.logical((p<=1)*(p>=0))
# u gets only the probabilities in [0, 1]
u <- p[!outside.range]
dens <- lambdas[2]/(lambdas[3] * (p^(lambdas[3] -1)) +
lambdas[4] * ((1 - p)^(lambdas[4] -1)))
dens
}
# ------------------------------------------------------------------------------
.qdgl.fmkl <-
function(p, lambdas)
{
# Check the values are OK)
if(!.gl.check.lambda(lambdas, param="fmkl", vect=TRUE)) {
stop(paste("The parameter values", paste(lambdas, collapse=" "),
"\ndo not produce a proper distribution with the FMKL",
"parameterisation - see \ndocumentation for gl.check.lambda"))
}
outside.range <- !as.logical((p<=1)*(p>=0))
# u gets only the probabilities in [0, 1]
u <- p[!outside.range]
# The density is given by 1/Q'(u)
dens <- lambdas[2]/(p^(lambdas[3] - 1) + (1 - p)^(lambdas[4] - 1))
dens
}
# ------------------------------------------------------------------------------
.qdgl.fm5 <-
function(p, lambdas)
{
# Check the values are OK)
if(!.gl.check.lambda(lambdas, param="fm5", vect=TRUE)) {
stop(paste("The parameter values", paste(lambdas, collapse=" "),
"\ndo not produce a proper distribution with the FM5",
"parameterisation - see \ndocumentation for gl.check.lambda"))
}
outside.range <- !as.logical((p<=1)*(p>=0))
# u gets only the probabilities in [0, 1]
u <- p[!outside.range]
# The density is given by 1/Q'(u)
dens <- lambdas[2]/((1-lambdas[5])*(u^(lambdas[3] - 1)) +
(1+lambdas[5])*((1 - u)^(lambdas[4] - 1)) )
dens
}
################################################################################
.rgld <-
function(n, lambda1=0, lambda2=NULL, lambda3=NULL, lambda4=NULL,
param="fmkl", lambda5=NULL)
{
# Check the parameters
lambdas <-
.gl.parameter.tidy(lambda1, lambda2, lambda3, lambda4, param, lambda5)
# Check the values are OK
if(!.gl.check.lambda(lambdas, param=param, vect=TRUE)) {
stop(paste("The parameter values", lambdas,
"\ndo not produce a proper distribution for the", param,
"parameterisation \n - see documentation for gl.check.lambda"))
}
# Produce the uniform data
p <- runif(n)
# convert to gl
res <- .qgld(p, lambda1=lambdas, param=param)
res
}
################################################################################
.gl.parameter.tidy <-
function(lambda1, lambda2=NULL, lambda3=NULL, lambda4=NULL, param="fmkl",
lambda5=NULL)
{
# Don't allow characters in lambda5 - common error with parameterisation stuff
if(is.character(lambda5)) {stop(paste("lambda5=", lambda5,
"It should be a number between -1 and 1"))}
# Don't allow numbers in parameterisation -
# included as a warning here, so the main one is a stop.
if(!is.character(param)) {warning(paste("param=", param,
"It shouldn't be a number, ",
"it should be a string describing the parameterisation"))
}
if(length(lambda1) > 1) {
# using a vector for the parameters.
# Check that there aren't values in the individual lambda arguments
if (!(is.null(lambda2) & is.null(lambda3)& is.null(lambda4) &
is.null(lambda5)) )
{ stop("Call includes vector version of the lambda parameters",
"as well as the \nscalar version")
}
if ((length(lambda1) < 4) | (length(lambda1) > 5 ) )
{ stop(paste("argument lambda1 has length", length(lambda1),
"\nThis should be 1 (lambda parameters given as seperate ",
"arguments), 4 (vector argument \n for RS or FMKL ",
"parameterisation) or 5 (vector argument for fm5",
"parameterisation"))
}
if (length(lambda1)== 5)
{ if (param != "fm5") {
stop(paste("argument lambda1 has length", length(lambda1),
"which is not valid for the", param, "\nparameterisation"))
}
# else --- fm5, in vector form, ready for gl.check.lambda
}
if (length(lambda1)== 4) {
if (param == "fm5" )
{ stop(paste("argument lambda1 has length 4, which is not",
" valid for the fm5 \nparameterisation"))
}
# else --- 4 parameter versions in vector form, ready for gl.check.lambda
}
} else {
# single parameter arguments - check they are there,
# then collect them together
if (is.null(lambda2)) { stop("No value for lambda2") }
if (is.null(lambda3)) { stop("No value for lambda3") }
if (is.null(lambda4)) { stop("No value for lambda4") }
if ((is.null(lambda5)) & param=="fm5" ) {
stop("No value for lambda5") }
if (!(is.null(lambda5)) & param!="fm5") {
stop(paste("lambda5=", lambda5, " but there is no lambda 5 for the\n",
param, "parameterisation")) }
if (param != "fm5") { # A 4 parameter version
lambda1 <- c(lambda1, lambda2, lambda3, lambda4)
} else { # fm5
lambda1 <- c(lambda1, lambda2, lambda3, lambda4, lambda5)
}
}
# There is now an error if there is the wrong number of parameters, and
# lambda1 returned as a vector with 4 or 5 elements
# as.double is needed to remove data.frame attributes if lambda1 was
# extracted from a data.frame
as.double(lambda1)
}
# ------------------------------------------------------------------------------
.gl.check.lambda <-
function(lambdas, lambda2=NULL, lambda3=NULL, lambda4=NULL,
param="fmkl", lambda5=NULL, vect=FALSE)
{
# Checks to see that the lambda values given are allowed.
# There is a function called .gl.parameter.tidy that does the tidying
# around of parameters. It return a single vector, which contains the
# parameters.
# If you call this after .gl.parameter.tidy, let it know with the vect=TRUE
# argument
# If vect=TRUE, we don't need to tidy
if (vect) {
if (!is.null(lambda3)) {
warning("lambda3 should be null because you claim the",
" parameters are in a vector")
}
} else {
lambdas <-
.gl.parameter.tidy(lambdas, lambda2, lambda3, lambda4, param, lambda5)
}
if(param=="fm5"){lambda5 = lambdas[5]}
lambda4 = lambdas[4]
lambda3 = lambdas[3]
lambda2 = lambdas[2]
lambda1 = lambdas[1]
# I did have a check for finite lambdas, but that caused a
# problem with data frames,
# so I removed it - still need to include the limit results
param <- switch(param,
# Different tests apply for each parameterisation
freimer=, # allows for alternate expressions
frm=, # allows for alternate expressions
FMKL=,
fmkl={
if (lambda2<=0) {return(FALSE)}
else {return(TRUE)}
},
ramberg=, # Ramberg & Schmeiser
ram=,
RS=,
rs={
if (lambda3*lambda4>0) { # regions 3 and 4
# all values of lambda 3 and lambda 4 OK
# check lambda 2
if ((lambda3>0)&(lambda4>0)) { # region 3 - l2 >0
if (lambda2<=0) {
ret <- FALSE
} else {
ret <- TRUE
}
}
if ((lambda3<0)&(lambda4<0)) { # region 4 - l2 <0
if (lambda2>=0) {
ret <- FALSE
} else {
ret <- TRUE
}
}
} else { # other quadrants - lambda 2 must be negative, and lambda3
# lambda 4 values need checking.
if (lambda2>=0) {return(FALSE)}
# Rectangular regions where RS is not defined
if ((lambda3>0)&(lambda3<1)&(lambda4<0)) {return(FALSE)}
if ((lambda4>0)&(lambda4<1)&(lambda3<0)) {return(FALSE)}
# Different here because there are a
# number of ways in which the parameters can fail.
#
# Curved regions where RS is not defined
# change to shorter var names
lc <- lambda3
ld <- lambda4
if ((lambda3>-1)&(lambda3<0)&(lambda4>1)) { # region 5 or not?
if ( ((1-lc)^(1-lc)*(ld-1)^(ld-1))/((ld-lc)^(ld-lc)) > -lc/ld )
{return(FALSE)}
else {return(TRUE)}
}
# Second curved region
if ((lambda4>-1)&(lambda4<0)&(lambda3>1)) { # region 6 or not?
if ( ((1-ld)^(1-ld)*(lc-1)^(lc-1))/((lc-ld)^(lc-ld)) > -ld/lc )
{return(FALSE)}
else {return(TRUE)}
}
# There may be some limit results that mean these are not correct, but
# I'll check that later
# This is not the place where the possible l3, l4 zero values should appear
if (lambda3 == 0) {
warning('lambda 3 = 0 with RS parameterisation - possible problem')
if (lambda4 == 0) {return(FALSE)}
else {return(TRUE)}
}
if (lambda4 == 0) {
warning('lambda 5 = 0 with RS parameterisation - possible problem')
if (lambda4 == 0) {return(FALSE)}
else {return(TRUE)}
}
# If we get here, then the parameters are OK.
ret <- TRUE
}
},
fm5={
# make lambda5 - in here so it doesn't stuff up the
# other parameterisations
lambda5 <- lambdas[5]
if (lambda2<=0) {ret <- FALSE}
else { #Check lambda5 - should be between -1 and 1, but
# I haven't checked this against a piece of paper
if ((lambda5 >= -1) & (lambda5 <= 1)) {
ret <- TRUE
} else {
ret <- FALSE
}
}
},
stop("Error when checking validity of parameters.\n",
" Parameterisation must be fmkl, rs or fm5")
) # closes "switch"
ret
}
################################################################################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.