# R/covfunctions.R In TLMoments: Calculate TL-Moments and Convert Them to Distribution Parameters

#### Defines functions CovParamtoQuanGEV_02GEV_01GEV_00CovTLtoParaparametricPWMCovPWMCov_GEVvCovLambdaToTau

```# This file contains internal functions used to calculate covariances
# of estimations of PWMs, TLMoments, Parameters, and Quantiles.

# Generate matrix for conversion of lambda-cov to tau-cov
CovLambdaToTau  <- function(l) {

K <- length(l)

if (K > 3) {
A <- rbind(
c(-l[2]/l[1]^2, rep(0, K-2)),
c(1/l[1], -l[-c(1,2)]/l[2]^2),
cbind(0, diag(rep(1/l[2], K-2)))
)
} else if (K == 3) {
A <- rbind(
c(-l[2]/l[1]^2, rep(0, K-2)),
c(1/l[1], -l[-c(1,2)]/l[2]^2),
cbind(0, 1/l[2])
)
} else {
A <- rbind(-l[2]/l[1]^2, 1/l[1])
}

A
}

v <- function(r, delta, scale, shape) {
G <- function(x, a, b = a) hypergeo::hypergeo(-b, -2*a, 1-a, -x)

if (delta == 0) {
out <- scale^2 * (-shape)^(-2) * (r+1)^(2*shape) * ( gamma(1-2*shape) * G(r/(r+1), shape) - gamma(1-shape)^2 )
} else if (delta == 1) {
out <- .5 * scale^2 * (-shape)^(-2) * ( (r+2)^(2*shape) * gamma(1-2*shape) * G(r/(r+2), shape) + (r+1)^shape * ( (r+1)^shape - 2*(r+2)^shape ) * gamma(1-shape)^2 )
} else {
s <- delta
out <- .5 * scale^2 * (-shape)^(-2) * gamma(1-2*shape) * (
(2*r+s+1)^(2*shape) * G(-(r+s)/(2*r+s+1), shape, -1) -
(2*r+s+2)^(2*shape) * G(-(r+s+1)/(2*r+s+2), shape, -1) +
(2*r+s+1)^(2*shape) * G(-(r)/(2*r+s+1), shape, -1) -
(2*r+s+2)^(2*shape) * G(-(r+1)/(2*r+s+2), shape, -1)
)
}
Re(out)
}

# @return Symmetric Matrix of dimensions determined by \code{betas}.
PWMCov_GEV <- function(betas = 0:3, scale, shape) {
if (abs(shape) <= 1e-06) shape <- 1e-6
mirs <- min(betas)
mars <- max(betas)
out <- matrix(NA, nrow = mars+1, ncol = mars+1)
for (r in mirs:mars) for (s in r:mars) {
if (r == s) {
out[r+1, s+1] <- v(r, s-r, scale, shape)
} else {
out[r+1, s+1] <- out[s+1, r+1] <- v(r, s-r, scale, shape)
}
}
out <- out[betas+1, betas+1, drop = FALSE]

colnames(out) <- rownames(out) <- paste0("beta", betas)
out
}

parametricPWMCov <- function(distribution, order = 0:3, ...) {
p <- list(...)
if (distribution == "gev") {
if (!all(c("scale", "shape") %in% names(p))) stop("scale and shape must be given. ")

return(PWMCov_GEV(order, p[["scale"]], p[["shape"]]))
}

stop("No functions available for given distribution. ")
}

# Matrices to calculate Cov of parameter estimations out of TL(0,0)- or TL(0,1)-Covs.
# Need L2 and L3 with corresponding trimmings.
CovTLtoPara <- function(distr, l2, l3, leftrim, rightrim) {
if (distr != "gev") stop("Only GEV for now. ")

if (distr == "gev") {
if (leftrim == 0 & rightrim == 0) {
GEV_00(l2, l3)
} else if (leftrim == 0 & rightrim == 1) {
GEV_01(l2, l3)
} else if (leftrim == 0 & rightrim == 2) {
GEV_02(l2, l3)
}
} else {
stop("Unsupported distr or trimming. ")
}
}

GEV_00 <- function(l2, l3) {

a <- -7.859
b <- -2.9554
kappa <- log(2)/log(3)

zeta <- 1/(l3+3*l2)
eta <- 2*l2*zeta - kappa
theta <- 2*zeta - 6*l2*zeta^2

# shape
g1 <- 0
g2 <- -2*l3*(2*b*kappa*l3-a*l3+6*b*kappa*l2-4*b*l2-3*a*l2)*zeta^3
g3 <- 2*l2*(2*b*kappa*l3-a*l3+6*b*kappa*l2-4*b*l2-3*a*l2)*zeta^3
r3 <- c(g1, g2, g3)

# scale
g2 <- (log(2)*(-b*eta^2-a*eta)*2^(b*eta^2+a*eta)*l2*(2*b*eta*theta+a*theta))/((1-2^(b*eta^2+a*eta))^2*gamma(-b*eta^2-a*eta+1))-((-b*eta^2-a*eta)*digamma(-b*eta^2-a*eta+1)*l2*(-2*b*eta*theta-a*theta))/((1-2^(b*eta^2+a*eta))*gamma(-b*eta^2-a*eta+1))+(l2*(-2*b*eta*theta-a*theta))/((1-2^(b*eta^2+a*eta))*gamma(-b*eta^2-a*eta+1))+(-b*eta^2-a*eta)/((1-2^(b*eta^2+a*eta))*gamma(-b*eta^2-a*eta+1))
g3 <- -((-b*eta^2-a*eta)*digamma(-b*eta^2-a*eta+1)*l2*(4*b*eta*l2*zeta^2+2*a*l2*zeta^2))/((1-2^(b*eta^2+a*eta))*gamma(-b*eta^2-a*eta+1))+(l2*(4*b*eta*l2*zeta^2+2*a*l2*zeta^2))/((1-2^(b*eta^2+a*eta))*gamma(-b*eta^2-a*eta+1))+(log(2)*(-b*eta^2-a*eta)*2^(b*eta^2+a*eta)*l2*(-4*b*eta*l2*zeta^2-2*a*l2*zeta^2))/((1-2^(b*eta^2+a*eta))^2*gamma(-b*eta^2-a*eta+1))
r2 <- c(0, g2, g3)

# loc
g2 <- -(log(2)*(-b*eta^2-a*eta)*2^(b*eta^2+a*eta)*(gamma(-b*eta^2-a*eta+1)-1)*l2*(2*b*eta*theta+a*theta))/((b*eta^2+a*eta)*(1-2^(b*eta^2+a*eta))^2*gamma(-b*eta^2-a*eta+1))+((-b*eta^2-a*eta)*(gamma(-b*eta^2-a*eta+1)-1)*l2*(2*b*eta*theta+a*theta))/((b*eta^2+a*eta)^2*(1-2^(b*eta^2+a*eta))*gamma(-b*eta^2-a*eta+1))+
((-b*eta^2-a*eta)*digamma(-b*eta^2-a*eta+1)*(gamma(-b*eta^2-a*eta+1)-1)*l2*(-2*b*eta*theta-a*theta))/((b*eta^2+a*eta)*(1-2^(b*eta^2+a*eta))*gamma(-b*eta^2-a*eta+1))-((gamma(-b*eta^2-a*eta+1)-1)*l2*(-2*b*eta*theta-a*theta))/((b*eta^2+a*eta)*(1-2^(b*eta^2+a*eta))*gamma(-b*eta^2-a*eta+1))-((-b*eta^2-a*eta)*digamma(-b*eta^2-a*eta+1)*l2*(-2*b*eta*theta-a*theta))/((b*eta^2+a*eta)*(1-2^(b*eta^2+a*eta)))-
((-b*eta^2-a*eta)*(gamma(-b*eta^2-a*eta+1)-1))/((b*eta^2+a*eta)*(1-2^(b*eta^2+a*eta))*gamma(-b*eta^2-a*eta+1))

g3 <- ((-b*eta^2-a*eta)*digamma(-b*eta^2-a*eta+1)*(gamma(-b*eta^2-a*eta+1)-1)*l2*(4*b*eta*l2*zeta^2+2*a*l2*zeta^2))/((b*eta^2+a*eta)*(1-2^(b*eta^2+a*eta))*gamma(-b*eta^2-a*eta+1))-((gamma(-b*eta^2-a*eta+1)-1)*l2*(4*b*eta*l2*zeta^2+2*a*l2*zeta^2))/((b*eta^2+a*eta)*(1-2^(b*eta^2+a*eta))*gamma(-b*eta^2-a*eta+1))-
((-b*eta^2-a*eta)*digamma(-b*eta^2-a*eta+1)*l2*(4*b*eta*l2*zeta^2+2*a*l2*zeta^2))/((b*eta^2+a*eta)*(1-2^(b*eta^2+a*eta)))-(log(2)*(-b*eta^2-a*eta)*2^(b*eta^2+a*eta)*(gamma(-b*eta^2-a*eta+1)-1)*l2*(-4*b*eta*l2*zeta^2-2*a*l2*zeta^2))/((b*eta^2+a*eta)*(1-2^(b*eta^2+a*eta))^2*gamma(-b*eta^2-a*eta+1))+
((-b*eta^2-a*eta)*(gamma(-b*eta^2-a*eta+1)-1)*l2*(-4*b*eta*l2*zeta^2-2*a*l2*zeta^2))/((b*eta^2+a*eta)^2*(1-2^(b*eta^2+a*eta))*gamma(-b*eta^2-a*eta+1))

r1 <- c(1, g2, g3)

rbind(r1, r2, r3)
}

GEV_01 <- function(l2, l3) {
kappa <- (2*log(2)-log(3))/(3*log(3)-2*log(4))
a <- -8.567394
b <- 0.675969

zeta <- 10/(9*(l3+2*l2))
eta <- 10*l2/(9*(l3+2*l2))
theta <- 20*l2/(9*(l3+2*l2)^2)

r3 <- (a + 2 * b * eta - kappa) * 10/(9*(l3+2*l2)^2) * c(
0, l3, -l2
)

g1 <- 0
g2 <- -(2*l2*(log(3)*3^(a*(eta-kappa)+b*(eta-kappa)^2)*(2*b*(eta-kappa)*(zeta-theta)+a*(zeta-theta))-log(2)*2^(a*(eta-kappa)+b*(eta-kappa)^2+1)*(2*b*(eta-kappa)*(zeta-theta)+a*(zeta-theta))))/(3*gamma(-a*(eta-kappa)-b*(eta-kappa)^2)*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1)^2)-(2*digamma(-a*(eta-kappa)-b*(eta-kappa)^2)*l2*(-2*b*(eta-kappa)*(zeta-theta)-a*(zeta-theta)))/(3*gamma(-a*(eta-kappa)-b*(eta-kappa)^2)*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1))+    2/(3*gamma(-a*(eta-kappa)-b*(eta-kappa)^2)*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1))
g3 <- -(2*l2*(log(3)*3^(a*(eta-kappa)+b*(eta-kappa)^2)*(-(20*b*(eta-kappa)*l2)/(9*(l3+2*l2)^2)-(10*a*l2)/(9*(l3+2*l2)^2))-log(2)*2^(a*(eta-kappa)+b*(eta-kappa)^2+1)*(-(20*b*(eta-kappa)*l2)/(9*(l3+2*l2)^2)-(10*a*l2)/(9*(l3+2*l2)^2))))/(3*gamma(-a*(eta-kappa)-b*(eta-kappa)^2)*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1)^2)  -(2*digamma(-a*(eta-kappa)-b*(eta-kappa)^2)*l2*((20*b*(eta-kappa)*l2)/(9*(l3+2*l2)^2)+(10*a*l2)/(9*(l3+2*l2)^2)))/(3*gamma(-a*(eta-kappa)-b*(eta-kappa)^2)*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1))

r2 <- c(g1, g2, g3)

g1 <- 1
g2 <- (2*(2^(a*(eta-kappa)+b*(eta-kappa)^2)-2)*l2*(log(3)*3^(a*(eta-kappa)+b*(eta-kappa)^2)*(2*b*(eta-kappa)*(zeta-theta)+a*(zeta-theta))-log(2)*2^(a*(eta-kappa)+b*(eta-kappa)^2+1)*(2*b*(eta-kappa)*(zeta-theta)+a*(zeta-theta)))  )/(3*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1)^2)-  (2*l2*(log(3)*3^(a*(eta-kappa)+b*(eta-kappa)^2)*(2*b*(eta-kappa)*(zeta-theta)+a*(zeta-theta))-log(2)*2^(a*(eta-kappa)+b*(eta-kappa)^2+1)*(2*b*(eta-kappa)*(zeta-theta)+a*(zeta-theta))))/(3*gamma(-a*(eta-kappa)-b*(eta-kappa)^2)*(a*(eta-kappa)+b*(eta-kappa)^2)*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1)^2)-  (log(2)*2^(a*(eta-kappa)+b*(eta-kappa)^2+1)*l2*(2*b*(eta-kappa)*(zeta-theta)+a*(zeta-theta)))/(3*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1))-  (2*l2*(2*b*(eta-kappa)*(zeta-theta)+a*(zeta-theta)))/(3*gamma(-a*(eta-kappa)-b*(eta-kappa)^2)*(a*(eta-kappa)+b*(eta-kappa)^2)^2*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1))-  (2*digamma(-a*(eta-kappa)-b*(eta-kappa)^2)*l2*(-2*b*(eta-kappa)*(zeta-theta)-a*(zeta-theta)))/(3*gamma(-a*(eta-kappa)-b*(eta-kappa)^2)*(a*(eta-kappa)+b*(eta-kappa)^2)*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1))-  (2*(2^(a*(eta-kappa)+b*(eta-kappa)^2)-2))/(3*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1))+    2/(3*gamma(-a*(eta-kappa)-b*(eta-kappa)^2)*(a*(eta-kappa)+b*(eta-kappa)^2)*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1))
g3 <- (2*(2^(a*(eta-kappa)+b*(eta-kappa)^2)-2)*l2*(log(3)*3^(a*(eta-kappa)+b*(eta-kappa)^2)*(-(20*b*(eta-kappa)*l2)/(9*(l3+2*l2)^2)-(10*a*l2)/(9*(l3+2*l2)^2))-log(2)*                                                       2^(a*(eta-kappa)+b*(eta-kappa)^2+1)*(-(20*b*(eta-kappa)*l2)/(9*(l3+2*l2)^2)-(10*a*l2)/(9*(l3+2*l2)^2))))/(3*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1)^2)-  (2*l2*(log(3)*3^(a*(eta-kappa)+b*(eta-kappa)^2)*(-(20*b*(eta-kappa)*l2)/(9*(l3+2*l2)^2)-(10*a*l2)/(9*(l3+2*l2)^2))-log(2)*2^(a*(eta-kappa)+b*(eta-kappa)^2+1)*(-(20*b*(eta-kappa)*l2)/(9*(l3+2*l2)^2)-(10*a*l2)/(9*(l3+2*l2)^2))))/(3*gamma(-a*(eta-kappa)-b*(eta-kappa)^2)*(a*(eta-kappa)+b*(eta-kappa)^2)*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1)^2)  -(2*digamma(-a*(eta-kappa)-b*(eta-kappa)^2)*l2*((20*b*(eta-kappa)*l2)/(9*(l3+2*l2)^2)+(10*a*l2)/(9*(l3+2*l2)^2)))/(3*gamma(-a*(eta-kappa)-b*(eta-kappa)^2)*(a*(eta-kappa)+b*(eta-kappa)^2)*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1))-  (log(2)*2^(a*(eta-kappa)+b*(eta-kappa)^2+1)*l2*(-(20*b*(eta-kappa)*l2)/(9*(l3+2*l2)^2)-(10*a*l2)/(9*(l3+2*l2)^2)))/(3*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1))-  (2*l2*(-(20*b*(eta-kappa)*l2)/(9*(l3+2*l2)^2)-(10*a*l2)/(9*(l3+2*l2)^2)))/(3*gamma(-a*(eta-kappa)-b*(eta-kappa)^2)*(a*(eta-kappa)+b*(eta-kappa)^2)^2*(-2^(a*(eta-kappa)+b*(eta-kappa)^2+1)+3^(a*(eta-kappa)+b*(eta-kappa)^2)+1))

r1 <- c(g1, g2, g3)

rbind(r1, r2, r3)
}

GEV_02 <- function(l2, l3) {
kappa <- (3*log(5)-8*log(4)+6*log(3)) / (log(4)-3*log(3)+3*log(2)) - 2
a <- 2.468959
b <- -1.130074
c <- 0.635912

zeta <- 6/5 * l3/l2
eta <- zeta - kappa
mu <- c*eta^3 + b*eta^2 + a*eta
theta <- 18*c*eta^2/(5*l2) + 12*b*eta/(5*l2) + 6*a/(5*l2)

r3 <- theta * c(
0,
-l3/l2,
1
)

g1 <- 0
g2 <- (((-l3*3^(mu+1)+l3*4^mu+3*l3*2^mu-l3)*digamma(-mu)+log(3)*l3*3^(mu+1)-log(4)*l3*4^mu-3*log(2)*l3*2^mu)*theta+3^(mu+1)-4^mu-3*2^mu+1)/((2*3^(2*mu+2)+9*2^(2*mu+1)+4^mu*(3*2^(mu+2)-4*3^(mu+1)-4)-3*2^(mu+2)+3^mu*(12-9*2^(mu+2))+2*4^(2*mu)+2)*gamma(-mu))
g3 <- -(((-l2*3^(mu+1)+l2*4^mu+3*l2*2^mu-l2)*digamma(-mu)+log(3)*l2*3^(mu+1)-log(4)*l2*4^mu-3*log(2)*l2*2^mu)*theta)/((2*3^(2*mu+2)+9*2^(2*mu+1)+4^mu*(3*2^(mu+2)-4*3^(mu+1)-4)-3*2^(mu+2)+3^mu*(12-9*2^(mu+2))+2*4^(2*mu)+2)*gamma(-mu))

r2 <- c(g1, g2, g3)

g1 <- 1
g2 <- -(((3*log(2)*l3*mu^2*2^(mu+1)+((log(4)-log(3))*l3*mu^2*3^mu+(3*log(2)-3*log(4))*l3*mu^2*2^mu+3*log(4)*l3*mu^2)*4^mu+((6*log(3)-6*log(2))*l3*mu^2*2^mu-8*log(3)*l3*mu^2)*3^mu)*gamma(-mu)+(l3*mu*3^(mu+1)-l3*mu*4^mu-3*l3*mu*2^mu+l3*mu)*digamma(-mu)+(log(4)*l3*mu+l3)*4^mu+(-3*log(3)*l3*mu-3*l3)*3^mu+(3*log(2)*l3*mu+3*l3)*2^mu-l3)*theta+(-mu^2*3^(2*mu+1)+3^mu*(3*mu^2*2^(mu+2)-10*mu^2)+3*mu^2*2^(mu+2)-9*mu^2*2^(2*mu)+(mu^2*3^mu-3*mu^2*2^mu+3*mu^2)*4^mu-3*mu^2)*gamma(-mu)-mu*3^(mu+1)+mu*4^mu+3*mu*2^mu-mu)/((2*mu^2*3^(2*mu+2)+9*mu^2*2^(2*mu+1)+4^mu*(3*mu^2*2^(mu+2)-4*mu^2*3^(mu+1)-4*mu^2)+3^mu*(12*mu^2-9*mu^2*2^(mu+2))-3*mu^2*2^(mu+2)+2*mu^2*4^(2*mu)+2*mu^2)*gamma(-mu))
g3 <-  (((3*log(2)*l2*mu^2*2^(mu+1)+((log(4)-log(3))*l2*mu^2*3^mu+(3*log(2)-3*log(4))*l2*mu^2*2^mu+3*log(4)*l2*mu^2)*4^mu+((6*log(3)-6*log(2))*l2*mu^2*2^mu-8*log(3)*l2*mu^2)*3^mu)*gamma(-mu)+(l2*mu*3^(mu+1)-l2*mu*4^mu-3*l2*mu*2^mu+l2*mu)*digamma(-mu)+(log(4)*l2*mu+l2)*4^mu+(-3*log(3)*l2*mu-3*l2)*3^mu+(3*log(2)*l2*mu+3*l2)*2^mu-l2)*theta)/((2*mu^2*3^(2*mu+2)+9*mu^2*2^(2*mu+1)+4^mu*(3*mu^2*2^(mu+2)-4*mu^2*3^(mu+1)-4*mu^2)+3^mu*(12*mu^2-9*mu^2*2^(mu+2))-3*mu^2*2^(mu+2)+2*mu^2*4^(2*mu)+2*mu^2)*gamma(-mu))

r1 <- c(g1, g2, g3)

rbind(r1, r2, r3)
}

# Matrices to convert parameter estimations to quantile estimations
CovParamtoQuan <- function(distr, param, p) {
if (distr == "gev") {
g <- param["shape"]
sigma <- param["scale"]
xi <- -log(p)
return(
cbind(
1,
1/g * (xi^(-g) - 1), #-((-log(p))^g - 1)/(g * (-log(p))^g),
sigma/g * (g^(-1) - xi^(-g) * (log(xi) + g^(-1))) #-((g*log(-log(p))-(-log(p))^g+1)*sigma)/(g^2*(-log(p))^g)
)
)
} else if (distr == "gumbel") {
return(
cbind(
1,
-log(-log(p))
)
)
} else {
stop("distr currently unsupported. ")
}
}
```

## Try the TLMoments package in your browser

Any scripts or data that you put into this service are public.

TLMoments documentation built on March 27, 2022, 5:07 p.m.