Nothing
## Check validity of GLD parameters
.is.egld = function(lmd){
out = TRUE
if(is.null(lmd)) out=FALSE
if(any(is.na(lmd))) out=FALSE
if(length(lmd)!=4) out=FALSE
if(lmd[2]<=0||lmd[3]<=0||lmd[4]<=0) out=FALSE
return(out)
}
##########################################################################
### EGLD basic functions
##########################################################################
pegld = function(x,lambda){
if(!.is.egld(lambda)) stop("Invalid EGLD parameters.")
pbeta((x-lambda[1])/lambda[2],lambda[3],lambda[4])
}
degld = function(x,lambda){
if(!.is.egld(lambda)) stop("Invalid EGLD parameters.")
xt = (x-lambda[1])/lambda[2]
dbeta(xt,lambda[3],lambda[4])/lambda[2]
}
qegld = function(p,lambda){
## x is a probability
if(!.is.egld(lambda)) stop("Invalid EGLD parameters.")
qbeta(p,lambda[3],lambda[4])*lambda[2]+lambda[1]
}
regld = function(n, lambda){
if(!.is.egld(lambda)) stop("Invalid EGLD parameters.")
lambda[1]+lambda[2]*rbeta(n,lambda[3],lambda[4])
}
fit.egld = function(x,xmin=NULL,xmax=NULL)
{
name <- deparse(substitute(x))
if(!is.numeric(x)) stop("\n\n\t 'x' must be numeric values!\n\n\n");
x = x[!is.na(x)]; n = length(x);
if(n <6) stop("\n\n\t Sample size is too small!\n\n\n");
if(n < 30) warning("\n Sample size might be too small to fit moment matching GLD!\n");
s = sd(x); gpoints = seq(min(x)-s,max(x)+s,length=100);
out = .bincount(x); os = out$y; nbins = length(out$y); xbin = out$x #
if(nbins<3) stop("Data has less than 3 bins. Double-check the raw data.")
if(is.null(xmin)){
if(is.null(xmax)){
type=0;xmin=min(x);xmax=max(x);
}else{
type=1;xmin=min(x);
}
}else{
if(is.null(xmax)){
type=2;xmax=max(x);
}else{
type=3;
}
}
y = (x-xmin)/(xmax-xmin)
s2 = var(y); mu=mean(y)
b3 = (1./mu-1.-s2/mu/mu)*mu^3./s2;
b4 = (1./mu-1.)*b3;
p = c(xmin,xmax,b3,b4)
out = .Fortran(.F_FitGBD, as.double(x), as.integer(length(x)),
as.integer(type), llk=as.double(0),p=as.double(p));
betas = out$p
gpoints = seq(betas[1],betas[1]+betas[2],length=100);
y = degld(gpoints,betas)
DF = nbins - 5
Fx = pegld(xbin[-1], betas); Fx = c(Fx[1], diff(Fx))
es = n*Fx; es = c(es,n-sum(es))
D2 = sum((os-es)^2/es)
pv = ifelse(DF>0,pchisq(D2,DF,lower.tail=FALSE),NA)
return(structure(list(x = gpoints,y = y,n = length(x),para=betas,
chisq = D2, p.value=pv, method="ml",
call = match.call(), data.name = name),
class = "egld"))
}
print.egld <- function (x, digits = NULL, ...)
{
tmp="";k=length(x$para)
for(i in 1:(k-1))
tmp = paste(tmp,formatC(x$para[i]),",",sep='')
tmp = paste(tmp,formatC(x$para[k]),sep='')
cat("\nCall:\n\t", deparse(x$call), "\n\nData: ", x$data.name,
" (", x$n, " obs.);", "\tMethod: = ", x$method,"\n",
"Parameters: \t'para' = (",tmp,
")\n\n", sep = "")
cat("\nGoodness-of-fit: \n\tChi-Square = ",x$chisq,
",\tp-value = ", format(x$p.value,6), "\n",sep = "")
print(summary(as.data.frame(x[c("x", "y")])), digits = digits,
...)
invisible(x)
}
plot.egld <-
function (x, main = NULL, xlab = NULL, ylab = "Density", type = "l",
zero.line = TRUE, ...)
{
tmp="";k=length(x$para)
for(i in 1:(k-1))
tmp = paste(tmp,formatC(x$para[i]),",",sep='')
tmp = paste(tmp,formatC(x$para[k]),sep='')
if (is.null(xlab))
xlab <- paste("N =", x$n, " Parameters =(",tmp,")")
if (is.null(main))
main <- deparse(x$call)
sele = x$y>0
plot.default(x$x[sele],x$y[sele], type = type,
main = main, xlab = xlab, ylab = ylab, ...)
if (zero.line)
abline(h = 0, lwd = 0.1, col = "gray")
invisible(NULL)
}
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.