Nothing
R/loggammarob.R
In robustloggamma: Robust Estimation of the Generalized log Gamma Model
Defines functions
loggammarob
loggammarob.control
loggammarob.QTau
loggammarob.WQTau
loggammarob.WL
loggammarob.ML
loggammarob.oneWL
loggammarob.test
print.loggammarob
summary.loggammarob
print.summary.loggammarob
Documented in
loggammarob
loggammarob.control
loggammarob.test
print.loggammarob
print.summary.loggammarob
summary.loggammarob
#############################################################
# loggammarob function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio.agostinelli@unitn.it
# Date: February, 5, 2015
# Version: 0.2
# Copyright (C) 2015 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob <- function(x, start=NULL, weights=rep(1, length(x)), method=c("oneWL", "WQTau", "WL", "QTau", "ML"), control, ...) {
method <- match.arg(method)
x <- na.omit(x)
if (!is.null(na <- attr(x, "na.action")))
weights <- weights[-na]
or <- order(x)
if (any(or!=1:length(x))) {
x <- sort(x)
warning("data 'x' are sorted")
}
weights <- weights[or]
if (missing(control))
control <- if (missing(method))
loggammarob.control(...)
else
loggammarob.control(method = method, ...)
cl <- match.call()
if (!missing(control) && !missing(method) && method != control$method) {
warning("Methods argument set by method is different from method in control\n", "Using method = ", method)
control$method <- method
}
## x.orig <- x
## mx <- median(x)
## sx <- mad(x)
## x <- (x - mx)/sx
if (control$method=="QTau")
result <- loggammarob.QTau(x, weights, control)
else if (method=="WQTau")
result <- loggammarob.WQTau(x, weights, control)
else if (method=="WL")
result <- loggammarob.WL(x, start, weights, control)
else if (method=="oneWL")
result <- loggammarob.oneWL(x, start, weights, control)
else if (method=="ML")
result <- loggammarob.ML(x, start, weights, control)
## result$mu <- result$mu*sx+mx
## result$sigma <- result$sigma*sx
result$eta <- Exp.response(result$mu, result$sigma, result$lambda)
## result$data <- x.orig
result$data <- x
result$method <- method
result$call <- cl
class(result) <- "loggammarob"
return(result)
}
#############################################################
# loggammarob.control function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio.agostinelli@unitn.it
# Date: November, 27, 2015
# Version: 0.1-2
# Copyright (C) 2015 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.control <- function(method="oneWL", tuning.rho=1.547647, tuning.psi=6.08, lower=-7, upper=7, n=201, max.it=750, refine.tol=1e-6, solve.tol=1e-7, nResample=100, bw=0.3, smooth=NULL, raf=c("NED","GKL","PWD","HD","SCHI2"), tau=1, subdivisions=1000, lambda.step=TRUE, sigma.step=TRUE, step=1, minw=0.04, nexp=1000, reparam=NULL, bootstrap=FALSE, bootstrap.lambda=NULL, qthreshold=0.9, nTML=2000, xmax=1e100, iter=1, pcut=0.997, compute.rd=FALSE, eps.outlier = function(nobs) 0.1 / nobs) {
#tuning.rho, tuning.chi # c1, c2
#lower,upper # optimization interval
#n # optimization grid
#max.it # maxit in refinement
#refine.tol # relative tolerance in refinement
#nResample # number of subsamples
## Parameters for Fully Itereated Weighted Likelihood and for One Step Weighted Likelihood
#bw # bandwidth
#smooth # smooth parameter, if not NULL then bw = smooth*\hat{sd}
#raf # residual adjustment function c("SCHI2", "NED", "HD", "GKL")
#tau # parameter for the GKL and PD family
#subdivisions
#lambda.step
#sigma.step
#step #step length in the oneWL function
#Used only in the oneWL function
#minw #below this value the weights are set to zero.
#nexp #number of quantile points used to approximate the Expected Jacobian.
#reparam #list for reparametrization of the sigma parameter.
#bootstrap #to bootstrap the oneWL.
#bootstrap.lambda#to bootstrap the oneWL.
raf <- match.arg(raf)
if (lower >= upper)
stop("'lower' must be smaller than 'upper'")
if (length(step) != 1 & length(step) != 3)
stop("'step' must be a scalar or a vector of length 3")
if (length(n) > 1 || n < 0)
stop("'n' must be a positive scalar")
n <- round(n)
n <- ifelse(n < 1, 1, n)
if (length(max.it) > 1 || max.it < 0)
stop("'max.it' must be a positive scalar")
max.it <- round(max.it)
max.it <- ifelse(max.it < 1, 1, max.it)
if (length(refine.tol) > 1 || refine.tol < 0)
stop("'refine.tol' must be a positive scalar")
if (length(bw) > 1 || bw < 0)
stop("'bw' must be a positive scalar")
if (!is.null(smooth) && (length(bw) > 1 || bw < 0))
stop("'smooth' must be NULL or a positive scalar")
if (length(subdivisions) > 1 || subdivisions < 0)
stop("'subdivisions' must be a positive scalar")
subdivisions <- round(subdivisions)
subdivisions <- ifelse(subdivisions < 1, 1, subdivisions)
if (!is.logical(lambda.step))
stop("'lambda.step' must be a logical")
if (!is.logical(sigma.step))
stop("'sigma.step' must be a logical")
if (length(minw) > 1 || minw > 1 | minw < 0)
stop("'minw' must be a scalar in [0,1]")
if (length(nexp) > 1 || nexp < 0)
stop("'nexp' must be a positive scalar")
nexp <- round(nexp)
nexp <- ifelse(nexp < 1, 1, nexp)
############ EXAMPLE OF REPARAM ################################
#reparam.loggamma <- list(gam=function(sigma) sqrt(sigma),
# gaminv=function(gam) gam^2,
# delta=function(sigma) 2*sqrt(sigma))
if (!is.null(reparam)) {
badreparam <- !is.list(reparam) || !is.function(reparam$gam) || !is.function(reparam$gaminv) || !is.function(reparam$delta)
if (badreparam)
stop("'reparam' must be a list of 3 functions named 'gam', 'gaminv' and 'delta'")
}
if (!is.logical(bootstrap))
stop("'bootstrap' must be a logical")
if (qthreshold < 0.5 | qthreshold > 1)
stop("'qthreshold' must be in the interval [0.5,1]")
nTML <- round(nTML)
if (nTML < 1)
stop("nTML must be greater than 1")
if (pcut < 0.5 | pcut > 1)
stop("'pcut' must be in the interval [0.5,1]")
res <- list(method=method, tuning.rho=tuning.rho, tuning.psi=tuning.psi, lower=lower, upper=upper, n=n, max.it=max.it, refine.tol=refine.tol, solve.tol=solve.tol, nResample=nResample, bw=bw, smooth=smooth, raf=raf, tau=tau, subdivisions=subdivisions, lambda.step=lambda.step, sigma.step=sigma.step,step=step,minw=minw, nexp=nexp, reparam=reparam, bootstrap=bootstrap, bootstrap.lambda=bootstrap.lambda, qthreshold=qthreshold, nTML=nTML, xmax=xmax, iter=iter, pcut=pcut, compute.rd=compute.rd, eps.outlier=eps.outlier)
return(res)
}
#############################################################
# loggammarob.QTau function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.QTau <- function(x, w=rep(1, length(x)), control=loggammarob.control()) {
lgrid <- seq(control$lower, control$upper, length.out=control$n)
if (!control$bootstrap) {
resQ <- WQtau(ri=x,w=w,lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
} else {
resQ <- WQtauboot(ri=x,w=w,lambda=control$bootstrap.lambda,step=0.25,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
}
result <- list(mu=resQ$mu2, sigma=resQ$sig2, lambda=resQ$lam2, weights=w, iterations=NULL, error=NULL)
return(result)
}
#############################################################
# loggammarob.WQTau function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.WQTau <- function(x, w=rep(1, length(x)), control=loggammarob.control()) {
lgrid <- seq(control$lower, control$upper, length.out=control$n)
resQ <- WQtau(ri=x,w=w,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
pp <- ppoints(length(x))
ql <- qloggamma(p=pp, lambda=resQ$lam1)
dl <- dloggamma(x=ql, lambda=resQ$lam1)
vl <- pp*(1-pp)/dl^2
wl <- 1/sqrt(vl)
resWQ <- WQtau(ri=x,w=wl,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
result <- list(mu=resWQ$mu2, sigma=resWQ$sig2, lambda=resWQ$lam2, weights=wl, iterations=NULL, error=NULL, QTau=list(mu=resQ$mu2, sigma=resQ$sig2, lambda=resQ$lam2, weights=w, iterations=NULL, error=NULL))
return(result)
}
#############################################################
# loggammarob.WL function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.WL <- function(x, start=NULL, w=rep(1, length(x)), control=loggammarob.control()) {
if (is.null(start)) {
lgrid <- seq(control$lower, control$upper, length.out=control$n)
resQ <- WQtau(ri=x,w=w,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
pp <- ppoints(length(x))
ql <- qloggamma(p=pp, lambda=resQ$lam1)
dl <- dloggamma(x=ql, lambda=resQ$lam1)
vl <- pp*(1-pp)/dl^2
wl <- 1/sqrt(vl)
resWQ <- WQtau(ri=x,w=wl,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
start <- c(resWQ$mu2, resWQ$sig2, resWQ$lam2)
QTau <- list(mu=resQ$mu2, sigma=resQ$sig2, lambda=resQ$lam2, weights=w, iterations=NULL, error=NULL)
WQTau <- list(mu=resWQ$mu2, sigma=resWQ$sig2, lambda=resWQ$lam2, weights=wl, iterations=NULL, error=NULL)
} else {
if (!is.numeric(start))
stop("'start' must be a numeric vector")
if (length(start)!=3)
stop("'start' must be a vector of length 3: mu, sigma2, lambda")
QTau <- NULL
WQTau <- NULL
}
if (!is.null(control$smooth))
control$bw <- control$smooth*sqrt(start[2])
resWL <- Disparity.WML.loggamma(y=x,mu0=start[1],sig0=start[2],lam0=start[3],lam.low=control$lower,lam.sup=control$upper,tol=control$refine.tol,maxit=control$max.it,lstep=control$lambda.step,sigstep=control$sigma.step,bw=control$bw,raf=control$raf,tau=control$tau,nmod=control$subdivisions,minw=control$minw)
result <- list(mu=resWL$mu, sigma=resWL$sig, lambda=resWL$lam, weights=resWL$weights, iterations=resWL$nit, error=NULL, QTau=QTau, WQTau=WQTau)
return(result)
}
#############################################################
# loggammarob.ML function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.ML <- function(x, start=NULL, w=rep(1, length(x)), control=loggammarob.control()) {
if (is.null(start)) {
lgrid <- seq(control$lower, control$upper, length.out=control$n)
resQ <- WQtau(ri=x,w=w,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
start <- c(resQ$mu2, resQ$sig2, resQ$lam2)
} else {
if (!is.numeric(start))
stop("'start' must be a numeric vector")
if (length(start)!=3)
stop("'start' must be a vector of length 3: mu, sigma2, lambda")
}
resML <- WML.loggamma(y=x,wi=w,mu0=start[1],sig0=start[2],lam0=start[3],lam.low=control$lower,lam.sup=control$upper,tol=control$refine.tol,maxit=control$max.it,lstep=control$lambda.step,sigstep=control$sigma.step)
result <- list(mu=resML$mu, sigma=resML$sig, lambda=resML$lam, weights=w, iterations=resML$nit, error=NULL)
return(result)
}
#############################################################
# loggammarob.oneWL function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.oneWL <- function(x, start=NULL, w=rep(1, length(x)), control=loggammarob.control()) {
if (is.null(start)) {
lgrid <- seq(control$lower, control$upper, length.out=control$n)
if (!control$bootstrap) {
resQ <- WQtau(ri=x,w=w,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
pp <- ppoints(length(x))
ql <- qloggamma(p=pp, lambda=resQ$lam1)
dl <- dloggamma(x=ql, lambda=resQ$lam1)
vl <- pp*(1-pp)/dl^2
wl <- 1/sqrt(vl)
} else {
resQ <- resWQ <- WQtauboot(ri=x,w=w,lambda=control$bootstrap.lambda,step=0.25,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
wl <- w
}
if (!control$bootstrap)
resWQ <- WQtau(ri=x,w=wl,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
start <- c(resWQ$mu2, resWQ$sig2, resWQ$lam2)
QTau <- list(mu=resQ$mu2, sigma=resQ$sig2, lambda=resQ$lam2, weights=w, iterations=NULL, error=NULL)
WQTau <- list(mu=resWQ$mu2, sigma=resWQ$sig2, lambda=resWQ$lam2, weights=wl, iterations=NULL, error=NULL)
} else {
if (!is.numeric(start))
stop("'start' must be a numeric vector")
if (length(start)!=3)
stop("'start' must be a vector of length 3: mu, sigma2, lambda")
QTau <- NULL
WQTau <- NULL
}
if (!is.null(control$smooth))
control$bw <- control$smooth*sqrt(start[2])
if (is.null(control$reparam))
resOWL <- WMLone.loggamma(yi.sorted=x,mu0=start[1],sig0=start[2],lam0=start[3],bw=control$bw,raf=control$raf,tau=control$tau,nmod=control$subdivisions,step=control$step,minw=control$minw,nexp=control$nexp)
else
resOWL <- WMLone.reparam.loggamma(yi.sorted=x,mu0=start[1],sig0=start[2],lam0=start[3],bw=control$bw,raf=control$raf,tau=control$tau,nmod=control$subdivisions,step=control$step,minw=control$minw,nexp=control$nexp,reparam=control$reparam)
result <- list(mu=resOWL$mu, sigma=resOWL$sig, lambda=resOWL$lam, weights=resOWL$weights, iterations=1, error=resOWL$error, step=control$step, QTau=QTau, WQTau=WQTau)
return(result)
}
#############################################################
# loggammarob.test function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.test <- function(x, mu=NULL, sigma=NULL, lambda=NULL, eta=NULL, type="Wald", conf.level = 0.95, prob=0.00001) {
## x an object from loggammarob
if (!(x$method %in% c("WL", "oneWL", "ML")))
stop("Inference is available only for 'WL', 'oneWL' and 'ML'")
if (is.null(mu) & is.null(sigma) & is.null(lambda) & is.null(eta)) {
mu <- 0
sigma <- 1
lambda <- 0
}
alternative <- "two.sided"
dname <- names(x$data)
type <- match.arg(type)
components <- c(!is.null(mu), !is.null(sigma), !is.null(lambda))
if (is.null(eta)) {
estimate <- c(x$mu, x$sigma, x$lambda)[components]
null.value <- c(mu, sigma, lambda)
if (type=="Wald") {
param <- estimate - null.value
df <- length(param)
acov <- ACOV.ML.loggamma(sigma=x$sigma,lambda=x$lambda, prob=prob)
if (df==1) {
hatJ <- 1/diag(acov$cov)[components]
} else {
hatJ <- acov$Fisher.Info
hatJ <- matrix(hatJ[as.logical(components%*%t(components))], nrow=length(param), byrow=TRUE)
}
wstat <- drop(sum(x$weights)*t(param)%*%hatJ%*%param)
pval <- pchisq(q=wstat, df=df, ncp=0, lower.tail = FALSE, log.p = FALSE)
estcov <- diag(acov$cov)[components]
if (df==1)
cint <- estimate+c(-1,1)*qnorm(1-(1-conf.level)/2)*sqrt(estcov)/sqrt(sum(x$weights))
else
cint <- c(NA,NA)
}
names(wstat) <- "ww"
names(df) <- "df"
names(null.value) <- c("mean", "scale", "shape")[components]
} else {
if (any(components))
warning("test with 'eta' together with other parameters is not implemented, only the test on 'eta' is performed")
eps <- 0.0001
estimate <- x$eta
null.value <- eta
if (type=="Wald") {
param <- estimate - null.value
if (x$lambda < eps)
warning("test and confidence intervals for 'eta' is not implemented for negative lambda")
hatvar <- AVAR.ML.eta.loggamma(mu=x$mu, sigma=x$sigma, lambda=x$lambda, prob=prob, eps=eps)
wstat <- sum(x$weights)*param^2/hatvar
df <- length(param)
pval <- pchisq(q=wstat, df=df, ncp=0, lower.tail = FALSE, log.p = FALSE)
cint <- estimate+c(-1,1)*qnorm(1-(1-conf.level)/2)*sqrt(hatvar)/sqrt(sum(x$weights))
}
names(wstat) <- "ww"
names(df) <- "df"
names(null.value) <- "mean(exp(X))"
}
type <- paste(ifelse(x$method=="ML","Classic","Weighted"), type, "Test based on", x$method, sep=" ")
attr(cint,"conf.level") <- conf.level
rval <- list(statistic=wstat, parameter=df, p.value=pval,
conf.int=cint, estimate=estimate, null.value=null.value,
alternative=alternative,
method=type, data.name=dname)
class(rval) <- "htest"
return(rval)
}
#############################################################
# print.loggammarob function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
######## PRINT
print.loggammarob <- function(x, digits = max(3, getOption("digits") - 3), ...) {
cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
cat("Location: ", format(x$mu, digits=digits))
cat(" Scale: ", format(x$sigma, digits=digits))
cat(" Shape: ", format(x$lambda, digits=digits))
cat(" E(exp(X)): ", format(x$eta, digits=digits))
cat("\n")
invisible(x)
}
#############################################################
# summary.loggammarob function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
########### SUMMARY
summary.loggammarob <- function(object, p=NULL, conf.level=0.95, prob=0.00001, ...) {
if ((object$method %in% c("WL", "oneWL", "ML"))) {
p <- p[ifelse(p > 0 & p < 1, TRUE, FALSE)]
nrep <- 4+length(p)
if (length(conf.level) == 1)
conf.level <- rep(conf.level, nrep)
if (length(conf.level) != nrep)
stop("'conf.level' must have length equal to 1 or 4+length(p)")
z <- qnorm(1-(1-conf.level)/2)
if (is.finite(object$eta)) {
asd <- sqrt(diag(ACOV.ML.loggamma(sigma=object$sigma,lambda=object$lambda,prob=prob)$cov))
asdeta <- sqrt(AVAR.ML.eta.loggamma(mu=object$mu, sigma=object$sigma, lambda=object$lambda, prob=prob, eps=0.0001, npoints=100000))
} else {
asd <- rep(NA,3)
asdeta <- NA
}
object$muconf.int <- object$mu+c(-1,1)*z[1]*asd[1]/sqrt(sum(object$weights))
object$sigmaconf.int <- object$sigma+c(-1,1)*z[2]*asd[2]/sqrt(sum(object$weights))
object$lambdaconf.int <- object$lambda+c(-1,1)*z[3]*asd[3]/sqrt(sum(object$weights))
object$etaconf.int <- object$eta+c(-1,1)*z[4]*asdeta/sqrt(sum(object$weights))
attr(object$muconf.int, "conf.level") <- conf.level[1]
attr(object$sigmaconf.int, "conf.level") <- conf.level[2]
attr(object$lambdaconf.int, "conf.level") <- conf.level[3]
attr(object$etaconf.int, "conf.level") <- conf.level[4]
object$muse <- asd[1]/sqrt(sum(object$weights))
object$sigmase <- asd[2]/sqrt(sum(object$weights))
object$lambdase <- asd[3]/sqrt(sum(object$weights))
object$etase <- asdeta/sqrt(sum(object$weights))
object$p <- p
if (!is.null(p)) {
object$q <- asdq <- rep(NA, length(p))
object$qconf.int <- matrix(NA, nrow=length(p), ncol=2)
for (i in 1:length(p)) {
object$q[i] <- qloggamma(p=p[i], mu=object$mu, sigma=object$sigma, lambda=object$lambda)
if (is.finite(object$eta))
asdq[i] <- sqrt(AVAR.ML.quantile.loggamma(p=p[i], mu=object$mu, sigma=object$sigma, lambda=object$lambda, prob=prob))
else
asdq[i] <- NA
object$qconf.int[i,] <- object$q[i]+c(-1,1)*z[4+i]*asdq[i]/sqrt(sum(object$weights))
}
##### object$qconf.int <- drop(qconf.int)
attr(object$qconf.int, "conf.level") <- conf.level[4+(1:length(p))]
object$qse <- asdq/sqrt(sum(object$weights))
}
}
object$call <- match.call()
class(object) <- "summary.loggammarob"
return(object)
}
#############################################################
# print.summary.loggammarob function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
#PRINT.SUMMARY
print.summary.loggammarob <- function(x, digits = max(3, getOption("digits") - 3), ...) {
cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
cat("Location: ", format(x$mu, digits=digits))
if ((x$method %in% c("WL", "oneWL", "ML"))) {
cat(" s.e. ", format(x$muse, digits=digits), "\n")
cat("(", format(x$muconf.int[1L], digits=digits), ", ",
format(x$muconf.int[2L], digits=digits), ") \n")
cat(format(100 * attr(x$muconf.int, "conf.level")),
"percent confidence interval\n\n")
} else {
cat("\n")
}
cat("Scale: ", format(x$sigma, digits=digits))
if ((x$method %in% c("WL", "oneWL", "ML"))) {
cat(" s.e. ", format(x$sigmase, digits=digits), "\n")
cat("(", format(x$sigmaconf.int[1L], digits=digits), ", ",
format(x$sigmaconf.int[2L], digits=digits), ") \n")
cat(format(100 * attr(x$sigmaconf.int, "conf.level")),
"percent confidence interval\n\n")
} else {
cat("\n")
}
cat("Shape: ", format(x$lambda, digits=digits))
if ((x$method %in% c("WL", "oneWL", "ML"))) {
cat(" s.e. ", format(x$lambdase, digits=digits), "\n")
cat("(", format(x$lambdaconf.int[1L], digits=digits), ", ",
format(x$lambdaconf.int[2L], digits=digits), ") \n")
cat(format(100 * attr(x$lambdaconf.int, "conf.level")),
"percent confidence interval\n\n")
} else {
cat("\n")
}
cat("Mean(exp(X)): ", format(x$eta, digits=digits))
if ((x$method %in% c("WL", "oneWL", "ML"))) {
cat(" s.e. ", format(x$etase, digits=digits), "\n")
cat("(", format(x$etaconf.int[1L], digits=digits), ", ",
format(x$etaconf.int[2L], digits=digits), ") \n")
cat(format(100 * attr(x$etaconf.int, "conf.level")),
"percent confidence interval\n\n\n")
} else {
cat("\n")
}
######### Quantiles estimation and confidence interval
if (!is.null(x$p)) {
for (i in 1:length(x$p)) {
cat("Quantile of order ", format(x$p[i], digits=digits), ": ", format(x$q[i], digits=digits))
if ((x$method %in% c("WL", "oneWL", "ML"))) {
cat(" s.e. ", format(x$qse[i], digits=digits), "\n")
cat("(", format(x$qconf.int[i,1L], digits=digits), ", ",
format(x$qconf.int[i,2L], digits=digits), ") \n")
cat(format(100 * attr(x$qconf.int, "conf.level")[i]),
"percent confidence interval\n\n\n")
} else {
cat("\n")
}
}
}
if ((x$method %in% c("WL", "oneWL", "ML")))
cat("\n\n")
else
cat("\n\nConfidence intervals and Standard Errors are available only for 'WL', 'oneWL' and 'ML'\n\n\n")
if (x$method %in% c("WL", "oneWL", "WQTau"))
summarizeRobWeights(x$weights, digits = digits, ...)
invisible(x)
}
Try the robustloggamma package in your browser
Any scripts or data that you put into this service are public.
robustloggamma documentation built on May 1, 2019, 9:20 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 loggammarob loggammarob.control loggammarob.QTau loggammarob.WQTau loggammarob.WL loggammarob.ML loggammarob.oneWL loggammarob.test print.loggammarob summary.loggammarob print.summary.loggammarob
Documented in loggammarob loggammarob.control loggammarob.test print.loggammarob print.summary.loggammarob summary.loggammarob
#############################################################
# loggammarob function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio.agostinelli@unitn.it
# Date: February, 5, 2015
# Version: 0.2
# Copyright (C) 2015 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob <- function(x, start=NULL, weights=rep(1, length(x)), method=c("oneWL", "WQTau", "WL", "QTau", "ML"), control, ...) {
method <- match.arg(method)
x <- na.omit(x)
if (!is.null(na <- attr(x, "na.action")))
weights <- weights[-na]
or <- order(x)
if (any(or!=1:length(x))) {
x <- sort(x)
warning("data 'x' are sorted")
}
weights <- weights[or]
if (missing(control))
control <- if (missing(method))
loggammarob.control(...)
else
loggammarob.control(method = method, ...)
cl <- match.call()
if (!missing(control) && !missing(method) && method != control$method) {
warning("Methods argument set by method is different from method in control\n", "Using method = ", method)
control$method <- method
}
## x.orig <- x
## mx <- median(x)
## sx <- mad(x)
## x <- (x - mx)/sx
if (control$method=="QTau")
result <- loggammarob.QTau(x, weights, control)
else if (method=="WQTau")
result <- loggammarob.WQTau(x, weights, control)
else if (method=="WL")
result <- loggammarob.WL(x, start, weights, control)
else if (method=="oneWL")
result <- loggammarob.oneWL(x, start, weights, control)
else if (method=="ML")
result <- loggammarob.ML(x, start, weights, control)
## result$mu <- result$mu*sx+mx
## result$sigma <- result$sigma*sx
result$eta <- Exp.response(result$mu, result$sigma, result$lambda)
## result$data <- x.orig
result$data <- x
result$method <- method
result$call <- cl
class(result) <- "loggammarob"
return(result)
}
#############################################################
# loggammarob.control function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio.agostinelli@unitn.it
# Date: November, 27, 2015
# Version: 0.1-2
# Copyright (C) 2015 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.control <- function(method="oneWL", tuning.rho=1.547647, tuning.psi=6.08, lower=-7, upper=7, n=201, max.it=750, refine.tol=1e-6, solve.tol=1e-7, nResample=100, bw=0.3, smooth=NULL, raf=c("NED","GKL","PWD","HD","SCHI2"), tau=1, subdivisions=1000, lambda.step=TRUE, sigma.step=TRUE, step=1, minw=0.04, nexp=1000, reparam=NULL, bootstrap=FALSE, bootstrap.lambda=NULL, qthreshold=0.9, nTML=2000, xmax=1e100, iter=1, pcut=0.997, compute.rd=FALSE, eps.outlier = function(nobs) 0.1 / nobs) {
#tuning.rho, tuning.chi # c1, c2
#lower,upper # optimization interval
#n # optimization grid
#max.it # maxit in refinement
#refine.tol # relative tolerance in refinement
#nResample # number of subsamples
## Parameters for Fully Itereated Weighted Likelihood and for One Step Weighted Likelihood
#bw # bandwidth
#smooth # smooth parameter, if not NULL then bw = smooth*\hat{sd}
#raf # residual adjustment function c("SCHI2", "NED", "HD", "GKL")
#tau # parameter for the GKL and PD family
#subdivisions
#lambda.step
#sigma.step
#step #step length in the oneWL function
#Used only in the oneWL function
#minw #below this value the weights are set to zero.
#nexp #number of quantile points used to approximate the Expected Jacobian.
#reparam #list for reparametrization of the sigma parameter.
#bootstrap #to bootstrap the oneWL.
#bootstrap.lambda#to bootstrap the oneWL.
raf <- match.arg(raf)
if (lower >= upper)
stop("'lower' must be smaller than 'upper'")
if (length(step) != 1 & length(step) != 3)
stop("'step' must be a scalar or a vector of length 3")
if (length(n) > 1 || n < 0)
stop("'n' must be a positive scalar")
n <- round(n)
n <- ifelse(n < 1, 1, n)
if (length(max.it) > 1 || max.it < 0)
stop("'max.it' must be a positive scalar")
max.it <- round(max.it)
max.it <- ifelse(max.it < 1, 1, max.it)
if (length(refine.tol) > 1 || refine.tol < 0)
stop("'refine.tol' must be a positive scalar")
if (length(bw) > 1 || bw < 0)
stop("'bw' must be a positive scalar")
if (!is.null(smooth) && (length(bw) > 1 || bw < 0))
stop("'smooth' must be NULL or a positive scalar")
if (length(subdivisions) > 1 || subdivisions < 0)
stop("'subdivisions' must be a positive scalar")
subdivisions <- round(subdivisions)
subdivisions <- ifelse(subdivisions < 1, 1, subdivisions)
if (!is.logical(lambda.step))
stop("'lambda.step' must be a logical")
if (!is.logical(sigma.step))
stop("'sigma.step' must be a logical")
if (length(minw) > 1 || minw > 1 | minw < 0)
stop("'minw' must be a scalar in [0,1]")
if (length(nexp) > 1 || nexp < 0)
stop("'nexp' must be a positive scalar")
nexp <- round(nexp)
nexp <- ifelse(nexp < 1, 1, nexp)
############ EXAMPLE OF REPARAM ################################
#reparam.loggamma <- list(gam=function(sigma) sqrt(sigma),
# gaminv=function(gam) gam^2,
# delta=function(sigma) 2*sqrt(sigma))
if (!is.null(reparam)) {
badreparam <- !is.list(reparam) || !is.function(reparam$gam) || !is.function(reparam$gaminv) || !is.function(reparam$delta)
if (badreparam)
stop("'reparam' must be a list of 3 functions named 'gam', 'gaminv' and 'delta'")
}
if (!is.logical(bootstrap))
stop("'bootstrap' must be a logical")
if (qthreshold < 0.5 | qthreshold > 1)
stop("'qthreshold' must be in the interval [0.5,1]")
nTML <- round(nTML)
if (nTML < 1)
stop("nTML must be greater than 1")
if (pcut < 0.5 | pcut > 1)
stop("'pcut' must be in the interval [0.5,1]")
res <- list(method=method, tuning.rho=tuning.rho, tuning.psi=tuning.psi, lower=lower, upper=upper, n=n, max.it=max.it, refine.tol=refine.tol, solve.tol=solve.tol, nResample=nResample, bw=bw, smooth=smooth, raf=raf, tau=tau, subdivisions=subdivisions, lambda.step=lambda.step, sigma.step=sigma.step,step=step,minw=minw, nexp=nexp, reparam=reparam, bootstrap=bootstrap, bootstrap.lambda=bootstrap.lambda, qthreshold=qthreshold, nTML=nTML, xmax=xmax, iter=iter, pcut=pcut, compute.rd=compute.rd, eps.outlier=eps.outlier)
return(res)
}
#############################################################
# loggammarob.QTau function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.QTau <- function(x, w=rep(1, length(x)), control=loggammarob.control()) {
lgrid <- seq(control$lower, control$upper, length.out=control$n)
if (!control$bootstrap) {
resQ <- WQtau(ri=x,w=w,lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
} else {
resQ <- WQtauboot(ri=x,w=w,lambda=control$bootstrap.lambda,step=0.25,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
}
result <- list(mu=resQ$mu2, sigma=resQ$sig2, lambda=resQ$lam2, weights=w, iterations=NULL, error=NULL)
return(result)
}
#############################################################
# loggammarob.WQTau function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.WQTau <- function(x, w=rep(1, length(x)), control=loggammarob.control()) {
lgrid <- seq(control$lower, control$upper, length.out=control$n)
resQ <- WQtau(ri=x,w=w,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
pp <- ppoints(length(x))
ql <- qloggamma(p=pp, lambda=resQ$lam1)
dl <- dloggamma(x=ql, lambda=resQ$lam1)
vl <- pp*(1-pp)/dl^2
wl <- 1/sqrt(vl)
resWQ <- WQtau(ri=x,w=wl,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
result <- list(mu=resWQ$mu2, sigma=resWQ$sig2, lambda=resWQ$lam2, weights=wl, iterations=NULL, error=NULL, QTau=list(mu=resQ$mu2, sigma=resQ$sig2, lambda=resQ$lam2, weights=w, iterations=NULL, error=NULL))
return(result)
}
#############################################################
# loggammarob.WL function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.WL <- function(x, start=NULL, w=rep(1, length(x)), control=loggammarob.control()) {
if (is.null(start)) {
lgrid <- seq(control$lower, control$upper, length.out=control$n)
resQ <- WQtau(ri=x,w=w,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
pp <- ppoints(length(x))
ql <- qloggamma(p=pp, lambda=resQ$lam1)
dl <- dloggamma(x=ql, lambda=resQ$lam1)
vl <- pp*(1-pp)/dl^2
wl <- 1/sqrt(vl)
resWQ <- WQtau(ri=x,w=wl,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
start <- c(resWQ$mu2, resWQ$sig2, resWQ$lam2)
QTau <- list(mu=resQ$mu2, sigma=resQ$sig2, lambda=resQ$lam2, weights=w, iterations=NULL, error=NULL)
WQTau <- list(mu=resWQ$mu2, sigma=resWQ$sig2, lambda=resWQ$lam2, weights=wl, iterations=NULL, error=NULL)
} else {
if (!is.numeric(start))
stop("'start' must be a numeric vector")
if (length(start)!=3)
stop("'start' must be a vector of length 3: mu, sigma2, lambda")
QTau <- NULL
WQTau <- NULL
}
if (!is.null(control$smooth))
control$bw <- control$smooth*sqrt(start[2])
resWL <- Disparity.WML.loggamma(y=x,mu0=start[1],sig0=start[2],lam0=start[3],lam.low=control$lower,lam.sup=control$upper,tol=control$refine.tol,maxit=control$max.it,lstep=control$lambda.step,sigstep=control$sigma.step,bw=control$bw,raf=control$raf,tau=control$tau,nmod=control$subdivisions,minw=control$minw)
result <- list(mu=resWL$mu, sigma=resWL$sig, lambda=resWL$lam, weights=resWL$weights, iterations=resWL$nit, error=NULL, QTau=QTau, WQTau=WQTau)
return(result)
}
#############################################################
# loggammarob.ML function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.ML <- function(x, start=NULL, w=rep(1, length(x)), control=loggammarob.control()) {
if (is.null(start)) {
lgrid <- seq(control$lower, control$upper, length.out=control$n)
resQ <- WQtau(ri=x,w=w,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
start <- c(resQ$mu2, resQ$sig2, resQ$lam2)
} else {
if (!is.numeric(start))
stop("'start' must be a numeric vector")
if (length(start)!=3)
stop("'start' must be a vector of length 3: mu, sigma2, lambda")
}
resML <- WML.loggamma(y=x,wi=w,mu0=start[1],sig0=start[2],lam0=start[3],lam.low=control$lower,lam.sup=control$upper,tol=control$refine.tol,maxit=control$max.it,lstep=control$lambda.step,sigstep=control$sigma.step)
result <- list(mu=resML$mu, sigma=resML$sig, lambda=resML$lam, weights=w, iterations=resML$nit, error=NULL)
return(result)
}
#############################################################
# loggammarob.oneWL function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.oneWL <- function(x, start=NULL, w=rep(1, length(x)), control=loggammarob.control()) {
if (is.null(start)) {
lgrid <- seq(control$lower, control$upper, length.out=control$n)
if (!control$bootstrap) {
resQ <- WQtau(ri=x,w=w,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
pp <- ppoints(length(x))
ql <- qloggamma(p=pp, lambda=resQ$lam1)
dl <- dloggamma(x=ql, lambda=resQ$lam1)
vl <- pp*(1-pp)/dl^2
wl <- 1/sqrt(vl)
} else {
resQ <- resWQ <- WQtauboot(ri=x,w=w,lambda=control$bootstrap.lambda,step=0.25,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
wl <- w
}
if (!control$bootstrap)
resWQ <- WQtau(ri=x,w=wl,lgrid=lgrid,c1=control$tuning.rho,c2=control$tuning.psi,N=control$nResample,maxit=control$max.it,tolr=control$refine.tol)
start <- c(resWQ$mu2, resWQ$sig2, resWQ$lam2)
QTau <- list(mu=resQ$mu2, sigma=resQ$sig2, lambda=resQ$lam2, weights=w, iterations=NULL, error=NULL)
WQTau <- list(mu=resWQ$mu2, sigma=resWQ$sig2, lambda=resWQ$lam2, weights=wl, iterations=NULL, error=NULL)
} else {
if (!is.numeric(start))
stop("'start' must be a numeric vector")
if (length(start)!=3)
stop("'start' must be a vector of length 3: mu, sigma2, lambda")
QTau <- NULL
WQTau <- NULL
}
if (!is.null(control$smooth))
control$bw <- control$smooth*sqrt(start[2])
if (is.null(control$reparam))
resOWL <- WMLone.loggamma(yi.sorted=x,mu0=start[1],sig0=start[2],lam0=start[3],bw=control$bw,raf=control$raf,tau=control$tau,nmod=control$subdivisions,step=control$step,minw=control$minw,nexp=control$nexp)
else
resOWL <- WMLone.reparam.loggamma(yi.sorted=x,mu0=start[1],sig0=start[2],lam0=start[3],bw=control$bw,raf=control$raf,tau=control$tau,nmod=control$subdivisions,step=control$step,minw=control$minw,nexp=control$nexp,reparam=control$reparam)
result <- list(mu=resOWL$mu, sigma=resOWL$sig, lambda=resOWL$lam, weights=resOWL$weights, iterations=1, error=resOWL$error, step=control$step, QTau=QTau, WQTau=WQTau)
return(result)
}
#############################################################
# loggammarob.test function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
loggammarob.test <- function(x, mu=NULL, sigma=NULL, lambda=NULL, eta=NULL, type="Wald", conf.level = 0.95, prob=0.00001) {
## x an object from loggammarob
if (!(x$method %in% c("WL", "oneWL", "ML")))
stop("Inference is available only for 'WL', 'oneWL' and 'ML'")
if (is.null(mu) & is.null(sigma) & is.null(lambda) & is.null(eta)) {
mu <- 0
sigma <- 1
lambda <- 0
}
alternative <- "two.sided"
dname <- names(x$data)
type <- match.arg(type)
components <- c(!is.null(mu), !is.null(sigma), !is.null(lambda))
if (is.null(eta)) {
estimate <- c(x$mu, x$sigma, x$lambda)[components]
null.value <- c(mu, sigma, lambda)
if (type=="Wald") {
param <- estimate - null.value
df <- length(param)
acov <- ACOV.ML.loggamma(sigma=x$sigma,lambda=x$lambda, prob=prob)
if (df==1) {
hatJ <- 1/diag(acov$cov)[components]
} else {
hatJ <- acov$Fisher.Info
hatJ <- matrix(hatJ[as.logical(components%*%t(components))], nrow=length(param), byrow=TRUE)
}
wstat <- drop(sum(x$weights)*t(param)%*%hatJ%*%param)
pval <- pchisq(q=wstat, df=df, ncp=0, lower.tail = FALSE, log.p = FALSE)
estcov <- diag(acov$cov)[components]
if (df==1)
cint <- estimate+c(-1,1)*qnorm(1-(1-conf.level)/2)*sqrt(estcov)/sqrt(sum(x$weights))
else
cint <- c(NA,NA)
}
names(wstat) <- "ww"
names(df) <- "df"
names(null.value) <- c("mean", "scale", "shape")[components]
} else {
if (any(components))
warning("test with 'eta' together with other parameters is not implemented, only the test on 'eta' is performed")
eps <- 0.0001
estimate <- x$eta
null.value <- eta
if (type=="Wald") {
param <- estimate - null.value
if (x$lambda < eps)
warning("test and confidence intervals for 'eta' is not implemented for negative lambda")
hatvar <- AVAR.ML.eta.loggamma(mu=x$mu, sigma=x$sigma, lambda=x$lambda, prob=prob, eps=eps)
wstat <- sum(x$weights)*param^2/hatvar
df <- length(param)
pval <- pchisq(q=wstat, df=df, ncp=0, lower.tail = FALSE, log.p = FALSE)
cint <- estimate+c(-1,1)*qnorm(1-(1-conf.level)/2)*sqrt(hatvar)/sqrt(sum(x$weights))
}
names(wstat) <- "ww"
names(df) <- "df"
names(null.value) <- "mean(exp(X))"
}
type <- paste(ifelse(x$method=="ML","Classic","Weighted"), type, "Test based on", x$method, sep=" ")
attr(cint,"conf.level") <- conf.level
rval <- list(statistic=wstat, parameter=df, p.value=pval,
conf.int=cint, estimate=estimate, null.value=null.value,
alternative=alternative,
method=type, data.name=dname)
class(rval) <- "htest"
return(rval)
}
#############################################################
# print.loggammarob function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
######## PRINT
print.loggammarob <- function(x, digits = max(3, getOption("digits") - 3), ...) {
cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
cat("Location: ", format(x$mu, digits=digits))
cat(" Scale: ", format(x$sigma, digits=digits))
cat(" Shape: ", format(x$lambda, digits=digits))
cat(" E(exp(X)): ", format(x$eta, digits=digits))
cat("\n")
invisible(x)
}
#############################################################
# summary.loggammarob function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
########### SUMMARY
summary.loggammarob <- function(object, p=NULL, conf.level=0.95, prob=0.00001, ...) {
if ((object$method %in% c("WL", "oneWL", "ML"))) {
p <- p[ifelse(p > 0 & p < 1, TRUE, FALSE)]
nrep <- 4+length(p)
if (length(conf.level) == 1)
conf.level <- rep(conf.level, nrep)
if (length(conf.level) != nrep)
stop("'conf.level' must have length equal to 1 or 4+length(p)")
z <- qnorm(1-(1-conf.level)/2)
if (is.finite(object$eta)) {
asd <- sqrt(diag(ACOV.ML.loggamma(sigma=object$sigma,lambda=object$lambda,prob=prob)$cov))
asdeta <- sqrt(AVAR.ML.eta.loggamma(mu=object$mu, sigma=object$sigma, lambda=object$lambda, prob=prob, eps=0.0001, npoints=100000))
} else {
asd <- rep(NA,3)
asdeta <- NA
}
object$muconf.int <- object$mu+c(-1,1)*z[1]*asd[1]/sqrt(sum(object$weights))
object$sigmaconf.int <- object$sigma+c(-1,1)*z[2]*asd[2]/sqrt(sum(object$weights))
object$lambdaconf.int <- object$lambda+c(-1,1)*z[3]*asd[3]/sqrt(sum(object$weights))
object$etaconf.int <- object$eta+c(-1,1)*z[4]*asdeta/sqrt(sum(object$weights))
attr(object$muconf.int, "conf.level") <- conf.level[1]
attr(object$sigmaconf.int, "conf.level") <- conf.level[2]
attr(object$lambdaconf.int, "conf.level") <- conf.level[3]
attr(object$etaconf.int, "conf.level") <- conf.level[4]
object$muse <- asd[1]/sqrt(sum(object$weights))
object$sigmase <- asd[2]/sqrt(sum(object$weights))
object$lambdase <- asd[3]/sqrt(sum(object$weights))
object$etase <- asdeta/sqrt(sum(object$weights))
object$p <- p
if (!is.null(p)) {
object$q <- asdq <- rep(NA, length(p))
object$qconf.int <- matrix(NA, nrow=length(p), ncol=2)
for (i in 1:length(p)) {
object$q[i] <- qloggamma(p=p[i], mu=object$mu, sigma=object$sigma, lambda=object$lambda)
if (is.finite(object$eta))
asdq[i] <- sqrt(AVAR.ML.quantile.loggamma(p=p[i], mu=object$mu, sigma=object$sigma, lambda=object$lambda, prob=prob))
else
asdq[i] <- NA
object$qconf.int[i,] <- object$q[i]+c(-1,1)*z[4+i]*asdq[i]/sqrt(sum(object$weights))
}
##### object$qconf.int <- drop(qconf.int)
attr(object$qconf.int, "conf.level") <- conf.level[4+(1:length(p))]
object$qse <- asdq/sqrt(sum(object$weights))
}
}
object$call <- match.call()
class(object) <- "summary.loggammarob"
return(object)
}
#############################################################
# print.summary.loggammarob function
# Author: C. Agostinelli, A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
# Maintainer e-mail: claudio@unive.it
# Date: January, 1, 2013
# Version: 0.1
# Copyright (C) 2013 C. Agostinelli A. Marazzi,
# V.J. Yohai and A. Randriamiharisoa
#############################################################
#PRINT.SUMMARY
print.summary.loggammarob <- function(x, digits = max(3, getOption("digits") - 3), ...) {
cat("\nCall:\n",deparse(x$call),"\n\n",sep="")
cat("Location: ", format(x$mu, digits=digits))
if ((x$method %in% c("WL", "oneWL", "ML"))) {
cat(" s.e. ", format(x$muse, digits=digits), "\n")
cat("(", format(x$muconf.int[1L], digits=digits), ", ",
format(x$muconf.int[2L], digits=digits), ") \n")
cat(format(100 * attr(x$muconf.int, "conf.level")),
"percent confidence interval\n\n")
} else {
cat("\n")
}
cat("Scale: ", format(x$sigma, digits=digits))
if ((x$method %in% c("WL", "oneWL", "ML"))) {
cat(" s.e. ", format(x$sigmase, digits=digits), "\n")
cat("(", format(x$sigmaconf.int[1L], digits=digits), ", ",
format(x$sigmaconf.int[2L], digits=digits), ") \n")
cat(format(100 * attr(x$sigmaconf.int, "conf.level")),
"percent confidence interval\n\n")
} else {
cat("\n")
}
cat("Shape: ", format(x$lambda, digits=digits))
if ((x$method %in% c("WL", "oneWL", "ML"))) {
cat(" s.e. ", format(x$lambdase, digits=digits), "\n")
cat("(", format(x$lambdaconf.int[1L], digits=digits), ", ",
format(x$lambdaconf.int[2L], digits=digits), ") \n")
cat(format(100 * attr(x$lambdaconf.int, "conf.level")),
"percent confidence interval\n\n")
} else {
cat("\n")
}
cat("Mean(exp(X)): ", format(x$eta, digits=digits))
if ((x$method %in% c("WL", "oneWL", "ML"))) {
cat(" s.e. ", format(x$etase, digits=digits), "\n")
cat("(", format(x$etaconf.int[1L], digits=digits), ", ",
format(x$etaconf.int[2L], digits=digits), ") \n")
cat(format(100 * attr(x$etaconf.int, "conf.level")),
"percent confidence interval\n\n\n")
} else {
cat("\n")
}
######### Quantiles estimation and confidence interval
if (!is.null(x$p)) {
for (i in 1:length(x$p)) {
cat("Quantile of order ", format(x$p[i], digits=digits), ": ", format(x$q[i], digits=digits))
if ((x$method %in% c("WL", "oneWL", "ML"))) {
cat(" s.e. ", format(x$qse[i], digits=digits), "\n")
cat("(", format(x$qconf.int[i,1L], digits=digits), ", ",
format(x$qconf.int[i,2L], digits=digits), ") \n")
cat(format(100 * attr(x$qconf.int, "conf.level")[i]),
"percent confidence interval\n\n\n")
} else {
cat("\n")
}
}
}
if ((x$method %in% c("WL", "oneWL", "ML")))
cat("\n\n")
else
cat("\n\nConfidence intervals and Standard Errors are available only for 'WL', 'oneWL' and 'ML'\n\n\n")
if (x$method %in% c("WL", "oneWL", "WQTau"))
summarizeRobWeights(x$weights, digits = digits, ...)
invisible(x)
}
Try the robustloggamma 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.