Nothing
R/family.qreg.R
In VGAM: Vector Generalized Linear and Additive Models
Defines functions
laplace
rlaplace
qlaplace
plaplace
dlaplace
fff
fff.control
rpolono
ppolono
triangle
triangle.control
ptriangle
qtriangle
rtriangle
dtriangle
hyperg
rbenini
qbenini
pbenini
dbenini
benini1
amlexponential
amlexponential.deviance
amlbinomial
amlbinomial.deviance
amlpoisson
amlpoisson.deviance
amlnormal
amlnormal.deviance
Wr2
Wr1
lmscreg.control
lms.yjn
lms.yjn2
lms.yjn2.control
gleg.weight.yjn.13
gleg.weight.yjn.12
gleg.weight.yjn.11
glag.weight.yjn.13
glag.weight.yjn.12
glag.weight.yjn.11
gh.weight.yjn.13
gh.weight.yjn.12
gh.weight.yjn.11
dpsi.dlambda.yjn
yeo.johnson
dyj.dy.yeojohnson
dy.dpsi.yeojohnson
lms.bcg
lms.bcn
lms.yjn.control
qlms.bcn
dlms.bcn
Documented in
amlbinomial
amlexponential
amlnormal
amlpoisson
benini1
dbenini
dlaplace
dlms.bcn
dpsi.dlambda.yjn
dtriangle
fff
fff.control
gh.weight.yjn.11
gh.weight.yjn.12
gh.weight.yjn.13
glag.weight.yjn.11
glag.weight.yjn.12
glag.weight.yjn.13
gleg.weight.yjn.11
gleg.weight.yjn.12
gleg.weight.yjn.13
hyperg
hyperg
laplace
lms.bcg
lms.bcn
lmscreg.control
lms.yjn
lms.yjn2
lms.yjn.control
pbenini
plaplace
ppolono
ptriangle
qbenini
qlaplace
qlms.bcn
qtriangle
rbenini
rlaplace
rpolono
rtriangle
Wr1
Wr2
yeo.johnson
# These functions are
# Copyright (C) 1998-2024 T.W. Yee, University of Auckland.
# All rights reserved.
dlms.bcn <- function(x, lambda = 1, mu = 0, sigma = 1,
tol0 = 0.001, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
zedd <- ((x/mu)^lambda - 1) / (lambda * sigma)
log.dz.dy <- (lambda - 1) * log(x/mu) - log(mu * sigma)
is.eff.0 <- abs(lambda) < tol0
if (any(is.eff.0)) {
zedd[is.eff.0] <- log(x[is.eff.0] /
mu[is.eff.0]) / sigma[is.eff.0]
log.dz.dy[is.eff.0] <- -log(x[is.eff.0] * sigma[is.eff.0])
}
logden <- dnorm(zedd, log = TRUE) + log.dz.dy
if (log.arg) logden else exp(logden)
}
qlms.bcn <- function(p, lambda = 1, mu = 0, sigma = 1) {
answer <- mu * (1 + lambda * sigma * qnorm(p))^(1/lambda)
answer
}
lms.bcn.control <-
lms.bcg.control <-
lms.yjn.control <- function(trace = TRUE, ...)
list(trace = trace)
lms.bcn <- function(percentiles = c(25, 50, 75),
zero = c("lambda", "sigma"),
llambda = "identitylink",
lmu = "identitylink",
lsigma = "loglink",
idf.mu = 4,
idf.sigma = 2,
ilambda = 1,
isigma = NULL,
tol0 = 0.001) {
llambda <- as.list(substitute(llambda))
elambda <- link2list(llambda)
llambda <- attr(elambda, "function.name")
lmu <- as.list(substitute(lmu))
emu <- link2list(lmu)
lmu <- attr(emu, "function.name")
lsigma <- as.list(substitute(lsigma))
esigma <- link2list(lsigma)
lsigma <- attr(esigma, "function.name")
if (!is.Numeric(tol0, positive = TRUE, length.arg = 1))
stop("bad input for argument 'tol0'")
if (!is.Numeric(ilambda))
stop("bad input for argument 'ilambda'")
if (length(isigma) &&
!is.Numeric(isigma, positive = TRUE))
stop("bad input for argument 'isigma'")
new("vglmff",
blurb = c("LMS ",
"quantile",
" regression (Box-Cox transformation to normality)\n",
"Links: ",
namesof("lambda", link = llambda, earg = elambda), ", ",
namesof("mu", link = lmu, earg = emu), ", ",
namesof("sigma", link = lsigma, earg = esigma)),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 3)
}), list( .zero = zero))),
infos = eval(substitute(function(...) {
list(M1 = 3,
Q1 = 1,
expected = TRUE,
multipleResponses = FALSE,
parameters.names = c("lambda", "mu", "sigma"),
llambda = .llambda ,
lmu = .lmu ,
lsigma = .lsigma ,
percentiles = .percentiles , # For the original fit only
true.mu = FALSE, # quantiles
zero = .zero )
}, list( .zero = zero,
.percentiles = percentiles,
.llambda = llambda, .lmu = lmu, .lsigma = lsigma ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
Is.positive.y = TRUE,
ncol.w.max = 1, ncol.y.max = 1)
predictors.names <-
c(namesof("lambda", .llambda, .elambda, short= TRUE),
namesof("mu", .lmu, .emu, short= TRUE),
namesof("sigma", .lsigma, .esigma, short= TRUE))
extra$percentiles <- .percentiles
if (!length(etastart)) {
Fit5 <- vsmooth.spline(x = x[, min(ncol(x), 2)],
y = y, w = w, df = .idf.mu )
fv.init <- c(predict(Fit5, x = x[, min(ncol(x), 2)])$y)
lambda.init <- if (is.Numeric( .ilambda )) .ilambda else 1.0
sigma.init <- if (is.null(.isigma)) {
myratio <- ((y/fv.init)^lambda.init - 1) / lambda.init
if (is.Numeric( .idf.sigma )) {
fit600 <- vsmooth.spline(x = x[, min(ncol(x), 2)],
y = myratio^2,
w = w, df = .idf.sigma)
sqrt(c(abs(predict(fit600, x = x[, min(ncol(x), 2)])$y)))
} else {
sqrt(var(myratio))
}
} else {
.isigma
}
etastart <-
cbind(theta2eta(lambda.init, .llambda , .elambda ),
theta2eta(fv.init, .lmu , .emu ),
theta2eta(sigma.init, .lsigma , .esigma ))
}
}), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.idf.mu = idf.mu,
.idf.sigma = idf.sigma,
.ilambda = ilambda, .isigma = isigma,
.percentiles = percentiles ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
pcent <- extra$percentiles
eta[, 1] <- eta2theta(eta[, 1], .llambda , .elambda )
eta[, 2] <- eta2theta(eta[, 2], .lmu , .emu )
eta[, 3] <- eta2theta(eta[, 3], .lsigma , .esigma )
qtplot.lms.bcn(percentiles = pcent, eta = eta)
}, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma
))),
last = eval(substitute(expression({
misc$links <- c(lambda = .llambda ,
mu = .lmu , sigma = .lsigma )
misc$earg <- list(lambda = .elambda ,
mu = .emu , sigma = .esigma )
misc$tol0 <- .tol0
misc$percentiles <- .percentiles # These are argument values
if (control$cdf) {
post$cdf <-
cdf.lms.bcn(y,
eta0 = matrix(c(lambda, mymu, sigma), ncol = 3,
dimnames = list(dimnames(x)[[1]], NULL)))
}
}), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.percentiles = percentiles,
.tol0 = tol0 ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
lambda <- eta2theta(eta[, 1], .llambda , .elambda )
muvec <- eta2theta(eta[, 2], .lmu , .emu )
sigma <- eta2theta(eta[, 3], .lsigma , .esigma )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- dlms.bcn(y, lam = lambda, mu = mu, sig = sigma,
tol0 = .tol0 , log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.tol0 = tol0 ))),
vfamily = c("lms.bcn", "lmscreg"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta2theta(eta[, 2], .lmu , earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
okay1 <- all(is.finite(mymu )) &&
all(is.finite(sigma )) && all(0 < sigma) &&
all(is.finite(lambda))
okay1
}, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.tol0 = tol0 ))),
deriv = eval(substitute(expression({
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta2theta(eta[, 2], .lmu , earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
zedd <- ((y / mymu)^lambda - 1) / (lambda * sigma)
z2m1 <- zedd * zedd - 1
dl.dlambda <- zedd * (zedd - log(y/mymu) / sigma) / lambda -
z2m1 * log(y/mymu)
dl.dmu <- zedd / (mymu * sigma) + z2m1 * lambda / mymu
dl.dsigma <- z2m1 / sigma
dlambda.deta <- dtheta.deta(lambda, .llambda , .elambda )
dmu.deta <- dtheta.deta(mymu, .lmu , .emu )
dsigma.deta <- dtheta.deta(sigma, .lsigma , .esigma )
c(w) * cbind(dl.dlambda * dlambda.deta,
dl.dmu * dmu.deta,
dl.dsigma * dsigma.deta)
}), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))),
weight = eval(substitute(expression({
wz <- matrix(NA_real_, n, 6)
wz[,iam(1, 1, M)] <- (7 * sigma^2 / 4) * dlambda.deta^2
wz[,iam(2, 2, M)] <- (1 + 2*(lambda*sigma)^2)/(mymu*sigma)^2 *
dmu.deta^2
wz[,iam(3, 3, M)] <- (2 / sigma^2) * dsigma.deta^2
wz[,iam(1, 2, M)] <- (-1 / (2 * mymu)) * dlambda.deta * dmu.deta
wz[,iam(1, 3, M)] <- (lambda * sigma) * dlambda.deta * dsigma.deta
wz[,iam(2, 3, M)] <- (2*lambda/(mymu * sigma)) *
dmu.deta * dsigma.deta
c(w) * wz
}),
list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))))
} # lms.bcn
lms.bcg <- function(percentiles = c(25, 50, 75),
zero = c("lambda", "sigma"),
llambda = "identitylink",
lmu = "identitylink",
lsigma = "loglink",
idf.mu = 4,
idf.sigma = 2,
ilambda = 1,
isigma = NULL) {
llambda <- as.list(substitute(llambda))
elambda <- link2list(llambda)
llambda <- attr(elambda, "function.name")
lmu <- as.list(substitute(lmu))
emu <- link2list(lmu)
lmu <- attr(emu, "function.name")
lsigma <- as.list(substitute(lsigma))
esigma <- link2list(lsigma)
lsigma <- attr(esigma, "function.name")
if (!is.Numeric(ilambda))
stop("bad input for argument 'ilambda'")
if (length(isigma) && !is.Numeric(isigma, positive = TRUE))
stop("bad input for argument 'isigma'")
new("vglmff",
blurb = c("LMS Quantile Regression ",
"(Box-Cox transformation to a Gamma distribution)\n",
"Links: ",
namesof("lambda", link = llambda, elambda), ", ",
namesof("mu", link = lmu, emu), ", ",
namesof("sigma", link = lsigma, esigma)),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 3)
}), list(.zero = zero))),
infos = eval(substitute(function(...) {
list(M1 = 3,
Q1 = 1,
expected = TRUE,
multipleResponses = FALSE,
parameters.names = c("lambda", "mu", "sigma"),
llambda = .llambda ,
lmu = .lmu ,
lsigma = .lsigma ,
percentiles = .percentiles , # For the original fit only
true.mu = FALSE, # quantiles
zero = .zero )
}, list( .zero = zero,
.percentiles = percentiles,
.llambda = llambda, .lmu = lmu, .lsigma = lsigma ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
Is.positive.y = TRUE,
ncol.w.max = 1, ncol.y.max = 1)
predictors.names <- c(
namesof("lambda", .llambda, .elambda, short = TRUE),
namesof("mu", .lmu, .emu, short = TRUE),
namesof("sigma", .lsigma, .esigma, short = TRUE))
extra$percentiles <- .percentiles
if (!length(etastart)) {
Fit5 <- vsmooth.spline(x = x[, min(ncol(x), 2)],
y = y, w = w, df = .idf.mu )
fv.init <- c(predict(Fit5, x = x[, min(ncol(x), 2)])$y)
lambda.init <- if (is.Numeric( .ilambda )) .ilambda else 1.0
sigma.init <- if (is.null( .isigma )) {
myratio <- ((y/fv.init)^lambda.init-1) / lambda.init
if (is.numeric( .idf.sigma ) &&
is.finite( .idf.sigma )) {
fit600 <- vsmooth.spline(x = x[, min(ncol(x), 2)],
y = (myratio)^2,
w = w, df = .idf.sigma )
sqrt(c(abs(predict(fit600, x = x[, min(ncol(x), 2)])$y)))
} else {
sqrt(var(myratio))
}
} else .isigma
etastart <-
cbind(theta2eta(lambda.init, .llambda , .elambda ),
theta2eta(fv.init, .lmu , .emu ),
theta2eta(sigma.init, .lsigma , .esigma ))
}
}), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.idf.mu = idf.mu,
.idf.sigma = idf.sigma,
.ilambda = ilambda, .isigma = isigma,
.percentiles = percentiles ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
pcent <- extra$percentiles
eta[, 1] <- eta2theta(eta[, 1], .llambda , .elambda )
eta[, 2] <- eta2theta(eta[, 2], .lmu , .emu )
eta[, 3] <- eta2theta(eta[, 3], .lsigma , .esigma )
qtplot.lms.bcg(percentiles = pcent, eta = eta)
}, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))),
last = eval(substitute(expression({
misc$link <- c(lambda = .llambda , mu = .lmu , sigma = .lsigma )
misc$earg <- list(lambda = .elambda , mu = .emu , sigma = .esigma )
misc$percentiles <- .percentiles # These are argument values
if (control$cdf) {
post$cdf <- cdf.lms.bcg(y, eta0 = matrix(c(lambda, mymu, sigma),
ncol = 3, dimnames = list(dimnames(x)[[1]], NULL)))
}
}), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.percentiles = percentiles ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mu <- eta2theta(eta[, 2], .lmu , earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
Gee <- (y / mu)^lambda
theta <- 1 / (sigma * lambda)^2
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * (log(abs(lambda)) + theta * (log(theta) +
log(Gee)-Gee) - lgamma(theta) - log(y))
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))),
vfamily = c("lms.bcg", "lmscreg"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta2theta(eta[, 2], .lmu , earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
okay1 <- all(is.finite(mymu )) &&
all(is.finite(sigma )) && all(0 < sigma) &&
all(is.finite(lambda))
okay1
}, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))),
deriv = eval(substitute(expression({
lambda <- eta2theta(eta[, 1], .llambda, earg = .elambda )
mymu <- eta2theta(eta[, 2], .lmu, earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma, earg = .esigma )
Gee <- (y / mymu)^lambda
theta <- 1 / (sigma * lambda)^2
dd <- digamma(theta)
dl.dlambda <- (1 + 2 * theta * (dd + Gee -1 -log(theta) -
0.5 * (Gee + 1) * log(Gee))) / lambda
dl.dmu <- lambda * theta * (Gee-1) / mymu
dl.dsigma <- 2*theta*(dd + Gee - log(theta * Gee)-1) / sigma
dlambda.deta <- dtheta.deta(lambda, link = .llambda , .elambda )
dmu.deta <- dtheta.deta(mymu, link = .lmu , .emu )
dsigma.deta <- dtheta.deta(sigma, link = .lsigma , .esigma )
cbind(dl.dlambda * dlambda.deta,
dl.dmu * dmu.deta,
dl.dsigma * dsigma.deta) * w
}), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))),
weight = eval(substitute(expression({
tritheta <- trigamma(theta)
wz <- matrix(0, n, 6)
if (TRUE) {
part2 <- dd + 2/theta - 2*log(theta)
wz[,iam(1, 1, M)] <- ((1 + theta*(tritheta*(1+4*theta) -
4*(1+1/theta) - log(theta)*(2/theta -
log(theta)) + dd*part2)) / lambda^2) *
dlambda.deta^2
} else {
temp <- mean( Gee*(log(Gee))^2 )
wz[,iam(1, 1, M)] <- ((4 * theta * (theta *
tritheta-1) - 1 +
theta*temp) / lambda^2) *
dlambda.deta^2
}
wz[,iam(2, 2, M)] <- dmu.deta^2 / (mymu * sigma)^2
wz[,iam(3, 3, M)] <- (4 * theta * (theta * tritheta -
1) / sigma^2) *
dsigma.deta^2
wz[,iam(1, 2, M)] <- (-theta * (dd + 1 / theta -
log(theta)) / mymu) *
dlambda.deta * dmu.deta
wz[,iam(1, 3, M)] <- 2 * theta^1.5 * (2 * theta * tritheta - 2 -
1 / theta) * dlambda.deta * dsigma.deta
c(w) * wz
}),
list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))))
} # lms.bcg
dy.dpsi.yeojohnson <- function(psi, lambda) {
L <- max(length(psi), length(lambda))
if (length(psi) != L) psi <- rep_len(psi, L)
if (length(lambda) != L) lambda <- rep_len(lambda, L)
ifelse(psi > 0, (1 + psi * lambda)^(1/lambda - 1),
(1 - (2-lambda) * psi)^((lambda - 1) / (2-
lambda)))
}
dyj.dy.yeojohnson <- function(y, lambda) {
L <- max(length(y), length(lambda))
if (length(y) != L) y <- rep_len(y, L)
if (length(lambda) != L) lambda <- rep_len(lambda, L)
ifelse(y>0, (1 + y)^(lambda - 1), (1 - y)^(1 - lambda))
}
yeo.johnson <-
function(y, lambda, derivative = 0,
epsilon = sqrt(.Machine$double.eps),
inverse = FALSE) {
if (!is.Numeric(derivative, length.arg = 1,
integer.valued = TRUE) ||
derivative < 0)
stop("argument 'derivative' must be a non-negative integer")
ans <- y
if (!is.Numeric(epsilon, length.arg = 1, positive = TRUE))
stop("argument 'epsilon' must be a single positive number")
L <- max(length(lambda), length(y))
if (length(y) != L) y <- rep_len(y, L)
if (length(lambda) != L) lambda <- rep_len(lambda, L)
if (inverse) {
if (derivative != 0)
stop("argument 'derivative' must 0 when inverse = TRUE")
if (any(index <- y >= 0 & abs(lambda ) > epsilon))
ans[index] <- (y[index]*lambda[index] +
1)^(1/lambda[index]) - 1
if (any(index <- y >= 0 & abs(lambda ) <= epsilon))
ans[index] <- expm1(y[index])
if (any(index <- y < 0 & abs(lambda-2) > epsilon))
ans[index] <- 1 - (-(2-lambda[index]) *
y[index]+1)^(1/(2-lambda[index]))
if (any(index <- y < 0 & abs(lambda-2) <= epsilon))
ans[index] <- -expm1(-y[index])
return(ans)
}
if (derivative == 0) {
if (any(index <- y >= 0 & abs(lambda ) > epsilon))
ans[index] <- ((y[index]+1)^(lambda[index]) -
1) / lambda[index]
if (any(index <- y >= 0 & abs(lambda ) <= epsilon))
ans[index] <- log1p(y[index])
if (any(index <- y < 0 & abs(lambda-2) > epsilon))
ans[index] <- -((-y[index]+1)^(2-lambda[index]) - 1)/(2 -
lambda[index])
if (any(index <- y < 0 & abs(lambda-2) <= epsilon))
ans[index] <- -log1p(-y[index])
} else {
psi <- Recall(y = y, lambda = lambda,
derivative = derivative - 1,
epsilon = epsilon, inverse = inverse)
if (any(index <- y >= 0 & abs(lambda ) > epsilon))
ans[index] <- ( (y[index]+1)^(lambda[index]) *
(log1p(y[index]))^(derivative) - derivative *
psi[index] ) / lambda[index]
if (any(index <- y >= 0 & abs(lambda ) <= epsilon))
ans[index] <- (log1p(y[index]))^(derivative +
1) / (derivative + 1)
if (any(index <- y < 0 & abs(lambda-2) > epsilon))
ans[index] <- -( (-y[index]+1)^(2-lambda[index]) *
(-log1p(-y[index]))^(derivative) - derivative *
psi[index] ) / (2-lambda[index])
if (any(index <- y < 0 & abs(lambda-2) <= epsilon))
ans[index] <- (-log1p(-y[index]))^(derivative +
1) / (derivative + 1)
}
ans
} # yeo.johnson
dpsi.dlambda.yjn <- function(psi, lambda, mymu, sigma,
derivative = 0, smallno = 1.0e-8) {
if (!is.Numeric(derivative, length.arg = 1,
integer.valued = TRUE) ||
derivative < 0)
stop("argument 'derivative' must be a non-negative integer")
if (!is.Numeric(smallno, length.arg = 1, positive = TRUE))
stop("argument 'smallno' must be a single positive number")
L <- max(length(psi), length(lambda), length(mymu),
length(sigma))
if (length(psi) != L) psi <- rep_len(psi, L)
if (length(lambda) != L) lambda <- rep_len(lambda, L)
if (length(mymu) != L) mymu <- rep_len(mymu, L)
if (length(sigma) != L) sigma <- rep_len(sigma, L)
answer <- matrix(NA_real_, L, derivative+1)
CC <- psi >= 0
BB <- ifelse(CC, lambda, -2+lambda)
AA <- psi * BB
temp8 <- if (derivative > 0) {
answer[,1:derivative] <-
Recall(psi = psi, lambda = lambda, mymu = mymu,
sigma = sigma,
derivative = derivative-1, smallno = smallno)
answer[,derivative] * derivative
} else {
0
}
answer[, 1+derivative] <- ((AA+1) *
(log1p(AA)/BB)^derivative - temp8) / BB
pos <- (CC & abs(lambda) <= smallno) |
(!CC & abs(lambda-2) <= smallno)
if (any(pos))
answer[pos,1+derivative] =
(answer[pos, 1]^(1+derivative))/(derivative+1)
answer
}
gh.weight.yjn.11 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
((derivmat[, 2]/sigma)^2 +
sqrt(2) * z * derivmat[, 3] / sigma) / sqrt(pi)
} else {
# Long-winded way
psi <- mymu + sqrt(2) * sigma * z
(1 / sqrt(pi)) *
(dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 1)[, 2]^2 +
(psi - mymu) *
dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 2)[, 3]) / sigma^2
}
}
gh.weight.yjn.12 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
(-derivmat[, 2]) / (sqrt(pi) * sigma^2)
} else {
psi <- mymu + sqrt(2) * sigma * z
(1 / sqrt(pi)) * (-dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 1)[, 2]) / sigma^2
}
}
gh.weight.yjn.13 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
sqrt(8 / pi) * (-derivmat[, 2]) * z / sigma^2
} else {
psi <- mymu + sqrt(2) * sigma * z
(1 / sqrt(pi)) *
(-2 * dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 1)[, 2]) *
(psi - mymu) / sigma^3
}
}
glag.weight.yjn.11 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
derivmat[, 4] * (derivmat[, 2]^2 + sqrt(2) *
sigma * z * derivmat[, 3])
} else {
psi <- mymu + sqrt(2) * sigma * z
discontinuity <- -mymu / (sqrt(2) * sigma)
(1 / (2 * sqrt((z-discontinuity^2)^2 + discontinuity^2))) *
(1 / sqrt(pi)) *
(dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 1)[, 2]^2 +
(psi - mymu) *
dpsi.dlambda.yjn(psi, lambda, mymu,
sigma, derivative = 2)[, 3]) / sigma^2
}
}
glag.weight.yjn.12 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
discontinuity <- -mymu / (sqrt(2) * sigma)
if (length(derivmat)) {
derivmat[, 4] * (-derivmat[, 2])
} else {
psi <- mymu + sqrt(2) * sigma * z
(1 / (2 * sqrt((z-discontinuity^2)^2 + discontinuity^2))) *
(1 / sqrt(pi)) *
(- dpsi.dlambda.yjn(psi, lambda, mymu,
sigma, derivative = 1)[, 2]) / sigma^2
}
}
glag.weight.yjn.13 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
derivmat[, 4] * (-derivmat[, 2]) * sqrt(8) * z
} else {
psi <- mymu + sqrt(2) * sigma * z
discontinuity <- -mymu / (sqrt(2) * sigma)
(1 / (2 * sqrt((z-discontinuity^2)^2 + discontinuity^2))) *
(1 / sqrt(pi)) *
(-2 * dpsi.dlambda.yjn(psi, lambda, mymu,
sigma, derivative = 1)[, 2]) *
(psi - mymu) / sigma^3
}
}
gleg.weight.yjn.11 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
derivmat[, 4] * (derivmat[, 2]^2 + sqrt(2) *
sigma*z* derivmat[, 3])
} else {
psi <- mymu + sqrt(2) * sigma * z
(exp(-z^2) / sqrt(pi)) *
(dpsi.dlambda.yjn(psi, lambda, mymu,
sigma, derivative = 1)[, 2]^2 +
(psi - mymu) *
dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 2)[, 3]) / sigma^2
}
}
gleg.weight.yjn.12 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
derivmat[, 4] * (- derivmat[, 2])
} else {
psi <- mymu + sqrt(2) * sigma * z
(exp(-z^2) / sqrt(pi)) *
(- dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 1)[, 2]) / sigma^2
}
}
gleg.weight.yjn.13 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
derivmat[, 4] * (-derivmat[, 2]) * sqrt(8) * z
} else {
psi <- mymu + sqrt(2) * sigma * z
(exp(-z^2) / sqrt(pi)) *
(-2 * dpsi.dlambda.yjn(psi, lambda, mymu,
sigma, derivative = 1)[, 2]) *
(psi - mymu) / sigma^3
}
}
lms.yjn2.control <- function(save.weights = TRUE, ...) {
list(save.weights = save.weights)
}
lms.yjn2 <- function(percentiles = c(25, 50, 75),
zero = c("lambda", "sigma"),
llambda = "identitylink",
lmu = "identitylink",
lsigma = "loglink",
idf.mu = 4,
idf.sigma = 2,
ilambda = 1.0,
isigma = NULL,
yoffset = NULL,
nsimEIM = 250) {
llambda <- as.list(substitute(llambda))
elambda <- link2list(llambda)
llambda <- attr(elambda, "function.name")
lmu <- as.list(substitute(lmu))
emu <- link2list(lmu)
lmu <- attr(emu, "function.name")
lsigma <- as.list(substitute(lsigma))
esigma <- link2list(lsigma)
lsigma <- attr(esigma, "function.name")
if (!is.Numeric(ilambda))
stop("bad input for argument 'ilambda'")
if (length(isigma) &&
!is.Numeric(isigma, positive = TRUE))
stop("bad input for argument 'isigma'")
new("vglmff",
blurb = c("LMS Quantile Regression (Yeo-Johnson transformation",
" to normality)\n",
"Links: ",
namesof("lambda", link = llambda, elambda), ", ",
namesof("mu", link = lmu, emu ), ", ",
namesof("sigma", link = lsigma, esigma )),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 3)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 3,
Q1 = 1,
expected = TRUE,
multipleResponses = FALSE,
parameters.names = c("lambda", "mu", "sigma"),
llambda = .llambda ,
lmu = .lmu ,
lsigma = .lsigma ,
percentiles = .percentiles , # For the original fit only
true.mu = FALSE, # quantiles
zero = .zero )
}, list( .zero = zero,
.percentiles = percentiles,
.llambda = llambda, .lmu = lmu, .lsigma = lsigma ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
ncol.w.max = 1, ncol.y.max = 1)
extra$percentiles <- .percentiles
predictors.names <-
c(namesof("lambda", .llambda, earg = .elambda, short= TRUE),
namesof("mu", .lmu, earg = .emu, short= TRUE),
namesof("sigma", .lsigma, earg = .esigma, short= TRUE))
y.save <- y
yoff <- if (is.Numeric( .yoffset)) .yoffset else -median(y)
extra$yoffset <- yoff
y <- y + yoff
if (!length(etastart)) {
lambda.init <- if (is.Numeric( .ilambda )) .ilambda else 1.
y.tx <- yeo.johnson(y, lambda.init)
fv.init =
if (smoothok <-
(length(unique(sort(x[, min(ncol(x), 2)]))) > 7)) {
fit700 <- vsmooth.spline(x = x[, min(ncol(x), 2)],
y = y.tx, w = w, df = .idf.mu )
c(predict(fit700, x = x[, min(ncol(x), 2)])$y)
} else {
rep_len(weighted.mean(y, w), n)
}
sigma.init <- if (!is.Numeric(.isigma)) {
if (is.Numeric( .idf.sigma) && smoothok) {
fit710 = vsmooth.spline(x = x[, min(ncol(x), 2)],
y = (y.tx - fv.init)^2,
w = w, df = .idf.sigma)
sqrt(c(abs(predict(fit710,
x = x[, min(ncol(x), 2)])$y)))
} else {
sqrt( sum( w * (y.tx - fv.init)^2 ) / sum(w) )
}
} else
.isigma
etastart <- matrix(0, n, 3)
etastart[, 1] <- theta2eta(lambda.init, .llambda, .elambda)
etastart[, 2] <- theta2eta(fv.init, .lmu, .emu)
etastart[, 3] <- theta2eta(sigma.init, .lsigma, .esigma)
}
}), list(.llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.ilambda = ilambda, .isigma = isigma,
.idf.mu = idf.mu,
.idf.sigma = idf.sigma,
.yoffset = yoffset,
.percentiles = percentiles ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
pcent <- extra$percentiles
eta[, 1] <- eta2theta(eta[, 1], .llambda, earg = .elambda)
eta[, 3] <- eta2theta(eta[, 3], .lsigma, earg = .esigma)
qtplot.lms.yjn(percentiles = pcent, eta = eta,
yoffset = extra$yoff)
}, list( .esigma = esigma, .elambda = elambda,
.llambda = llambda,
.lsigma = lsigma ))),
last = eval(substitute(expression({
misc$link <- c(lambda = .llambda, mu = .lmu, sigma = .lsigma)
misc$earg <- list(lambda = .elambda, mu = .emu, sigma = .esigma)
misc$nsimEIM <- .nsimEIM
misc$percentiles <- .percentiles # These are argument values
misc[["yoffset"]] <- extra$yoffset
y <- y.save # Restore back the value; to be attached to object
if (control$cdf) {
post$cdf <- cdf.lms.yjn(y + misc$yoffset,
eta0=matrix(c(lambda,mymu,sigma),
ncol=3, dimnames = list(dimnames(x)[[1]], NULL)))
}
}), list(.percentiles = percentiles,
.elambda = elambda, .emu = emu, .esigma = esigma,
.nsimEIM=nsimEIM,
.llambda = llambda, .lmu = lmu, .lsigma = lsigma ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mu <- eta2theta(eta[, 2], .lmu , earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
psi <- yeo.johnson(y, lambda)
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * (-log(sigma) - 0.5 * ((psi-mu)/sigma)^2 +
(lambda-1) * sign(y) * log1p(abs(y)))
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .elambda = elambda, .emu = emu, .esigma = esigma,
.llambda = llambda, .lmu = lmu,
.lsigma = lsigma ))),
vfamily = c("lms.yjn2", "lmscreg"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta2theta(eta[, 2], .lmu , earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
okay1 <- all(is.finite(mymu )) &&
all(is.finite(sigma )) && all(0 < sigma) &&
all(is.finite(lambda))
okay1
}, list( .elambda = elambda, .emu = emu, .esigma = esigma,
.llambda = llambda, .lmu = lmu,
.lsigma = lsigma ))),
deriv = eval(substitute(expression({
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta2theta(eta[, 2], .lmu , .emu )
sigma <- eta2theta(eta[, 3], .lsigma , .esigma )
dlambda.deta <- dtheta.deta(lambda, link = .llambda, .elambda)
dmu.deta <- dtheta.deta(mymu, link = .lmu, .emu)
dsigma.deta <- dtheta.deta(sigma, link = .lsigma, .esigma)
psi <- yeo.johnson(y, lambda)
d1 <- yeo.johnson(y, lambda, deriv = 1)
AA <- (psi - mymu) / sigma
dl.dlambda <- -AA * d1 /sigma + sign(y) * log1p(abs(y))
dl.dmu <- AA / sigma
dl.dsigma <- (AA^2 -1) / sigma
dthetas.detas <- cbind(dlambda.deta, dmu.deta, dsigma.deta)
c(w) * cbind(dl.dlambda, dl.dmu, dl.dsigma) * dthetas.detas
}), list( .elambda = elambda, .emu = emu, .esigma = esigma,
.llambda = llambda, .lmu = lmu,
.lsigma = lsigma ))),
weight = eval(substitute(expression({
run.varcov <- 0
ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)
for (ii in 1:( .nsimEIM )) {
psi <- rnorm(n, mymu, sigma)
ysim <- yeo.johnson(y = psi, lam = lambda, inverse = TRUE)
d1 <- yeo.johnson(ysim, lambda, deriv = 1)
AA <- (psi - mymu) / sigma
dl.dlambda <- -AA * d1 /sigma + sign(ysim) * log1p(abs(ysim))
dl.dmu <- AA / sigma
dl.dsigma <- (AA^2 -1) / sigma
rm(ysim)
temp3 <- cbind(dl.dlambda, dl.dmu, dl.dsigma)
run.varcov <- ((ii-1) * run.varcov +
temp3[,ind1$row.index]*temp3[,ind1$col.index]) / ii
}
if (intercept.only)
run.varcov <- matrix(colMeans(run.varcov),
n, ncol(run.varcov), byrow = TRUE)
wz <- run.varcov * dthetas.detas[,ind1$row] *
dthetas.detas[,ind1$col]
dimnames(wz) <- list(rownames(wz), NULL) # Remove the colnames
c(w) * wz
}), list(.lsigma = lsigma,
.esigma = esigma, .elambda = elambda,
.nsimEIM=nsimEIM,
.llambda = llambda))))
} # lms.yjn2
lms.yjn <-
function(percentiles = c(25, 50, 75),
zero = c("lambda", "sigma"),
llambda = "identitylink",
lsigma = "loglink",
idf.mu = 4,
idf.sigma = 2,
ilambda = 1.0,
isigma = NULL,
rule = c(10, 5),
yoffset = NULL,
diagW = FALSE, iters.diagW = 6) {
llambda <- as.list(substitute(llambda))
elambda <- link2list(llambda)
llambda <- attr(elambda, "function.name")
lsigma <- as.list(substitute(lsigma))
esigma <- link2list(lsigma)
lsigma <- attr(esigma, "function.name")
rule <- rule[1] # No of pts (common) 4 all the quadrature schemes
if (rule != 5 && rule != 10)
stop("only rule=5 or 10 is supported")
new("vglmff",
blurb = c("LMS Quantile Regression ",
"(Yeo-Johnson transformation to normality)\n",
"Links: ",
namesof("lambda", link = llambda, earg = elambda), ", ",
namesof("mu", link = "identitylink", list()), ", ",
namesof("sigma", link = lsigma, earg = esigma)),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 3)
}), list(.zero = zero))),
infos = eval(substitute(function(...) {
list(M1 = 3,
Q1 = 1,
expected = TRUE,
multipleResponses = FALSE,
parameters.names = c("lambda", "mu", "sigma"),
llambda = .llambda ,
lmu = "identitylink",
lsigma = .lsigma ,
percentiles = .percentiles , # For the original fit only
true.mu = FALSE, # quantiles
zero = .zero )
}, list( .zero = zero,
.percentiles = percentiles,
.llambda = llambda, .lsigma = lsigma ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
ncol.w.max = 1, ncol.y.max = 1)
predictors.names <-
c(namesof("lambda", .llambda, earg = .elambda , short = TRUE),
"mu",
namesof("sigma", .lsigma, earg = .esigma , short = TRUE))
extra$percentiles <- .percentiles
y.save <- y
yoff <- if (is.Numeric( .yoffset )) .yoffset else -median(y)
extra$yoffset <- yoff
y <- y + yoff
if (!length(etastart)) {
lambda.init <- if (is.Numeric( .ilambda )) .ilambda else 1.0
y.tx <- yeo.johnson(y, lambda.init)
if (smoothok <-
(length(unique(sort(x[, min(ncol(x), 2)]))) > 7)) {
fit700 <- vsmooth.spline(x = x[, min(ncol(x), 2)],
y = y.tx, w = w, df = .idf.mu )
fv.init <- c(predict(fit700, x = x[, min(ncol(x), 2)])$y)
} else {
fv.init <- rep_len(weighted.mean(y, w), n)
}
sigma.init <- if (!is.Numeric( .isigma )) {
if (is.Numeric( .idf.sigma) &&
smoothok) {
fit710 = vsmooth.spline(x = x[, min(ncol(x), 2)],
y = (y.tx - fv.init)^2,
w = w, df = .idf.sigma)
sqrt(c(abs(predict(fit710,
x = x[, min(ncol(x), 2)])$y)))
} else {
sqrt( sum( w * (y.tx - fv.init)^2 ) / sum(w) )
}
} else
.isigma
etastart <-
cbind(theta2eta(lambda.init, .llambda , earg = .elambda ),
fv.init,
theta2eta(sigma.init, .lsigma , earg = .esigma ))
}
}), list( .llambda = llambda, .lsigma = lsigma,
.elambda = elambda, .esigma = esigma,
.ilambda = ilambda, .isigma = isigma,
.idf.mu = idf.mu,
.idf.sigma = idf.sigma,
.yoffset = yoffset,
.percentiles = percentiles ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
pcent <- extra$percentiles
eta[, 1] <- eta2theta(eta[, 1], .llambda , earg = .elambda )
eta[, 3] <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
qtplot.lms.yjn(percentiles = pcent,
eta = eta, yoffset = extra$yoff)
}, list(.percentiles = percentiles,
.esigma = esigma,
.elambda = elambda,
.llambda = llambda,
.lsigma = lsigma))),
last = eval(substitute(expression({
misc$link <- c(lambda = .llambda ,
mu = "identitylink",
sigma = .lsigma )
misc$earg <- list(lambda = .elambda ,
mu = list(theta = NULL),
sigma = .esigma )
misc$percentiles <- .percentiles # These are argument values
misc$true.mu <- FALSE # $fitted is not a true mu
misc[["yoffset"]] <- extra$yoff
y <- y.save # Restore back the value; to be attached to object
if (control$cdf) {
post$cdf =
cdf.lms.yjn(y + misc$yoffset,
eta0 = matrix(c(lambda,mymu,sigma),
ncol = 3,
dimnames = list(dimnames(x)[[1]], NULL)))
}
}), list( .percentiles = percentiles,
.elambda = elambda, .esigma = esigma,
.llambda = llambda, .lsigma = lsigma))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL, summation = TRUE) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mu <- eta[, 2]
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
psi <- yeo.johnson(y, lambda)
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * (-log(sigma) - 0.5 * ((psi-mu)/sigma)^2 +
(lambda-1) * sign(y) * log1p(abs(y)))
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .esigma = esigma, .elambda = elambda,
.lsigma = lsigma, .llambda = llambda))),
vfamily = c("lms.yjn", "lmscreg"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta[, 2]
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
okay1 <- all(is.finite(mymu )) &&
all(is.finite(sigma )) && all(0 < sigma) &&
all(is.finite(lambda))
okay1
}, list( .esigma = esigma, .elambda = elambda,
.lsigma = lsigma, .llambda = llambda))),
deriv = eval(substitute(expression({
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta[, 2]
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
psi <- yeo.johnson(y, lambda)
d1 <- yeo.johnson(y, lambda, deriv = 1)
AA <- (psi - mymu) / sigma
dl.dlambda <- -AA * d1 /sigma + sign(y) * log1p(abs(y))
dl.dmu <- AA / sigma
dl.dsigma <- (AA^2 -1) / sigma
dlambda.deta <- dtheta.deta(lambda, link = .llambda, .elambda )
dsigma.deta <- dtheta.deta(sigma, link = .lsigma, .esigma )
cbind(dl.dlambda * dlambda.deta,
dl.dmu,
dl.dsigma * dsigma.deta) * c(w)
}), list( .esigma = esigma, .elambda = elambda,
.lsigma = lsigma, .llambda = llambda ))),
weight = eval(substitute(expression({
wz <- matrix(0, n, 6)
wz[,iam(2, 2, M)] <- 1 / sigma^2
wz[,iam(3, 3, M)] <- 2 * wz[,iam(2, 2, M)] # 2 / sigma^2
if ( .rule == 10) {
glag.abs = c(0.13779347054,0.729454549503,
1.80834290174,3.40143369785,
5.55249614006,8.33015274676,
11.8437858379,16.2792578314,
21.996585812, 29.9206970123)
glag.wts = c(0.308441115765, 0.401119929155, 0.218068287612,
0.0620874560987, 0.00950151697517, 0.000753008388588,
2.82592334963e-5,
4.24931398502e-7, 1.83956482398e-9, 9.91182721958e-13)
} else {
glag.abs = c(0.2635603197180449, 1.4134030591060496,
3.5964257710396850,
7.0858100058570503, 12.6408008442729685)
glag.wts = c(5.217556105826727e-01, 3.986668110832433e-01,
7.594244968176882e-02,
3.611758679927785e-03, 2.336997238583738e-05)
}
if ( .rule == 10) {
sgh.abs = c(0.03873852801690856, 0.19823332465268367,
0.46520116404433082,
0.81686197962535023, 1.23454146277833154,
1.70679833036403172,
2.22994030591819214, 2.80910399394755972,
3.46387269067033854,
4.25536209637269280)
sgh.wts = c(9.855210713854302e-02, 2.086780884700499e-01,
2.520517066468666e-01,
1.986843323208932e-01,9.719839905023238e-02,
2.702440190640464e-02,
3.804646170194185e-03, 2.288859354675587e-04,
4.345336765471935e-06,
1.247734096219375e-08)
} else {
sgh.abs = c(0.1002421519682381, 0.4828139660462573,
1.0609498215257607,
1.7797294185202606, 2.6697603560875995)
sgh.wts = c(0.2484061520284881475,0.3923310666523834311,
0.2114181930760276606,
0.0332466603513424663, 0.0008248533445158026)
}
if ( .rule == 10) {
gleg.abs = c(-0.973906528517, -0.865063366689, -0.679409568299,
-0.433395394129, -0.148874338982)
gleg.abs = c(gleg.abs, rev(-gleg.abs))
gleg.wts = c(0.0666713443087, 0.149451349151, 0.219086362516,
0.26926671931, 0.295524224715)
gleg.wts = c(gleg.wts, rev(gleg.wts))
} else {
gleg.abs = c(-0.9061798459386643,-0.5384693101056820, 0,
0.5384693101056828, 0.9061798459386635)
gleg.wts = c(0.2369268850561853,0.4786286704993680,
0.5688888888888889,
0.4786286704993661, 0.2369268850561916)
}
discontinuity = -mymu/(sqrt(2)*sigma)
LL <- pmin(discontinuity, 0)
UU <- pmax(discontinuity, 0)
if (FALSE) {
AA <- (UU-LL)/2
for (kk in seq_along(gleg.wts)) {
temp1 <- AA * gleg.wts[kk]
abscissae <- (UU+LL)/2 + AA * gleg.abs[kk]
psi <- mymu + sqrt(2) * sigma * abscissae
temp9 <- dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 2)
temp9 <- cbind(temp9,
exp(-abscissae^2) / (sqrt(pi) * sigma^2))
wz[,iam(1, 1, M)] <- wz[,iam(1, 1, M)] + temp1 *
gleg.weight.yjn.11(abscissae, lambda, mymu, sigma, temp9)
wz[,iam(1, 2, M)] <- wz[,iam(1, 2, M)] + temp1 *
gleg.weight.yjn.12(abscissae, lambda, mymu, sigma, temp9)
wz[,iam(1, 3, M)] <- wz[,iam(1, 3, M)] + temp1 *
gleg.weight.yjn.13(abscissae, lambda, mymu, sigma, temp9)
}
} else {
temp9 <- .Fortran("yjngintf", as.double(LL),
as.double(UU),
as.double(gleg.abs), as.double(gleg.wts),
as.integer(n),
as.integer(length(gleg.abs)), as.double(lambda),
as.double(mymu), as.double(sigma),
answer = double(3*n),
eps = as.double(1.0e-5))$ans
dim(temp9) <- c(3, n)
wz[,iam(1, 1, M)] <- temp9[1,]
wz[,iam(1, 2, M)] <- temp9[2,]
wz[,iam(1, 3, M)] <- temp9[3,]
}
for (kk in seq_along(sgh.wts)) {
abscissae <- sign(-discontinuity) * sgh.abs[kk]
psi <- mymu + sqrt(2) * sigma * abscissae # abscissae = z
temp9 <- dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 2)
wz[,iam(1, 1, M)] <- wz[,iam(1, 1, M)] + sgh.wts[kk] *
gh.weight.yjn.11(abscissae, lambda, mymu, sigma, temp9)
wz[,iam(1, 2, M)] <- wz[,iam(1, 2, M)] + sgh.wts[kk] *
gh.weight.yjn.12(abscissae, lambda, mymu, sigma, temp9)
wz[,iam(1, 3, M)] <- wz[,iam(1, 3, M)] + sgh.wts[kk] *
gh.weight.yjn.13(abscissae, lambda, mymu, sigma, temp9)
}
temp1 <- exp(-discontinuity^2)
for (kk in seq_along(glag.wts)) {
abscissae <- sign(discontinuity) * sqrt(glag.abs[kk]) +
discontinuity^2
psi <- mymu + sqrt(2) * sigma * abscissae
temp9 <- dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 2)
temp9 <- cbind(temp9,
1 / (2 * sqrt((abscissae-discontinuity^2)^2 +
discontinuity^2) *
sqrt(pi) * sigma^2))
temp7 <- temp1 * glag.wts[kk]
wz[,iam(1, 1, M)] <- wz[,iam(1, 1, M)] + temp7 *
glag.weight.yjn.11(abscissae, lambda, mymu, sigma, temp9)
wz[,iam(1, 2, M)] <- wz[,iam(1, 2, M)] + temp7 *
glag.weight.yjn.12(abscissae, lambda, mymu, sigma, temp9)
wz[,iam(1, 3, M)] <- wz[,iam(1, 3, M)] + temp7 *
glag.weight.yjn.13(abscissae, lambda, mymu, sigma, temp9)
}
wz[, iam(1, 1, M)] <- wz[, iam(1, 1, M)] * dlambda.deta^2
wz[, iam(1, 2, M)] <- wz[, iam(1, 2, M)] * dlambda.deta
wz[, iam(1, 3, M)] <- wz[, iam(1, 3, M)] *
dsigma.deta * dlambda.deta
if ( .diagW && iter <= .iters.diagW) {
wz[,iam(1, 2, M)] <- wz[, iam(1, 3, M)] <- 0
}
wz[, iam(2, 3, M)] <- wz[, iam(2, 3, M)] * dsigma.deta
wz[, iam(3, 3, M)] <- wz[, iam(3, 3, M)] * dsigma.deta^2
c(w) * wz
}), list( .lsigma = lsigma, .llambda = llambda,
.esigma = esigma, .elambda = elambda,
.rule = rule, .diagW = diagW,
.iters.diagW = iters.diagW ))))
} # lms.yjn
lmscreg.control <-
function(cdf = TRUE, at.arg = NULL, x0 = NULL, ...) {
if (!is.logical(cdf)) {
warning("'cdf' is not logical; using TRUE instead")
cdf <- TRUE
}
list(cdf = cdf, at.arg = at.arg, x0 = x0)
}
Wr1 <- function(r, w) ifelse(r <= 0, 1, w)
Wr2 <- function(r, w) (r <= 0) * 1 + (r > 0) * w
amlnormal.deviance <-
function(mu, y, w, residuals = FALSE,
eta, extra = NULL) {
M <- length(extra$w.aml)
if (M > 1) y <- matrix(y, extra$n, extra$M)
devi <- cbind((y - mu)^2)
if (residuals) {
stop("not sure here")
wz <- VGAM.weights.function(w = w, M = extra$M, n = extra$n)
return((y - mu) * sqrt(wz) * matrix(extra$w.aml,extra$n,extra$M))
} else {
all.deviances <- numeric(M)
myresid <- matrix(y,extra$n,extra$M) - cbind(mu)
for (ii in 1:M)
all.deviances[ii] <- sum(c(w) * devi[, ii] *
Wr1(myresid[, ii],
w = extra$w.aml[ii]))
}
if (is.logical(extra$individual) && extra$individual)
all.deviances else sum(all.deviances)
} # amlnormal.deviance
amlnormal <-
function(w.aml = 1, parallel = FALSE,
lexpectile = "identitylink",
iexpectile = NULL,
imethod = 1, digw = 4) {
if (!is.Numeric(w.aml, positive = TRUE))
stop("argument 'w.aml' must be a vector of positive values")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 3)
stop("argument 'imethod' must be 1, 2 or 3")
lexpectile <- as.list(substitute(lexpectile))
eexpectile <- link2list(lexpectile)
lexpectile <- attr(eexpectile, "function.name")
if (length(iexpectile) && !is.Numeric(iexpectile))
stop("bad input for argument 'iexpectile'")
new("vglmff",
blurb = c("Asymmetric least squares quantile regression\n\n",
"Links: ",
namesof("expectile", link = lexpectile, eexpectile)),
constraints = eval(substitute(expression({
constraints <- cm.VGAM(matrix(1, M, 1), x = x,
bool = .parallel ,
constraints = constraints)
}), list( .parallel = parallel ))),
deviance = function(mu, y, w, residuals = FALSE, eta,
extra = NULL) {
amlnormal.deviance(mu = mu, y = y, w = w, residuals = residuals,
eta = eta, extra = extra)
},
initialize = eval(substitute(expression({
extra$w.aml <- .w.aml
temp5 <-
w.y.check(w = w, y = y,
ncol.w.max = 1, ncol.y.max = 1,
out.wy = TRUE,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
extra$M <- M <- length(extra$w.aml) # Recycle if necessary
extra$n <- n
extra$y.names <- y.names <-
paste0("w.aml = ", round(extra$w.aml, digits = .digw ))
predictors.names <- c(namesof(
paste0("expectile(",y.names,")"), .lexpectile ,
earg = .eexpectile, tag = FALSE))
if (!length(etastart)) {
mean.init <-
if ( .imethod == 1)
rep_len(median(y), n) else
if ( .imethod == 2 || .imethod == 3)
rep_len(weighted.mean(y, w), n) else {
junk <- lm.wfit(x = x, y = c(y), w = c(w))
junk$fitted
}
if ( .imethod == 3)
mean.init <- abs(mean.init) + 0.01
if (length( .iexpectile))
mean.init <- matrix( .iexpectile, n, M, byrow = TRUE)
etastart <-
matrix(theta2eta(mean.init, .lexpectile,
earg = .eexpectile), n, M)
}
}), list( .lexpectile = lexpectile, .eexpectile = eexpectile,
.iexpectile = iexpectile,
.imethod = imethod, .digw = digw, .w.aml = w.aml ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
ans <- eta <- as.matrix(eta)
for (ii in 1:ncol(eta))
ans[, ii] <- eta2theta(eta[, ii], .lexpectile , .eexpectile )
dimnames(ans) <- list(dimnames(eta)[[1]], extra$y.names)
ans
}, list( .lexpectile = lexpectile, .eexpectile = eexpectile ))),
last = eval(substitute(expression({
misc$link <- rep_len(.lexpectile , M)
names(misc$link) <- extra$y.names
misc$earg <- vector("list", M)
for (ilocal in 1:M)
misc$earg[[ilocal]] <- list(theta = NULL)
names(misc$earg) <- names(misc$link)
misc$parallel <- .parallel
misc$expected <- TRUE
extra$percentile <- numeric(M) # These are estimates (empirical)
misc$multipleResponses <- TRUE
for (ii in 1:M) {
use.w <- if (M > 1 && NCOL(w) == M) w[, ii] else w
extra$percentile[ii] <- 100 *
weighted.mean(myresid[, ii] <= 0, use.w)
}
names(extra$percentile) <- names(misc$link)
extra$individual <- TRUE
if (!(M > 1 && NCOL(w) == M)) {
extra$deviance <-
amlnormal.deviance(mu = mu, y = y, w = w,
residuals = FALSE, eta = eta,
extra = extra)
names(extra$deviance) <- extra$y.names
}
}), list( .lexpectile = lexpectile,
.eexpectile = eexpectile, .parallel = parallel ))),
vfamily = c("amlnormal"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
mymu <- eta2theta(eta, .lexpectile , earg = .eexpectile )
okay1 <- all(is.finite(mymu))
okay1
}, list( .lexpectile = lexpectile,
.eexpectile = eexpectile ))),
deriv = eval(substitute(expression({
mymu <- eta2theta(eta, .lexpectile , earg = .eexpectile )
dexpectile.deta <- dtheta.deta(mymu, .lexpectile ,
earg = .eexpectile )
myresid <- matrix(y,extra$n,extra$M) - cbind(mu)
wor1 <- Wr2(myresid, w = matrix(extra$w.aml, extra$n, extra$M,
byrow = TRUE))
c(w) * myresid * wor1 * dexpectile.deta
}), list( .lexpectile = lexpectile,
.eexpectile = eexpectile ))),
weight = eval(substitute(expression({
wz <- c(w) * wor1 * dexpectile.deta^2
wz
}), list( .lexpectile = lexpectile,
.eexpectile = eexpectile ))))
} # amlnormal
amlpoisson.deviance <-
function(mu, y, w, residuals = FALSE, eta,
extra = NULL) {
M <- length(extra$w.aml)
if (M > 1) y <- matrix(y,extra$n,extra$M)
nz <- y > 0
devi <- cbind(-(y - mu))
devi[nz] <- devi[nz] + y[nz] * log(y[nz]/mu[nz])
if (residuals) {
stop("not sure here")
return(sign(y - mu) * sqrt(2 * abs(devi) * w) *
matrix(extra$w,extra$n,extra$M))
} else {
all.deviances <- numeric(M)
myresid <- matrix(y,extra$n,extra$M) - cbind(mu)
for (ii in 1:M) all.deviances[ii] <- 2 *
sum(c(w) * devi[, ii] *
Wr1(myresid[, ii],
w=extra$w.aml[ii]))
}
if (is.logical(extra$individual) && extra$individual)
all.deviances else sum(all.deviances)
} # amlpoisson.deviance
amlpoisson <-
function(w.aml = 1, parallel = FALSE, imethod = 1,
digw = 4, link = "loglink") {
if (!is.Numeric(w.aml, positive = TRUE))
stop("'w.aml' must be a vector of positive values")
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
new("vglmff",
blurb = c("Poisson expectile regression by",
" asymmetric maximum likelihood estimation\n\n",
"Link: ", namesof("expectile", link, earg = earg)),
constraints = eval(substitute(expression({
constraints <- cm.VGAM(matrix(1, M, 1), x = x,
bool = .parallel ,
constraints = constraints)
}), list( .parallel = parallel ))),
deviance = function(mu, y, w, residuals = FALSE,
eta, extra = NULL) {
amlpoisson.deviance(mu = mu, y = y, w = w,
residuals = residuals,
eta = eta, extra = extra)
},
initialize = eval(substitute(expression({
extra$w.aml <- .w.aml
temp5 <-
w.y.check(w = w, y = y,
ncol.w.max = 1, ncol.y.max = 1,
out.wy = TRUE,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
extra$M <- M <- length(extra$w.aml) # Recycle if necessary
extra$n <- n
extra$y.names <- y.names <-
paste0("w.aml = ", round(extra$w.aml, digits = .digw ))
extra$individual <- FALSE
predictors.names <-
c(namesof(paste0("expectile(", y.names, ")"),
.link , earg = .earg , tag = FALSE))
if (!length(etastart)) {
mean.init <- if ( .imethod == 2)
rep_len(median(y), n) else
if ( .imethod == 1)
rep_len(weighted.mean(y, w), n) else {
junk = lm.wfit(x = x, y = c(y), w = c(w))
abs(junk$fitted)
}
etastart <-
matrix(theta2eta(mean.init, .link , earg = .earg ), n, M)
}
}), list( .link = link, .earg = earg, .imethod = imethod,
.digw = digw, .w.aml = w.aml ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
mu.ans <- eta <- as.matrix(eta)
for (ii in 1:ncol(eta))
mu.ans[, ii] <- eta2theta(eta[, ii], .link , earg = .earg )
dimnames(mu.ans) <- list(dimnames(eta)[[1]], extra$y.names)
mu.ans
}, list( .link = link, .earg = earg ))),
last = eval(substitute(expression({
misc$multipleResponses <- TRUE
misc$expected <- TRUE
misc$parallel <- .parallel
misc$link <- rep_len( .link , M)
names(misc$link) <- extra$y.names
misc$earg <- vector("list", M)
for (ilocal in 1:M)
misc$earg[[ilocal]] <- list(theta = NULL)
names(misc$earg) <- names(misc$link)
extra$percentile <- numeric(M) # These are estimates (empirical)
for (ii in 1:M)
extra$percentile[ii] <- 100 *
weighted.mean(myresid[, ii] <= 0, w)
names(extra$percentile) <- names(misc$link)
extra$individual <- TRUE
extra$deviance <- amlpoisson.deviance(mu = mu, y = y, w = w,
residuals = FALSE,
eta = eta, extra = extra)
names(extra$deviance) <- extra$y.names
}), list( .link = link, .earg = earg, .parallel = parallel ))),
linkfun = eval(substitute(function(mu, extra = NULL) {
theta2eta(mu, link = .link , earg = .earg )
}, list( .link = link, .earg = earg ))),
vfamily = c("amlpoisson"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
mymu <- eta2theta(eta, .link , earg = .earg )
okay1 <- all(is.finite(mymu)) && all(0 < mymu)
okay1
}, list( .link = link, .earg = earg ))),
deriv = eval(substitute(expression({
mymu <- eta2theta(eta, .link , earg = .earg )
dexpectile.deta <- dtheta.deta(mymu, .link , earg = .earg )
myresid <- matrix(y,extra$n,extra$M) - cbind(mu)
wor1 <- Wr2(myresid, w = matrix(extra$w.aml, extra$n, extra$M,
byrow = TRUE))
c(w) * myresid * wor1 * (dexpectile.deta / mymu)
}), list( .link = link, .earg = earg ))),
weight = eval(substitute(expression({
use.mu <- mymu
use.mu[use.mu < .Machine$double.eps^(3/4)] <-
.Machine$double.eps^(3/4)
wz <- c(w) * wor1 * use.mu * (dexpectile.deta / mymu)^2
wz
}), list( .link = link, .earg = earg ))))
} # amlpoisson
amlbinomial.deviance <-
function(mu, y, w, residuals = FALSE,
eta, extra = NULL) {
M <- length(extra$w.aml)
if (M > 1) y <- matrix(y,extra$n,extra$M)
devy <- y
nz <- y != 0
devy[nz] <- y[nz] * log(y[nz])
nz <- (1 - y) != 0
devy[nz] <- devy[nz] + (1 - y[nz]) * log1p(-y[nz])
devmu <- y * log(mu) + (1 - y) * log1p(-mu)
if (any(small <- mu * (1 - mu) < .Machine$double.eps)) {
warning("fitted values close to 0 or 1")
smu <- mu[small]
sy <- y[small]
smu <- ifelse(smu < .Machine$double.eps,
.Machine$double.eps, smu)
onemsmu <- ifelse((1 - smu) < .Machine$double.eps,
.Machine$double.eps, 1 - smu)
devmu[small] <- sy * log(smu) + (1 - sy) * log(onemsmu)
}
devi <- 2 * (devy - devmu)
if (residuals) {
stop("not sure here")
return(sign(y - mu) * sqrt(abs(devi) * w))
} else {
all.deviances <- numeric(M)
myresid <- matrix(y,extra$n,extra$M) -
matrix(mu,extra$n,extra$M)
for (ii in 1:M) all.deviances[ii] <- sum(c(w) * devi[, ii] *
Wr1(myresid[, ii], w=extra$w.aml[ii]))
}
if (is.logical(extra$individual) && extra$individual)
all.deviances else sum(all.deviances)
} # amlbinomial.deviance
amlbinomial <-
function(w.aml = 1, parallel = FALSE, digw = 4,
link = "logitlink") {
if (!is.Numeric(w.aml, positive = TRUE))
stop("'w.aml' must be a vector of positive values")
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
new("vglmff",
blurb = c("Logistic expectile regression by ",
"asymmetric maximum likelihood estimation\n\n",
"Link: ", namesof("expectile", link, earg = earg)),
constraints = eval(substitute(expression({
constraints <- cm.VGAM(matrix(1, M, 1), x = x,
bool = .parallel ,
constraints = constraints)
}), list( .parallel = parallel ))),
deviance = function(mu, y, w, residuals = FALSE,
eta, extra = NULL) {
amlbinomial.deviance(mu = mu, y = y, w = w,
residuals = residuals,
eta = eta, extra = extra)
},
initialize = eval(substitute(expression({
{
NCOL <- function (x)
if (is.array(x) && length(dim(x)) > 1 ||
is.data.frame(x)) ncol(x) else as.integer(1)
if (NCOL(y) == 1) {
if (is.factor(y)) y <- y != levels(y)[1]
nn <- rep_len(1, n)
if (!all(y >= 0 & y <= 1))
stop("response values must be in [0, 1]")
if (!length(mustart) && !length(etastart))
mustart <- (0.5 + w * y) / (1 + w)
no.successes <- w * y
if (any(abs(no.successes -
round(no.successes)) > 0.001))
stop("Number of successes must be integer-valued")
} else if (NCOL(y) == 2) {
if (any(abs(y - round(y)) > 0.001))
stop("Count data must be integer-valued")
nn <- y[, 1] + y[, 2]
y <- ifelse(nn > 0, y[, 1]/nn, 0)
w <- w * nn
if (!length(mustart) && !length(etastart))
mustart <- (0.5 + nn * y) / (1 + nn)
} else
stop("Response not of the right form")
}
extra$w.aml <- .w.aml
if (ncol(y <- cbind(y)) != 1)
stop("response must be a vector or a 1-column matrix")
extra$M <- M <- length(extra$w.aml) # Recycle if necessary
extra$n <- n
extra$y.names <- y.names <-
paste0("w.aml = ", round(extra$w.aml, digits = .digw ),
sep = "")
extra$individual <- FALSE
predictors.names <-
c(namesof(paste0("expectile(", y.names, ")"),
.link , earg = .earg , tag = FALSE))
if (!length(etastart)) {
etastart <- matrix(theta2eta(mustart, .link , .earg ),
n, M)
mustart <- NULL
}
}), list( .link = link, .earg = earg,
.digw = digw, .w.aml = w.aml ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
mu.ans <- eta <- as.matrix(eta)
for (ii in 1:ncol(eta))
mu.ans[, ii] <- eta2theta(eta[, ii], .link , earg = .earg )
dimnames(mu.ans) <- list(dimnames(eta)[[1]], extra$y.names)
mu.ans
}, list( .link = link, .earg = earg ))),
last = eval(substitute(expression({
misc$link <- rep_len(.link , M)
names(misc$link) <- extra$y.names
misc$earg <- vector("list", M)
for (ilocal in 1:M)
misc$earg[[ilocal]] <- list(theta = NULL)
names(misc$earg) <- names(misc$link)
misc$parallel <- .parallel
misc$expected <- TRUE
extra$percentile <- numeric(M) # These are estimates (empirical)
for (ii in 1:M)
extra$percentile[ii] <- 100 *
weighted.mean(myresid[, ii] <= 0, w)
names(extra$percentile) <- names(misc$link)
extra$individual <- TRUE
extra$deviance <- amlbinomial.deviance(mu = mu, y = y, w = w,
residuals = FALSE, eta = eta, extra = extra)
names(extra$deviance) <- extra$y.names
}), list( .link = link, .earg = earg, .parallel = parallel ))),
linkfun = eval(substitute(function(mu, extra = NULL) {
theta2eta(mu, link = .link , earg = .earg )
}, list( .link = link, .earg = earg ))),
vfamily = c("amlbinomial"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
mymu <- eta2theta(eta, .link , earg = .earg )
okay1 <- all(is.finite(mymu)) && all(0 < mymu & mymu < 1)
okay1
}, list( .link = link, .earg = earg ))),
deriv = eval(substitute(expression({
mymu <- eta2theta(eta, .link , earg = .earg )
use.mu <- mymu
use.mu[use.mu < .Machine$double.eps^(3/4)] <-
.Machine$double.eps^(3/4)
dexpectile.deta <- dtheta.deta(use.mu, .link , earg = .earg )
myresid <- matrix(y,extra$n,extra$M) - cbind(mu)
wor1 <- Wr2(myresid, w = matrix(extra$w.aml, extra$n, extra$M,
byrow = TRUE))
c(w) * myresid * wor1 * (dexpectile.deta / (use.mu * (1-use.mu)))
}), list( .link = link, .earg = earg ))),
weight = eval(substitute(expression({
wz <- c(w) * wor1 * (dexpectile.deta^2 / (use.mu * (1 - use.mu)))
wz
}), list( .link = link, .earg = earg))))
} # amlbinomial
amlexponential.deviance <-
function(mu, y, w, residuals = FALSE,
eta, extra = NULL) {
M <- length(extra$w.aml)
if (M > 1) y <- matrix(y,extra$n,extra$M)
devy <- cbind(-log(y) - 1)
devi <- cbind(-log(mu) - y / mu)
if (residuals) {
stop("not sure here")
return(sign(y - mu) * sqrt(2 * abs(devi) * w) *
matrix(extra$w,extra$n,extra$M))
} else {
all.deviances <- numeric(M)
myresid <- matrix(y,extra$n,extra$M) - cbind(mu)
for (ii in 1:M) all.deviances[ii] = 2 * sum(c(w) *
(devy[, ii] - devi[, ii]) *
Wr1(myresid[, ii], w=extra$w.aml[ii]))
}
if (is.logical(extra$individual) && extra$individual)
all.deviances else sum(all.deviances)
} # amlexponential.deviance
amlexponential <-
function(w.aml = 1, parallel = FALSE, imethod = 1,
digw = 4, link = "loglink") {
if (!is.Numeric(w.aml, positive = TRUE))
stop("'w.aml' must be a vector of positive values")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 3)
stop("argument 'imethod' must be 1, 2 or 3")
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
y.names <- paste0("w.aml = ", round(w.aml, digits = digw))
predictors.names <- c(namesof(
paste0("expectile(", y.names,")"), link, earg = earg))
predictors.names <- paste(predictors.names, collapse = ", ")
new("vglmff",
blurb = c("Exponential expectile regression by",
" asymmetric maximum likelihood estimation\n\n",
"Link: ", predictors.names),
constraints = eval(substitute(expression({
constraints <- cm.VGAM(matrix(1, M, 1), x = x,
bool = .parallel ,
constraints = constraints)
}), list( .parallel = parallel ))),
deviance = function(mu, y, w, residuals = FALSE,
eta, extra = NULL) {
amlexponential.deviance(mu = mu, y = y, w = w,
residuals = residuals,
eta = eta, extra = extra)
},
initialize = eval(substitute(expression({
extra$w.aml <- .w.aml
temp5 <-
w.y.check(w = w, y = y,
Is.positive.y = TRUE,
ncol.w.max = 1, ncol.y.max = 1,
out.wy = TRUE,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
extra$M <- M <- length(extra$w.aml) # Recycle if necessary
extra$n <- n
extra$y.names <- y.names <-
paste0("w.aml = ", round(extra$w.aml, digits = .digw ))
extra$individual = FALSE
predictors.names <- c(namesof(
paste0("expectile(", y.names, ")"),
.link , earg = .earg , tag = FALSE))
if (!length(etastart)) {
mean.init <- if ( .imethod == 1)
rep_len(median(y), n) else
if ( .imethod == 2)
rep_len(weighted.mean(y, w), n) else {
1 / (y + 1)
}
etastart <- matrix(theta2eta(mean.init, .link , .earg ),
n, M)
}
}), list( .link = link, .earg = earg, .imethod = imethod,
.digw = digw, .w.aml = w.aml ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
mu.ans <- eta <- as.matrix(eta)
for (ii in 1:ncol(eta))
mu.ans[, ii] <- eta2theta(eta[, ii], .link , earg = .earg )
dimnames(mu.ans) <- list(dimnames(eta)[[1]], extra$y.names)
mu.ans
}, list( .link = link, .earg = earg ))),
last = eval(substitute(expression({
misc$multipleResponses <- TRUE
misc$expected <- TRUE
misc$parallel <- .parallel
misc$link <- rep_len( .link , M)
names(misc$link) <- extra$y.names
misc$earg <- vector("list", M)
for (ilocal in 1:M)
misc$earg[[ilocal]] <- list(theta = NULL)
names(misc$earg) <- names(misc$link)
extra$percentile <- numeric(M) # These are estimates (empirical)
for (ii in 1:M)
extra$percentile[ii] <- 100 *
weighted.mean(myresid[, ii] <= 0, w)
names(extra$percentile) <- names(misc$link)
extra$individual <- TRUE
extra$deviance =
amlexponential.deviance(mu = mu, y = y, w = w,
residuals = FALSE,
eta = eta, extra = extra)
names(extra$deviance) <- extra$y.names
}), list( .link = link, .earg = earg, .parallel = parallel ))),
linkfun = eval(substitute(function(mu, extra = NULL) {
theta2eta(mu, link = .link , earg = .earg )
}, list( .link = link, .earg = earg ))),
vfamily = c("amlexponential"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
mymu <- eta2theta(eta, .link , earg = .earg )
okay1 <- all(is.finite(mymu)) && all(0 < mymu)
okay1
}, list( .link = link, .earg = earg ))),
deriv = eval(substitute(expression({
mymu <- eta2theta(eta, .link , earg = .earg )
bigy <- matrix(y,extra$n,extra$M)
dl.dmu <- (bigy - mymu) / mymu^2
dmu.deta <- dtheta.deta(mymu, .link , earg = .earg )
myresid <- bigy - cbind(mymu)
wor1 <- Wr2(myresid, w = matrix(extra$w.aml, extra$n, extra$M,
byrow = TRUE))
c(w) * wor1 * dl.dmu * dmu.deta
}), list( .link = link, .earg = earg ))),
weight = eval(substitute(expression({
ned2l.dmu2 <- 1 / mymu^2
wz <- c(w) * wor1 * ned2l.dmu2 * dmu.deta^2
wz
}), list( .link = link, .earg = earg ))))
} # amlexponential
benini1 <-
function(y0 = stop("argument 'y0' must be specified"),
lshape = "loglink",
ishape = NULL, imethod = 1, zero = NULL,
parallel = FALSE,
type.fitted = c("percentiles", "Qlink"),
percentiles = 50) {
type.fitted <- match.arg(type.fitted,
c("percentiles", "Qlink"))[1]
lshape <- as.list(substitute(lshape))
eshape <- link2list(lshape)
lshape <- attr(eshape, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
if (!is.Numeric(y0, positive = TRUE, length.arg = 1))
stop("bad input for argument 'y0'")
new("vglmff",
blurb = c("1-parameter Benini distribution\n\n",
"Link: ",
namesof("shape", lshape, earg = eshape),
"\n", "\n",
"Median: qbenini(p = 0.5, y0, shape)"),
constraints = eval(substitute(expression({
constraints <- cm.VGAM(matrix(1, M, 1), x = x,
bool = .parallel ,
constraints, apply.int = FALSE)
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 1)
}), list( .parallel = parallel,
.zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 1,
Q1 = 1,
expected = TRUE,
multipleResponses = TRUE,
parameters.names = c("shape"),
parallel = .parallel ,
percentiles = .percentiles ,
type.fitted = .type.fitted ,
lshape = .lshape ,
eshape = .eshape ,
zero = .zero )
}, list( .parallel = parallel,
.zero = zero,
.percentiles = percentiles ,
.type.fitted = type.fitted,
.eshape = eshape,
.lshape = lshape))),
initialize = eval(substitute(expression({
temp5 <-
w.y.check(w = w, y = y,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
ncoly <- ncol(y)
M1 <- 1
M <- M1 * ncoly
extra$ncoly <- ncoly
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
extra$percentiles <- .percentiles
extra$M1 <- M1
mynames1 <- paste0("shape", if (ncoly > 1) 1:ncoly else "")
predictors.names <-
namesof(mynames1, .lshape , earg = .eshape , tag = FALSE)
extra$y0 <- .y0 # Of unit length; 20181205; to make things easy.
if (any(y <= extra$y0))
stop("some values of the response are > argument 'y0' values")
if (!length(etastart)) {
probs.y <- (1:3) / 4
qofy <- quantile(rep(y, times = w), probs = probs.y)
if ( .imethod == 1) {
shape.init <- mean(-log1p(-probs.y) / (log(qofy))^2)
} else {
shape.init <- median(-log1p(-probs.y) / (log(qofy))^2)
}
shape.init <- matrix(if (length( .ishape )) .ishape else shape.init,
n, ncoly, byrow = TRUE)
etastart <- cbind(theta2eta(shape.init, .lshape , earg = .eshape ))
}
}), list( .imethod = imethod,
.ishape = ishape,
.lshape = lshape, .eshape = eshape,
.percentiles = percentiles,
.type.fitted = type.fitted,
.y0 = y0 ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
type.fitted <-
if (length(extra$type.fitted)) {
extra$type.fitted
} else {
warning("cannot find 'type.fitted'. Returning the 'median'.")
extra$percentiles <- 50 # Overwrite whatever was there
"percentiles"
}
type.fitted <- match.arg(type.fitted,
c("percentiles", "Qlink"))[1]
if (type.fitted == "Qlink") {
eta2theta(eta, link = "loglink")
} else {
shape <- eta2theta(eta, .lshape , earg = .eshape )
pcent <- extra$percentiles
perc.mat <- matrix(pcent, NROW(eta), length(pcent),
byrow = TRUE) / 100
fv <-
switch(type.fitted,
"percentiles" = qbenini(perc.mat,
y0 = extra$y0,
shape = matrix(shape, nrow(perc.mat), ncol(perc.mat))))
if (type.fitted == "percentiles")
fv <- label.cols.y(fv, colnames.y = extra$colnames.y,
NOS = NCOL(eta), percentiles = pcent,
one.on.one = FALSE)
fv
}
}, list( .lshape = lshape, .eshape = eshape ))),
last = eval(substitute(expression({
M1 <- extra$M1
misc$link <- c(rep_len( .lshape , ncoly))
names(misc$link) <- mynames1
misc$earg <- vector("list", M)
names(misc$earg) <- mynames1
for (ii in 1:ncoly) {
misc$earg[[ii]] <- .eshape
}
extra$y0 <- .y0
}), list( .lshape = lshape,
.eshape = eshape, .y0 = y0 ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL,
summation = TRUE) {
shape <- eta2theta(eta, .lshape , earg = .eshape )
y0 <- extra$y0
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * dbenini(y, y0 = y0, sh = shape, log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lshape = lshape, .eshape = eshape ))),
vfamily = c("benini1"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
shape <- eta2theta(eta, .lshape , earg = .eshape )
okay1 <- all(is.finite(shape)) && all(0 < shape)
okay1
}, list( .lshape = lshape, .eshape = eshape ))),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
extra <- object@extra
shape <- eta2theta(eta, .lshape , earg = .eshape )
y0 <- extra$y0
rbenini(nsim * length(shape), y0 = y0, shape = shape)
}, list( .lshape = lshape, .eshape = eshape ))),
deriv = eval(substitute(expression({
shape <- eta2theta(eta, .lshape , earg = .eshape )
y0 <- extra$y0
dl.dshape <- 1/shape - (log(y/y0))^2
dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )
c(w) * dl.dshape * dshape.deta
}), list( .lshape = lshape, .eshape = eshape ))),
weight = eval(substitute(expression({
ned2l.dshape2 <- 1 / shape^2
wz <- ned2l.dshape2 * dshape.deta^2
c(w) * wz
}), list( .lshape = lshape, .eshape = eshape ))))
} # benini1
dbenini <- function(x, y0, shape, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
N <- max(length(x), length(shape), length(y0))
if (length(x) != N) x <- rep_len(x, N)
if (length(shape) != N) shape <- rep_len(shape, N)
if (length(y0) != N) y0 <- rep_len(y0, N)
logdensity <- rep_len(log(0), N)
xok <- (x > y0)
tempxok <- log(x[xok]/y0[xok])
logdensity[xok] <- log(2*shape[xok]) - shape[xok] * tempxok^2 +
log(tempxok) - log(x[xok])
logdensity[is.infinite(x)] <- log(0) # 20141209 KaiH
if (log.arg) logdensity else exp(logdensity)
}
pbenini <-
function(q, y0, shape, lower.tail = TRUE, log.p = FALSE) {
if (!is.Numeric(q))
stop("bad input for argument 'q'")
if (!is.Numeric(shape, positive = TRUE))
stop("bad input for argument 'shape'")
if (!is.Numeric(y0, positive = TRUE))
stop("bad input for argument 'y0'")
if (!is.logical(lower.tail) || length(lower.tail ) != 1)
stop("bad input for argument 'lower.tail'")
if (!is.logical(log.p) || length(log.p) != 1)
stop("bad input for argument 'log.p'")
N <- max(length(q), length(shape), length(y0))
if (length(q) != N) q <- rep_len(q, N)
if (length(shape) != N) shape <- rep_len(shape, N)
if (length(y0) != N) y0 <- rep_len(y0, N)
ans <- y0 * 0
ok <- q > y0
if (lower.tail) {
if (log.p) {
ans[ok] <- log(-expm1(-shape[ok] * (log(q[ok]/y0[ok]))^2))
ans[q <= y0 ] <- -Inf
} else {
ans[ok] <- -expm1(-shape[ok] * (log(q[ok]/y0[ok]))^2)
}
} else {
if (log.p) {
ans[ok] <- -shape[ok] * (log(q[ok]/y0[ok]))^2
ans[q <= y0] <- 0
} else {
ans[ok] <- exp(-shape[ok] * (log(q[ok]/y0[ok]))^2)
ans[q <= y0] <- 1
}
}
ans
}
qbenini <- function(p, y0, shape, lower.tail = TRUE, log.p = FALSE) {
if (!is.logical(lower.tail) || length(lower.tail ) != 1)
stop("bad input for argument 'lower.tail'")
if (!is.logical(log.p) || length(log.p) != 1)
stop("bad input for argument 'log.p'")
if (lower.tail) {
if (log.p) {
ln.p <- p
ans <- y0 * exp(sqrt(-log(-expm1(ln.p)) / shape))
} else {
ans <- y0 * exp(sqrt(-log1p(-p) / shape))
}
} else {
if (log.p) {
ln.p <- p
ans <- y0 * exp(sqrt(-ln.p / shape))
} else {
ans <- y0 * exp(sqrt(-log(p) / shape))
}
}
ans[y0 <= 0] <- NaN
ans
}
rbenini <- function(n, y0, shape) {
y0 * exp(sqrt(-log(runif(n)) / shape))
}
hyperg <- function(N = NULL, D = NULL,
lprob = "logitlink",
iprob = NULL) {
inputN <- is.Numeric(N, positive = TRUE)
inputD <- is.Numeric(D, positive = TRUE)
if (inputD && inputN)
stop("only one of 'N' and 'D' is to be inputted")
if (!inputD && !inputN)
stop("one of 'N' and 'D' needs to be inputted")
lprob <- as.list(substitute(lprob))
earg <- link2list(lprob)
lprob <- attr(earg, "function.name")
new("vglmff",
blurb = c("Hypergeometric distribution\n\n",
"Link: ",
namesof("prob", lprob, earg = earg), "\n",
"Mean: D/N\n"),
initialize = eval(substitute(expression({
NCOL <- function (x)
if (is.array(x) && length(dim(x)) > 1 ||
is.data.frame(x)) ncol(x) else as.integer(1)
if (NCOL(y) == 1) {
if (is.factor(y)) y <- y != levels(y)[1]
nn <- rep_len(1, n)
if (!all(y >= 0 & y <= 1))
stop("response values must be in [0, 1]")
mustart <- (0.5 + w * y) / (1 + w)
no.successes <- w * y
if (any(abs(no.successes - round(no.successes)) > 0.001))
stop("Number of successes must be integer-valued")
} else if (NCOL(y) == 2) {
if (any(abs(y - round(y)) > 0.001))
stop("Count data must be integer-valued")
nn <- y[, 1] + y[, 2]
y <- ifelse(nn > 0, y[, 1]/nn, 0)
w <- w * nn
mustart <- (0.5 + nn * y) / (1 + nn)
mustart[mustart >= 1] <- 0.95
} else
stop("Response not of the right form")
predictors.names <-
namesof("prob", .lprob , earg = .earg , tag = FALSE)
extra$Nvector <- .N
extra$Dvector <- .D
extra$Nunknown <- length(extra$Nvector) == 0
if (!length(etastart)) {
init.prob <- if (length( .iprob))
rep_len( .iprob, n) else
mustart
etastart <- matrix(init.prob, n, NCOL(y))
}
}), list( .lprob = lprob, .earg = earg, .N = N, .D = D,
.iprob = iprob ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
eta2theta(eta, .lprob , earg = .earg )
}, list( .lprob = lprob, .earg = earg ))),
last = eval(substitute(expression({
misc$link <- c("prob" = .lprob)
misc$earg <- list("prob" = .earg )
misc$Dvector <- .D
misc$Nvector <- .N
}), list( .N = N, .D = D, .lprob = lprob, .earg = earg ))),
linkfun = eval(substitute(function(mu, extra = NULL) {
theta2eta(mu, .lprob, earg = .earg )
}, list( .lprob = lprob, .earg = earg ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL,
summation = TRUE) {
N <- extra$Nvector
Dvec <- extra$Dvector
prob <- mu
yvec <- w * y
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
if (extra$Nunknown) {
tmp12 <- Dvec * (1-prob) / prob
(lgamma(1+tmp12) + lgamma(1+Dvec/prob-w) -
lgamma(1+tmp12-w+yvec) - lgamma(1+Dvec/prob))
} else {
(lgamma(1+N*prob) + lgamma(1+N*(1-prob)) -
lgamma(1+N*prob-yvec) -
lgamma(1+N*(1-prob) -w + yvec))
}
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lprob = lprob, .earg = earg ))),
vfamily = c("hyperg"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
prob <- eta2theta(eta, .lprob , earg = .earg )
okay1 <- all(is.finite(prob)) && all(0 < prob & prob < 1)
okay1
}, list( .lprob = lprob, .earg = earg ))),
deriv = eval(substitute(expression({
prob <- mu # equivalently, eta2theta(eta, .lprob, .earg )
dprob.deta <- dtheta.deta(prob, .lprob, earg = .earg )
Dvec <- extra$Dvector
Nvec <- extra$Nvector
yvec <- w * y
if (extra$Nunknown) {
tmp72 <- -Dvec / prob^2
tmp12 <- Dvec * (1-prob) / prob
dl.dprob <- tmp72 * (digamma(1 + tmp12) +
digamma(1 + Dvec/prob -w) -
digamma(1 + tmp12-w+yvec) - digamma(1 + Dvec/prob))
} else {
dl.dprob <- Nvec * (digamma(1+Nvec*prob) -
digamma(1+Nvec*(1-prob)) -
digamma(1+Nvec*prob-yvec) +
digamma(1+Nvec*(1-prob)-w+yvec))
}
c(w) * dl.dprob * dprob.deta
}), list( .lprob = lprob, .earg = earg ))),
weight = eval(substitute(expression({
if (extra$Nunknown) {
tmp722 <- tmp72^2
tmp13 <- 2*Dvec / prob^3
d2l.dprob2 <- tmp722 * (trigamma(1 + tmp12) +
trigamma(1 + Dvec/prob - w) -
trigamma(1 + tmp12 - w + yvec) -
trigamma(1 + Dvec/prob)) +
tmp13 * (digamma(1 + tmp12) +
digamma(1 + Dvec/prob - w) -
digamma(1 + tmp12 - w + yvec) -
digamma(1 + Dvec/prob))
} else {
d2l.dprob2 <- Nvec^2 * (trigamma(1+Nvec*prob) +
trigamma(1+Nvec*(1-prob)) -
trigamma(1+Nvec*prob-yvec) -
trigamma(1+Nvec*(1-prob)-w+yvec))
}
d2prob.deta2 <- d2theta.deta2(prob, .lprob , earg = .earg )
wz <- -(dprob.deta^2) * d2l.dprob2
wz <- c(w) * wz
wz[wz < .Machine$double.eps] <- .Machine$double.eps
wz
}), list( .lprob = lprob, .earg = earg ))))
} # hyperg
dtriangle <- function(x, theta, lower = 0, upper = 1, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
N <- max(length(x), length(theta), length(lower), length(upper))
if (length(x) != N) x <- rep_len(x, N)
if (length(theta) != N) theta <- rep_len(theta, N)
if (length(lower) != N) lower <- rep_len(lower, N)
if (length(upper) != N) upper <- rep_len(upper, N)
denom1 <- ((upper-lower)*(theta-lower))
denom2 <- ((upper-lower)*(upper-theta))
logdensity <- rep_len(log(0), N)
xok.neg <- (lower < x) & (x <= theta)
xok.pos <- (theta <= x) & (x < upper)
logdensity[xok.neg] =
log(2 * (x[xok.neg] - lower[xok.neg]) / denom1[xok.neg])
logdensity[xok.pos] =
log(2 * (upper[xok.pos] - x[xok.pos]) / denom2[xok.pos])
logdensity[lower >= upper] <- NaN
logdensity[lower > theta] <- NaN
logdensity[upper < theta] <- NaN
if (log.arg) logdensity else exp(logdensity)
}
rtriangle <- function(n, theta, lower = 0, upper = 1) {
use.n <- if ((length.n <- length(n)) > 1) length.n else
if (!is.Numeric(n, integer.valued = TRUE,
length.arg = 1, positive = TRUE))
stop("bad input for argument 'n'") else n
if (!is.Numeric(theta))
stop("bad input for argument 'theta'")
if (!is.Numeric(lower))
stop("bad input for argument 'lower'")
if (!is.Numeric(upper))
stop("bad input for argument 'upper'")
if (!all(lower < theta & theta < upper))
stop("lower < theta < upper values are required")
N <- use.n
lower <- rep_len(lower, N)
upper <- rep_len(upper, N)
theta <- rep_len(theta, N)
t1 <- sqrt(runif(n))
t2 <- sqrt(runif(n))
ifelse(runif(n) < (theta - lower) / (upper - lower),
lower + (theta - lower) * t1,
upper - (upper - theta) * t2)
}
qtriangle <- function(p, theta, lower = 0, upper = 1,
lower.tail = TRUE, log.p = FALSE) {
if (!is.logical(lower.tail) || length(lower.tail ) != 1)
stop("bad input for argument 'lower.tail'")
if (!is.logical(log.p) || length(log.p) != 1)
stop("bad input for argument 'log.p'")
N <- max(length(p), length(theta), length(lower), length(upper))
if (length(p) != N) p <- rep_len(p, N)
if (length(theta) != N) theta <- rep_len(theta, N)
if (length(lower) != N) lower <- rep_len(lower, N)
if (length(upper) != N) upper <- rep_len(upper, N)
ans <- NA_real_ * p
if (lower.tail) {
if (log.p) {
Neg <- (exp(ln.p) <= (theta - lower) / (upper - lower))
temp1 <- exp(ln.p) * (upper - lower) * (theta - lower)
Pos <- (exp(ln.p) >= (theta - lower) / (upper - lower))
pstar <- (exp(ln.p) - (theta - lower) / (upper - lower)) /
((upper - theta) / (upper - lower))
} else {
Neg <- (p <= (theta - lower) / (upper - lower))
temp1 <- p * (upper - lower) * (theta - lower)
Pos <- (p >= (theta - lower) / (upper - lower))
pstar <- (p - (theta - lower) / (upper - lower)) /
((upper - theta) / (upper - lower))
}
} else {
if (log.p) {
ln.p <- p
Neg <- (exp(ln.p) >= (upper- theta) / (upper - lower))
temp1 <- -expm1(ln.p) * (upper - lower) * (theta - lower)
Pos <- (exp(ln.p) <= (upper- theta) / (upper - lower))
pstar <- (-expm1(ln.p) - (theta - lower) / (upper - lower)) /
((upper - theta) / (upper - lower))
} else {
Neg <- (p >= (upper- theta) / (upper - lower))
temp1 <- (1 - p) * (upper - lower) * (theta - lower)
Pos <- (p <= (upper- theta) / (upper - lower))
pstar <- ((upper- theta) / (upper - lower) - p) /
((upper - theta) / (upper - lower))
}
}
ans[ Neg] <- lower[ Neg] + sqrt(temp1[ Neg])
if (any(Pos)) {
qstar <- cbind(1 - sqrt(1-pstar), 1 + sqrt(1-pstar))
qstar <- qstar[Pos,, drop = FALSE]
qstar <- ifelse(qstar[, 1] >= 0 & qstar[, 1] <= 1,
qstar[, 1],
qstar[, 2])
ans[Pos] <- theta[Pos] + qstar * (upper - theta)[Pos]
}
ans[theta < lower | theta > upper] <- NaN
ans
}
ptriangle <- function(q, theta, lower = 0, upper = 1,
lower.tail = TRUE, log.p = FALSE) {
N <- max(length(q), length(theta), length(lower), length(upper))
if (length(q) != N) q <- rep_len(q, N)
if (length(theta) != N) theta <- rep_len(theta, N)
if (length(lower) != N) lower <- rep_len(lower, N)
if (length(upper) != N) upper <- rep_len(upper, N)
if (!is.logical(lower.tail) || length(lower.tail ) != 1)
stop("bad input for argument 'lower.tail'")
if (!is.logical(log.p) || length(log.p) != 1)
stop("bad input for argument 'log.p'")
ans <- q * 0
qstar <- (q - lower)^2 / ((upper - lower) * (theta - lower))
Neg <- (lower <= q & q <= theta)
ans[Neg] <- if (lower.tail) {
if (log.p) {
(log(qstar))[Neg]
} else {
qstar[Neg]
}
} else {
if (log.p) {
(log1p(-qstar))[Neg]
} else {
1 - qstar[Neg]
}
}
Pos <- (theta <= q & q <= upper)
qstar <- (q - theta) / (upper-theta)
if (lower.tail) {
if (log.p) {
ans[Pos] <- log(((theta-lower)/(upper-lower))[Pos] +
(qstar * (2-qstar) *
(upper-theta) / (upper - lower))[Pos])
ans[q <= lower] <- -Inf
ans[q >= upper] <- 0
} else {
ans[Pos] <- ((theta-lower)/(upper-lower))[Pos] +
(qstar * (2-qstar) *
(upper-theta) / (upper - lower))[Pos]
ans[q <= lower] <- 0
ans[q >= upper] <- 1
}
} else {
if (log.p) {
ans[Pos] <- log(((upper - theta)/(upper-lower))[Pos] +
(qstar * (2-qstar) *
(upper-theta) / (upper - lower))[Pos])
ans[q <= lower] <- 0
ans[q >= upper] <- -Inf
} else {
ans[Pos] <- ((upper - theta)/(upper-lower))[Pos] +
(qstar * (2-qstar) *
(upper-theta) / (upper - lower))[Pos]
ans[q <= lower] <- 1
ans[q >= upper] <- 0
}
}
ans[theta < lower | theta > upper] <- NaN
ans
}
triangle.control <- function(stepsize = 0.33, maxit = 100, ...) {
list(stepsize = stepsize, maxit = maxit)
}
triangle <-
function(lower = 0, upper = 1,
link = extlogitlink(min = 0, max = 1),
itheta = NULL) {
if (!is.Numeric(lower))
stop("bad input for argument 'lower'")
if (!is.Numeric(upper))
stop("bad input for argument 'upper'")
if (!all(lower < upper))
stop("lower < upper values are required")
if (length(itheta) && !is.Numeric(itheta))
stop("bad input for 'itheta'")
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
if (length(earg$min) && any(earg$min != lower))
stop("argument 'lower' does not match the 'link'")
if (length(earg$max) && any(earg$max != upper))
stop("argument 'upper' does not match the 'link'")
new("vglmff",
blurb = c("Triangle distribution\n\n",
"Link: ",
namesof("theta", link, earg = earg)),
infos = eval(substitute(function(...) {
list(M1 = 1,
Q1 = 1,
parameters.names = c("theta"),
link = .link )
}, list( .link = link ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
ncol.w.max = 1,
ncol.y.max = 1)
extra$lower <- rep_len( .lower , n)
extra$upper <- rep_len( .upper , n)
if (any(y <= extra$lower | y >= extra$upper))
stop("some y values in [lower,upper] detected")
predictors.names <-
namesof("theta", .link , earg = .earg , tag = FALSE)
if (!length(etastart)) {
Theta.init <- if (length( .itheta )) .itheta else {
weighted.mean(y, w)
}
Theta.init <- rep_len(Theta.init, n)
etastart <- theta2eta(Theta.init, .link , earg = .earg )
}
}), list( .link = link, .earg = earg, .itheta=itheta,
.upper = upper, .lower = lower ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
Theta <- eta2theta(eta, .link , earg = .earg )
lower <- extra$lower
upper <- extra$upper
mu1 <- (lower + upper + Theta) / 3
mu1
}, list( .link = link, .earg = earg ))),
last = eval(substitute(expression({
misc$link <- c(theta = .link )
misc$earg <- list(theta = .earg )
misc$expected <- TRUE
}), list( .link = link, .earg = earg ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL,
summation = TRUE) {
Theta <- eta2theta(eta, .link , earg = .earg )
lower <- extra$lower
upper <- extra$upper
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * dtriangle(y, theta = Theta, lower = lower,
upper = upper, log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .link = link, .earg = earg ))),
vfamily = c("triangle"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
Theta <- eta2theta(eta, .link , earg = .earg )
okay1 <- all(is.finite(Theta)) &&
all(extra$lower < Theta & Theta < extra$upper)
okay1
}, list( .link = link, .earg = earg ))),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
extra <- object@extra
Theta <- eta2theta(eta, .link , earg = .earg )
lower <- extra$lower
upper <- extra$upper
rtriangle(nsim * length(Theta),
theta = Theta, lower = lower, upper = upper)
}, list( .link = link, .earg = earg ))),
deriv = eval(substitute(expression({
Theta <- eta2theta(eta, .link , earg = .earg )
dTheta.deta <- dtheta.deta(Theta, .link , earg = .earg )
pos <- y > Theta
neg <- y < Theta
lower <- extra$lower
upper <- extra$upper
dl.dTheta <- 0 * y
dl.dTheta[neg] <- -1 / (Theta[neg]-lower[neg])
dl.dTheta[pos] <- 1 / (upper[pos]-Theta[pos])
c(w) * dl.dTheta * dTheta.deta
}), list( .link = link, .earg = earg ))),
weight = eval(substitute(expression({
var.dl.dTheta <- 1 / ((Theta - lower) * (upper - Theta))
wz <- var.dl.dTheta * dTheta.deta^2
c(w) * wz
}), list( .link = link, .earg = earg ))))
}
dpolono <-
function (x, meanlog = 0, sdlog = 1, bigx = 170, ...) {
mapply(function(x, meanlog, sdlog, ...) {
if (abs(x) > floor(x)) { # zero prob for -ve or non-integer
0
} else
if (x == Inf) { # 20141215 KaiH
0
} else
if (x > bigx) {
z <- (log(x) - meanlog) / sdlog
(1 + (z^2 + log(x) - meanlog - 1) / (2 * x * sdlog^2)) *
exp(-0.5 * z^2) / (sqrt(2 * pi) * sdlog * x)
} else
integrate( function(t) exp(t * x - exp(t) -
0.5 * ((t - meanlog) / sdlog)^2),
lower = -Inf,
upper = Inf, ...)$value / (sqrt(2 * pi) *
sdlog * exp(lgamma(x + 1.0)))
}, x, meanlog, sdlog, ...)
}
ppolono <-
function(q, meanlog = 0, sdlog = 1,
isOne = 1 - sqrt( .Machine$double.eps ), ...) {
.cumprob <- rep_len(0, length(q))
.cumprob[q == Inf] <- 1 # special case
q <- floor(q)
ii <- -1
while (any(xActive <- ((.cumprob < isOne) & (q > ii))))
.cumprob[xActive] <- .cumprob[xActive] +
dpolono(ii <- (ii+1), meanlog, sdlog, ...)
.cumprob
}
rpolono <- function(n, meanlog = 0, sdlog = 1) {
lambda <- rlnorm(n = n, meanlog = meanlog, sdlog = sdlog)
rpois(n = n, lambda = lambda)
}
fff.control <- function(save.weights = TRUE, ...) {
list(save.weights = save.weights)
}
fff <-
function(link = "loglink",
idf1 = NULL, idf2 = NULL, nsimEIM = 100, # ncp = 0,
imethod = 1, zero = NULL) {
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
if (!is.Numeric(nsimEIM, length.arg = 1,
integer.valued = TRUE) ||
nsimEIM <= 10)
stop("argument 'nsimEIM' should be an integer greater than 10")
ncp <- 0
if (any(ncp != 0))
warning("not sure about ncp != 0 wrt dl/dtheta")
new("vglmff",
blurb = c("F-distribution\n\n",
"Links: ",
namesof("df1", link, earg = earg), ", ",
namesof("df2", link, earg = earg),
"\n", "\n",
"Mean: df2/(df2-2) provided df2>2 and ncp = 0", "\n",
"Variance: ",
"2*df2^2*(df1+df2-2)/(df1*(df2-2)^2*(df2-4)) ",
"provided df2>4 and ncp = 0"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 2)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
multipleResponses = FALSE,
parameters.names = c("df1", "df2"),
zero = .zero )
}, list( .zero = zero ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
ncol.w.max = 1,
ncol.y.max = 1)
predictors.names <-
c(namesof("df1", .link , earg = .earg , tag = FALSE),
namesof("df2", .link , earg = .earg , tag = FALSE))
if (!length(etastart)) {
if ( .imethod == 1) {
df2.init <- b <- 2*mean(y) / (mean(y)-1)
df1.init <- 2*b^2*(b-2)/(var(y)*(b-2)^2 * (b-4) - 2*b^2)
if (df2.init < 4) df2.init <- 5
if (df1.init < 2) df1.init <- 3
} else {
df2.init <- b <- 2*median(y) / (median(y)-1)
summy <- summary(y)
var.est <- summy[5] - summy[2]
df1.init <- 2*b^2*(b-2)/(var.est*(b-2)^2 * (b-4) - 2*b^2)
}
df1.init <- if (length( .idf1 ))
rep_len( .idf1 , n) else
rep_len(df1.init, n)
df2.init <- if (length( .idf2 ))
rep_len( .idf2 , n) else
rep_len(1, n)
etastart <- cbind(theta2eta(df1.init, .link , earg = .earg ),
theta2eta(df2.init, .link , earg = .earg ))
}
}), list( .imethod = imethod, .idf1 = idf1, .earg = earg,
.idf2 = idf2, .link = link ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
df2 <- eta2theta(eta[, 2], .link , earg = .earg )
ans <- df2 * NA
ans[df2 > 2] <- df2[df2 > 2] / (df2[df2 > 2] - 2)
ans
}, list( .link = link, .earg = earg ))),
last = eval(substitute(expression({
misc$link <- c(df1 = .link , df2 = .link )
misc$earg <- list(df1 = .earg , df2 = .earg )
misc$nsimEIM <- .nsimEIM
misc$ncp <- .ncp
}), list( .link = link, .earg = earg,
.ncp = ncp,
.nsimEIM = nsimEIM ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
df1 <- eta2theta(eta[, 1], .link , earg = .earg )
df2 <- eta2theta(eta[, 2], .link , earg = .earg )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * df(x = y, df1 = df1, df2 = df2,
ncp = .ncp , log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .link = link, .earg = earg, .ncp = ncp ))),
vfamily = c("fff"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
df1 <- eta2theta(eta[, 1], .link , earg = .earg )
df2 <- eta2theta(eta[, 2], .link , earg = .earg )
okay1 <- all(is.finite(df1)) && all(0 < df1) &&
all(is.finite(df2)) && all(0 < df2)
okay1
}, list( .link = link, .earg = earg, .ncp = ncp ))),
deriv = eval(substitute(expression({
df1 <- eta2theta(eta[, 1], .link , earg = .earg )
df2 <- eta2theta(eta[, 2], .link , earg = .earg )
dl.ddf1 <- 0.5*digamma(0.5*(df1+df2)) +
0.5 + 0.5*log(df1/df2) +
0.5*log(y) - 0.5*digamma(0.5*df1) -
0.5*(df1+df2)*(y/df2) / (1 + df1*y/df2) -
0.5*log1p(df1*y/df2)
dl.ddf2 <- 0.5*digamma(0.5*(df1+df2)) - 0.5*df1/df2 -
0.5*digamma(0.5*df2) -
0.5*(df1+df2) * (-df1*y/df2^2) / (1 + df1*y/df2) -
0.5*log1p(df1*y/df2)
ddf1.deta <- dtheta.deta(df1, .link , earg = .earg )
ddf2.deta <- dtheta.deta(df2, .link , earg = .earg )
dthetas.detas <- cbind(ddf1.deta, ddf2.deta)
c(w) * dthetas.detas * cbind(dl.ddf1, dl.ddf2)
}), list( .link = link, .earg = earg ))),
weight = eval(substitute(expression({
run.varcov <- 0
ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)
for (ii in 1:( .nsimEIM )) {
ysim <- rf(n = n, df1=df1, df2=df2)
dl.ddf1 <- 0.5*digamma(0.5*(df1+df2)) +
0.5 + 0.5*log(df1/df2) +
0.5*log(ysim) - 0.5*digamma(0.5*df1) -
0.5*(df1+df2)*(ysim/df2) / (1 + df1*ysim/df2) -
0.5*log1p(df1*ysim/df2)
dl.ddf2 <- 0.5*digamma(0.5*(df1+df2)) -
0.5*df1/df2 -
0.5*digamma(0.5*df2) -
0.5*(df1+df2) *
(-df1*ysim/df2^2)/(1 + df1*ysim/df2) -
0.5*log1p(df1*ysim/df2)
rm(ysim)
temp3 <- cbind(dl.ddf1, dl.ddf2)
run.varcov <- ((ii-1) * run.varcov +
temp3[,ind1$row.index]*
temp3[,ind1$col.index]) / ii
}
wz <- if (intercept.only)
matrix(colMeans(run.varcov),
n, ncol(run.varcov), byrow = TRUE) else run.varcov
wz <- c(w) * wz * dthetas.detas[, ind1$row] *
dthetas.detas[, ind1$col]
wz
}), list( .link = link, .earg = earg, .nsimEIM = nsimEIM,
.ncp = ncp ))))
} # fff
dlaplace <-
function(x, location = 0, scale = 1, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
logdensity <- (-abs(x-location)/scale) - log(2*scale)
if (log.arg) logdensity else exp(logdensity)
} # dlaplace
plaplace <-
function(q, location = 0, scale = 1,
lower.tail = TRUE, log.p =FALSE) {
zedd <- (q - location) / scale
if (!is.logical(lower.tail) || length(lower.tail ) != 1)
stop("bad input for argument 'lower.tail'")
if (!is.logical(log.p) || length(log.p) != 1)
stop("bad input for argument 'log.p'")
L <- max(length(q), length(location), length(scale))
if (length(q) != L) q <- rep_len(q, L)
if (length(location) != L) location <- rep_len(location, L)
if (length(scale) != L) scale <- rep_len(scale, L)
if (lower.tail) {
if (log.p) {
ans <- ifelse(q < location, log(0.5) + zedd,
log1p(- 0.5 * exp(-zedd)))
} else {
ans <- ifelse(q < location, 0.5 * exp(zedd),
1 - 0.5 * exp(-zedd))
}
} else {
if (log.p) {
ans <- ifelse(q < location, log1p(- 0.5 * exp(zedd)),
log(0.5) - zedd)
} else {
ans <- ifelse(q < location, 1 - 0.5 *
exp(zedd), 0.5 * exp(-zedd))
}
}
ans[scale <= 0] <- NaN
ans
} # plaplace
qlaplace <-
function(p, location = 0, scale = 1,
lower.tail = TRUE, log.p = FALSE) {
if (!is.logical(lower.tail) || length(lower.tail ) != 1)
stop("bad input for argument 'lower.tail'")
if (!is.logical(log.p) || length(log.p) != 1)
stop("bad input for argument 'log.p'")
L <- max(length(p), length(location), length(scale))
if (length(p) != L) p <- rep_len(p, L)
if (length(location) != L) location <- rep_len(location, L)
if (length(scale) != L) scale <- rep_len(scale, L)
if (lower.tail) {
if (log.p) {
ln.p <- p
ans <- location - sign(exp(ln.p)-0.5) * scale *
log(2 * ifelse(exp(ln.p) < 0.5,
exp(ln.p), -expm1(ln.p)))
} else {
ans <- location - sign(p-0.5) * scale *
log(2 * ifelse(p < 0.5, p, 1-p))
}
} else {
if (log.p) {
ln.p <- p
ans <- location - sign(0.5 - exp(ln.p)) * scale *
log(2 * ifelse(-expm1(ln.p) < 0.5,
-expm1(ln.p), exp(ln.p)))
# ans[ln.p > 0] <- NaN
} else {
ans <- location - sign(0.5 - p) * scale *
log(2 * ifelse(p > 0.5, 1 - p, p))
}
}
ans[scale <= 0] <- NaN
ans
} # qlaplace
rlaplace <-
function(n, location = 0, scale = 1) {
use.n <- if ((length.n <- length(n)) > 1) length.n else
if (!is.Numeric(n, integer.valued = TRUE,
length.arg = 1, positive = TRUE))
stop("bad input for argument 'n'") else n
if (!is.Numeric(scale, positive = TRUE))
stop("'scale' must be positive")
location <- rep_len(location, use.n)
scale <- rep_len(scale, use.n)
rrrr <- runif(use.n)
location - sign(rrrr - 0.5) * scale *
(log(2) + ifelse(rrrr < 0.5, log(rrrr), log1p(-rrrr)))
} # rlaplace
laplace <-
function(llocation = "identitylink",
lscale = "loglink",
ilocation = NULL, iscale = NULL,
imethod = 1,
zero = "scale") {
llocat <- as.list(substitute(llocation))
elocat <- link2list(llocat)
llocat <- attr(elocat, "function.name")
ilocat <- ilocation
lscale <- as.list(substitute(lscale))
escale <- link2list(lscale)
lscale <- attr(escale, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 3)
stop("argument 'imethod' must be 1 or 2 or 3")
if (length(iscale) &&
!is.Numeric(iscale, positive = TRUE))
stop("bad input for argument 'iscale'")
new("vglmff",
blurb = c("Two-parameter Laplace distribution\n\n",
"Links: ",
namesof("location", llocat, earg = elocat), ", ",
namesof("scale", lscale, earg = escale),
"\n", "\n",
"Mean: location", "\n",
"Variance: 2*scale^2"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 2)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
multipleResponses = FALSE,
parameters.names = c("location", "scale"),
summary.pvalues = FALSE,
zero = .zero )
}, list( .zero = zero ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
ncol.w.max = 1,
ncol.y.max = 1)
predictors.names <-
c(namesof("location", .llocat , earg = .elocat, tag = FALSE),
namesof("scale", .lscale , earg = .escale, tag = FALSE))
if (!length(etastart)) {
if ( .imethod == 1) {
locat.init <- median(y)
scale.init <- sqrt(var(y) / 2)
} else if ( .imethod == 2) {
locat.init <- weighted.mean(y, w)
scale.init <- sqrt(var(y) / 2)
} else {
locat.init <- median(y)
scale.init <- sqrt(sum(c(w)*abs(y-median(y ))) / (sum(w) *2))
}
locat.init <- if (length( .ilocat ))
rep_len( .ilocat , n) else
rep_len(locat.init, n)
scale.init <- if (length( .iscale ))
rep_len( .iscale , n) else
rep_len(scale.init, n)
etastart <-
cbind(theta2eta(locat.init, .llocat , earg = .elocat ),
theta2eta(scale.init, .lscale , earg = .escale ))
}
}), list( .imethod = imethod,
.elocat = elocat, .escale = escale,
.llocat = llocat, .lscale = lscale,
.ilocat = ilocat, .iscale = iscale ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
eta2theta(eta[, 1], .llocat , earg = .elocat )
}, list( .elocat = elocat, .llocat = llocat ))),
last = eval(substitute(expression({
misc$link <- c(location = .llocat , scale = .lscale )
misc$earg <- list(location = .elocat , scale = .escale )
misc$expected <- TRUE
misc$RegCondOK <- FALSE # Save this for later
}), list( .escale = escale, .lscale = lscale,
.elocat = elocat, .llocat = llocat ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
locat <- eta2theta(eta[, 1], .llocat , earg = .elocat )
Scale <- eta2theta(eta[, 2], .lscale , earg = .escale )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * dlaplace(x = y, locat = locat,
scale = Scale, log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .escale = escale, .lscale = lscale,
.elocat = elocat, .llocat = llocat ))),
vfamily = c("laplace"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
Locat <- eta2theta(eta[, 1], .llocat , earg = .elocat )
Scale <- eta2theta(eta[, 2], .lscale , earg = .escale )
okay1 <- all(is.finite(Locat)) &&
all(is.finite(Scale)) && all(0 < Scale)
okay1
}, list( .escale = escale, .lscale = lscale,
.elocat = elocat, .llocat = llocat ))),
deriv = eval(substitute(expression({
Locat <- eta2theta(eta[, 1], .llocat , earg = .elocat )
Scale <- eta2theta(eta[, 2], .lscale , earg = .escale )
zedd <- abs(y-Locat) / Scale
dl.dLocat <- sign(y - Locat) / Scale
dl.dscale <- zedd / Scale - 1 / Scale
dLocat.deta <- dtheta.deta(Locat, .llocat , earg = .elocat )
dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale )
c(w) * cbind(dl.dLocat * dLocat.deta,
dl.dscale * dscale.deta)
}), list( .escale = escale, .lscale = lscale,
.elocat = elocat, .llocat = llocat ))),
weight = eval(substitute(expression({
d2l.dLocat2 <- d2l.dscale2 <- 1 / Scale^2
wz <- matrix(0, nrow = n, ncol = M) # diagonal
wz[,iam(1, 1, M)] <- d2l.dLocat2 * dLocat.deta^2
wz[,iam(2, 2, M)] <- d2l.dscale2 * dscale.deta^2
c(w) * wz
}), list( .escale = escale, .lscale = lscale,
.elocat = elocat, .llocat = llocat ))))
} # laplace
Try the VGAM package in your browser
Any scripts or data that you put into this service are public.
VGAM documentation built on Sept. 18, 2024, 9:09 a.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 laplace rlaplace qlaplace plaplace dlaplace fff fff.control rpolono ppolono triangle triangle.control ptriangle qtriangle rtriangle dtriangle hyperg rbenini qbenini pbenini dbenini benini1 amlexponential amlexponential.deviance amlbinomial amlbinomial.deviance amlpoisson amlpoisson.deviance amlnormal amlnormal.deviance Wr2 Wr1 lmscreg.control lms.yjn lms.yjn2 lms.yjn2.control gleg.weight.yjn.13 gleg.weight.yjn.12 gleg.weight.yjn.11 glag.weight.yjn.13 glag.weight.yjn.12 glag.weight.yjn.11 gh.weight.yjn.13 gh.weight.yjn.12 gh.weight.yjn.11 dpsi.dlambda.yjn yeo.johnson dyj.dy.yeojohnson dy.dpsi.yeojohnson lms.bcg lms.bcn lms.yjn.control qlms.bcn dlms.bcn
Documented in amlbinomial amlexponential amlnormal amlpoisson benini1 dbenini dlaplace dlms.bcn dpsi.dlambda.yjn dtriangle fff fff.control gh.weight.yjn.11 gh.weight.yjn.12 gh.weight.yjn.13 glag.weight.yjn.11 glag.weight.yjn.12 glag.weight.yjn.13 gleg.weight.yjn.11 gleg.weight.yjn.12 gleg.weight.yjn.13 hyperg hyperg laplace lms.bcg lms.bcn lmscreg.control lms.yjn lms.yjn2 lms.yjn.control pbenini plaplace ppolono ptriangle qbenini qlaplace qlms.bcn qtriangle rbenini rlaplace rpolono rtriangle Wr1 Wr2 yeo.johnson
# These functions are
# Copyright (C) 1998-2024 T.W. Yee, University of Auckland.
# All rights reserved.
dlms.bcn <- function(x, lambda = 1, mu = 0, sigma = 1,
tol0 = 0.001, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
zedd <- ((x/mu)^lambda - 1) / (lambda * sigma)
log.dz.dy <- (lambda - 1) * log(x/mu) - log(mu * sigma)
is.eff.0 <- abs(lambda) < tol0
if (any(is.eff.0)) {
zedd[is.eff.0] <- log(x[is.eff.0] /
mu[is.eff.0]) / sigma[is.eff.0]
log.dz.dy[is.eff.0] <- -log(x[is.eff.0] * sigma[is.eff.0])
}
logden <- dnorm(zedd, log = TRUE) + log.dz.dy
if (log.arg) logden else exp(logden)
}
qlms.bcn <- function(p, lambda = 1, mu = 0, sigma = 1) {
answer <- mu * (1 + lambda * sigma * qnorm(p))^(1/lambda)
answer
}
lms.bcn.control <-
lms.bcg.control <-
lms.yjn.control <- function(trace = TRUE, ...)
list(trace = trace)
lms.bcn <- function(percentiles = c(25, 50, 75),
zero = c("lambda", "sigma"),
llambda = "identitylink",
lmu = "identitylink",
lsigma = "loglink",
idf.mu = 4,
idf.sigma = 2,
ilambda = 1,
isigma = NULL,
tol0 = 0.001) {
llambda <- as.list(substitute(llambda))
elambda <- link2list(llambda)
llambda <- attr(elambda, "function.name")
lmu <- as.list(substitute(lmu))
emu <- link2list(lmu)
lmu <- attr(emu, "function.name")
lsigma <- as.list(substitute(lsigma))
esigma <- link2list(lsigma)
lsigma <- attr(esigma, "function.name")
if (!is.Numeric(tol0, positive = TRUE, length.arg = 1))
stop("bad input for argument 'tol0'")
if (!is.Numeric(ilambda))
stop("bad input for argument 'ilambda'")
if (length(isigma) &&
!is.Numeric(isigma, positive = TRUE))
stop("bad input for argument 'isigma'")
new("vglmff",
blurb = c("LMS ",
"quantile",
" regression (Box-Cox transformation to normality)\n",
"Links: ",
namesof("lambda", link = llambda, earg = elambda), ", ",
namesof("mu", link = lmu, earg = emu), ", ",
namesof("sigma", link = lsigma, earg = esigma)),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 3)
}), list( .zero = zero))),
infos = eval(substitute(function(...) {
list(M1 = 3,
Q1 = 1,
expected = TRUE,
multipleResponses = FALSE,
parameters.names = c("lambda", "mu", "sigma"),
llambda = .llambda ,
lmu = .lmu ,
lsigma = .lsigma ,
percentiles = .percentiles , # For the original fit only
true.mu = FALSE, # quantiles
zero = .zero )
}, list( .zero = zero,
.percentiles = percentiles,
.llambda = llambda, .lmu = lmu, .lsigma = lsigma ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
Is.positive.y = TRUE,
ncol.w.max = 1, ncol.y.max = 1)
predictors.names <-
c(namesof("lambda", .llambda, .elambda, short= TRUE),
namesof("mu", .lmu, .emu, short= TRUE),
namesof("sigma", .lsigma, .esigma, short= TRUE))
extra$percentiles <- .percentiles
if (!length(etastart)) {
Fit5 <- vsmooth.spline(x = x[, min(ncol(x), 2)],
y = y, w = w, df = .idf.mu )
fv.init <- c(predict(Fit5, x = x[, min(ncol(x), 2)])$y)
lambda.init <- if (is.Numeric( .ilambda )) .ilambda else 1.0
sigma.init <- if (is.null(.isigma)) {
myratio <- ((y/fv.init)^lambda.init - 1) / lambda.init
if (is.Numeric( .idf.sigma )) {
fit600 <- vsmooth.spline(x = x[, min(ncol(x), 2)],
y = myratio^2,
w = w, df = .idf.sigma)
sqrt(c(abs(predict(fit600, x = x[, min(ncol(x), 2)])$y)))
} else {
sqrt(var(myratio))
}
} else {
.isigma
}
etastart <-
cbind(theta2eta(lambda.init, .llambda , .elambda ),
theta2eta(fv.init, .lmu , .emu ),
theta2eta(sigma.init, .lsigma , .esigma ))
}
}), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.idf.mu = idf.mu,
.idf.sigma = idf.sigma,
.ilambda = ilambda, .isigma = isigma,
.percentiles = percentiles ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
pcent <- extra$percentiles
eta[, 1] <- eta2theta(eta[, 1], .llambda , .elambda )
eta[, 2] <- eta2theta(eta[, 2], .lmu , .emu )
eta[, 3] <- eta2theta(eta[, 3], .lsigma , .esigma )
qtplot.lms.bcn(percentiles = pcent, eta = eta)
}, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma
))),
last = eval(substitute(expression({
misc$links <- c(lambda = .llambda ,
mu = .lmu , sigma = .lsigma )
misc$earg <- list(lambda = .elambda ,
mu = .emu , sigma = .esigma )
misc$tol0 <- .tol0
misc$percentiles <- .percentiles # These are argument values
if (control$cdf) {
post$cdf <-
cdf.lms.bcn(y,
eta0 = matrix(c(lambda, mymu, sigma), ncol = 3,
dimnames = list(dimnames(x)[[1]], NULL)))
}
}), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.percentiles = percentiles,
.tol0 = tol0 ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
lambda <- eta2theta(eta[, 1], .llambda , .elambda )
muvec <- eta2theta(eta[, 2], .lmu , .emu )
sigma <- eta2theta(eta[, 3], .lsigma , .esigma )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- dlms.bcn(y, lam = lambda, mu = mu, sig = sigma,
tol0 = .tol0 , log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.tol0 = tol0 ))),
vfamily = c("lms.bcn", "lmscreg"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta2theta(eta[, 2], .lmu , earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
okay1 <- all(is.finite(mymu )) &&
all(is.finite(sigma )) && all(0 < sigma) &&
all(is.finite(lambda))
okay1
}, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.tol0 = tol0 ))),
deriv = eval(substitute(expression({
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta2theta(eta[, 2], .lmu , earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
zedd <- ((y / mymu)^lambda - 1) / (lambda * sigma)
z2m1 <- zedd * zedd - 1
dl.dlambda <- zedd * (zedd - log(y/mymu) / sigma) / lambda -
z2m1 * log(y/mymu)
dl.dmu <- zedd / (mymu * sigma) + z2m1 * lambda / mymu
dl.dsigma <- z2m1 / sigma
dlambda.deta <- dtheta.deta(lambda, .llambda , .elambda )
dmu.deta <- dtheta.deta(mymu, .lmu , .emu )
dsigma.deta <- dtheta.deta(sigma, .lsigma , .esigma )
c(w) * cbind(dl.dlambda * dlambda.deta,
dl.dmu * dmu.deta,
dl.dsigma * dsigma.deta)
}), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))),
weight = eval(substitute(expression({
wz <- matrix(NA_real_, n, 6)
wz[,iam(1, 1, M)] <- (7 * sigma^2 / 4) * dlambda.deta^2
wz[,iam(2, 2, M)] <- (1 + 2*(lambda*sigma)^2)/(mymu*sigma)^2 *
dmu.deta^2
wz[,iam(3, 3, M)] <- (2 / sigma^2) * dsigma.deta^2
wz[,iam(1, 2, M)] <- (-1 / (2 * mymu)) * dlambda.deta * dmu.deta
wz[,iam(1, 3, M)] <- (lambda * sigma) * dlambda.deta * dsigma.deta
wz[,iam(2, 3, M)] <- (2*lambda/(mymu * sigma)) *
dmu.deta * dsigma.deta
c(w) * wz
}),
list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))))
} # lms.bcn
lms.bcg <- function(percentiles = c(25, 50, 75),
zero = c("lambda", "sigma"),
llambda = "identitylink",
lmu = "identitylink",
lsigma = "loglink",
idf.mu = 4,
idf.sigma = 2,
ilambda = 1,
isigma = NULL) {
llambda <- as.list(substitute(llambda))
elambda <- link2list(llambda)
llambda <- attr(elambda, "function.name")
lmu <- as.list(substitute(lmu))
emu <- link2list(lmu)
lmu <- attr(emu, "function.name")
lsigma <- as.list(substitute(lsigma))
esigma <- link2list(lsigma)
lsigma <- attr(esigma, "function.name")
if (!is.Numeric(ilambda))
stop("bad input for argument 'ilambda'")
if (length(isigma) && !is.Numeric(isigma, positive = TRUE))
stop("bad input for argument 'isigma'")
new("vglmff",
blurb = c("LMS Quantile Regression ",
"(Box-Cox transformation to a Gamma distribution)\n",
"Links: ",
namesof("lambda", link = llambda, elambda), ", ",
namesof("mu", link = lmu, emu), ", ",
namesof("sigma", link = lsigma, esigma)),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 3)
}), list(.zero = zero))),
infos = eval(substitute(function(...) {
list(M1 = 3,
Q1 = 1,
expected = TRUE,
multipleResponses = FALSE,
parameters.names = c("lambda", "mu", "sigma"),
llambda = .llambda ,
lmu = .lmu ,
lsigma = .lsigma ,
percentiles = .percentiles , # For the original fit only
true.mu = FALSE, # quantiles
zero = .zero )
}, list( .zero = zero,
.percentiles = percentiles,
.llambda = llambda, .lmu = lmu, .lsigma = lsigma ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
Is.positive.y = TRUE,
ncol.w.max = 1, ncol.y.max = 1)
predictors.names <- c(
namesof("lambda", .llambda, .elambda, short = TRUE),
namesof("mu", .lmu, .emu, short = TRUE),
namesof("sigma", .lsigma, .esigma, short = TRUE))
extra$percentiles <- .percentiles
if (!length(etastart)) {
Fit5 <- vsmooth.spline(x = x[, min(ncol(x), 2)],
y = y, w = w, df = .idf.mu )
fv.init <- c(predict(Fit5, x = x[, min(ncol(x), 2)])$y)
lambda.init <- if (is.Numeric( .ilambda )) .ilambda else 1.0
sigma.init <- if (is.null( .isigma )) {
myratio <- ((y/fv.init)^lambda.init-1) / lambda.init
if (is.numeric( .idf.sigma ) &&
is.finite( .idf.sigma )) {
fit600 <- vsmooth.spline(x = x[, min(ncol(x), 2)],
y = (myratio)^2,
w = w, df = .idf.sigma )
sqrt(c(abs(predict(fit600, x = x[, min(ncol(x), 2)])$y)))
} else {
sqrt(var(myratio))
}
} else .isigma
etastart <-
cbind(theta2eta(lambda.init, .llambda , .elambda ),
theta2eta(fv.init, .lmu , .emu ),
theta2eta(sigma.init, .lsigma , .esigma ))
}
}), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.idf.mu = idf.mu,
.idf.sigma = idf.sigma,
.ilambda = ilambda, .isigma = isigma,
.percentiles = percentiles ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
pcent <- extra$percentiles
eta[, 1] <- eta2theta(eta[, 1], .llambda , .elambda )
eta[, 2] <- eta2theta(eta[, 2], .lmu , .emu )
eta[, 3] <- eta2theta(eta[, 3], .lsigma , .esigma )
qtplot.lms.bcg(percentiles = pcent, eta = eta)
}, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))),
last = eval(substitute(expression({
misc$link <- c(lambda = .llambda , mu = .lmu , sigma = .lsigma )
misc$earg <- list(lambda = .elambda , mu = .emu , sigma = .esigma )
misc$percentiles <- .percentiles # These are argument values
if (control$cdf) {
post$cdf <- cdf.lms.bcg(y, eta0 = matrix(c(lambda, mymu, sigma),
ncol = 3, dimnames = list(dimnames(x)[[1]], NULL)))
}
}), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.percentiles = percentiles ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mu <- eta2theta(eta[, 2], .lmu , earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
Gee <- (y / mu)^lambda
theta <- 1 / (sigma * lambda)^2
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * (log(abs(lambda)) + theta * (log(theta) +
log(Gee)-Gee) - lgamma(theta) - log(y))
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))),
vfamily = c("lms.bcg", "lmscreg"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta2theta(eta[, 2], .lmu , earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
okay1 <- all(is.finite(mymu )) &&
all(is.finite(sigma )) && all(0 < sigma) &&
all(is.finite(lambda))
okay1
}, list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))),
deriv = eval(substitute(expression({
lambda <- eta2theta(eta[, 1], .llambda, earg = .elambda )
mymu <- eta2theta(eta[, 2], .lmu, earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma, earg = .esigma )
Gee <- (y / mymu)^lambda
theta <- 1 / (sigma * lambda)^2
dd <- digamma(theta)
dl.dlambda <- (1 + 2 * theta * (dd + Gee -1 -log(theta) -
0.5 * (Gee + 1) * log(Gee))) / lambda
dl.dmu <- lambda * theta * (Gee-1) / mymu
dl.dsigma <- 2*theta*(dd + Gee - log(theta * Gee)-1) / sigma
dlambda.deta <- dtheta.deta(lambda, link = .llambda , .elambda )
dmu.deta <- dtheta.deta(mymu, link = .lmu , .emu )
dsigma.deta <- dtheta.deta(sigma, link = .lsigma , .esigma )
cbind(dl.dlambda * dlambda.deta,
dl.dmu * dmu.deta,
dl.dsigma * dsigma.deta) * w
}), list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))),
weight = eval(substitute(expression({
tritheta <- trigamma(theta)
wz <- matrix(0, n, 6)
if (TRUE) {
part2 <- dd + 2/theta - 2*log(theta)
wz[,iam(1, 1, M)] <- ((1 + theta*(tritheta*(1+4*theta) -
4*(1+1/theta) - log(theta)*(2/theta -
log(theta)) + dd*part2)) / lambda^2) *
dlambda.deta^2
} else {
temp <- mean( Gee*(log(Gee))^2 )
wz[,iam(1, 1, M)] <- ((4 * theta * (theta *
tritheta-1) - 1 +
theta*temp) / lambda^2) *
dlambda.deta^2
}
wz[,iam(2, 2, M)] <- dmu.deta^2 / (mymu * sigma)^2
wz[,iam(3, 3, M)] <- (4 * theta * (theta * tritheta -
1) / sigma^2) *
dsigma.deta^2
wz[,iam(1, 2, M)] <- (-theta * (dd + 1 / theta -
log(theta)) / mymu) *
dlambda.deta * dmu.deta
wz[,iam(1, 3, M)] <- 2 * theta^1.5 * (2 * theta * tritheta - 2 -
1 / theta) * dlambda.deta * dsigma.deta
c(w) * wz
}),
list( .llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma ))))
} # lms.bcg
dy.dpsi.yeojohnson <- function(psi, lambda) {
L <- max(length(psi), length(lambda))
if (length(psi) != L) psi <- rep_len(psi, L)
if (length(lambda) != L) lambda <- rep_len(lambda, L)
ifelse(psi > 0, (1 + psi * lambda)^(1/lambda - 1),
(1 - (2-lambda) * psi)^((lambda - 1) / (2-
lambda)))
}
dyj.dy.yeojohnson <- function(y, lambda) {
L <- max(length(y), length(lambda))
if (length(y) != L) y <- rep_len(y, L)
if (length(lambda) != L) lambda <- rep_len(lambda, L)
ifelse(y>0, (1 + y)^(lambda - 1), (1 - y)^(1 - lambda))
}
yeo.johnson <-
function(y, lambda, derivative = 0,
epsilon = sqrt(.Machine$double.eps),
inverse = FALSE) {
if (!is.Numeric(derivative, length.arg = 1,
integer.valued = TRUE) ||
derivative < 0)
stop("argument 'derivative' must be a non-negative integer")
ans <- y
if (!is.Numeric(epsilon, length.arg = 1, positive = TRUE))
stop("argument 'epsilon' must be a single positive number")
L <- max(length(lambda), length(y))
if (length(y) != L) y <- rep_len(y, L)
if (length(lambda) != L) lambda <- rep_len(lambda, L)
if (inverse) {
if (derivative != 0)
stop("argument 'derivative' must 0 when inverse = TRUE")
if (any(index <- y >= 0 & abs(lambda ) > epsilon))
ans[index] <- (y[index]*lambda[index] +
1)^(1/lambda[index]) - 1
if (any(index <- y >= 0 & abs(lambda ) <= epsilon))
ans[index] <- expm1(y[index])
if (any(index <- y < 0 & abs(lambda-2) > epsilon))
ans[index] <- 1 - (-(2-lambda[index]) *
y[index]+1)^(1/(2-lambda[index]))
if (any(index <- y < 0 & abs(lambda-2) <= epsilon))
ans[index] <- -expm1(-y[index])
return(ans)
}
if (derivative == 0) {
if (any(index <- y >= 0 & abs(lambda ) > epsilon))
ans[index] <- ((y[index]+1)^(lambda[index]) -
1) / lambda[index]
if (any(index <- y >= 0 & abs(lambda ) <= epsilon))
ans[index] <- log1p(y[index])
if (any(index <- y < 0 & abs(lambda-2) > epsilon))
ans[index] <- -((-y[index]+1)^(2-lambda[index]) - 1)/(2 -
lambda[index])
if (any(index <- y < 0 & abs(lambda-2) <= epsilon))
ans[index] <- -log1p(-y[index])
} else {
psi <- Recall(y = y, lambda = lambda,
derivative = derivative - 1,
epsilon = epsilon, inverse = inverse)
if (any(index <- y >= 0 & abs(lambda ) > epsilon))
ans[index] <- ( (y[index]+1)^(lambda[index]) *
(log1p(y[index]))^(derivative) - derivative *
psi[index] ) / lambda[index]
if (any(index <- y >= 0 & abs(lambda ) <= epsilon))
ans[index] <- (log1p(y[index]))^(derivative +
1) / (derivative + 1)
if (any(index <- y < 0 & abs(lambda-2) > epsilon))
ans[index] <- -( (-y[index]+1)^(2-lambda[index]) *
(-log1p(-y[index]))^(derivative) - derivative *
psi[index] ) / (2-lambda[index])
if (any(index <- y < 0 & abs(lambda-2) <= epsilon))
ans[index] <- (-log1p(-y[index]))^(derivative +
1) / (derivative + 1)
}
ans
} # yeo.johnson
dpsi.dlambda.yjn <- function(psi, lambda, mymu, sigma,
derivative = 0, smallno = 1.0e-8) {
if (!is.Numeric(derivative, length.arg = 1,
integer.valued = TRUE) ||
derivative < 0)
stop("argument 'derivative' must be a non-negative integer")
if (!is.Numeric(smallno, length.arg = 1, positive = TRUE))
stop("argument 'smallno' must be a single positive number")
L <- max(length(psi), length(lambda), length(mymu),
length(sigma))
if (length(psi) != L) psi <- rep_len(psi, L)
if (length(lambda) != L) lambda <- rep_len(lambda, L)
if (length(mymu) != L) mymu <- rep_len(mymu, L)
if (length(sigma) != L) sigma <- rep_len(sigma, L)
answer <- matrix(NA_real_, L, derivative+1)
CC <- psi >= 0
BB <- ifelse(CC, lambda, -2+lambda)
AA <- psi * BB
temp8 <- if (derivative > 0) {
answer[,1:derivative] <-
Recall(psi = psi, lambda = lambda, mymu = mymu,
sigma = sigma,
derivative = derivative-1, smallno = smallno)
answer[,derivative] * derivative
} else {
0
}
answer[, 1+derivative] <- ((AA+1) *
(log1p(AA)/BB)^derivative - temp8) / BB
pos <- (CC & abs(lambda) <= smallno) |
(!CC & abs(lambda-2) <= smallno)
if (any(pos))
answer[pos,1+derivative] =
(answer[pos, 1]^(1+derivative))/(derivative+1)
answer
}
gh.weight.yjn.11 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
((derivmat[, 2]/sigma)^2 +
sqrt(2) * z * derivmat[, 3] / sigma) / sqrt(pi)
} else {
# Long-winded way
psi <- mymu + sqrt(2) * sigma * z
(1 / sqrt(pi)) *
(dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 1)[, 2]^2 +
(psi - mymu) *
dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 2)[, 3]) / sigma^2
}
}
gh.weight.yjn.12 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
(-derivmat[, 2]) / (sqrt(pi) * sigma^2)
} else {
psi <- mymu + sqrt(2) * sigma * z
(1 / sqrt(pi)) * (-dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 1)[, 2]) / sigma^2
}
}
gh.weight.yjn.13 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
sqrt(8 / pi) * (-derivmat[, 2]) * z / sigma^2
} else {
psi <- mymu + sqrt(2) * sigma * z
(1 / sqrt(pi)) *
(-2 * dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 1)[, 2]) *
(psi - mymu) / sigma^3
}
}
glag.weight.yjn.11 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
derivmat[, 4] * (derivmat[, 2]^2 + sqrt(2) *
sigma * z * derivmat[, 3])
} else {
psi <- mymu + sqrt(2) * sigma * z
discontinuity <- -mymu / (sqrt(2) * sigma)
(1 / (2 * sqrt((z-discontinuity^2)^2 + discontinuity^2))) *
(1 / sqrt(pi)) *
(dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 1)[, 2]^2 +
(psi - mymu) *
dpsi.dlambda.yjn(psi, lambda, mymu,
sigma, derivative = 2)[, 3]) / sigma^2
}
}
glag.weight.yjn.12 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
discontinuity <- -mymu / (sqrt(2) * sigma)
if (length(derivmat)) {
derivmat[, 4] * (-derivmat[, 2])
} else {
psi <- mymu + sqrt(2) * sigma * z
(1 / (2 * sqrt((z-discontinuity^2)^2 + discontinuity^2))) *
(1 / sqrt(pi)) *
(- dpsi.dlambda.yjn(psi, lambda, mymu,
sigma, derivative = 1)[, 2]) / sigma^2
}
}
glag.weight.yjn.13 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
derivmat[, 4] * (-derivmat[, 2]) * sqrt(8) * z
} else {
psi <- mymu + sqrt(2) * sigma * z
discontinuity <- -mymu / (sqrt(2) * sigma)
(1 / (2 * sqrt((z-discontinuity^2)^2 + discontinuity^2))) *
(1 / sqrt(pi)) *
(-2 * dpsi.dlambda.yjn(psi, lambda, mymu,
sigma, derivative = 1)[, 2]) *
(psi - mymu) / sigma^3
}
}
gleg.weight.yjn.11 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
derivmat[, 4] * (derivmat[, 2]^2 + sqrt(2) *
sigma*z* derivmat[, 3])
} else {
psi <- mymu + sqrt(2) * sigma * z
(exp(-z^2) / sqrt(pi)) *
(dpsi.dlambda.yjn(psi, lambda, mymu,
sigma, derivative = 1)[, 2]^2 +
(psi - mymu) *
dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 2)[, 3]) / sigma^2
}
}
gleg.weight.yjn.12 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
derivmat[, 4] * (- derivmat[, 2])
} else {
psi <- mymu + sqrt(2) * sigma * z
(exp(-z^2) / sqrt(pi)) *
(- dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 1)[, 2]) / sigma^2
}
}
gleg.weight.yjn.13 <-
function(z, lambda, mymu, sigma, derivmat = NULL) {
if (length(derivmat)) {
derivmat[, 4] * (-derivmat[, 2]) * sqrt(8) * z
} else {
psi <- mymu + sqrt(2) * sigma * z
(exp(-z^2) / sqrt(pi)) *
(-2 * dpsi.dlambda.yjn(psi, lambda, mymu,
sigma, derivative = 1)[, 2]) *
(psi - mymu) / sigma^3
}
}
lms.yjn2.control <- function(save.weights = TRUE, ...) {
list(save.weights = save.weights)
}
lms.yjn2 <- function(percentiles = c(25, 50, 75),
zero = c("lambda", "sigma"),
llambda = "identitylink",
lmu = "identitylink",
lsigma = "loglink",
idf.mu = 4,
idf.sigma = 2,
ilambda = 1.0,
isigma = NULL,
yoffset = NULL,
nsimEIM = 250) {
llambda <- as.list(substitute(llambda))
elambda <- link2list(llambda)
llambda <- attr(elambda, "function.name")
lmu <- as.list(substitute(lmu))
emu <- link2list(lmu)
lmu <- attr(emu, "function.name")
lsigma <- as.list(substitute(lsigma))
esigma <- link2list(lsigma)
lsigma <- attr(esigma, "function.name")
if (!is.Numeric(ilambda))
stop("bad input for argument 'ilambda'")
if (length(isigma) &&
!is.Numeric(isigma, positive = TRUE))
stop("bad input for argument 'isigma'")
new("vglmff",
blurb = c("LMS Quantile Regression (Yeo-Johnson transformation",
" to normality)\n",
"Links: ",
namesof("lambda", link = llambda, elambda), ", ",
namesof("mu", link = lmu, emu ), ", ",
namesof("sigma", link = lsigma, esigma )),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 3)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 3,
Q1 = 1,
expected = TRUE,
multipleResponses = FALSE,
parameters.names = c("lambda", "mu", "sigma"),
llambda = .llambda ,
lmu = .lmu ,
lsigma = .lsigma ,
percentiles = .percentiles , # For the original fit only
true.mu = FALSE, # quantiles
zero = .zero )
}, list( .zero = zero,
.percentiles = percentiles,
.llambda = llambda, .lmu = lmu, .lsigma = lsigma ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
ncol.w.max = 1, ncol.y.max = 1)
extra$percentiles <- .percentiles
predictors.names <-
c(namesof("lambda", .llambda, earg = .elambda, short= TRUE),
namesof("mu", .lmu, earg = .emu, short= TRUE),
namesof("sigma", .lsigma, earg = .esigma, short= TRUE))
y.save <- y
yoff <- if (is.Numeric( .yoffset)) .yoffset else -median(y)
extra$yoffset <- yoff
y <- y + yoff
if (!length(etastart)) {
lambda.init <- if (is.Numeric( .ilambda )) .ilambda else 1.
y.tx <- yeo.johnson(y, lambda.init)
fv.init =
if (smoothok <-
(length(unique(sort(x[, min(ncol(x), 2)]))) > 7)) {
fit700 <- vsmooth.spline(x = x[, min(ncol(x), 2)],
y = y.tx, w = w, df = .idf.mu )
c(predict(fit700, x = x[, min(ncol(x), 2)])$y)
} else {
rep_len(weighted.mean(y, w), n)
}
sigma.init <- if (!is.Numeric(.isigma)) {
if (is.Numeric( .idf.sigma) && smoothok) {
fit710 = vsmooth.spline(x = x[, min(ncol(x), 2)],
y = (y.tx - fv.init)^2,
w = w, df = .idf.sigma)
sqrt(c(abs(predict(fit710,
x = x[, min(ncol(x), 2)])$y)))
} else {
sqrt( sum( w * (y.tx - fv.init)^2 ) / sum(w) )
}
} else
.isigma
etastart <- matrix(0, n, 3)
etastart[, 1] <- theta2eta(lambda.init, .llambda, .elambda)
etastart[, 2] <- theta2eta(fv.init, .lmu, .emu)
etastart[, 3] <- theta2eta(sigma.init, .lsigma, .esigma)
}
}), list(.llambda = llambda, .lmu = lmu, .lsigma = lsigma,
.elambda = elambda, .emu = emu, .esigma = esigma,
.ilambda = ilambda, .isigma = isigma,
.idf.mu = idf.mu,
.idf.sigma = idf.sigma,
.yoffset = yoffset,
.percentiles = percentiles ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
pcent <- extra$percentiles
eta[, 1] <- eta2theta(eta[, 1], .llambda, earg = .elambda)
eta[, 3] <- eta2theta(eta[, 3], .lsigma, earg = .esigma)
qtplot.lms.yjn(percentiles = pcent, eta = eta,
yoffset = extra$yoff)
}, list( .esigma = esigma, .elambda = elambda,
.llambda = llambda,
.lsigma = lsigma ))),
last = eval(substitute(expression({
misc$link <- c(lambda = .llambda, mu = .lmu, sigma = .lsigma)
misc$earg <- list(lambda = .elambda, mu = .emu, sigma = .esigma)
misc$nsimEIM <- .nsimEIM
misc$percentiles <- .percentiles # These are argument values
misc[["yoffset"]] <- extra$yoffset
y <- y.save # Restore back the value; to be attached to object
if (control$cdf) {
post$cdf <- cdf.lms.yjn(y + misc$yoffset,
eta0=matrix(c(lambda,mymu,sigma),
ncol=3, dimnames = list(dimnames(x)[[1]], NULL)))
}
}), list(.percentiles = percentiles,
.elambda = elambda, .emu = emu, .esigma = esigma,
.nsimEIM=nsimEIM,
.llambda = llambda, .lmu = lmu, .lsigma = lsigma ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mu <- eta2theta(eta[, 2], .lmu , earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
psi <- yeo.johnson(y, lambda)
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * (-log(sigma) - 0.5 * ((psi-mu)/sigma)^2 +
(lambda-1) * sign(y) * log1p(abs(y)))
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .elambda = elambda, .emu = emu, .esigma = esigma,
.llambda = llambda, .lmu = lmu,
.lsigma = lsigma ))),
vfamily = c("lms.yjn2", "lmscreg"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta2theta(eta[, 2], .lmu , earg = .emu )
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
okay1 <- all(is.finite(mymu )) &&
all(is.finite(sigma )) && all(0 < sigma) &&
all(is.finite(lambda))
okay1
}, list( .elambda = elambda, .emu = emu, .esigma = esigma,
.llambda = llambda, .lmu = lmu,
.lsigma = lsigma ))),
deriv = eval(substitute(expression({
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta2theta(eta[, 2], .lmu , .emu )
sigma <- eta2theta(eta[, 3], .lsigma , .esigma )
dlambda.deta <- dtheta.deta(lambda, link = .llambda, .elambda)
dmu.deta <- dtheta.deta(mymu, link = .lmu, .emu)
dsigma.deta <- dtheta.deta(sigma, link = .lsigma, .esigma)
psi <- yeo.johnson(y, lambda)
d1 <- yeo.johnson(y, lambda, deriv = 1)
AA <- (psi - mymu) / sigma
dl.dlambda <- -AA * d1 /sigma + sign(y) * log1p(abs(y))
dl.dmu <- AA / sigma
dl.dsigma <- (AA^2 -1) / sigma
dthetas.detas <- cbind(dlambda.deta, dmu.deta, dsigma.deta)
c(w) * cbind(dl.dlambda, dl.dmu, dl.dsigma) * dthetas.detas
}), list( .elambda = elambda, .emu = emu, .esigma = esigma,
.llambda = llambda, .lmu = lmu,
.lsigma = lsigma ))),
weight = eval(substitute(expression({
run.varcov <- 0
ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)
for (ii in 1:( .nsimEIM )) {
psi <- rnorm(n, mymu, sigma)
ysim <- yeo.johnson(y = psi, lam = lambda, inverse = TRUE)
d1 <- yeo.johnson(ysim, lambda, deriv = 1)
AA <- (psi - mymu) / sigma
dl.dlambda <- -AA * d1 /sigma + sign(ysim) * log1p(abs(ysim))
dl.dmu <- AA / sigma
dl.dsigma <- (AA^2 -1) / sigma
rm(ysim)
temp3 <- cbind(dl.dlambda, dl.dmu, dl.dsigma)
run.varcov <- ((ii-1) * run.varcov +
temp3[,ind1$row.index]*temp3[,ind1$col.index]) / ii
}
if (intercept.only)
run.varcov <- matrix(colMeans(run.varcov),
n, ncol(run.varcov), byrow = TRUE)
wz <- run.varcov * dthetas.detas[,ind1$row] *
dthetas.detas[,ind1$col]
dimnames(wz) <- list(rownames(wz), NULL) # Remove the colnames
c(w) * wz
}), list(.lsigma = lsigma,
.esigma = esigma, .elambda = elambda,
.nsimEIM=nsimEIM,
.llambda = llambda))))
} # lms.yjn2
lms.yjn <-
function(percentiles = c(25, 50, 75),
zero = c("lambda", "sigma"),
llambda = "identitylink",
lsigma = "loglink",
idf.mu = 4,
idf.sigma = 2,
ilambda = 1.0,
isigma = NULL,
rule = c(10, 5),
yoffset = NULL,
diagW = FALSE, iters.diagW = 6) {
llambda <- as.list(substitute(llambda))
elambda <- link2list(llambda)
llambda <- attr(elambda, "function.name")
lsigma <- as.list(substitute(lsigma))
esigma <- link2list(lsigma)
lsigma <- attr(esigma, "function.name")
rule <- rule[1] # No of pts (common) 4 all the quadrature schemes
if (rule != 5 && rule != 10)
stop("only rule=5 or 10 is supported")
new("vglmff",
blurb = c("LMS Quantile Regression ",
"(Yeo-Johnson transformation to normality)\n",
"Links: ",
namesof("lambda", link = llambda, earg = elambda), ", ",
namesof("mu", link = "identitylink", list()), ", ",
namesof("sigma", link = lsigma, earg = esigma)),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 3)
}), list(.zero = zero))),
infos = eval(substitute(function(...) {
list(M1 = 3,
Q1 = 1,
expected = TRUE,
multipleResponses = FALSE,
parameters.names = c("lambda", "mu", "sigma"),
llambda = .llambda ,
lmu = "identitylink",
lsigma = .lsigma ,
percentiles = .percentiles , # For the original fit only
true.mu = FALSE, # quantiles
zero = .zero )
}, list( .zero = zero,
.percentiles = percentiles,
.llambda = llambda, .lsigma = lsigma ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
ncol.w.max = 1, ncol.y.max = 1)
predictors.names <-
c(namesof("lambda", .llambda, earg = .elambda , short = TRUE),
"mu",
namesof("sigma", .lsigma, earg = .esigma , short = TRUE))
extra$percentiles <- .percentiles
y.save <- y
yoff <- if (is.Numeric( .yoffset )) .yoffset else -median(y)
extra$yoffset <- yoff
y <- y + yoff
if (!length(etastart)) {
lambda.init <- if (is.Numeric( .ilambda )) .ilambda else 1.0
y.tx <- yeo.johnson(y, lambda.init)
if (smoothok <-
(length(unique(sort(x[, min(ncol(x), 2)]))) > 7)) {
fit700 <- vsmooth.spline(x = x[, min(ncol(x), 2)],
y = y.tx, w = w, df = .idf.mu )
fv.init <- c(predict(fit700, x = x[, min(ncol(x), 2)])$y)
} else {
fv.init <- rep_len(weighted.mean(y, w), n)
}
sigma.init <- if (!is.Numeric( .isigma )) {
if (is.Numeric( .idf.sigma) &&
smoothok) {
fit710 = vsmooth.spline(x = x[, min(ncol(x), 2)],
y = (y.tx - fv.init)^2,
w = w, df = .idf.sigma)
sqrt(c(abs(predict(fit710,
x = x[, min(ncol(x), 2)])$y)))
} else {
sqrt( sum( w * (y.tx - fv.init)^2 ) / sum(w) )
}
} else
.isigma
etastart <-
cbind(theta2eta(lambda.init, .llambda , earg = .elambda ),
fv.init,
theta2eta(sigma.init, .lsigma , earg = .esigma ))
}
}), list( .llambda = llambda, .lsigma = lsigma,
.elambda = elambda, .esigma = esigma,
.ilambda = ilambda, .isigma = isigma,
.idf.mu = idf.mu,
.idf.sigma = idf.sigma,
.yoffset = yoffset,
.percentiles = percentiles ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
pcent <- extra$percentiles
eta[, 1] <- eta2theta(eta[, 1], .llambda , earg = .elambda )
eta[, 3] <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
qtplot.lms.yjn(percentiles = pcent,
eta = eta, yoffset = extra$yoff)
}, list(.percentiles = percentiles,
.esigma = esigma,
.elambda = elambda,
.llambda = llambda,
.lsigma = lsigma))),
last = eval(substitute(expression({
misc$link <- c(lambda = .llambda ,
mu = "identitylink",
sigma = .lsigma )
misc$earg <- list(lambda = .elambda ,
mu = list(theta = NULL),
sigma = .esigma )
misc$percentiles <- .percentiles # These are argument values
misc$true.mu <- FALSE # $fitted is not a true mu
misc[["yoffset"]] <- extra$yoff
y <- y.save # Restore back the value; to be attached to object
if (control$cdf) {
post$cdf =
cdf.lms.yjn(y + misc$yoffset,
eta0 = matrix(c(lambda,mymu,sigma),
ncol = 3,
dimnames = list(dimnames(x)[[1]], NULL)))
}
}), list( .percentiles = percentiles,
.elambda = elambda, .esigma = esigma,
.llambda = llambda, .lsigma = lsigma))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL, summation = TRUE) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mu <- eta[, 2]
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
psi <- yeo.johnson(y, lambda)
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * (-log(sigma) - 0.5 * ((psi-mu)/sigma)^2 +
(lambda-1) * sign(y) * log1p(abs(y)))
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .esigma = esigma, .elambda = elambda,
.lsigma = lsigma, .llambda = llambda))),
vfamily = c("lms.yjn", "lmscreg"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta[, 2]
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
okay1 <- all(is.finite(mymu )) &&
all(is.finite(sigma )) && all(0 < sigma) &&
all(is.finite(lambda))
okay1
}, list( .esigma = esigma, .elambda = elambda,
.lsigma = lsigma, .llambda = llambda))),
deriv = eval(substitute(expression({
lambda <- eta2theta(eta[, 1], .llambda , earg = .elambda )
mymu <- eta[, 2]
sigma <- eta2theta(eta[, 3], .lsigma , earg = .esigma )
psi <- yeo.johnson(y, lambda)
d1 <- yeo.johnson(y, lambda, deriv = 1)
AA <- (psi - mymu) / sigma
dl.dlambda <- -AA * d1 /sigma + sign(y) * log1p(abs(y))
dl.dmu <- AA / sigma
dl.dsigma <- (AA^2 -1) / sigma
dlambda.deta <- dtheta.deta(lambda, link = .llambda, .elambda )
dsigma.deta <- dtheta.deta(sigma, link = .lsigma, .esigma )
cbind(dl.dlambda * dlambda.deta,
dl.dmu,
dl.dsigma * dsigma.deta) * c(w)
}), list( .esigma = esigma, .elambda = elambda,
.lsigma = lsigma, .llambda = llambda ))),
weight = eval(substitute(expression({
wz <- matrix(0, n, 6)
wz[,iam(2, 2, M)] <- 1 / sigma^2
wz[,iam(3, 3, M)] <- 2 * wz[,iam(2, 2, M)] # 2 / sigma^2
if ( .rule == 10) {
glag.abs = c(0.13779347054,0.729454549503,
1.80834290174,3.40143369785,
5.55249614006,8.33015274676,
11.8437858379,16.2792578314,
21.996585812, 29.9206970123)
glag.wts = c(0.308441115765, 0.401119929155, 0.218068287612,
0.0620874560987, 0.00950151697517, 0.000753008388588,
2.82592334963e-5,
4.24931398502e-7, 1.83956482398e-9, 9.91182721958e-13)
} else {
glag.abs = c(0.2635603197180449, 1.4134030591060496,
3.5964257710396850,
7.0858100058570503, 12.6408008442729685)
glag.wts = c(5.217556105826727e-01, 3.986668110832433e-01,
7.594244968176882e-02,
3.611758679927785e-03, 2.336997238583738e-05)
}
if ( .rule == 10) {
sgh.abs = c(0.03873852801690856, 0.19823332465268367,
0.46520116404433082,
0.81686197962535023, 1.23454146277833154,
1.70679833036403172,
2.22994030591819214, 2.80910399394755972,
3.46387269067033854,
4.25536209637269280)
sgh.wts = c(9.855210713854302e-02, 2.086780884700499e-01,
2.520517066468666e-01,
1.986843323208932e-01,9.719839905023238e-02,
2.702440190640464e-02,
3.804646170194185e-03, 2.288859354675587e-04,
4.345336765471935e-06,
1.247734096219375e-08)
} else {
sgh.abs = c(0.1002421519682381, 0.4828139660462573,
1.0609498215257607,
1.7797294185202606, 2.6697603560875995)
sgh.wts = c(0.2484061520284881475,0.3923310666523834311,
0.2114181930760276606,
0.0332466603513424663, 0.0008248533445158026)
}
if ( .rule == 10) {
gleg.abs = c(-0.973906528517, -0.865063366689, -0.679409568299,
-0.433395394129, -0.148874338982)
gleg.abs = c(gleg.abs, rev(-gleg.abs))
gleg.wts = c(0.0666713443087, 0.149451349151, 0.219086362516,
0.26926671931, 0.295524224715)
gleg.wts = c(gleg.wts, rev(gleg.wts))
} else {
gleg.abs = c(-0.9061798459386643,-0.5384693101056820, 0,
0.5384693101056828, 0.9061798459386635)
gleg.wts = c(0.2369268850561853,0.4786286704993680,
0.5688888888888889,
0.4786286704993661, 0.2369268850561916)
}
discontinuity = -mymu/(sqrt(2)*sigma)
LL <- pmin(discontinuity, 0)
UU <- pmax(discontinuity, 0)
if (FALSE) {
AA <- (UU-LL)/2
for (kk in seq_along(gleg.wts)) {
temp1 <- AA * gleg.wts[kk]
abscissae <- (UU+LL)/2 + AA * gleg.abs[kk]
psi <- mymu + sqrt(2) * sigma * abscissae
temp9 <- dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 2)
temp9 <- cbind(temp9,
exp(-abscissae^2) / (sqrt(pi) * sigma^2))
wz[,iam(1, 1, M)] <- wz[,iam(1, 1, M)] + temp1 *
gleg.weight.yjn.11(abscissae, lambda, mymu, sigma, temp9)
wz[,iam(1, 2, M)] <- wz[,iam(1, 2, M)] + temp1 *
gleg.weight.yjn.12(abscissae, lambda, mymu, sigma, temp9)
wz[,iam(1, 3, M)] <- wz[,iam(1, 3, M)] + temp1 *
gleg.weight.yjn.13(abscissae, lambda, mymu, sigma, temp9)
}
} else {
temp9 <- .Fortran("yjngintf", as.double(LL),
as.double(UU),
as.double(gleg.abs), as.double(gleg.wts),
as.integer(n),
as.integer(length(gleg.abs)), as.double(lambda),
as.double(mymu), as.double(sigma),
answer = double(3*n),
eps = as.double(1.0e-5))$ans
dim(temp9) <- c(3, n)
wz[,iam(1, 1, M)] <- temp9[1,]
wz[,iam(1, 2, M)] <- temp9[2,]
wz[,iam(1, 3, M)] <- temp9[3,]
}
for (kk in seq_along(sgh.wts)) {
abscissae <- sign(-discontinuity) * sgh.abs[kk]
psi <- mymu + sqrt(2) * sigma * abscissae # abscissae = z
temp9 <- dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 2)
wz[,iam(1, 1, M)] <- wz[,iam(1, 1, M)] + sgh.wts[kk] *
gh.weight.yjn.11(abscissae, lambda, mymu, sigma, temp9)
wz[,iam(1, 2, M)] <- wz[,iam(1, 2, M)] + sgh.wts[kk] *
gh.weight.yjn.12(abscissae, lambda, mymu, sigma, temp9)
wz[,iam(1, 3, M)] <- wz[,iam(1, 3, M)] + sgh.wts[kk] *
gh.weight.yjn.13(abscissae, lambda, mymu, sigma, temp9)
}
temp1 <- exp(-discontinuity^2)
for (kk in seq_along(glag.wts)) {
abscissae <- sign(discontinuity) * sqrt(glag.abs[kk]) +
discontinuity^2
psi <- mymu + sqrt(2) * sigma * abscissae
temp9 <- dpsi.dlambda.yjn(psi, lambda, mymu, sigma,
derivative = 2)
temp9 <- cbind(temp9,
1 / (2 * sqrt((abscissae-discontinuity^2)^2 +
discontinuity^2) *
sqrt(pi) * sigma^2))
temp7 <- temp1 * glag.wts[kk]
wz[,iam(1, 1, M)] <- wz[,iam(1, 1, M)] + temp7 *
glag.weight.yjn.11(abscissae, lambda, mymu, sigma, temp9)
wz[,iam(1, 2, M)] <- wz[,iam(1, 2, M)] + temp7 *
glag.weight.yjn.12(abscissae, lambda, mymu, sigma, temp9)
wz[,iam(1, 3, M)] <- wz[,iam(1, 3, M)] + temp7 *
glag.weight.yjn.13(abscissae, lambda, mymu, sigma, temp9)
}
wz[, iam(1, 1, M)] <- wz[, iam(1, 1, M)] * dlambda.deta^2
wz[, iam(1, 2, M)] <- wz[, iam(1, 2, M)] * dlambda.deta
wz[, iam(1, 3, M)] <- wz[, iam(1, 3, M)] *
dsigma.deta * dlambda.deta
if ( .diagW && iter <= .iters.diagW) {
wz[,iam(1, 2, M)] <- wz[, iam(1, 3, M)] <- 0
}
wz[, iam(2, 3, M)] <- wz[, iam(2, 3, M)] * dsigma.deta
wz[, iam(3, 3, M)] <- wz[, iam(3, 3, M)] * dsigma.deta^2
c(w) * wz
}), list( .lsigma = lsigma, .llambda = llambda,
.esigma = esigma, .elambda = elambda,
.rule = rule, .diagW = diagW,
.iters.diagW = iters.diagW ))))
} # lms.yjn
lmscreg.control <-
function(cdf = TRUE, at.arg = NULL, x0 = NULL, ...) {
if (!is.logical(cdf)) {
warning("'cdf' is not logical; using TRUE instead")
cdf <- TRUE
}
list(cdf = cdf, at.arg = at.arg, x0 = x0)
}
Wr1 <- function(r, w) ifelse(r <= 0, 1, w)
Wr2 <- function(r, w) (r <= 0) * 1 + (r > 0) * w
amlnormal.deviance <-
function(mu, y, w, residuals = FALSE,
eta, extra = NULL) {
M <- length(extra$w.aml)
if (M > 1) y <- matrix(y, extra$n, extra$M)
devi <- cbind((y - mu)^2)
if (residuals) {
stop("not sure here")
wz <- VGAM.weights.function(w = w, M = extra$M, n = extra$n)
return((y - mu) * sqrt(wz) * matrix(extra$w.aml,extra$n,extra$M))
} else {
all.deviances <- numeric(M)
myresid <- matrix(y,extra$n,extra$M) - cbind(mu)
for (ii in 1:M)
all.deviances[ii] <- sum(c(w) * devi[, ii] *
Wr1(myresid[, ii],
w = extra$w.aml[ii]))
}
if (is.logical(extra$individual) && extra$individual)
all.deviances else sum(all.deviances)
} # amlnormal.deviance
amlnormal <-
function(w.aml = 1, parallel = FALSE,
lexpectile = "identitylink",
iexpectile = NULL,
imethod = 1, digw = 4) {
if (!is.Numeric(w.aml, positive = TRUE))
stop("argument 'w.aml' must be a vector of positive values")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 3)
stop("argument 'imethod' must be 1, 2 or 3")
lexpectile <- as.list(substitute(lexpectile))
eexpectile <- link2list(lexpectile)
lexpectile <- attr(eexpectile, "function.name")
if (length(iexpectile) && !is.Numeric(iexpectile))
stop("bad input for argument 'iexpectile'")
new("vglmff",
blurb = c("Asymmetric least squares quantile regression\n\n",
"Links: ",
namesof("expectile", link = lexpectile, eexpectile)),
constraints = eval(substitute(expression({
constraints <- cm.VGAM(matrix(1, M, 1), x = x,
bool = .parallel ,
constraints = constraints)
}), list( .parallel = parallel ))),
deviance = function(mu, y, w, residuals = FALSE, eta,
extra = NULL) {
amlnormal.deviance(mu = mu, y = y, w = w, residuals = residuals,
eta = eta, extra = extra)
},
initialize = eval(substitute(expression({
extra$w.aml <- .w.aml
temp5 <-
w.y.check(w = w, y = y,
ncol.w.max = 1, ncol.y.max = 1,
out.wy = TRUE,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
extra$M <- M <- length(extra$w.aml) # Recycle if necessary
extra$n <- n
extra$y.names <- y.names <-
paste0("w.aml = ", round(extra$w.aml, digits = .digw ))
predictors.names <- c(namesof(
paste0("expectile(",y.names,")"), .lexpectile ,
earg = .eexpectile, tag = FALSE))
if (!length(etastart)) {
mean.init <-
if ( .imethod == 1)
rep_len(median(y), n) else
if ( .imethod == 2 || .imethod == 3)
rep_len(weighted.mean(y, w), n) else {
junk <- lm.wfit(x = x, y = c(y), w = c(w))
junk$fitted
}
if ( .imethod == 3)
mean.init <- abs(mean.init) + 0.01
if (length( .iexpectile))
mean.init <- matrix( .iexpectile, n, M, byrow = TRUE)
etastart <-
matrix(theta2eta(mean.init, .lexpectile,
earg = .eexpectile), n, M)
}
}), list( .lexpectile = lexpectile, .eexpectile = eexpectile,
.iexpectile = iexpectile,
.imethod = imethod, .digw = digw, .w.aml = w.aml ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
ans <- eta <- as.matrix(eta)
for (ii in 1:ncol(eta))
ans[, ii] <- eta2theta(eta[, ii], .lexpectile , .eexpectile )
dimnames(ans) <- list(dimnames(eta)[[1]], extra$y.names)
ans
}, list( .lexpectile = lexpectile, .eexpectile = eexpectile ))),
last = eval(substitute(expression({
misc$link <- rep_len(.lexpectile , M)
names(misc$link) <- extra$y.names
misc$earg <- vector("list", M)
for (ilocal in 1:M)
misc$earg[[ilocal]] <- list(theta = NULL)
names(misc$earg) <- names(misc$link)
misc$parallel <- .parallel
misc$expected <- TRUE
extra$percentile <- numeric(M) # These are estimates (empirical)
misc$multipleResponses <- TRUE
for (ii in 1:M) {
use.w <- if (M > 1 && NCOL(w) == M) w[, ii] else w
extra$percentile[ii] <- 100 *
weighted.mean(myresid[, ii] <= 0, use.w)
}
names(extra$percentile) <- names(misc$link)
extra$individual <- TRUE
if (!(M > 1 && NCOL(w) == M)) {
extra$deviance <-
amlnormal.deviance(mu = mu, y = y, w = w,
residuals = FALSE, eta = eta,
extra = extra)
names(extra$deviance) <- extra$y.names
}
}), list( .lexpectile = lexpectile,
.eexpectile = eexpectile, .parallel = parallel ))),
vfamily = c("amlnormal"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
mymu <- eta2theta(eta, .lexpectile , earg = .eexpectile )
okay1 <- all(is.finite(mymu))
okay1
}, list( .lexpectile = lexpectile,
.eexpectile = eexpectile ))),
deriv = eval(substitute(expression({
mymu <- eta2theta(eta, .lexpectile , earg = .eexpectile )
dexpectile.deta <- dtheta.deta(mymu, .lexpectile ,
earg = .eexpectile )
myresid <- matrix(y,extra$n,extra$M) - cbind(mu)
wor1 <- Wr2(myresid, w = matrix(extra$w.aml, extra$n, extra$M,
byrow = TRUE))
c(w) * myresid * wor1 * dexpectile.deta
}), list( .lexpectile = lexpectile,
.eexpectile = eexpectile ))),
weight = eval(substitute(expression({
wz <- c(w) * wor1 * dexpectile.deta^2
wz
}), list( .lexpectile = lexpectile,
.eexpectile = eexpectile ))))
} # amlnormal
amlpoisson.deviance <-
function(mu, y, w, residuals = FALSE, eta,
extra = NULL) {
M <- length(extra$w.aml)
if (M > 1) y <- matrix(y,extra$n,extra$M)
nz <- y > 0
devi <- cbind(-(y - mu))
devi[nz] <- devi[nz] + y[nz] * log(y[nz]/mu[nz])
if (residuals) {
stop("not sure here")
return(sign(y - mu) * sqrt(2 * abs(devi) * w) *
matrix(extra$w,extra$n,extra$M))
} else {
all.deviances <- numeric(M)
myresid <- matrix(y,extra$n,extra$M) - cbind(mu)
for (ii in 1:M) all.deviances[ii] <- 2 *
sum(c(w) * devi[, ii] *
Wr1(myresid[, ii],
w=extra$w.aml[ii]))
}
if (is.logical(extra$individual) && extra$individual)
all.deviances else sum(all.deviances)
} # amlpoisson.deviance
amlpoisson <-
function(w.aml = 1, parallel = FALSE, imethod = 1,
digw = 4, link = "loglink") {
if (!is.Numeric(w.aml, positive = TRUE))
stop("'w.aml' must be a vector of positive values")
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
new("vglmff",
blurb = c("Poisson expectile regression by",
" asymmetric maximum likelihood estimation\n\n",
"Link: ", namesof("expectile", link, earg = earg)),
constraints = eval(substitute(expression({
constraints <- cm.VGAM(matrix(1, M, 1), x = x,
bool = .parallel ,
constraints = constraints)
}), list( .parallel = parallel ))),
deviance = function(mu, y, w, residuals = FALSE,
eta, extra = NULL) {
amlpoisson.deviance(mu = mu, y = y, w = w,
residuals = residuals,
eta = eta, extra = extra)
},
initialize = eval(substitute(expression({
extra$w.aml <- .w.aml
temp5 <-
w.y.check(w = w, y = y,
ncol.w.max = 1, ncol.y.max = 1,
out.wy = TRUE,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
extra$M <- M <- length(extra$w.aml) # Recycle if necessary
extra$n <- n
extra$y.names <- y.names <-
paste0("w.aml = ", round(extra$w.aml, digits = .digw ))
extra$individual <- FALSE
predictors.names <-
c(namesof(paste0("expectile(", y.names, ")"),
.link , earg = .earg , tag = FALSE))
if (!length(etastart)) {
mean.init <- if ( .imethod == 2)
rep_len(median(y), n) else
if ( .imethod == 1)
rep_len(weighted.mean(y, w), n) else {
junk = lm.wfit(x = x, y = c(y), w = c(w))
abs(junk$fitted)
}
etastart <-
matrix(theta2eta(mean.init, .link , earg = .earg ), n, M)
}
}), list( .link = link, .earg = earg, .imethod = imethod,
.digw = digw, .w.aml = w.aml ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
mu.ans <- eta <- as.matrix(eta)
for (ii in 1:ncol(eta))
mu.ans[, ii] <- eta2theta(eta[, ii], .link , earg = .earg )
dimnames(mu.ans) <- list(dimnames(eta)[[1]], extra$y.names)
mu.ans
}, list( .link = link, .earg = earg ))),
last = eval(substitute(expression({
misc$multipleResponses <- TRUE
misc$expected <- TRUE
misc$parallel <- .parallel
misc$link <- rep_len( .link , M)
names(misc$link) <- extra$y.names
misc$earg <- vector("list", M)
for (ilocal in 1:M)
misc$earg[[ilocal]] <- list(theta = NULL)
names(misc$earg) <- names(misc$link)
extra$percentile <- numeric(M) # These are estimates (empirical)
for (ii in 1:M)
extra$percentile[ii] <- 100 *
weighted.mean(myresid[, ii] <= 0, w)
names(extra$percentile) <- names(misc$link)
extra$individual <- TRUE
extra$deviance <- amlpoisson.deviance(mu = mu, y = y, w = w,
residuals = FALSE,
eta = eta, extra = extra)
names(extra$deviance) <- extra$y.names
}), list( .link = link, .earg = earg, .parallel = parallel ))),
linkfun = eval(substitute(function(mu, extra = NULL) {
theta2eta(mu, link = .link , earg = .earg )
}, list( .link = link, .earg = earg ))),
vfamily = c("amlpoisson"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
mymu <- eta2theta(eta, .link , earg = .earg )
okay1 <- all(is.finite(mymu)) && all(0 < mymu)
okay1
}, list( .link = link, .earg = earg ))),
deriv = eval(substitute(expression({
mymu <- eta2theta(eta, .link , earg = .earg )
dexpectile.deta <- dtheta.deta(mymu, .link , earg = .earg )
myresid <- matrix(y,extra$n,extra$M) - cbind(mu)
wor1 <- Wr2(myresid, w = matrix(extra$w.aml, extra$n, extra$M,
byrow = TRUE))
c(w) * myresid * wor1 * (dexpectile.deta / mymu)
}), list( .link = link, .earg = earg ))),
weight = eval(substitute(expression({
use.mu <- mymu
use.mu[use.mu < .Machine$double.eps^(3/4)] <-
.Machine$double.eps^(3/4)
wz <- c(w) * wor1 * use.mu * (dexpectile.deta / mymu)^2
wz
}), list( .link = link, .earg = earg ))))
} # amlpoisson
amlbinomial.deviance <-
function(mu, y, w, residuals = FALSE,
eta, extra = NULL) {
M <- length(extra$w.aml)
if (M > 1) y <- matrix(y,extra$n,extra$M)
devy <- y
nz <- y != 0
devy[nz] <- y[nz] * log(y[nz])
nz <- (1 - y) != 0
devy[nz] <- devy[nz] + (1 - y[nz]) * log1p(-y[nz])
devmu <- y * log(mu) + (1 - y) * log1p(-mu)
if (any(small <- mu * (1 - mu) < .Machine$double.eps)) {
warning("fitted values close to 0 or 1")
smu <- mu[small]
sy <- y[small]
smu <- ifelse(smu < .Machine$double.eps,
.Machine$double.eps, smu)
onemsmu <- ifelse((1 - smu) < .Machine$double.eps,
.Machine$double.eps, 1 - smu)
devmu[small] <- sy * log(smu) + (1 - sy) * log(onemsmu)
}
devi <- 2 * (devy - devmu)
if (residuals) {
stop("not sure here")
return(sign(y - mu) * sqrt(abs(devi) * w))
} else {
all.deviances <- numeric(M)
myresid <- matrix(y,extra$n,extra$M) -
matrix(mu,extra$n,extra$M)
for (ii in 1:M) all.deviances[ii] <- sum(c(w) * devi[, ii] *
Wr1(myresid[, ii], w=extra$w.aml[ii]))
}
if (is.logical(extra$individual) && extra$individual)
all.deviances else sum(all.deviances)
} # amlbinomial.deviance
amlbinomial <-
function(w.aml = 1, parallel = FALSE, digw = 4,
link = "logitlink") {
if (!is.Numeric(w.aml, positive = TRUE))
stop("'w.aml' must be a vector of positive values")
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
new("vglmff",
blurb = c("Logistic expectile regression by ",
"asymmetric maximum likelihood estimation\n\n",
"Link: ", namesof("expectile", link, earg = earg)),
constraints = eval(substitute(expression({
constraints <- cm.VGAM(matrix(1, M, 1), x = x,
bool = .parallel ,
constraints = constraints)
}), list( .parallel = parallel ))),
deviance = function(mu, y, w, residuals = FALSE,
eta, extra = NULL) {
amlbinomial.deviance(mu = mu, y = y, w = w,
residuals = residuals,
eta = eta, extra = extra)
},
initialize = eval(substitute(expression({
{
NCOL <- function (x)
if (is.array(x) && length(dim(x)) > 1 ||
is.data.frame(x)) ncol(x) else as.integer(1)
if (NCOL(y) == 1) {
if (is.factor(y)) y <- y != levels(y)[1]
nn <- rep_len(1, n)
if (!all(y >= 0 & y <= 1))
stop("response values must be in [0, 1]")
if (!length(mustart) && !length(etastart))
mustart <- (0.5 + w * y) / (1 + w)
no.successes <- w * y
if (any(abs(no.successes -
round(no.successes)) > 0.001))
stop("Number of successes must be integer-valued")
} else if (NCOL(y) == 2) {
if (any(abs(y - round(y)) > 0.001))
stop("Count data must be integer-valued")
nn <- y[, 1] + y[, 2]
y <- ifelse(nn > 0, y[, 1]/nn, 0)
w <- w * nn
if (!length(mustart) && !length(etastart))
mustart <- (0.5 + nn * y) / (1 + nn)
} else
stop("Response not of the right form")
}
extra$w.aml <- .w.aml
if (ncol(y <- cbind(y)) != 1)
stop("response must be a vector or a 1-column matrix")
extra$M <- M <- length(extra$w.aml) # Recycle if necessary
extra$n <- n
extra$y.names <- y.names <-
paste0("w.aml = ", round(extra$w.aml, digits = .digw ),
sep = "")
extra$individual <- FALSE
predictors.names <-
c(namesof(paste0("expectile(", y.names, ")"),
.link , earg = .earg , tag = FALSE))
if (!length(etastart)) {
etastart <- matrix(theta2eta(mustart, .link , .earg ),
n, M)
mustart <- NULL
}
}), list( .link = link, .earg = earg,
.digw = digw, .w.aml = w.aml ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
mu.ans <- eta <- as.matrix(eta)
for (ii in 1:ncol(eta))
mu.ans[, ii] <- eta2theta(eta[, ii], .link , earg = .earg )
dimnames(mu.ans) <- list(dimnames(eta)[[1]], extra$y.names)
mu.ans
}, list( .link = link, .earg = earg ))),
last = eval(substitute(expression({
misc$link <- rep_len(.link , M)
names(misc$link) <- extra$y.names
misc$earg <- vector("list", M)
for (ilocal in 1:M)
misc$earg[[ilocal]] <- list(theta = NULL)
names(misc$earg) <- names(misc$link)
misc$parallel <- .parallel
misc$expected <- TRUE
extra$percentile <- numeric(M) # These are estimates (empirical)
for (ii in 1:M)
extra$percentile[ii] <- 100 *
weighted.mean(myresid[, ii] <= 0, w)
names(extra$percentile) <- names(misc$link)
extra$individual <- TRUE
extra$deviance <- amlbinomial.deviance(mu = mu, y = y, w = w,
residuals = FALSE, eta = eta, extra = extra)
names(extra$deviance) <- extra$y.names
}), list( .link = link, .earg = earg, .parallel = parallel ))),
linkfun = eval(substitute(function(mu, extra = NULL) {
theta2eta(mu, link = .link , earg = .earg )
}, list( .link = link, .earg = earg ))),
vfamily = c("amlbinomial"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
mymu <- eta2theta(eta, .link , earg = .earg )
okay1 <- all(is.finite(mymu)) && all(0 < mymu & mymu < 1)
okay1
}, list( .link = link, .earg = earg ))),
deriv = eval(substitute(expression({
mymu <- eta2theta(eta, .link , earg = .earg )
use.mu <- mymu
use.mu[use.mu < .Machine$double.eps^(3/4)] <-
.Machine$double.eps^(3/4)
dexpectile.deta <- dtheta.deta(use.mu, .link , earg = .earg )
myresid <- matrix(y,extra$n,extra$M) - cbind(mu)
wor1 <- Wr2(myresid, w = matrix(extra$w.aml, extra$n, extra$M,
byrow = TRUE))
c(w) * myresid * wor1 * (dexpectile.deta / (use.mu * (1-use.mu)))
}), list( .link = link, .earg = earg ))),
weight = eval(substitute(expression({
wz <- c(w) * wor1 * (dexpectile.deta^2 / (use.mu * (1 - use.mu)))
wz
}), list( .link = link, .earg = earg))))
} # amlbinomial
amlexponential.deviance <-
function(mu, y, w, residuals = FALSE,
eta, extra = NULL) {
M <- length(extra$w.aml)
if (M > 1) y <- matrix(y,extra$n,extra$M)
devy <- cbind(-log(y) - 1)
devi <- cbind(-log(mu) - y / mu)
if (residuals) {
stop("not sure here")
return(sign(y - mu) * sqrt(2 * abs(devi) * w) *
matrix(extra$w,extra$n,extra$M))
} else {
all.deviances <- numeric(M)
myresid <- matrix(y,extra$n,extra$M) - cbind(mu)
for (ii in 1:M) all.deviances[ii] = 2 * sum(c(w) *
(devy[, ii] - devi[, ii]) *
Wr1(myresid[, ii], w=extra$w.aml[ii]))
}
if (is.logical(extra$individual) && extra$individual)
all.deviances else sum(all.deviances)
} # amlexponential.deviance
amlexponential <-
function(w.aml = 1, parallel = FALSE, imethod = 1,
digw = 4, link = "loglink") {
if (!is.Numeric(w.aml, positive = TRUE))
stop("'w.aml' must be a vector of positive values")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 3)
stop("argument 'imethod' must be 1, 2 or 3")
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
y.names <- paste0("w.aml = ", round(w.aml, digits = digw))
predictors.names <- c(namesof(
paste0("expectile(", y.names,")"), link, earg = earg))
predictors.names <- paste(predictors.names, collapse = ", ")
new("vglmff",
blurb = c("Exponential expectile regression by",
" asymmetric maximum likelihood estimation\n\n",
"Link: ", predictors.names),
constraints = eval(substitute(expression({
constraints <- cm.VGAM(matrix(1, M, 1), x = x,
bool = .parallel ,
constraints = constraints)
}), list( .parallel = parallel ))),
deviance = function(mu, y, w, residuals = FALSE,
eta, extra = NULL) {
amlexponential.deviance(mu = mu, y = y, w = w,
residuals = residuals,
eta = eta, extra = extra)
},
initialize = eval(substitute(expression({
extra$w.aml <- .w.aml
temp5 <-
w.y.check(w = w, y = y,
Is.positive.y = TRUE,
ncol.w.max = 1, ncol.y.max = 1,
out.wy = TRUE,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
extra$M <- M <- length(extra$w.aml) # Recycle if necessary
extra$n <- n
extra$y.names <- y.names <-
paste0("w.aml = ", round(extra$w.aml, digits = .digw ))
extra$individual = FALSE
predictors.names <- c(namesof(
paste0("expectile(", y.names, ")"),
.link , earg = .earg , tag = FALSE))
if (!length(etastart)) {
mean.init <- if ( .imethod == 1)
rep_len(median(y), n) else
if ( .imethod == 2)
rep_len(weighted.mean(y, w), n) else {
1 / (y + 1)
}
etastart <- matrix(theta2eta(mean.init, .link , .earg ),
n, M)
}
}), list( .link = link, .earg = earg, .imethod = imethod,
.digw = digw, .w.aml = w.aml ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
mu.ans <- eta <- as.matrix(eta)
for (ii in 1:ncol(eta))
mu.ans[, ii] <- eta2theta(eta[, ii], .link , earg = .earg )
dimnames(mu.ans) <- list(dimnames(eta)[[1]], extra$y.names)
mu.ans
}, list( .link = link, .earg = earg ))),
last = eval(substitute(expression({
misc$multipleResponses <- TRUE
misc$expected <- TRUE
misc$parallel <- .parallel
misc$link <- rep_len( .link , M)
names(misc$link) <- extra$y.names
misc$earg <- vector("list", M)
for (ilocal in 1:M)
misc$earg[[ilocal]] <- list(theta = NULL)
names(misc$earg) <- names(misc$link)
extra$percentile <- numeric(M) # These are estimates (empirical)
for (ii in 1:M)
extra$percentile[ii] <- 100 *
weighted.mean(myresid[, ii] <= 0, w)
names(extra$percentile) <- names(misc$link)
extra$individual <- TRUE
extra$deviance =
amlexponential.deviance(mu = mu, y = y, w = w,
residuals = FALSE,
eta = eta, extra = extra)
names(extra$deviance) <- extra$y.names
}), list( .link = link, .earg = earg, .parallel = parallel ))),
linkfun = eval(substitute(function(mu, extra = NULL) {
theta2eta(mu, link = .link , earg = .earg )
}, list( .link = link, .earg = earg ))),
vfamily = c("amlexponential"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
mymu <- eta2theta(eta, .link , earg = .earg )
okay1 <- all(is.finite(mymu)) && all(0 < mymu)
okay1
}, list( .link = link, .earg = earg ))),
deriv = eval(substitute(expression({
mymu <- eta2theta(eta, .link , earg = .earg )
bigy <- matrix(y,extra$n,extra$M)
dl.dmu <- (bigy - mymu) / mymu^2
dmu.deta <- dtheta.deta(mymu, .link , earg = .earg )
myresid <- bigy - cbind(mymu)
wor1 <- Wr2(myresid, w = matrix(extra$w.aml, extra$n, extra$M,
byrow = TRUE))
c(w) * wor1 * dl.dmu * dmu.deta
}), list( .link = link, .earg = earg ))),
weight = eval(substitute(expression({
ned2l.dmu2 <- 1 / mymu^2
wz <- c(w) * wor1 * ned2l.dmu2 * dmu.deta^2
wz
}), list( .link = link, .earg = earg ))))
} # amlexponential
benini1 <-
function(y0 = stop("argument 'y0' must be specified"),
lshape = "loglink",
ishape = NULL, imethod = 1, zero = NULL,
parallel = FALSE,
type.fitted = c("percentiles", "Qlink"),
percentiles = 50) {
type.fitted <- match.arg(type.fitted,
c("percentiles", "Qlink"))[1]
lshape <- as.list(substitute(lshape))
eshape <- link2list(lshape)
lshape <- attr(eshape, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
if (!is.Numeric(y0, positive = TRUE, length.arg = 1))
stop("bad input for argument 'y0'")
new("vglmff",
blurb = c("1-parameter Benini distribution\n\n",
"Link: ",
namesof("shape", lshape, earg = eshape),
"\n", "\n",
"Median: qbenini(p = 0.5, y0, shape)"),
constraints = eval(substitute(expression({
constraints <- cm.VGAM(matrix(1, M, 1), x = x,
bool = .parallel ,
constraints, apply.int = FALSE)
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 1)
}), list( .parallel = parallel,
.zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 1,
Q1 = 1,
expected = TRUE,
multipleResponses = TRUE,
parameters.names = c("shape"),
parallel = .parallel ,
percentiles = .percentiles ,
type.fitted = .type.fitted ,
lshape = .lshape ,
eshape = .eshape ,
zero = .zero )
}, list( .parallel = parallel,
.zero = zero,
.percentiles = percentiles ,
.type.fitted = type.fitted,
.eshape = eshape,
.lshape = lshape))),
initialize = eval(substitute(expression({
temp5 <-
w.y.check(w = w, y = y,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
ncoly <- ncol(y)
M1 <- 1
M <- M1 * ncoly
extra$ncoly <- ncoly
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
extra$percentiles <- .percentiles
extra$M1 <- M1
mynames1 <- paste0("shape", if (ncoly > 1) 1:ncoly else "")
predictors.names <-
namesof(mynames1, .lshape , earg = .eshape , tag = FALSE)
extra$y0 <- .y0 # Of unit length; 20181205; to make things easy.
if (any(y <= extra$y0))
stop("some values of the response are > argument 'y0' values")
if (!length(etastart)) {
probs.y <- (1:3) / 4
qofy <- quantile(rep(y, times = w), probs = probs.y)
if ( .imethod == 1) {
shape.init <- mean(-log1p(-probs.y) / (log(qofy))^2)
} else {
shape.init <- median(-log1p(-probs.y) / (log(qofy))^2)
}
shape.init <- matrix(if (length( .ishape )) .ishape else shape.init,
n, ncoly, byrow = TRUE)
etastart <- cbind(theta2eta(shape.init, .lshape , earg = .eshape ))
}
}), list( .imethod = imethod,
.ishape = ishape,
.lshape = lshape, .eshape = eshape,
.percentiles = percentiles,
.type.fitted = type.fitted,
.y0 = y0 ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
type.fitted <-
if (length(extra$type.fitted)) {
extra$type.fitted
} else {
warning("cannot find 'type.fitted'. Returning the 'median'.")
extra$percentiles <- 50 # Overwrite whatever was there
"percentiles"
}
type.fitted <- match.arg(type.fitted,
c("percentiles", "Qlink"))[1]
if (type.fitted == "Qlink") {
eta2theta(eta, link = "loglink")
} else {
shape <- eta2theta(eta, .lshape , earg = .eshape )
pcent <- extra$percentiles
perc.mat <- matrix(pcent, NROW(eta), length(pcent),
byrow = TRUE) / 100
fv <-
switch(type.fitted,
"percentiles" = qbenini(perc.mat,
y0 = extra$y0,
shape = matrix(shape, nrow(perc.mat), ncol(perc.mat))))
if (type.fitted == "percentiles")
fv <- label.cols.y(fv, colnames.y = extra$colnames.y,
NOS = NCOL(eta), percentiles = pcent,
one.on.one = FALSE)
fv
}
}, list( .lshape = lshape, .eshape = eshape ))),
last = eval(substitute(expression({
M1 <- extra$M1
misc$link <- c(rep_len( .lshape , ncoly))
names(misc$link) <- mynames1
misc$earg <- vector("list", M)
names(misc$earg) <- mynames1
for (ii in 1:ncoly) {
misc$earg[[ii]] <- .eshape
}
extra$y0 <- .y0
}), list( .lshape = lshape,
.eshape = eshape, .y0 = y0 ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL,
summation = TRUE) {
shape <- eta2theta(eta, .lshape , earg = .eshape )
y0 <- extra$y0
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * dbenini(y, y0 = y0, sh = shape, log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lshape = lshape, .eshape = eshape ))),
vfamily = c("benini1"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
shape <- eta2theta(eta, .lshape , earg = .eshape )
okay1 <- all(is.finite(shape)) && all(0 < shape)
okay1
}, list( .lshape = lshape, .eshape = eshape ))),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
extra <- object@extra
shape <- eta2theta(eta, .lshape , earg = .eshape )
y0 <- extra$y0
rbenini(nsim * length(shape), y0 = y0, shape = shape)
}, list( .lshape = lshape, .eshape = eshape ))),
deriv = eval(substitute(expression({
shape <- eta2theta(eta, .lshape , earg = .eshape )
y0 <- extra$y0
dl.dshape <- 1/shape - (log(y/y0))^2
dshape.deta <- dtheta.deta(shape, .lshape , earg = .eshape )
c(w) * dl.dshape * dshape.deta
}), list( .lshape = lshape, .eshape = eshape ))),
weight = eval(substitute(expression({
ned2l.dshape2 <- 1 / shape^2
wz <- ned2l.dshape2 * dshape.deta^2
c(w) * wz
}), list( .lshape = lshape, .eshape = eshape ))))
} # benini1
dbenini <- function(x, y0, shape, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
N <- max(length(x), length(shape), length(y0))
if (length(x) != N) x <- rep_len(x, N)
if (length(shape) != N) shape <- rep_len(shape, N)
if (length(y0) != N) y0 <- rep_len(y0, N)
logdensity <- rep_len(log(0), N)
xok <- (x > y0)
tempxok <- log(x[xok]/y0[xok])
logdensity[xok] <- log(2*shape[xok]) - shape[xok] * tempxok^2 +
log(tempxok) - log(x[xok])
logdensity[is.infinite(x)] <- log(0) # 20141209 KaiH
if (log.arg) logdensity else exp(logdensity)
}
pbenini <-
function(q, y0, shape, lower.tail = TRUE, log.p = FALSE) {
if (!is.Numeric(q))
stop("bad input for argument 'q'")
if (!is.Numeric(shape, positive = TRUE))
stop("bad input for argument 'shape'")
if (!is.Numeric(y0, positive = TRUE))
stop("bad input for argument 'y0'")
if (!is.logical(lower.tail) || length(lower.tail ) != 1)
stop("bad input for argument 'lower.tail'")
if (!is.logical(log.p) || length(log.p) != 1)
stop("bad input for argument 'log.p'")
N <- max(length(q), length(shape), length(y0))
if (length(q) != N) q <- rep_len(q, N)
if (length(shape) != N) shape <- rep_len(shape, N)
if (length(y0) != N) y0 <- rep_len(y0, N)
ans <- y0 * 0
ok <- q > y0
if (lower.tail) {
if (log.p) {
ans[ok] <- log(-expm1(-shape[ok] * (log(q[ok]/y0[ok]))^2))
ans[q <= y0 ] <- -Inf
} else {
ans[ok] <- -expm1(-shape[ok] * (log(q[ok]/y0[ok]))^2)
}
} else {
if (log.p) {
ans[ok] <- -shape[ok] * (log(q[ok]/y0[ok]))^2
ans[q <= y0] <- 0
} else {
ans[ok] <- exp(-shape[ok] * (log(q[ok]/y0[ok]))^2)
ans[q <= y0] <- 1
}
}
ans
}
qbenini <- function(p, y0, shape, lower.tail = TRUE, log.p = FALSE) {
if (!is.logical(lower.tail) || length(lower.tail ) != 1)
stop("bad input for argument 'lower.tail'")
if (!is.logical(log.p) || length(log.p) != 1)
stop("bad input for argument 'log.p'")
if (lower.tail) {
if (log.p) {
ln.p <- p
ans <- y0 * exp(sqrt(-log(-expm1(ln.p)) / shape))
} else {
ans <- y0 * exp(sqrt(-log1p(-p) / shape))
}
} else {
if (log.p) {
ln.p <- p
ans <- y0 * exp(sqrt(-ln.p / shape))
} else {
ans <- y0 * exp(sqrt(-log(p) / shape))
}
}
ans[y0 <= 0] <- NaN
ans
}
rbenini <- function(n, y0, shape) {
y0 * exp(sqrt(-log(runif(n)) / shape))
}
hyperg <- function(N = NULL, D = NULL,
lprob = "logitlink",
iprob = NULL) {
inputN <- is.Numeric(N, positive = TRUE)
inputD <- is.Numeric(D, positive = TRUE)
if (inputD && inputN)
stop("only one of 'N' and 'D' is to be inputted")
if (!inputD && !inputN)
stop("one of 'N' and 'D' needs to be inputted")
lprob <- as.list(substitute(lprob))
earg <- link2list(lprob)
lprob <- attr(earg, "function.name")
new("vglmff",
blurb = c("Hypergeometric distribution\n\n",
"Link: ",
namesof("prob", lprob, earg = earg), "\n",
"Mean: D/N\n"),
initialize = eval(substitute(expression({
NCOL <- function (x)
if (is.array(x) && length(dim(x)) > 1 ||
is.data.frame(x)) ncol(x) else as.integer(1)
if (NCOL(y) == 1) {
if (is.factor(y)) y <- y != levels(y)[1]
nn <- rep_len(1, n)
if (!all(y >= 0 & y <= 1))
stop("response values must be in [0, 1]")
mustart <- (0.5 + w * y) / (1 + w)
no.successes <- w * y
if (any(abs(no.successes - round(no.successes)) > 0.001))
stop("Number of successes must be integer-valued")
} else if (NCOL(y) == 2) {
if (any(abs(y - round(y)) > 0.001))
stop("Count data must be integer-valued")
nn <- y[, 1] + y[, 2]
y <- ifelse(nn > 0, y[, 1]/nn, 0)
w <- w * nn
mustart <- (0.5 + nn * y) / (1 + nn)
mustart[mustart >= 1] <- 0.95
} else
stop("Response not of the right form")
predictors.names <-
namesof("prob", .lprob , earg = .earg , tag = FALSE)
extra$Nvector <- .N
extra$Dvector <- .D
extra$Nunknown <- length(extra$Nvector) == 0
if (!length(etastart)) {
init.prob <- if (length( .iprob))
rep_len( .iprob, n) else
mustart
etastart <- matrix(init.prob, n, NCOL(y))
}
}), list( .lprob = lprob, .earg = earg, .N = N, .D = D,
.iprob = iprob ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
eta2theta(eta, .lprob , earg = .earg )
}, list( .lprob = lprob, .earg = earg ))),
last = eval(substitute(expression({
misc$link <- c("prob" = .lprob)
misc$earg <- list("prob" = .earg )
misc$Dvector <- .D
misc$Nvector <- .N
}), list( .N = N, .D = D, .lprob = lprob, .earg = earg ))),
linkfun = eval(substitute(function(mu, extra = NULL) {
theta2eta(mu, .lprob, earg = .earg )
}, list( .lprob = lprob, .earg = earg ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL,
summation = TRUE) {
N <- extra$Nvector
Dvec <- extra$Dvector
prob <- mu
yvec <- w * y
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
if (extra$Nunknown) {
tmp12 <- Dvec * (1-prob) / prob
(lgamma(1+tmp12) + lgamma(1+Dvec/prob-w) -
lgamma(1+tmp12-w+yvec) - lgamma(1+Dvec/prob))
} else {
(lgamma(1+N*prob) + lgamma(1+N*(1-prob)) -
lgamma(1+N*prob-yvec) -
lgamma(1+N*(1-prob) -w + yvec))
}
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lprob = lprob, .earg = earg ))),
vfamily = c("hyperg"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
prob <- eta2theta(eta, .lprob , earg = .earg )
okay1 <- all(is.finite(prob)) && all(0 < prob & prob < 1)
okay1
}, list( .lprob = lprob, .earg = earg ))),
deriv = eval(substitute(expression({
prob <- mu # equivalently, eta2theta(eta, .lprob, .earg )
dprob.deta <- dtheta.deta(prob, .lprob, earg = .earg )
Dvec <- extra$Dvector
Nvec <- extra$Nvector
yvec <- w * y
if (extra$Nunknown) {
tmp72 <- -Dvec / prob^2
tmp12 <- Dvec * (1-prob) / prob
dl.dprob <- tmp72 * (digamma(1 + tmp12) +
digamma(1 + Dvec/prob -w) -
digamma(1 + tmp12-w+yvec) - digamma(1 + Dvec/prob))
} else {
dl.dprob <- Nvec * (digamma(1+Nvec*prob) -
digamma(1+Nvec*(1-prob)) -
digamma(1+Nvec*prob-yvec) +
digamma(1+Nvec*(1-prob)-w+yvec))
}
c(w) * dl.dprob * dprob.deta
}), list( .lprob = lprob, .earg = earg ))),
weight = eval(substitute(expression({
if (extra$Nunknown) {
tmp722 <- tmp72^2
tmp13 <- 2*Dvec / prob^3
d2l.dprob2 <- tmp722 * (trigamma(1 + tmp12) +
trigamma(1 + Dvec/prob - w) -
trigamma(1 + tmp12 - w + yvec) -
trigamma(1 + Dvec/prob)) +
tmp13 * (digamma(1 + tmp12) +
digamma(1 + Dvec/prob - w) -
digamma(1 + tmp12 - w + yvec) -
digamma(1 + Dvec/prob))
} else {
d2l.dprob2 <- Nvec^2 * (trigamma(1+Nvec*prob) +
trigamma(1+Nvec*(1-prob)) -
trigamma(1+Nvec*prob-yvec) -
trigamma(1+Nvec*(1-prob)-w+yvec))
}
d2prob.deta2 <- d2theta.deta2(prob, .lprob , earg = .earg )
wz <- -(dprob.deta^2) * d2l.dprob2
wz <- c(w) * wz
wz[wz < .Machine$double.eps] <- .Machine$double.eps
wz
}), list( .lprob = lprob, .earg = earg ))))
} # hyperg
dtriangle <- function(x, theta, lower = 0, upper = 1, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
N <- max(length(x), length(theta), length(lower), length(upper))
if (length(x) != N) x <- rep_len(x, N)
if (length(theta) != N) theta <- rep_len(theta, N)
if (length(lower) != N) lower <- rep_len(lower, N)
if (length(upper) != N) upper <- rep_len(upper, N)
denom1 <- ((upper-lower)*(theta-lower))
denom2 <- ((upper-lower)*(upper-theta))
logdensity <- rep_len(log(0), N)
xok.neg <- (lower < x) & (x <= theta)
xok.pos <- (theta <= x) & (x < upper)
logdensity[xok.neg] =
log(2 * (x[xok.neg] - lower[xok.neg]) / denom1[xok.neg])
logdensity[xok.pos] =
log(2 * (upper[xok.pos] - x[xok.pos]) / denom2[xok.pos])
logdensity[lower >= upper] <- NaN
logdensity[lower > theta] <- NaN
logdensity[upper < theta] <- NaN
if (log.arg) logdensity else exp(logdensity)
}
rtriangle <- function(n, theta, lower = 0, upper = 1) {
use.n <- if ((length.n <- length(n)) > 1) length.n else
if (!is.Numeric(n, integer.valued = TRUE,
length.arg = 1, positive = TRUE))
stop("bad input for argument 'n'") else n
if (!is.Numeric(theta))
stop("bad input for argument 'theta'")
if (!is.Numeric(lower))
stop("bad input for argument 'lower'")
if (!is.Numeric(upper))
stop("bad input for argument 'upper'")
if (!all(lower < theta & theta < upper))
stop("lower < theta < upper values are required")
N <- use.n
lower <- rep_len(lower, N)
upper <- rep_len(upper, N)
theta <- rep_len(theta, N)
t1 <- sqrt(runif(n))
t2 <- sqrt(runif(n))
ifelse(runif(n) < (theta - lower) / (upper - lower),
lower + (theta - lower) * t1,
upper - (upper - theta) * t2)
}
qtriangle <- function(p, theta, lower = 0, upper = 1,
lower.tail = TRUE, log.p = FALSE) {
if (!is.logical(lower.tail) || length(lower.tail ) != 1)
stop("bad input for argument 'lower.tail'")
if (!is.logical(log.p) || length(log.p) != 1)
stop("bad input for argument 'log.p'")
N <- max(length(p), length(theta), length(lower), length(upper))
if (length(p) != N) p <- rep_len(p, N)
if (length(theta) != N) theta <- rep_len(theta, N)
if (length(lower) != N) lower <- rep_len(lower, N)
if (length(upper) != N) upper <- rep_len(upper, N)
ans <- NA_real_ * p
if (lower.tail) {
if (log.p) {
Neg <- (exp(ln.p) <= (theta - lower) / (upper - lower))
temp1 <- exp(ln.p) * (upper - lower) * (theta - lower)
Pos <- (exp(ln.p) >= (theta - lower) / (upper - lower))
pstar <- (exp(ln.p) - (theta - lower) / (upper - lower)) /
((upper - theta) / (upper - lower))
} else {
Neg <- (p <= (theta - lower) / (upper - lower))
temp1 <- p * (upper - lower) * (theta - lower)
Pos <- (p >= (theta - lower) / (upper - lower))
pstar <- (p - (theta - lower) / (upper - lower)) /
((upper - theta) / (upper - lower))
}
} else {
if (log.p) {
ln.p <- p
Neg <- (exp(ln.p) >= (upper- theta) / (upper - lower))
temp1 <- -expm1(ln.p) * (upper - lower) * (theta - lower)
Pos <- (exp(ln.p) <= (upper- theta) / (upper - lower))
pstar <- (-expm1(ln.p) - (theta - lower) / (upper - lower)) /
((upper - theta) / (upper - lower))
} else {
Neg <- (p >= (upper- theta) / (upper - lower))
temp1 <- (1 - p) * (upper - lower) * (theta - lower)
Pos <- (p <= (upper- theta) / (upper - lower))
pstar <- ((upper- theta) / (upper - lower) - p) /
((upper - theta) / (upper - lower))
}
}
ans[ Neg] <- lower[ Neg] + sqrt(temp1[ Neg])
if (any(Pos)) {
qstar <- cbind(1 - sqrt(1-pstar), 1 + sqrt(1-pstar))
qstar <- qstar[Pos,, drop = FALSE]
qstar <- ifelse(qstar[, 1] >= 0 & qstar[, 1] <= 1,
qstar[, 1],
qstar[, 2])
ans[Pos] <- theta[Pos] + qstar * (upper - theta)[Pos]
}
ans[theta < lower | theta > upper] <- NaN
ans
}
ptriangle <- function(q, theta, lower = 0, upper = 1,
lower.tail = TRUE, log.p = FALSE) {
N <- max(length(q), length(theta), length(lower), length(upper))
if (length(q) != N) q <- rep_len(q, N)
if (length(theta) != N) theta <- rep_len(theta, N)
if (length(lower) != N) lower <- rep_len(lower, N)
if (length(upper) != N) upper <- rep_len(upper, N)
if (!is.logical(lower.tail) || length(lower.tail ) != 1)
stop("bad input for argument 'lower.tail'")
if (!is.logical(log.p) || length(log.p) != 1)
stop("bad input for argument 'log.p'")
ans <- q * 0
qstar <- (q - lower)^2 / ((upper - lower) * (theta - lower))
Neg <- (lower <= q & q <= theta)
ans[Neg] <- if (lower.tail) {
if (log.p) {
(log(qstar))[Neg]
} else {
qstar[Neg]
}
} else {
if (log.p) {
(log1p(-qstar))[Neg]
} else {
1 - qstar[Neg]
}
}
Pos <- (theta <= q & q <= upper)
qstar <- (q - theta) / (upper-theta)
if (lower.tail) {
if (log.p) {
ans[Pos] <- log(((theta-lower)/(upper-lower))[Pos] +
(qstar * (2-qstar) *
(upper-theta) / (upper - lower))[Pos])
ans[q <= lower] <- -Inf
ans[q >= upper] <- 0
} else {
ans[Pos] <- ((theta-lower)/(upper-lower))[Pos] +
(qstar * (2-qstar) *
(upper-theta) / (upper - lower))[Pos]
ans[q <= lower] <- 0
ans[q >= upper] <- 1
}
} else {
if (log.p) {
ans[Pos] <- log(((upper - theta)/(upper-lower))[Pos] +
(qstar * (2-qstar) *
(upper-theta) / (upper - lower))[Pos])
ans[q <= lower] <- 0
ans[q >= upper] <- -Inf
} else {
ans[Pos] <- ((upper - theta)/(upper-lower))[Pos] +
(qstar * (2-qstar) *
(upper-theta) / (upper - lower))[Pos]
ans[q <= lower] <- 1
ans[q >= upper] <- 0
}
}
ans[theta < lower | theta > upper] <- NaN
ans
}
triangle.control <- function(stepsize = 0.33, maxit = 100, ...) {
list(stepsize = stepsize, maxit = maxit)
}
triangle <-
function(lower = 0, upper = 1,
link = extlogitlink(min = 0, max = 1),
itheta = NULL) {
if (!is.Numeric(lower))
stop("bad input for argument 'lower'")
if (!is.Numeric(upper))
stop("bad input for argument 'upper'")
if (!all(lower < upper))
stop("lower < upper values are required")
if (length(itheta) && !is.Numeric(itheta))
stop("bad input for 'itheta'")
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
if (length(earg$min) && any(earg$min != lower))
stop("argument 'lower' does not match the 'link'")
if (length(earg$max) && any(earg$max != upper))
stop("argument 'upper' does not match the 'link'")
new("vglmff",
blurb = c("Triangle distribution\n\n",
"Link: ",
namesof("theta", link, earg = earg)),
infos = eval(substitute(function(...) {
list(M1 = 1,
Q1 = 1,
parameters.names = c("theta"),
link = .link )
}, list( .link = link ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
ncol.w.max = 1,
ncol.y.max = 1)
extra$lower <- rep_len( .lower , n)
extra$upper <- rep_len( .upper , n)
if (any(y <= extra$lower | y >= extra$upper))
stop("some y values in [lower,upper] detected")
predictors.names <-
namesof("theta", .link , earg = .earg , tag = FALSE)
if (!length(etastart)) {
Theta.init <- if (length( .itheta )) .itheta else {
weighted.mean(y, w)
}
Theta.init <- rep_len(Theta.init, n)
etastart <- theta2eta(Theta.init, .link , earg = .earg )
}
}), list( .link = link, .earg = earg, .itheta=itheta,
.upper = upper, .lower = lower ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
Theta <- eta2theta(eta, .link , earg = .earg )
lower <- extra$lower
upper <- extra$upper
mu1 <- (lower + upper + Theta) / 3
mu1
}, list( .link = link, .earg = earg ))),
last = eval(substitute(expression({
misc$link <- c(theta = .link )
misc$earg <- list(theta = .earg )
misc$expected <- TRUE
}), list( .link = link, .earg = earg ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta, extra = NULL,
summation = TRUE) {
Theta <- eta2theta(eta, .link , earg = .earg )
lower <- extra$lower
upper <- extra$upper
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * dtriangle(y, theta = Theta, lower = lower,
upper = upper, log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .link = link, .earg = earg ))),
vfamily = c("triangle"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
Theta <- eta2theta(eta, .link , earg = .earg )
okay1 <- all(is.finite(Theta)) &&
all(extra$lower < Theta & Theta < extra$upper)
okay1
}, list( .link = link, .earg = earg ))),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
extra <- object@extra
Theta <- eta2theta(eta, .link , earg = .earg )
lower <- extra$lower
upper <- extra$upper
rtriangle(nsim * length(Theta),
theta = Theta, lower = lower, upper = upper)
}, list( .link = link, .earg = earg ))),
deriv = eval(substitute(expression({
Theta <- eta2theta(eta, .link , earg = .earg )
dTheta.deta <- dtheta.deta(Theta, .link , earg = .earg )
pos <- y > Theta
neg <- y < Theta
lower <- extra$lower
upper <- extra$upper
dl.dTheta <- 0 * y
dl.dTheta[neg] <- -1 / (Theta[neg]-lower[neg])
dl.dTheta[pos] <- 1 / (upper[pos]-Theta[pos])
c(w) * dl.dTheta * dTheta.deta
}), list( .link = link, .earg = earg ))),
weight = eval(substitute(expression({
var.dl.dTheta <- 1 / ((Theta - lower) * (upper - Theta))
wz <- var.dl.dTheta * dTheta.deta^2
c(w) * wz
}), list( .link = link, .earg = earg ))))
}
dpolono <-
function (x, meanlog = 0, sdlog = 1, bigx = 170, ...) {
mapply(function(x, meanlog, sdlog, ...) {
if (abs(x) > floor(x)) { # zero prob for -ve or non-integer
0
} else
if (x == Inf) { # 20141215 KaiH
0
} else
if (x > bigx) {
z <- (log(x) - meanlog) / sdlog
(1 + (z^2 + log(x) - meanlog - 1) / (2 * x * sdlog^2)) *
exp(-0.5 * z^2) / (sqrt(2 * pi) * sdlog * x)
} else
integrate( function(t) exp(t * x - exp(t) -
0.5 * ((t - meanlog) / sdlog)^2),
lower = -Inf,
upper = Inf, ...)$value / (sqrt(2 * pi) *
sdlog * exp(lgamma(x + 1.0)))
}, x, meanlog, sdlog, ...)
}
ppolono <-
function(q, meanlog = 0, sdlog = 1,
isOne = 1 - sqrt( .Machine$double.eps ), ...) {
.cumprob <- rep_len(0, length(q))
.cumprob[q == Inf] <- 1 # special case
q <- floor(q)
ii <- -1
while (any(xActive <- ((.cumprob < isOne) & (q > ii))))
.cumprob[xActive] <- .cumprob[xActive] +
dpolono(ii <- (ii+1), meanlog, sdlog, ...)
.cumprob
}
rpolono <- function(n, meanlog = 0, sdlog = 1) {
lambda <- rlnorm(n = n, meanlog = meanlog, sdlog = sdlog)
rpois(n = n, lambda = lambda)
}
fff.control <- function(save.weights = TRUE, ...) {
list(save.weights = save.weights)
}
fff <-
function(link = "loglink",
idf1 = NULL, idf2 = NULL, nsimEIM = 100, # ncp = 0,
imethod = 1, zero = NULL) {
link <- as.list(substitute(link))
earg <- link2list(link)
link <- attr(earg, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
if (!is.Numeric(nsimEIM, length.arg = 1,
integer.valued = TRUE) ||
nsimEIM <= 10)
stop("argument 'nsimEIM' should be an integer greater than 10")
ncp <- 0
if (any(ncp != 0))
warning("not sure about ncp != 0 wrt dl/dtheta")
new("vglmff",
blurb = c("F-distribution\n\n",
"Links: ",
namesof("df1", link, earg = earg), ", ",
namesof("df2", link, earg = earg),
"\n", "\n",
"Mean: df2/(df2-2) provided df2>2 and ncp = 0", "\n",
"Variance: ",
"2*df2^2*(df1+df2-2)/(df1*(df2-2)^2*(df2-4)) ",
"provided df2>4 and ncp = 0"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 2)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
multipleResponses = FALSE,
parameters.names = c("df1", "df2"),
zero = .zero )
}, list( .zero = zero ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
ncol.w.max = 1,
ncol.y.max = 1)
predictors.names <-
c(namesof("df1", .link , earg = .earg , tag = FALSE),
namesof("df2", .link , earg = .earg , tag = FALSE))
if (!length(etastart)) {
if ( .imethod == 1) {
df2.init <- b <- 2*mean(y) / (mean(y)-1)
df1.init <- 2*b^2*(b-2)/(var(y)*(b-2)^2 * (b-4) - 2*b^2)
if (df2.init < 4) df2.init <- 5
if (df1.init < 2) df1.init <- 3
} else {
df2.init <- b <- 2*median(y) / (median(y)-1)
summy <- summary(y)
var.est <- summy[5] - summy[2]
df1.init <- 2*b^2*(b-2)/(var.est*(b-2)^2 * (b-4) - 2*b^2)
}
df1.init <- if (length( .idf1 ))
rep_len( .idf1 , n) else
rep_len(df1.init, n)
df2.init <- if (length( .idf2 ))
rep_len( .idf2 , n) else
rep_len(1, n)
etastart <- cbind(theta2eta(df1.init, .link , earg = .earg ),
theta2eta(df2.init, .link , earg = .earg ))
}
}), list( .imethod = imethod, .idf1 = idf1, .earg = earg,
.idf2 = idf2, .link = link ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
df2 <- eta2theta(eta[, 2], .link , earg = .earg )
ans <- df2 * NA
ans[df2 > 2] <- df2[df2 > 2] / (df2[df2 > 2] - 2)
ans
}, list( .link = link, .earg = earg ))),
last = eval(substitute(expression({
misc$link <- c(df1 = .link , df2 = .link )
misc$earg <- list(df1 = .earg , df2 = .earg )
misc$nsimEIM <- .nsimEIM
misc$ncp <- .ncp
}), list( .link = link, .earg = earg,
.ncp = ncp,
.nsimEIM = nsimEIM ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
df1 <- eta2theta(eta[, 1], .link , earg = .earg )
df2 <- eta2theta(eta[, 2], .link , earg = .earg )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * df(x = y, df1 = df1, df2 = df2,
ncp = .ncp , log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .link = link, .earg = earg, .ncp = ncp ))),
vfamily = c("fff"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
df1 <- eta2theta(eta[, 1], .link , earg = .earg )
df2 <- eta2theta(eta[, 2], .link , earg = .earg )
okay1 <- all(is.finite(df1)) && all(0 < df1) &&
all(is.finite(df2)) && all(0 < df2)
okay1
}, list( .link = link, .earg = earg, .ncp = ncp ))),
deriv = eval(substitute(expression({
df1 <- eta2theta(eta[, 1], .link , earg = .earg )
df2 <- eta2theta(eta[, 2], .link , earg = .earg )
dl.ddf1 <- 0.5*digamma(0.5*(df1+df2)) +
0.5 + 0.5*log(df1/df2) +
0.5*log(y) - 0.5*digamma(0.5*df1) -
0.5*(df1+df2)*(y/df2) / (1 + df1*y/df2) -
0.5*log1p(df1*y/df2)
dl.ddf2 <- 0.5*digamma(0.5*(df1+df2)) - 0.5*df1/df2 -
0.5*digamma(0.5*df2) -
0.5*(df1+df2) * (-df1*y/df2^2) / (1 + df1*y/df2) -
0.5*log1p(df1*y/df2)
ddf1.deta <- dtheta.deta(df1, .link , earg = .earg )
ddf2.deta <- dtheta.deta(df2, .link , earg = .earg )
dthetas.detas <- cbind(ddf1.deta, ddf2.deta)
c(w) * dthetas.detas * cbind(dl.ddf1, dl.ddf2)
}), list( .link = link, .earg = earg ))),
weight = eval(substitute(expression({
run.varcov <- 0
ind1 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)
for (ii in 1:( .nsimEIM )) {
ysim <- rf(n = n, df1=df1, df2=df2)
dl.ddf1 <- 0.5*digamma(0.5*(df1+df2)) +
0.5 + 0.5*log(df1/df2) +
0.5*log(ysim) - 0.5*digamma(0.5*df1) -
0.5*(df1+df2)*(ysim/df2) / (1 + df1*ysim/df2) -
0.5*log1p(df1*ysim/df2)
dl.ddf2 <- 0.5*digamma(0.5*(df1+df2)) -
0.5*df1/df2 -
0.5*digamma(0.5*df2) -
0.5*(df1+df2) *
(-df1*ysim/df2^2)/(1 + df1*ysim/df2) -
0.5*log1p(df1*ysim/df2)
rm(ysim)
temp3 <- cbind(dl.ddf1, dl.ddf2)
run.varcov <- ((ii-1) * run.varcov +
temp3[,ind1$row.index]*
temp3[,ind1$col.index]) / ii
}
wz <- if (intercept.only)
matrix(colMeans(run.varcov),
n, ncol(run.varcov), byrow = TRUE) else run.varcov
wz <- c(w) * wz * dthetas.detas[, ind1$row] *
dthetas.detas[, ind1$col]
wz
}), list( .link = link, .earg = earg, .nsimEIM = nsimEIM,
.ncp = ncp ))))
} # fff
dlaplace <-
function(x, location = 0, scale = 1, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
logdensity <- (-abs(x-location)/scale) - log(2*scale)
if (log.arg) logdensity else exp(logdensity)
} # dlaplace
plaplace <-
function(q, location = 0, scale = 1,
lower.tail = TRUE, log.p =FALSE) {
zedd <- (q - location) / scale
if (!is.logical(lower.tail) || length(lower.tail ) != 1)
stop("bad input for argument 'lower.tail'")
if (!is.logical(log.p) || length(log.p) != 1)
stop("bad input for argument 'log.p'")
L <- max(length(q), length(location), length(scale))
if (length(q) != L) q <- rep_len(q, L)
if (length(location) != L) location <- rep_len(location, L)
if (length(scale) != L) scale <- rep_len(scale, L)
if (lower.tail) {
if (log.p) {
ans <- ifelse(q < location, log(0.5) + zedd,
log1p(- 0.5 * exp(-zedd)))
} else {
ans <- ifelse(q < location, 0.5 * exp(zedd),
1 - 0.5 * exp(-zedd))
}
} else {
if (log.p) {
ans <- ifelse(q < location, log1p(- 0.5 * exp(zedd)),
log(0.5) - zedd)
} else {
ans <- ifelse(q < location, 1 - 0.5 *
exp(zedd), 0.5 * exp(-zedd))
}
}
ans[scale <= 0] <- NaN
ans
} # plaplace
qlaplace <-
function(p, location = 0, scale = 1,
lower.tail = TRUE, log.p = FALSE) {
if (!is.logical(lower.tail) || length(lower.tail ) != 1)
stop("bad input for argument 'lower.tail'")
if (!is.logical(log.p) || length(log.p) != 1)
stop("bad input for argument 'log.p'")
L <- max(length(p), length(location), length(scale))
if (length(p) != L) p <- rep_len(p, L)
if (length(location) != L) location <- rep_len(location, L)
if (length(scale) != L) scale <- rep_len(scale, L)
if (lower.tail) {
if (log.p) {
ln.p <- p
ans <- location - sign(exp(ln.p)-0.5) * scale *
log(2 * ifelse(exp(ln.p) < 0.5,
exp(ln.p), -expm1(ln.p)))
} else {
ans <- location - sign(p-0.5) * scale *
log(2 * ifelse(p < 0.5, p, 1-p))
}
} else {
if (log.p) {
ln.p <- p
ans <- location - sign(0.5 - exp(ln.p)) * scale *
log(2 * ifelse(-expm1(ln.p) < 0.5,
-expm1(ln.p), exp(ln.p)))
# ans[ln.p > 0] <- NaN
} else {
ans <- location - sign(0.5 - p) * scale *
log(2 * ifelse(p > 0.5, 1 - p, p))
}
}
ans[scale <= 0] <- NaN
ans
} # qlaplace
rlaplace <-
function(n, location = 0, scale = 1) {
use.n <- if ((length.n <- length(n)) > 1) length.n else
if (!is.Numeric(n, integer.valued = TRUE,
length.arg = 1, positive = TRUE))
stop("bad input for argument 'n'") else n
if (!is.Numeric(scale, positive = TRUE))
stop("'scale' must be positive")
location <- rep_len(location, use.n)
scale <- rep_len(scale, use.n)
rrrr <- runif(use.n)
location - sign(rrrr - 0.5) * scale *
(log(2) + ifelse(rrrr < 0.5, log(rrrr), log1p(-rrrr)))
} # rlaplace
laplace <-
function(llocation = "identitylink",
lscale = "loglink",
ilocation = NULL, iscale = NULL,
imethod = 1,
zero = "scale") {
llocat <- as.list(substitute(llocation))
elocat <- link2list(llocat)
llocat <- attr(elocat, "function.name")
ilocat <- ilocation
lscale <- as.list(substitute(lscale))
escale <- link2list(lscale)
lscale <- attr(escale, "function.name")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 3)
stop("argument 'imethod' must be 1 or 2 or 3")
if (length(iscale) &&
!is.Numeric(iscale, positive = TRUE))
stop("bad input for argument 'iscale'")
new("vglmff",
blurb = c("Two-parameter Laplace distribution\n\n",
"Links: ",
namesof("location", llocat, earg = elocat), ", ",
namesof("scale", lscale, earg = escale),
"\n", "\n",
"Mean: location", "\n",
"Variance: 2*scale^2"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M,
predictors.names = predictors.names,
M1 = 2)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
multipleResponses = FALSE,
parameters.names = c("location", "scale"),
summary.pvalues = FALSE,
zero = .zero )
}, list( .zero = zero ))),
initialize = eval(substitute(expression({
w.y.check(w = w, y = y,
ncol.w.max = 1,
ncol.y.max = 1)
predictors.names <-
c(namesof("location", .llocat , earg = .elocat, tag = FALSE),
namesof("scale", .lscale , earg = .escale, tag = FALSE))
if (!length(etastart)) {
if ( .imethod == 1) {
locat.init <- median(y)
scale.init <- sqrt(var(y) / 2)
} else if ( .imethod == 2) {
locat.init <- weighted.mean(y, w)
scale.init <- sqrt(var(y) / 2)
} else {
locat.init <- median(y)
scale.init <- sqrt(sum(c(w)*abs(y-median(y ))) / (sum(w) *2))
}
locat.init <- if (length( .ilocat ))
rep_len( .ilocat , n) else
rep_len(locat.init, n)
scale.init <- if (length( .iscale ))
rep_len( .iscale , n) else
rep_len(scale.init, n)
etastart <-
cbind(theta2eta(locat.init, .llocat , earg = .elocat ),
theta2eta(scale.init, .lscale , earg = .escale ))
}
}), list( .imethod = imethod,
.elocat = elocat, .escale = escale,
.llocat = llocat, .lscale = lscale,
.ilocat = ilocat, .iscale = iscale ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
eta2theta(eta[, 1], .llocat , earg = .elocat )
}, list( .elocat = elocat, .llocat = llocat ))),
last = eval(substitute(expression({
misc$link <- c(location = .llocat , scale = .lscale )
misc$earg <- list(location = .elocat , scale = .escale )
misc$expected <- TRUE
misc$RegCondOK <- FALSE # Save this for later
}), list( .escale = escale, .lscale = lscale,
.elocat = elocat, .llocat = llocat ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
locat <- eta2theta(eta[, 1], .llocat , earg = .elocat )
Scale <- eta2theta(eta[, 2], .lscale , earg = .escale )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * dlaplace(x = y, locat = locat,
scale = Scale, log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .escale = escale, .lscale = lscale,
.elocat = elocat, .llocat = llocat ))),
vfamily = c("laplace"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
Locat <- eta2theta(eta[, 1], .llocat , earg = .elocat )
Scale <- eta2theta(eta[, 2], .lscale , earg = .escale )
okay1 <- all(is.finite(Locat)) &&
all(is.finite(Scale)) && all(0 < Scale)
okay1
}, list( .escale = escale, .lscale = lscale,
.elocat = elocat, .llocat = llocat ))),
deriv = eval(substitute(expression({
Locat <- eta2theta(eta[, 1], .llocat , earg = .elocat )
Scale <- eta2theta(eta[, 2], .lscale , earg = .escale )
zedd <- abs(y-Locat) / Scale
dl.dLocat <- sign(y - Locat) / Scale
dl.dscale <- zedd / Scale - 1 / Scale
dLocat.deta <- dtheta.deta(Locat, .llocat , earg = .elocat )
dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale )
c(w) * cbind(dl.dLocat * dLocat.deta,
dl.dscale * dscale.deta)
}), list( .escale = escale, .lscale = lscale,
.elocat = elocat, .llocat = llocat ))),
weight = eval(substitute(expression({
d2l.dLocat2 <- d2l.dscale2 <- 1 / Scale^2
wz <- matrix(0, nrow = n, ncol = M) # diagonal
wz[,iam(1, 1, M)] <- d2l.dLocat2 * dLocat.deta^2
wz[,iam(2, 2, M)] <- d2l.dscale2 * dscale.deta^2
c(w) * wz
}), list( .escale = escale, .lscale = lscale,
.elocat = elocat, .llocat = llocat ))))
} # laplace
Try the VGAM 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.