Nothing
R/builtin-glGld.R
In fBasics: Rmetrics - Markets and Basic Statistics
Defines functions
.gl.check.lambda
.gl.parameter.tidy
.rgld
.qdgl.fm5
.qdgl.fmkl
.qdgl.rs
.qdgld
.qgl.rs
.qgl.fm5
.qgl.fmkl
.qgld
.pgld
.dgld
# 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
}
################################################################################
Try the fBasics package in your browser
Any scripts or data that you put into this service are public.
fBasics documentation built on Nov. 3, 2023, 5:10 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .gl.check.lambda .gl.parameter.tidy .rgld .qdgl.fm5 .qdgl.fmkl .qdgl.rs .qdgld .qgl.rs .qgl.fm5 .qgl.fmkl .qgld .pgld .dgld
# 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
}
################################################################################
Try the fBasics package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.