Nothing
#' @title
#' Converting TL-moments to distribution parameters
#' @description
#' Converts TL-moments (or PWMs) to distribution parameters. By now, conversions for gev, gumbel, gpd, and ln3
#' are available. Important trimming options are calculated through known formulas (see references for
#' some of them), other options are calculated through a numerical optimization.
#'
#' @param x object returned by TLMoments (or PWMs, in which case TLMoments with no trimming is used).
#' @param distr character object defining the distribution. Supported types are
#' "gev", "gum", "gpd", and "ln3".
#' @param ... additional arguments.
#'
#' @return numeric vector, matrix, list, or data.frame of parameter estimates with
#' class \code{parameters}. The object contains the following attributes: \itemize{
#' \item \code{distribution}: a character indicating the used distribution
#' \item \code{source}: a list with background information (used function, data, n, formula, trimmings; mainly for internal purposes)
#' }
#' The attributes are hidden in the print-function for a clearer presentation.
#'
#' @references Elamir, E. A. H. (2010). Optimal choices for trimming in trimmed L-moment method. Applied Mathematical Sciences, 4(58), 2881-2890.
#' @references Fischer, S., Fried, R., & Schumann, A. (2015). Examination for robustness of parametric estimators for flood statistics in the context of extraordinary extreme events. Hydrology and Earth System Sciences Discussions, 12, 8553-8576.
#' @references Hosking, J. R. (1990). L-moments: analysis and estimation of distributions using linear combinations of order statistics. Journal of the Royal Statistical Society. Series B (Methodological), 105-124.
#' @references Hosking, J. R. M. (2007). Some theory and practical uses of trimmed L-moments. Journal of Statistical Planning and Inference, 137(9), 3024-3039.
#' @seealso \code{\link{PWMs}}, \code{\link{TLMoments}}, \code{\link{quantiles}}, \code{\link{summary.parameters}}, \code{\link{as.parameters}}. Built-in distributions: \code{\link{pgev}}, \code{\link{pgum}}, \code{\link{pgpd}}, \code{\link{pln3}}.
#'
#' @examples
#' xmat <- matrix(rgev(100, shape = .2), nc = 4)
#' xvec <- xmat[, 3]
#' xlist <- lapply(1L:ncol(xmat), function(i) xmat[, i])
#' xdat <- data.frame(
#' station = rep(letters[1:2], each = 50),
#' season = rep(c("S", "W"), 50),
#' hq = as.vector(xmat)
#' )
#'
#' # TLMoments-objects or PWMs-objects can be used. However, in case of PWMs
#' # simply the TLMoments(., leftrim = 0, rightrim = 0)-variant is used.
#'
#' parameters(PWMs(xvec), "gev")
#' tlm <- TLMoments(xvec, leftrim = 0, rightrim = 0)
#' parameters(tlm, "gev")
#'
#' tlm <- TLMoments(xmat, leftrim = 1, rightrim = 1)
#' parameters(tlm, "gum")
#'
#' tlm <- TLMoments(xlist)
#' parameters(tlm, "gpd")
#'
#' tlm <- TLMoments(xdat, hq ~ station, leftrim = 0, rightrim = 2)
#' parameters(tlm, "gev")
#'
#' tlm <- TLMoments(xdat, hq ~ station + season, leftrim = 0, rightrim = 2)
#' parameters(tlm, "gev")
#'
#'
#' # If no explicit formula is implemented, it is tried to calculate
#' # parameters numerically. The attribute source$param.computation.method
#' # indicates if this is the case.
#'
#' tlm <- TLMoments(rgum(200, loc = 5, scale = 2), leftrim = 1, rightrim = 4)
#' parameters(tlm, "gum")
#'
#' tlm <- TLMoments(rgev(200, loc = 10, scale = 5, shape = .4), leftrim = 2, rightrim = 2)
#' parameters(tlm, "gev")
#'
#' tlm <- TLMoments(rln3(200, loc = 3, scale = 1.5, shape = 2), leftrim = 0, rightrim = 1)
#' parameters(tlm, "ln3")
#'
#' # Numerical calculation is A LOT slower:
#' \dontrun{
#' system.time(replicate(500,
#' parameters(TLMoments(rgum(100, loc = 5, scale = 2), 1, 1), "gum")
#' ))[3]
#' system.time(replicate(500,
#' parameters(TLMoments(rgum(100, loc = 5, scale = 2), 1, 2), "gum")
#' ))[3]
#' }
#'
#' # Using magrittr
#' library(magrittr)
#'
#' TLMoments(rgpd(500, loc = 10, scale = 3, shape = .3), rightrim = 0) %>%
#' parameters("gpd")
#'
#' TLMoments(rgpd(500, loc = 10, scale = 3, shape = .3), rightrim = 0) %>%
#' parameters("gpd", u = 10)
#'
#' TLMoments(rgpd(500, loc = 10, scale = 3, shape = .3), rightrim = 1) %>%
#' parameters("gpd")
#'
#' TLMoments(rgpd(500, loc = 10, scale = 3, shape = .3), rightrim = 2) %>%
#' parameters("gpd")
#'
#' @export
parameters <- function(x, distr, ...) {
if (!inherits(x, c("TLMoments", "PWMs"))) stop("tlm has to be of class TLMoments or PWMs")
valid.distr <- c("gev", "gum", "gpd", "ln3")
if (!(distr %in% valid.distr))
stop(paste0("distr has to be one of ", paste(valid.distr, collapse = ", "), "."))
# Check for sufficient data:
if (distr == "gev") {
if (!all(c(1, 2, 3) %in% attr(x, "order"))) stop("Parameter calculation of \"gev\" needs at least the first three TL-moments")
} else if (distr == "gum") {
if (!all(c(1, 2) %in% attr(x, "order"))) stop("Parameter calculation of \"gum\" needs at least the first two TL-moments")
} else if (distr == "gpd" && "u" %in% names(list(...))) {
if (!all(c(1, 2) %in% attr(x, "order"))) stop("Parameter calculation of \"gpd\" needs at least the first two TL-moments")
} else if (distr == "gpd") {
if (!all(c(1, 2, 3) %in% attr(x, "order"))) stop("Parameter calculation of \"gpd\" needs at least the first three TL-moments or a specified threshold value (u = ...)")
} else if (distr == "ln3") {
if (!all(c(1, 2, 3) %in% attr(x, "order"))) stop("Parameter calculation of \"ln3\" needs at least the first three TL-moments")
}
UseMethod("parameters")
}
#' @title returnParameters
#' @description Sets attributions to parameters objects and returns them. This function is for internal use.
#' @param out -
#' @param distribution -
#' @param ... -
#' @return An object of class parameters.
returnParameters <- function(out, distribution, ...) {
class <- class(out)
args <- list(...)
# If no func attribute is set, set to
if (!exists("func", args)) args$func <- "parameters"
# If more than one func attributes are given, concatenate them
if (sum(names(args) == "func") >= 2) {
newfunc <- as.vector(unlist(args[names(args) == "func"]))
args$func <- NULL
args$func <- newfunc
}
# Attributes of parameters
# distribution
# source: func
# data (if calculated)
# input (if not calculated)
# n (if calculated)
# formula (if data is data.frame)
# trimmings (if coming from TLMoments)
# lambdas (if coming from TLMoments)
# tl.order (if coming from TLMoments)
# class: "parameters", class
attr(out, "distribution") <- distribution
attr(out, "source") <- args
class(out) <- c("parameters", class)
out
}
#' @describeIn parameters parameters for PWMs-object
#' @method parameters PWMs
#' @export
parameters.PWMs <- function(x, distr, ...) {
parameters.TLMoments(TLMoments(x), distr = distr)
}
#' @describeIn parameters parameters for TLMoments-object
#' @method parameters TLMoments
#' @export
parameters.TLMoments <- function(x, distr, ...) {
UseMethod("parameters.TLMoments", x$lambdas)
}
#' @method parameters.TLMoments numeric
#' @export
parameters.TLMoments.numeric <- function(x, distr, ...) {
leftrim <- attr(x, "leftrim")
rightrim <- attr(x, "rightrim")
out <- param.est(x$lambdas, leftrim, rightrim, distr, ...)
param.computation.method <- attr(out, "param.computation.method")
attr(out, "param.computation.method") <- NULL
do.call(returnParameters, c(
list(out = out, distribution = distr),
func = "parameters",
lambdas = list(x$lambdas),
trimmings = list(c(attr(x, "leftrim"), attr(x, "rightrim"))),
tl.order = list(attr(x, "order")),
param.computation.method = param.computation.method,
attr(x, "source")
))
}
#' @method parameters.TLMoments matrix
#' @export
parameters.TLMoments.matrix <- function(x, distr, ...) {
leftrim <- attr(x, "leftrim")
rightrim <- attr(x, "rightrim")
out <- lapply(1L:ncol(x$lambdas), function(i) {
param.est(x$lambdas[, i], leftrim, rightrim, distr, ...)
})
param.computation.method <- attr(out[[1]], "param.computation.method")
out <- simplify2array(out)
do.call(returnParameters, c(
list(out = out, distribution = distr),
func = "parameters",
lambdas = list(x$lambdas),
trimmings = list(c(attr(x, "leftrim"), attr(x, "rightrim"))),
tl.order = list(attr(x, "order")),
param.computation.method = param.computation.method,
attr(x, "source")
))
}
#' @method parameters.TLMoments list
#' @export
parameters.TLMoments.list <- function(x, distr, ...) {
leftrim <- attr(x, "leftrim")
rightrim <- attr(x, "rightrim")
out <- lapply(seq_along(x$lambdas), function(i) {
param.est(x$lambdas[[i]], leftrim, rightrim, distr, ...)
})
param.computation.method <- attr(out[[1]], "param.computation.method")
# Delete attributes...
for (i in seq_along(out)) {
attr(out[[i]], "param.computation.method") <- NULL
}
do.call(returnParameters, c(
list(out = out, distribution = distr),
func = "parameters",
lambdas = list(x$lambdas),
trimmings = list(c(attr(x, "leftrim"), attr(x, "rightrim"))),
tl.order = list(attr(x, "order")),
param.computation.method = param.computation.method,
attr(x, "source")
))
}
#' @method parameters.TLMoments data.frame
#' @export
parameters.TLMoments.data.frame <- function(x, distr, ...) {
leftrim <- attr(x, "leftrim")
rightrim <- attr(x, "rightrim")
nam <- getFormulaSides(attr(x, "source")$formula)
lambdas <- t(x$lambdas[!(names(x$lambdas) %in% nam$rhs)])
out <- lapply(1L:ncol(lambdas), function(i) {
param.est(lambdas[, i], leftrim, rightrim, distr, ...)
})
param.computation.method <- attr(out[[1]], "param.computation.method")
out <- simplify2array(out)
out <- cbind(x$lambdas[nam$rhs], as.data.frame(t(out)))
do.call(returnParameters, c(
list(out = out, distribution = distr),
func = "parameters",
lambdas = list(x$lambdas),
trimmings = list(c(attr(x, "leftrim"), attr(x, "rightrim"))),
tl.order = list(attr(x, "order")),
param.computation.method = param.computation.method,
attr(x, "source")
))
}
#' @export
print.parameters <- function(x, ...) {
if (inherits(x, "data.frame")) {
print.data.frame(x)
return(invisible(x))
}
tmp <- x
attributes(tmp) <- NULL
dim(tmp) <- dim(x)
names(tmp) <- names(x)
dimnames(tmp) <- dimnames(x)
print(tmp)
invisible(x)
}
param.est <- function(lambdas, leftrim, rightrim, distr, ...) {
if (length(lambdas) < 2 && distr == "gpd" && "u" %in% names(list(...)))
stop("Insufficient L-moments to calculate parameters. ")
if (length(lambdas) < 3 && distr == "gpd" && !("u" %in% names(list(...))))
stop("Insufficient L-moments to calculate parameters. ")
if (length(lambdas) < 3 && distr %in% c("gev", "ln3"))
stop("Insufficient L-moments to calculate parameters. ")
if (length(lambdas) < 2 && distr == "gum")
stop("Insufficient L-moments to calculate parameters. ")
switch(
distr,
gev = gev.est(lambdas[1], lambdas[2], lambdas[3]/lambdas[2], leftrim, rightrim),
gum = gum.est(lambdas[1], lambdas[2], leftrim, rightrim),
gpd = gpd.est(lambdas[1], lambdas[2], lambdas[3]/lambdas[2], leftrim, rightrim, ...),
ln3 = ln3.est(lambdas[1], lambdas[2], lambdas[3]/lambdas[2], leftrim, rightrim),
stop("Not implemented. ")
)
}
#### GUM ####
gum.est <- function(l1, l2, s, t) {
switch(paste(s, t, sep = "-"),
`0-0` = gum.tl00.est(l1, l2),
`0-1` = gum.tl01.est(l1, l2),
`1-0` = gum.tl10.est(l1, l2),
`1-1` = gum.tl11.est(l1, l2),
gum.numerical.est(l1, l2, s, t))
}
# TL(0,0)=L
gum.tl00.est <- function(l1, l2) {
scale <- l2 / log(2)
loc <- l1 - 0.577 * scale
out <- setNames(c(loc, scale), c("loc", "scale"))
attr(out, "param.computation.method") <- "analytical"
return(out)
}
# TL(0,1), TL(1,0), TL(1,1) from Elamir, 2010: Optimal Choices for Trimming in Trimmed L-moment Method
gum.tl01.est <- function(l1, l2) {
scale <- l2 / 0.431
loc <- l1 + 0.116 * scale
out <- setNames(c(loc, scale), c("loc", "scale"))
attr(out, "param.computation.method") <- "analytical"
return(out)
}
gum.tl10.est <- function(l1, l2) {
scale <- l2 / 0.607
loc <- l1 - 1.269 * scale
out <- setNames(c(loc, scale), c("loc", "scale"))
attr(out, "param.computation.method") <- "analytical"
return(out)
}
gum.tl11.est <- function(l1, l2) {
scale <- l2 / 0.353
loc <- l1 - .459 * scale
out <- setNames(c(loc, scale), c("loc", "scale"))
attr(out, "param.computation.method") <- "analytical"
return(out)
}
gum.numerical.est <- function(l1, l2, s, t) {
#warning("Calculating numerical solution. ")
# scale
f1 <- function(scale) calcTLMom(2, s, t, qgum, loc = 0, scale = scale)[2] - l2
i <- 0; while (sign(f1(10^-i)) == sign(f1(10^i))) i <- i+1
u <- uniroot(f1, c(10^-i, 10^i))
scale <- u$root
# loc
f2 <- function(loc) calcTLMom(2, s, t, qgum, loc = loc, scale = scale)[1] - l1
i <- 0; while (sign(f2(-10^i)) == sign(f2(10^i))) i <- i+1
u <- uniroot(f2, c(-10^i, 10^i))
loc <- u$root
out <- setNames(c(loc, scale), c("loc", "scale"))
attr(out, "param.computation.method") <- "numerical"
return(out)
}
#### GEV ####
gev.est <- function(l1, l2, t3, s, t) {
switch(paste(s, t, sep = "-"),
`0-0` = gev.tl00.est(l1, l2, t3),
`0-1` = gev.tl01.est(l1, l2, t3),
`0-2` = gev.tl02.est(l1, l2, t3),
`1-1` = gev.tl11.est(l1, l2, t3),
`1-0` = gev.tl10.est(l1, l2, t3),
gev.numerical.est(l1, l2, t3, s, t))
}
# TL(0,0)=L der GEV
gev.tl00.est <- function(l1, l2, t3) {
z <- 2/(3+t3) - log(2)/log(3)
k <- 7.8590 * z + 2.9554 * z^2
# Hosking (1990)
gk <- gamma(1+k)
al <- l2*k/((1-2^-k)*gk)
xi <- l1 + al*(gk-1)/k
out <- setNames(c(xi, al, -k), c("loc", "scale", "shape"))
attr(out, "param.computation.method") <- "analytical"
return(out)
}
# TL(0,1) der GEV
gev.tl01.est <- function(l1, l2, t3) {
z <- 1/(2+t3) * 10/9 - (2*log(2)-log(3))/(3*log(3)-2*log(4))
k <- 8.567394*z - 0.675969*z^2
gk <- gamma(k)
V3 <- (1/3)^k
V2 <- (1/2)^k
al <- l2 * 2/3 * 1/gk * 1/(V3 - 2*V2 + 1)
xi <- l1 - al/k - al*gk*(V2 - 2)
out <- setNames(c(xi, al, -k), c("loc", "scale", "shape"))
attr(out, "param.computation.method") <- "analytical"
return(out)
}
# TL(0,2) der GEV
gev.tl02.est <- function(l1, l2, t3) {
z <- (t3+5/3) * 6/5 - (3*log(5)-8*log(4)+6*log(3))/(log(4)-3*log(3)+3*log(2))
k <- -2.468959*z + 1.130074*z^2 - 0.635912*z^3
gk <- gamma(k)
V4 <- (1/4)^k
V3 <- (1/3)^k
V2 <- (1/2)^k
al <- l2 * 2 * 1/gk * 1/(-4*V4 + 12*V3 - 12*V2 + 4)
xi <- l1 - al/k - al*gk*(-V3 + 3*V2 - 3)
out <- setNames(c(xi, al, -k), c("loc", "scale", "shape"))
attr(out, "param.computation.method") <- "analytical"
return(out)
}
# TL(1,0) der GEV
gev.tl10.est <- function(l1, l2, t3) {
z <- (t3-40/3) * 9/20 - (log(3)-log(4))/(log(2)-log(3))
k <- -9.066941*z - 3.374925*z^2 - 0.303208*z^3
gk <- gamma(k)
V3 <- (1/3)^k
V2 <- (1/2)^k
al <- l2 * 2 * 1/gk * 1/(3*V2 - 3*V3)
xi <- l1 - al/k + al*gk*V2
out <- setNames(c(xi, al, -k), c("loc", "scale", "shape"))
attr(out, "param.computation.method") <- "analytical"
return(out)
}
# TL(1,1) der GEV
gev.tl11.est <- function(l1, l2, t3) {
z <- t3 * 9/20 + (log(3)-2*log(4)+log(5))/(log(2)-2*log(3)+log(4))
k <- 25.3171*z - 91.5507*z^2 + 110.0626*z^3 - 46.5518*z^4
gk <- gamma(k)
V4 <- (1/4)^k
V3 <- (1/3)^k
V2 <- (1/2)^k
al <- l2 * 1/gk * 1/(3*V2 - 6*V3 + 3*V4)
xi <- l1 - al/k - al*gk*(-3*V2 + 2*V3)
out <- setNames(c(xi, al, -k), c("loc", "scale", "shape"))
attr(out, "param.computation.method") <- "analytical"
return(out)
}
gev.numerical.est <- function(l1, l2, t3, s, t) {
#warning("Calculating numerical solution. ")
# shape
f1 <- function(shape) {
a <- calcTLMom(3, s, t, qgev, loc = 0, scale = 1, shape = shape)
a[3]/a[2] - t3
}
i <- 0; while (sign(f1(-10^i)) == sign(f1(10^i))) i <- i+1
u <- uniroot(f1, c(-10^i, 10^i))
shape <- u$root
# scale (positiv)
f2 <- function(scale) calcTLMom(2, s, t, qgev, loc = 0, scale = scale, shape = shape)[2] - l2
i <- 0; while (sign(f2(10^-i)) == sign(f2(10^i))) i <- i+1
u <- uniroot(f2, c(10^-i, 10^i))
scale <- u$root
# loc
f3 <- function(loc) calcTLMom(1, s, t, qgev, loc = loc, scale = scale, shape = shape)[1] - l1
i <- 0; while (sign(f3(-10^i)) == sign(f3(10^i))) i <- i+1
u <- uniroot(f3, c(-10^i, 10^i))
loc <- u$root
out <- setNames(c(loc, scale, shape), c("loc", "scale", "shape"))
attr(out, "param.computation.method") <- "numerical"
return(out)
}
#### LN3 ####
ln3.est <- function(l1, l2, t3, s, t) {
switch(paste(s, t, sep = "-"),
`0-0` = ln3.tl00.est(l1, l2, t3),
ln3.numerical.est(l1, l2, t3, s, t))
}
ln3.tl00.est <- function(l1, l2, t3) {
erf <- function(x) 2 * pnorm(x * sqrt(2)) - 1
# 1. step, T3 <=> shape (positiv)
f1 <- function(s) 6/sqrt(pi) * 1/erf(s/2) * integrate(function(x) erf(x/sqrt(3) * exp(-x^2)), 0, s/2)$val - t3
i <- 1; while (sign(f1(10^-i)) == sign(f1(sqrt(exp(i))))) i <- i+1
u <- uniroot(f1, c(10^-i, sqrt(exp(i))), tol = .Machine$double.eps^0.5, check.conv = TRUE)
shape <- u$root
# 2. scale and loc
scale <- log(l2 / erf(shape/2)) - shape^2/2
loc <- l1 - exp(scale + shape^2/2)
out <- setNames(c(loc, scale, shape), c("loc", "scale", "shape"))
attr(out, "param.computation.method") <- "analytical"
return(out)
}
ln3.numerical.est <- function(l1, l2, t3, s, t) {
#warning("Calculating numerical solution. ")
# 1. step, T3 <=> shape (positiv)
f1 <- function(shape) {
a <- calcTLMom(3, s, t, qln3, loc = 0, scale = 1, shape = shape)
a[3]/a[2] - t3
}
i <- 1; while (sign(f1(10^-i)) == sign(f1(sqrt(exp(i))))) i <- i+1
u <- uniroot(f1, c(10^-i, sqrt(exp(i))))
shape <- u$root
# 2. step, L2 <=> scale
f2 <- function(scale) calcTLMom(2, s, t, qln3, loc = 0, scale = scale, shape = shape)[2] - l2
i <- 0; while (sign(f2(-10^i)) == sign(f2(10^i))) i <- i+1
u <- uniroot(f2, c(-10^i, 10^i))
scale <- u$root
# 3. step L1 <=> loc
f3 <- function(loc) calcTLMom(1, s, t, qln3, loc = loc, scale = scale, shape = shape)[1] - l1
i <- 0; while (sign(f3(-10^i)) == sign(f3(10^i))) i <- i+1
u <- uniroot(f3, c(-10^i, 10^i))
loc <- u$root
out <- setNames(c(loc, scale, shape), c("loc", "scale", "shape"))
attr(out, "param.computation.method") <- "numerical"
return(out)
}
#### GPD ####
gpd.est <- function(l1, l2, t3, s, t, ...) {
switch(paste(s, t, sep = "-"),
`0-0` = gpd.tl00.est(l1, l2, t3, ...),
`0-1` = gpd.tl01.est(l1, l2, t3, ...),
`0-2` = gpd.tl02.est(l1, l2, t3, ...),
gpd.numerical.est(l1, l2, t3, s, t, ...))
}
# GPD
gpd.tl00.est <- function(l1, l2, t3, u = NULL) {
if (!is.null(u) & is.numeric(u)) {
shape <- -(l1-u)/l2 + 2
scale <- (l1-u) * (-shape+1)
loc <- u
} else {
shape <- (3*t3 - 1)/(t3 + 1)
scale <- l2 * (shape-1) * (shape-2)
loc <- l1 + scale/(shape-1)
}
out <- setNames(c(loc, scale, shape), c("loc", "scale", "shape"))
attr(out, "param.computation.method") <- "analytical"
return(out)
}
gpd.tl01.est <- function(l1, l2, t3, u = NULL) {
if (!is.null(u) & is.numeric(u)) {
# from Hosking, 2007.
shape <- -3/2 * (l1-u)/l2 + 3
scale <- (l1-u) * (-shape+2)
loc <- u
} else {
# from Fischer, 2015.
shape <- (36*t3 - 8)/(9*t3 + 8)
scale <- 2/3 * l2 * (shape-2) * (shape-3)
loc <- l1 + scale/(shape-2)
}
out <- setNames(c(loc, scale, shape), c("loc", "scale", "shape"))
attr(out, "param.computation.method") <- "analytical"
return(out)
}
gpd.tl02.est <- function(l1, l2, t3, u = NULL) {
if (!is.null(u) & is.numeric(u)) {
# from Hosking, 2007.
shape <- -2 * (l1-u)/l2 + 4
scale <- (l1-u) * (-shape+3)
loc <- u
} else {
# from Fischer, 2015.
shape <- (30*t3 - 5)/(6*t3 + 5)
scale <- 1/2 * l2 * (shape-3) * (shape-4)
loc <- l1 + scale/(shape-3)
}
out <- setNames(c(loc, scale, shape), c("loc", "scale", "shape"))
attr(out, "param.computation.method") <- "analytical"
return(out)
}
gpd.numerical.est <- function(l1, l2, t3, s, t, u = NULL) {
#warning("Calculating numerical solution. ")
if (!is.null(u) & is.numeric(u)) {
# shape
f1u <- function(shape, u) {
a <- calcTLMom(2, s, t, qgpd, loc = u, scale = 1, shape = shape)
a[2] - l2
}
i <- 0; while (sign(f1u(-10^i, u)) == sign(f1u(10^i, u))) i <- i+1
ur <- uniroot(f1u, c(-10^i, 10^i), u = u)
shape <- ur$root
# scale
f2u <- function(scale, u) calcTLMom(1, s, t, qgpd, loc = u, scale = scale, shape = shape)[1] - l1
i <- 0; while (sign(f2u(10^-i, u)) == sign(f2u(10^i, u))) i <- i+1
ur <- uniroot(f2u, c(10^-i, 10^i), u = u)
scale <- ur$root
# loc
loc <- u
} else {
# shape
f1 <- function(shape) {
a <- calcTLMom(3, s, t, qgpd, loc = 0, scale = 1, shape = shape)
a[3]/a[2] - t3
}
i <- 0; while (sign(f1(-10^i)) == sign(f1(10^i))) i <- i+1
ur <- uniroot(f1, c(-10^i, 10^i))
shape <- ur$root
# scale
f2 <- function(scale) calcTLMom(2, s, t, qgpd, loc = 0, scale = scale, shape = shape)[2] - l2
i <- 0; while (sign(f2(10^-i)) == sign(f2(10^i))) i <- i+1
ur <- uniroot(f2, c(10^-i, 10^i))
scale <- ur$root
# loc
f3 <- function(loc) calcTLMom(1, s, t, qgpd, loc = loc, scale = scale, shape = shape)[1] - l1
i <- 0; while (sign(f3(-10^i)) == sign(f3(10^i))) i <- i+1
ur <- uniroot(f3, c(-10^i, 10^i))
loc <- ur$root
}
out <- setNames(c(loc, scale, shape), c("loc", "scale", "shape"))
attr(out, "param.computation.method") <- "numerical"
return(out)
}
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.