R/helper-functions.R

In pmoments: Partial Moments

#################################################################################
##
##   R package pmoments by Alexios Ghalanos Copyright (C) 2008
##   This file is part of the R package pmoments.
##
##   The R package pmoments is free software: you can redistribute it and/or modify
##   it under the terms of the GNU General Public License as published by
##   the Free Software Foundation, either version 3 of the License, or
##   (at your option) any later version.
##
##   The R package pmoments 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 General Public License for more details.
##
#################################################################################

epsilon<-.Machine$double.eps

.sumC = function(x)
{
	z = embed(x,2)
	return(z[,1]+z[,2])
}

# Professor Ignacio Grossmann suggested this on the gams faq
# (www.gams.com/docs/FAQ/MODELING.htm)
.max.linear<-function(x)
{
	(sqrt(x^2 + epsilon^2 ) + x ) / 2
}

.min.linear<-function(x)
{
	-(sqrt((-x)^2 + epsilon^2 ) - x ) / 2
}

#setMethod("as.matrix",signature(x="PME"),.extractMoments)

.extractMoments<-function(x,...)
{
	if(x@datatype=="single")
	{
		moment=x@moment
		z1=.extractsAnalytical(x)
		z2=.extractsSample(x)
		z3=.extractsBoot(x)
		z4=.extractsJack(x)
		pmtable=cbind(z1,z2,z3,z4)
		rownames(pmtable)<-paste("PM",moment,sep="")
		pmtable<-as.data.frame(pmtable)
	} else{
		datanames=x@datanames
		moment=x@moment
		z1=.extractmAnalytical(x)
		z2=.extractmSample(x)
		z3=.extractmBoot(x)
		z4=.extractmJack(x)
		nonemp<-x@distribution[x@distribution!="empirical"]
		if(!is.null(z1)) z1n<-dim(z1)[1] else z1n<-0
		if(!is.null(z2)) z2n<-dim(z2)[1] else z2n<-0
		if(!is.null(z3)) z3n<-dim(z3)[1] else z3n<-0
		if(!is.null(z4)) z4n<-dim(z4)[1] else z4n<-0
		
		dnames=c(rownames(z1),rownames(z2),rownames(z3),rownames(z4))

		dist=c(nonemp,rep("empirical",length(datanames)-length(nonemp)))
		Method=c(rep("Analytical",z1n),rep("Sample",z2n),rep("Boot",z3n),rep("Jack",z4n))
		pmtable<-data.frame(Names=dnames,Method=Method,Distribution=dist,Moment=as.numeric(rbind(z1,z2,z3,z4)))
		# sort in order it was received
		orignames<-x@datanames
		pmtable<-pmtable[match(pmtable[,1],x@datanames),]
	}
	
	return(new("data.frame",pmtable))
	
}

setMethod("as.data.frame",signature(x="PME",row.names="missing", optional="missing"),.extractMoments)

.extractsAnalytical<-function(object)
{
	method=object@method
	if(any(method=="analytical")){
	pm=object@pm
	method=object@method
	m=length(object@moment)
	d=which(method=="analytical")
	df=which(object@distribution!="empirical")
	fn=length(df)
	fm=matrix(NA,ncol=fn,nrow=m)
	for(i in 1:fn){
		fm[,i]=as.numeric(lapply(object@pm$analytical,FUN=function(x) x[[i]][1]))
	}
	colnames(fm)<-object@distribution[df]
	} else{
		fm=NULL
	}
	return(fm)
}

.extractsSample<-function(object)
{
	pm=object@pm
	method=object@method
	m=length(object@moment)
	if(any(method=="sample")){
		fs=matrix(NA,ncol=1,nrow=m)
		fs[,1]=as.numeric(lapply(pm$sample,FUN=function(x) x))
		colnames(fs)<-"sample"
	} else{
		fs=NULL
	}
	return(fs)
}

.extractsBoot<-function(object)
{
	pm=object@pm
	method=object@method
	m=length(object@moment)
	if(any(method=="boot")){
		fb=matrix(NA,ncol=1,nrow=m)
		fb[,1]=as.numeric(lapply(pm$boot,FUN=function(x) mean(x$thetastar)))
		colnames(fb)<-"boot"
	} else{
		fb=NULL
	}
	return(fb)
}

.extractsJack<-function(object)
{
	pm=object@pm
	method=object@method
	m=length(object@moment)
	if(any(method=="jack")){
		fk=matrix(NA,ncol=1,nrow=m)
		fk[,1]=as.numeric(lapply(pm$jack,FUN=function(x) mean(x$jack.values)))
		colnames(fk)<-"jack"
	} else{
		fk=NULL
	}
	return(fk)
}


.extractmAnalytical<-function(object)
{
	method=object@method
	if(any(method=="analytical")){
		pm=object@pm
		d=which(method=="analytical")
		fm=matrix(NA,ncol=1,nrow=length(d))
		rownames(fm)<-object@datanames[d]
		colnames(fm)<-"analytical"
		for(i in 1:length(d)){
			fm[i,1]=as.numeric(object@pm[[d[i]]])
		}
	} else{
		fm=NULL
	}
	return(fm)
}

.extractmSample<-function(object)
{
	method=object@method
	if(any(method=="sample")){
		pm=object@pm
		d=which(method=="sample")
		fm=matrix(NA,ncol=1,nrow=length(d))
		rownames(fm)<-object@datanames[d]
		colnames(fm)<-"sample"
		for(i in 1:length(d)){
			fm[i,1]=as.numeric(object@pm[[d[i]]])
		}
	} else{
		fm=NULL
	}
	return(fm)
}

.extractmBoot<-function(object)
{
	method=object@method
	if(any(method=="boot")){
		pm=object@pm
		d=which(method=="boot")
		fm=matrix(NA,ncol=1,nrow=length(d))
		rownames(fm)<-object@datanames[d]
		colnames(fm)<-"boot"
		for(i in 1:length(d)){
			fm[i,1]=as.numeric(mean(object@pm[[d[i]]]$thetastar))
		}
	} else{
		fm=NULL
	}
	return(fm)
}

.extractmJack<-function(object)
{
	method=object@method
	if(any(method=="jack")){
		pm=object@pm
		d=which(method=="jack")
		fm=matrix(NA,ncol=1,nrow=length(d))
		rownames(fm)<-object@datanames[d]
		colnames(fm)<-"jack"
		for(i in 1:length(d)){
			fm[i,1]=as.numeric(mean(object@pm[[d[i]]]$jack.values))
		}
	} else{
		fm=NULL
	}
	return(fm)
}


.colorgradient<-function(n = 50, low.col = 0.6, high.col=0.7, saturation = 0.8) {
	if (n < 2) stop("n must be greater than 2")
	n1 <- n%/%2
	n2 <- n - n1
	c(hsv(low.col, saturation, seq(1,0.5,length=n1)),
			hsv(high.col, saturation, seq(0.5,1,length=n2)))
}