context("Generalized gamma distribution")
tol <- 1e-06
x <- c(-1,1,2,3,4); mu <- 2.2; sigma <- 1.6; Q <- 0; P <- 1
delta <- (Q^2 + 2*P)^{1/2}
s1 <- 2 / (Q^2 + 2*P + Q*delta); s2 <- 2 / (Q^2 + 2*P - Q*delta)
x <- c(-1,1,2,3,4); shape <- 2.2; scale <- 1.6; k <- 1.9
test_that("Generalized gamma reductions: d",{
expect_equal(dgengamma(c(-1,1,2,3,4), mu=0, sigma=1, Q=1), # FIXME value for x=0 and add here
c(0, 0.367879441171442, 0.135335283236613, 0.0497870683678639, 0.0183156388887342))
expect_equal(dgengamma(c(-1,0,1,2,3,4), mu=0, sigma=1, Q=0),
dlnorm(c(-1,0,1,2,3,4), meanlog=0, sdlog=1))
expect_equal(dgengamma(c(-1,0,1,2,3,4), mu=0, sigma=1:3, Q=0),
dlnorm(c(-1,0,1,2,3,4), meanlog=0, sdlog=1:3))
expect_equal(dgengamma(c(-1,0,1,2,3,4), mu=0, sigma=1:3, Q=0, log=TRUE),
dlnorm(c(-1,0,1,2,3,4), meanlog=0, sdlog=1:3, log=TRUE))
expect_equal(dgengamma(c(1,2,3,4), mu=0.1, sigma=1.2, Q=1),
dweibull(c(1,2,3,4), shape=1/1.2, scale=exp(0.1))) # only defined for x>0 anyway
x <- c(1,2,3,4); mu <- 0.4; sigma <- 1.2
expect_equal(dgengamma(x, mu=mu, sigma=sigma, Q=sigma),
dgamma(x, shape=1/sigma^2, scale=exp(mu)*sigma^2))
# FIXME add limiting value for x=0
expect_equal(dgengamma.orig(x, shape=shape, scale=scale, k=k),
dgengamma(x, mu=log(scale) + log(k)/shape, sigma=1/(shape*sqrt(k)), Q=1/sqrt(k)))
})
test_that("Generalized gamma reductions: p",{
expect_equal(pgengamma(c(-1,0,1,2,3,4), mu=0, sigma=1, Q=1),
c(0, 0, 0.632120558828558, 0.864664716763387, 0.950212931632136, 0.981684361111266))
expect_equal(pgengamma(c(-1,0,1,2,3,4), mu=0, sigma=1, Q=0),
plnorm(c(-1,0,1,2,3,4), meanlog=0, sdlog=1))
expect_equal(pgengamma(c(-1,0,1,2,3,4), mu=0.1, sigma=1.2, Q=1),
pweibull(c(-1,0,1,2,3,4), shape=1/1.2, scale=exp(0.1)))
x <- c(1,2,3,4); mu <- 0.4; sigma <- 1.2
expect_equal(pgengamma(x, mu=mu, sigma=sigma, Q=sigma),
pgamma(x, shape=1/sigma^2, scale=exp(mu)*sigma^2))
})
test_that("Generalized gamma reductions: q",{
expect_equal(qgengamma(p=0.25, mu=0, sigma=1, Q=1), 0.287682072451781)
expect_equal(qgengamma(p=0.25, mu=0, sigma=1, Q=1), qgeneric(pgengamma, p=0.25, mu=0, sigma=1, Q=1))
expect_equal(qgengamma(p=0, mu=0, sigma=1, Q=1), 0)
expect_equal(qgengamma(p=0.25, mu=0, sigma=1, Q=0), qlnorm(p=0.25, meanlog=0, sdlog=1))
expect_equal(qgengamma(p=0.25, mu=0, sigma=1.2, Q=0), qlnorm(p=0.25, meanlog=0, sdlog=1.2))
expect_equal(qgengamma(p=0.25, mu=0.1, sigma=1.2, Q=1), qweibull(p=0.25, scale=exp(0.1), shape=1/1.2))
expect_equal(qgengamma(pgengamma(q=c(-2,-1,0,1,2,3,4), mu=0, sigma=1, Q=1), mu=0, sigma=1, Q=1),
c(0,0,0,1,2,3,4))
expect_equal(qgengamma(pgengamma(q=c(-2,-1,0,1,2,3,4), mu=0, sigma=1, Q=-1), mu=0, sigma=1, Q=-1),
c(0,0,0,1,2,3,4))
expect_equal(qlnorm(pgengamma(q=c(-2,-1,0,1,2,3,4), mu=0, sigma=1.2, Q=0), meanlog=0, sdlog=1.2),
c(0,0,0,1,2,3,4))
expect_equal(qlnorm(pgengamma(q=c(-2,-1,0,1,2,3,4), mu=0, sigma=1, Q=0, log.p=TRUE), meanlog=0, sdlog=1, log.p=TRUE),
c(0,0,0,1,2,3,4))
expect_equal(qgengamma(p=0.25, mu=0, sigma=1, Q=1, lower.tail=TRUE), qgeneric(pgengamma, p=0.25, mu=0, sigma=1, Q=1, lower.tail=TRUE))
})
test_that("Generalized gamma reductions: r",{
rgengamma(n=10, mu=0, sigma=1, Q=0)
set.seed(22061976)
x <- rgengamma(n=10, mu=0, sigma=1.1, Q=0)
set.seed(22061976)
y <- rlnorm(n=10, meanlog=0, sdlog=1.1)
expect_equal(x, y)
})
test_that("Gives errors",
{
expect_warning(dgengamma(1, 1, -2, 1), "Negative")
expect_warning(pgengamma(1, 1, -2, 1), "Negative")
expect_warning(qgengamma(0.1, 1, -2, 1), "Negative")
expect_warning(rgengamma(1, 1, -2, 1), "Negative")
})
test_that("Generalized gamma definition in Stata manual",{
Sgg <- function(t, mu, sigma, kappa){
IGF <- function(a, x){ pgamma(x, a) } # incomplete gamma function
gamma <- 1 / kappa^2
z <- (log(t) - mu)/sigma
z[kappa<0] <- -z[kappa<0] # Stata manual uses sign(0) = 1
u <- gamma*exp(abs(kappa)*z)
ifelse(kappa > 0,
1 - IGF(gamma, u),
ifelse(kappa==0,
1 - pnorm(z),
IGF(gamma, u)))
}
expect_equal(
Sgg(c(1,2,3),c(-1,2,0.2),c(1,2,1),c(-1, 0, 1)),
pgengamma(c(1,2,3),c(-1,2,0.2),c(1,2,1),c(-1, 0, 1), lower.tail=FALSE)
)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.