1# These functions are 2# Copyright (C) 1998-2021 T.W. Yee, University of Auckland. 3# All rights reserved. 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19dgumbelII <- function(x, scale = 1, shape, log = FALSE) { 20 21 22 if (!is.logical(log.arg <- log) || length(log) != 1) 23 stop("bad input for argument 'log'") 24 rm(log) 25 26 LLL <- max(length(x), length(shape), length(scale)) 27 if (length(x) != LLL) x <- rep_len(x, LLL) 28 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 29 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 30 31 32 ans <- x 33 index0 <- (x < 0) & is.finite(x) & !is.na(x) 34 35 ans[!index0] <- log(shape[!index0] / scale[!index0]) + 36 (shape[!index0] + 1) * log(scale[!index0] / x[!index0]) - 37 (x[!index0] / scale[!index0])^(-shape[!index0]) 38 ans[index0] <- log(0) 39 ans[x == Inf] <- log(0) 40 41 if (log.arg) { 42 } else { 43 ans <- exp(ans) 44 ans[index0] <- 0 45 ans[x == Inf] <- 0 46 } 47 ans[shape <= 0 | scale <= 0] <- NaN 48 ans 49} 50 51 52 53pgumbelII <- function(q, scale = 1, shape, 54 lower.tail = TRUE, log.p = FALSE) { 55 56 # 20150121 KaiH 57 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 58 stop("bad input for argument 'lower.tail'") 59 60 # 20150121 KaiH 61 if (!is.logical(log.p) || length(log.p) != 1) 62 stop("bad input for argument 'log.p'") 63 64 LLL <- max(length(q), length(shape), length(scale)) 65 if (length(q) != LLL) q <- rep_len(q, LLL) 66 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 67 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 68 69 # 20150121 KaiH 70 if (lower.tail) { 71 if (log.p) { 72 ans <- -(q / scale)^(-shape) 73 ans[q <= 0 ] <- -Inf 74 ans[q == Inf] <- 0 75 } else { 76 ans <- exp(-(q / scale)^(-shape)) 77 ans[q <= 0] <- 0 78 ans[q == Inf] <- 1 79 } 80 } else { 81 if (log.p) { 82 ans <- log(-expm1(-(q / scale)^(-shape))) 83 ans[q <= 0] <- 0 84 ans[q == Inf] <- -Inf 85 } else { 86 ans <- -expm1(-(q / scale)^(-shape)) 87 ans[q <= 0] <- 1 88 ans[q == Inf] <- 0 89 } 90 } 91 ans[shape <= 0 | scale <= 0] <- NaN 92 ans 93} 94 95 96 97 98 99 100 101 102 103 104qgumbelII <- function(p, scale = 1, shape, 105 lower.tail = TRUE, log.p = FALSE) { 106 107 108 109 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 110 stop("bad input for argument 'lower.tail'") 111 if (!is.logical(log.p) || length(log.p) != 1) 112 stop("bad input for argument 'log.p'") 113 114 115 116 LLL <- max(length(p), length(shape), length(scale)) 117 if (length(p) != LLL) p <- rep_len(p, LLL) 118 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 119 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 120 121 122 if (lower.tail) { 123 if (log.p) { 124 ln.p <- p 125 ans <- scale * (-ln.p)^(-1 / shape) 126 ans[ln.p > 0] <- NaN 127 } else { # Default 128 ans <- scale * (-log(p))^(-1 / shape) 129 ans[p < 0] <- NaN 130 ans[p == 0] <- 0 131 ans[p == 1] <- Inf 132 ans[p > 1] <- NaN 133 } 134 } else { 135 if (log.p) { 136 ln.p <- p 137 ans <- scale * (-log(-expm1(ln.p)))^(-1 / shape) 138 ans[ln.p > 0] <- NaN 139 } else { 140 ans <- scale * (-log1p(-p))^(-1 / shape) 141 ans[p < 0] <- NaN 142 ans[p == 0] <- Inf 143 ans[p == 1] <- 0 144 ans[p > 1] <- NaN 145 } 146 } 147 148 ans[shape <= 0 | scale <= 0] <- NaN 149 ans 150} 151 152 153rgumbelII <- function(n, scale = 1, shape) { 154 qgumbelII(runif(n), shape = shape, scale = scale) 155} 156 157 158 159 160 161 162 163 164 165 gumbelII <- 166 function(lscale = "loglink", lshape = "loglink", 167 iscale = NULL, ishape = NULL, 168 probs.y = c(0.2, 0.5, 0.8), 169 perc.out = NULL, # 50, 170 imethod = 1, zero = "shape", nowarning = FALSE) { 171 172 173 174 175 lshape <- as.list(substitute(lshape)) 176 eshape <- link2list(lshape) 177 lshape <- attr(eshape, "function.name") 178 179 lscale <- as.list(substitute(lscale)) 180 escale <- link2list(lscale) 181 lscale <- attr(escale, "function.name") 182 183 184 if (!is.Numeric(imethod, length.arg = 1, 185 integer.valued = TRUE, positive = TRUE) || 186 imethod > 2) 187 stop("argument 'imethod' must be 1 or 2") 188 if (!is.Numeric(probs.y, positive = TRUE) || 189 length(probs.y) < 2 || 190 max(probs.y) >= 1) 191 stop("bad input for argument 'probs.y'") 192 if (length(perc.out)) 193 if (!is.Numeric(perc.out, positive = TRUE) || 194 max(probs.y) >= 100) 195 stop("bad input for argument 'perc.out'") 196 197 198 if (length(ishape)) 199 if (!is.Numeric(ishape, positive = TRUE)) 200 stop("argument 'ishape' values must be positive") 201 if (length(iscale)) 202 if (!is.Numeric(iscale, positive = TRUE)) 203 stop("argument 'iscale' values must be positive") 204 205 206 new("vglmff", 207 blurb = c("Gumbel Type II distribution\n\n", 208 "Links: ", 209 namesof("scale", lscale, escale), ", ", 210 namesof("shape", lshape, eshape), "\n", 211 "Mean: scale^(1/shape) * gamma(1 - 1 / shape)\n", 212 "Variance: scale^(2/shape) * (gamma(1 - 2/shape) - ", 213 "gamma(1 + 1/shape)^2)"), 214 constraints = eval(substitute(expression({ 215 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 216 predictors.names = predictors.names, 217 M1 = 2) 218 219 }), list( .zero = zero ))), 220 221 infos = eval(substitute(function(...) { 222 list(M1 = 2, 223 Q1 = 1, 224 parameters.names = c("scale", "shape"), 225 perc.out = .perc.out , 226 zero = .zero ) 227 }, list( .zero = zero, 228 .perc.out = perc.out 229 ))), 230 231 initialize = eval(substitute(expression({ 232 233 temp5 <- 234 w.y.check(w = w, y = y, 235 Is.positive.y = TRUE, 236 ncol.w.max = Inf, 237 ncol.y.max = Inf, 238 out.wy = TRUE, 239 colsyperw = 1, 240 maximize = TRUE) 241 w <- temp5$w 242 y <- temp5$y 243 244 ncoly <- ncol(y) 245 M1 <- 2 246 extra$ncoly <- ncoly 247 extra$M1 <- M1 248 M <- M1 * ncoly 249 250 251 mynames1 <- param.names("scale", ncoly, skip1 = TRUE) 252 mynames2 <- param.names("shape", ncoly, skip1 = TRUE) 253 254 255 predictors.names <- 256 c(namesof(mynames1, .lscale , .escale , tag = FALSE), 257 namesof(mynames2, .lshape , .eshape , tag = FALSE))[ 258 interleave.VGAM(M, M1 = M1)] 259 260 261 Shape.init <- matrix(if (length( .ishape )) .ishape else 0 + NA, 262 n, ncoly, byrow = TRUE) 263 Scale.init <- matrix(if (length( .iscale )) .iscale else 0 + NA, 264 n, ncoly, byrow = TRUE) 265 266 if (!length(etastart)) { 267 if (!length( .ishape ) || 268 !length( .iscale )) { 269 for (ilocal in 1:ncoly) { 270 271 anyc <- FALSE # extra$leftcensored | extra$rightcensored 272 i11 <- if ( .imethod == 1) anyc else FALSE # can be all data 273 probs.y <- .probs.y 274 xvec <- log(-log(probs.y)) 275 fit0 <- lsfit(y = xvec, 276 x = log(quantile(y[!i11, ilocal], 277 probs = probs.y ))) 278 279 280 if (!is.Numeric(Shape.init[, ilocal])) 281 Shape.init[, ilocal] <- -fit0$coef["X"] 282 if (!is.Numeric(Scale.init[, ilocal])) 283 Scale.init[, ilocal] <- 284 exp(fit0$coef["Intercept"] / Shape.init[, ilocal]) 285 } # ilocal 286 287 etastart <- 288 cbind(theta2eta(Scale.init, .lscale , .escale ), 289 theta2eta(Shape.init, .lshape , .eshape ))[, 290 interleave.VGAM(M, M1 = M1)] 291 } 292 } 293 }), list( 294 .lscale = lscale, .lshape = lshape, 295 .escale = escale, .eshape = eshape, 296 .iscale = iscale, .ishape = ishape, 297 .probs.y = probs.y, 298 .imethod = imethod ) )), 299 linkinv = eval(substitute(function(eta, extra = NULL) { 300 Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale ) 301 Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape ) 302 Shape <- as.matrix(Shape) 303 304 if (length( .perc.out ) > 1 && ncol(Shape) > 1) 305 stop("argument 'perc.out' should be of length one since ", 306 "there are multiple responses") 307 308 if (!length( .perc.out )) { 309 return(Scale * gamma(1 - 1 / Shape)) 310 } 311 312 ans <- if (length( .perc.out ) > 1) { 313 qgumbelII(p = matrix( .perc.out / 100, length(Shape), 314 length( .perc.out ), byrow = TRUE), 315 shape = Shape, scale = Scale) 316 } else { 317 qgumbelII(p = .perc.out / 100, shape = Shape, scale = Scale) 318 } 319 colnames(ans) <- paste(as.character( .perc.out ), "%", sep = "") 320 ans 321 }, list( 322 .lscale = lscale, .lshape = lshape, 323 .escale = escale, .eshape = eshape, 324 .perc.out = perc.out ) )), 325 last = eval(substitute(expression({ 326 327 328 M1 <- extra$M1 329 misc$link <- 330 c(rep_len( .lscale , ncoly), 331 rep_len( .lshape , ncoly))[interleave.VGAM(M, M1 = M1)] 332 temp.names <- c(mynames1, mynames2)[interleave.VGAM(M, M1 = M1)] 333 names(misc$link) <- temp.names 334 335 misc$earg <- vector("list", M) 336 names(misc$earg) <- temp.names 337 for (ii in 1:ncoly) { 338 misc$earg[[M1*ii-1]] <- .escale 339 misc$earg[[M1*ii ]] <- .eshape 340 } 341 342 misc$M1 <- M1 343 misc$imethod <- .imethod 344 misc$expected <- TRUE 345 misc$multipleResponses <- TRUE 346 misc$perc.out <- .perc.out 347 misc$true.mu <- FALSE # @fitted is not a true mu 348 349 350 }), list( 351 .lscale = lscale, .lshape = lshape, 352 .escale = escale, .eshape = eshape, 353 .perc.out = perc.out, 354 .imethod = imethod ) )), 355 loglikelihood = eval(substitute( 356 function(mu, y, w, residuals = FALSE, eta, 357 extra = NULL, 358 summation = TRUE) { 359 Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale ) 360 Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape ) 361 if (residuals) { 362 stop("loglikelihood residuals not implemented yet") 363 } else { 364 ll.elts <- c(w) * dgumbelII(x = y, shape = Shape, 365 scale = Scale, log = TRUE) 366 if (summation) { 367 sum(ll.elts) 368 } else { 369 ll.elts 370 } 371 } 372 }, list( .lscale = lscale, .lshape = lshape, 373 .escale = escale, .eshape = eshape 374 ) )), 375 vfamily = c("gumbelII"), 376 validparams = eval(substitute(function(eta, y, extra = NULL) { 377 Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale ) 378 Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape ) 379 okay1 <- all(is.finite(Scale)) && all(0 < Scale) && 380 all(is.finite(Shape)) && all(0 < Shape) 381 okay1 382 }, list( .lscale = lscale, .lshape = lshape, 383 .escale = escale, .eshape = eshape) )), 384 385 386 simslot = eval(substitute( 387 function(object, nsim) { 388 389 pwts <- if (length(pwts <- object@prior.weights) > 0) 390 pwts else weights(object, type = "prior") 391 if (any(pwts != 1)) 392 warning("ignoring prior weights") 393 eta <- predict(object) 394 Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale ) 395 Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape ) 396 rgumbelII(nsim * length(Scale), shape = Shape, scale = Scale) 397 }, list( .lscale = lscale, .lshape = lshape, 398 .escale = escale, .eshape = eshape 399 ) )), 400 401 402 403 deriv = eval(substitute(expression({ 404 M1 <- 2 405 Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale ) 406 Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape ) 407 408 dl.dshape <- 1 / Shape + log(Scale / y) - 409 log(Scale / y) * (Scale / y)^Shape 410 dl.dscale <- Shape / Scale - (Shape / y) * (Scale / y)^(Shape - 1) 411 412 413 dscale.deta <- dtheta.deta(Scale, .lscale , .escale ) 414 dshape.deta <- dtheta.deta(Shape, .lshape , .eshape ) 415 416 myderiv <- c(w) * cbind(dl.dscale, dl.dshape) * 417 cbind(dscale.deta, dshape.deta) 418 myderiv[, interleave.VGAM(M, M1 = M1)] 419 }), list( .lscale = lscale, .lshape = lshape, 420 .escale = escale, .eshape = eshape 421 ) )), 422 weight = eval(substitute(expression({ 423 EulerM <- -digamma(1.0) 424 425 426 ned2l.dshape2 <- (1 + trigamma(2) + digamma(2)^2) / Shape^2 427 ned2l.dscale2 <- (Shape / Scale)^2 428 ned2l.dshapescale <- digamma(2) / Scale 429 430 431 wz <- array(c(c(w) * ned2l.dscale2 * dscale.deta^2, 432 c(w) * ned2l.dshape2 * dshape.deta^2, 433 c(w) * ned2l.dshapescale * dscale.deta * dshape.deta), 434 dim = c(n, M / M1, 3)) 435 wz <- arwz2wz(wz, M = M, M1 = M1) 436 wz 437 }), list( .lscale = lscale, .lshape = lshape )))) 438} 439 440 441 442 443 444dmbeard <- function(x, shape, scale = 1, rho, epsilon, log = FALSE) { 445 446 447 if (!is.logical(log.arg <- log) || length(log) != 1) 448 stop("bad input for argument 'log'") 449 rm(log) 450 451 LLL <- max(length(x), length(shape), length(scale), 452 length(rho), length(epsilon)) 453 if (length(x) != LLL) x <- rep_len(x, LLL) 454 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 455 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 456 if (length(rho) != LLL) rho <- rep_len(rho, LLL) 457 if (length(epsilon) != LLL) epsilon <- rep_len(epsilon, LLL) 458 459 460 index0 <- (x < 0) 461 462 ans <- log(epsilon * exp(-x * scale) + shape) + 463 (-epsilon * x - 464 ((rho * epsilon - 1) / (rho * scale)) * 465 (log1p(rho * shape) - 466 log(exp(-x * scale) + rho * shape) - scale * x)) - 467 log(exp(-x * scale) + shape * rho) 468 469 ans[index0] <- log(0) 470 ans[x == Inf] <- log(0) 471 472 if (log.arg) { 473 } else { 474 ans <- exp(ans) 475 ans[index0] <- 0 476 ans[x == Inf] <- 0 477 } 478 ans[shape <= 0 | scale <= 0 | rho <= 0 | epsilon <= 0] <- NaN 479 ans 480} 481 482 483pmbeard <- function(q, shape, scale = 1, rho, epsilon) { 484 485 LLL <- max(length(q), length(shape), length(scale), 486 length(rho), length(epsilon)) 487 if (length(q) != LLL) q <- rep_len(q, LLL) 488 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 489 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 490 if (length(rho) != LLL) rho <- rep_len(rho, LLL) 491 if (length(epsilon) != LLL) epsilon <- rep_len(epsilon, LLL) 492 493 494 ans <- -expm1(-epsilon * q - 495 ((rho * epsilon - 1) / (rho * scale)) * 496 (log1p(rho * shape) - 497 log(exp(-scale * q) + rho * shape) - scale * q)) 498 ans[(q <= 0)] <- 0 499 ans[shape <= 0 | scale <= 0 | rho <= 0 | epsilon <= 0] <- NaN 500 ans[q == Inf] <- 1 501 ans 502} 503 504 505 506 507 508 509 510dmperks <- function(x, scale = 1, shape, epsilon, log = FALSE) { 511 512 if (!is.logical(log.arg <- log) || length(log) != 1) 513 stop("bad input for argument 'log'") 514 rm(log) 515 516 LLL <- max(length(x), length(shape), length(scale), length(epsilon)) 517 if (length(x) != LLL) x <- rep_len(x, LLL) 518 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 519 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 520 if (length(epsilon) != LLL) epsilon <- rep_len(epsilon, LLL) 521 522 523 index0 <- (x < 0) 524 ans <- log(epsilon * exp(-x * scale) + shape) + 525 (-epsilon * x - 526 ((epsilon - 1) / scale) * 527 (log1p(shape) - 528 log(shape + exp(-x * scale)) -x * scale)) - 529 log(exp(-x * scale) + shape) 530 531 ans[index0] <- log(0) 532 ans[x == Inf] <- log(0) 533 if (log.arg) { 534 } else { 535 ans <- exp(ans) 536 ans[index0] <- 0 537 ans[x == Inf] <- 0 538 } 539 ans[shape <= 0 | scale <= 0 | epsilon <= 0] <- NaN 540 ans 541} 542 543 544 545pmperks <- function(q, scale = 1, shape, epsilon) { 546 547 LLL <- max(length(q), length(shape), length(scale)) 548 if (length(q) != LLL) q <- rep_len(q, LLL) 549 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 550 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 551 552 553 ans <- -expm1(-epsilon * q - 554 ((epsilon - 1) / scale) * 555 (log1p(shape) - 556 log(shape + exp(-q * scale)) - q * scale)) 557 558 ans[(q <= 0)] <- 0 559 ans[shape <= 0 | scale <= 0] <- NaN 560 ans[q == Inf] <- 1 561 ans 562} 563 564 565 566 567 568 569 570 571 572 573 574 575dbeard <- function(x, shape, scale = 1, rho, log = FALSE) { 576 577 warning("does not integrate to unity") 578 579 if (!is.logical(log.arg <- log) || length(log) != 1) 580 stop("bad input for argument 'log'") 581 rm(log) 582 583 LLL <- max(length(x), length(shape), length(scale), length(rho)) 584 if (length(x) != LLL) x <- rep_len(x, LLL) 585 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 586 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 587 if (length(rho) != LLL) rho <- rep_len(rho, LLL) 588 589 index0 <- (x < 0) 590 ans <- log(shape) - x * scale * (rho^(-1 / scale)) + 591 log(rho) + log(scale) + 592 (rho^(-1 / scale)) * log1p(shape * rho) - 593 (1 + rho^(-1 / scale)) * 594 log(shape * rho + exp(-x * scale)) 595 ans[index0] <- log(0) 596 ans[x == Inf] <- log(0) 597 598 599 if (log.arg) { 600 } else { 601 ans <- exp(ans) 602 ans[index0] <- 0 603 ans[x == Inf] <- 0 604 } 605 ans[shape <= 0 | scale <= 0 | rho <= 0] <- NaN 606 ans 607} 608 609 610 611 612 613 614dbeard <- function(x, shape, scale = 1, rho, log = FALSE) { 615 alpha <- shape 616 beta <- scale 617 618 warning("does not integrate to unity") 619 620 ret <- ifelse(x <= 0 | beta <= 0, NaN, 621 exp(alpha+beta*x)*(1+exp(alpha+rho))**(exp(-rho/beta))/ 622 (1+exp(alpha+rho+beta*x))**(1+exp(-rho/beta))) 623 ret 624} 625 626 627 628qbeard <- function(x, u = 0.5, alpha = 1, beta = 1,rho = 1) { 629 ret <- 630 ifelse(x <= 0 | u <= 0 | u >= 1 | 631 length(x) != length(u) | beta <= 0, 632 NaN, 633 (1/beta) * (log((u**(-beta*exp(rho))) * 634 (1+exp(alpha+rho+beta*x))-1)-alpha-rho)-x) 635 636 return(ret) 637} 638 639 640 641 642 643 644 645 646 647 648dperks <- function(x, scale = 1, shape, log = FALSE) { 649 650 if (!is.logical(log.arg <- log) || length(log) != 1) 651 stop("bad input for argument 'log'") 652 rm(log) 653 654 LLL <- max(length(x), length(shape), length(scale)) 655 if (length(x) != LLL) x <- rep_len(x, LLL) 656 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 657 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 658 659 index0 <- (x < 0) 660 ans <- log(shape) - x + 661 log1p(shape) / scale - 662 (1 + 1 / scale) * log(shape + exp(-x * scale)) 663 ans[index0] <- log(0) 664 ans[x == Inf] <- log(0) 665 666 if (log.arg) { 667 } else { 668 ans <- exp(ans) 669 ans[index0] <- 0 670 ans[x == Inf] <- 0 671 } 672 ans[shape <= 0 | scale <= 0] <- NaN 673 ans 674} 675 676 677 678pperks <- function(q, scale = 1, shape, 679 lower.tail = TRUE, log.p = FALSE) { 680 681 682 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 683 stop("bad input for argument 'lower.tail'") 684 685 if (!is.logical(log.p) || length(log.p) != 1) 686 stop("bad input for argument 'log.p'") 687 688 689 LLL <- max(length(q), length(shape), length(scale)) 690 if (length(q) != LLL) q <- rep_len(q, LLL) 691 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 692 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 693 694 logS <- -q + (log1p(shape) - 695 log(shape + exp(-q * scale))) / scale 696 697 698 if (lower.tail) { 699 if (log.p) { 700 ans <- log(-expm1(logS)) 701 ans[q <= 0 ] <- -Inf 702 ans[q == Inf] <- 0 703 } else { 704 ans <- -expm1(logS) 705 ans[q <= 0] <- 0 706 ans[q == Inf] <- 1 707 } 708 } else { 709 if (log.p) { 710 ans <- logS 711 ans[q <= 0] <- 0 712 ans[q == Inf] <- -Inf 713 } else { 714 ans <- exp(logS) 715 ans[q <= 0] <- 1 716 ans[q == Inf] <- 0 717 } 718 } 719 720 ans[shape <= 0 | scale <= 0] <- NaN 721 ans 722} 723 724 725 qperks <- 726 function(p, scale = 1, shape, lower.tail = TRUE, log.p = FALSE) { 727 728 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 729 stop("bad input for argument 'lower.tail'") 730 731 if (!is.logical(log.p) || length(log.p) != 1) 732 stop("bad input for argument 'log.p'") 733 734 LLL <- max(length(p), length(shape), length(scale)) 735 if (length(p) != LLL) p <- rep_len(p, LLL) 736 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 737 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 738 739 740 if (lower.tail) { 741 if (log.p) { 742 ln.p <- p 743 tmp <- scale * log(-expm1(ln.p)) 744 onemFb <- exp(tmp) 745 ans <- (log1p(shape - onemFb) - log(shape) - tmp) / scale 746 ans[ln.p > 0] <- NaN 747 } else { 748 tmp <- scale * log1p(-p) 749 onemFb <- exp(tmp) 750 ans <- (log1p(shape - onemFb) - log(shape) - tmp) / scale 751 ans[p < 0] <- NaN 752 ans[p == 0] <- 0 753 ans[p == 1] <- Inf 754 ans[p > 1] <- NaN 755 } 756 } else { 757 if (log.p) { 758 ln.p <- p 759 tmp <- scale * ln.p 760 onemFb <- exp(tmp) 761 ans <- (log1p(shape - onemFb) - log(shape) - tmp) / scale 762 ans[ln.p > 0] <- NaN 763 } else { 764 tmp <- scale * log(p) 765 onemFb <- exp(tmp) 766 ans <- (log1p(shape - onemFb) - log(shape) - tmp) / scale 767 ans[p < 0] <- NaN 768 ans[p == 0] <- Inf 769 ans[p == 1] <- 0 770 ans[p > 1] <- NaN 771 } 772 } 773 774 ans[shape <= 0 | scale <= 0] <- NaN 775 ans 776} 777 778 779 780rperks <- function(n, scale = 1, shape) { 781 qperks(runif(n), scale = scale, shape = shape) 782} 783 784 785 786 787 788perks.control <- function(save.weights = TRUE, ...) { 789 list(save.weights = save.weights) 790} 791 792 793 perks <- 794 function(lscale = "loglink", lshape = "loglink", 795 iscale = NULL, ishape = NULL, 796 gscale = exp(-5:5), gshape = exp(-5:5), 797 nsimEIM = 500, 798 oim.mean = FALSE, 799 zero = NULL, nowarning = FALSE) { 800 801 802 803 lshape <- as.list(substitute(lshape)) 804 eshape <- link2list(lshape) 805 lshape <- attr(eshape, "function.name") 806 807 lscale <- as.list(substitute(lscale)) 808 escale <- link2list(lscale) 809 lscale <- attr(escale, "function.name") 810 811 812 if (!is.Numeric(nsimEIM, length.arg = 1, 813 integer.valued = TRUE)) 814 stop("bad input for argument 'nsimEIM'") 815 if (nsimEIM <= 50) 816 warning("argument 'nsimEIM' should be an integer ", 817 "greater than 50, say") 818 819 820 if (length(ishape)) 821 if (!is.Numeric(ishape, positive = TRUE)) 822 stop("argument 'ishape' values must be positive") 823 if (length(iscale)) 824 if (!is.Numeric(iscale, positive = TRUE)) 825 stop("argument 'iscale' values must be positive") 826 827 828 829 830 if (!is.logical(oim.mean) || length(oim.mean) != 1) 831 stop("bad input for argument 'oim.mean'") 832 833 834 835 new("vglmff", 836 blurb = c("Perks' distribution\n\n", 837 "Links: ", 838 namesof("scale", lscale, escale), ", ", 839 namesof("shape", lshape, eshape), "\n", 840 "Median: qperks(p = 0.5, scale = scale, shape = shape)"), 841 842 constraints = eval(substitute(expression({ 843 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 844 predictors.names = predictors.names, 845 M1 = 2) 846 }), list( .zero = zero ))), 847 848 infos = eval(substitute(function(...) { 849 list(M1 = 2, 850 Q1 = 1, 851 nsimEIM = .nsimEIM , 852 parameters.names = c("scale", "shape"), 853 zero = .zero ) 854 }, list( .zero = zero, 855 .nsimEIM = nsimEIM ))), 856 initialize = eval(substitute(expression({ 857 858 temp5 <- 859 w.y.check(w = w, y = y, 860 Is.positive.y = TRUE, 861 ncol.w.max = Inf, 862 ncol.y.max = Inf, 863 out.wy = TRUE, 864 colsyperw = 1, 865 maximize = TRUE) 866 w <- temp5$w 867 y <- temp5$y 868 869 870 871 ncoly <- ncol(y) 872 M1 <- 2 873 extra$ncoly <- ncoly 874 extra$M1 <- M1 875 M <- M1 * ncoly 876 877 878 mynames1 <- param.names("scale", ncoly, skip1 = TRUE) 879 mynames2 <- param.names("shape", ncoly, skip1 = TRUE) 880 predictors.names <- 881 c(namesof(mynames1, .lscale , .escale , tag = FALSE), 882 namesof(mynames2, .lshape , .eshape , tag = FALSE))[ 883 interleave.VGAM(M, M1 = M1)] 884 885 886 887 if (!length(etastart)) { 888 889 matH <- matrix(if (length( .ishape )) .ishape else 0 + NA, 890 n, ncoly, byrow = TRUE) 891 matC <- matrix(if (length( .iscale )) .iscale else 0 + NA, 892 n, ncoly, byrow = TRUE) 893 894 shape.grid <- .gshape 895 scale.grid <- .gscale 896 897 for (spp. in 1:ncoly) { 898 yvec <- y[, spp.] 899 wvec <- w[, spp.] 900 901 902 perks.Loglikfun2 <- function(scaleval, shapeval, 903 y, x, w, extraargs) { 904 sum(c(w) * dperks(x = y, shape = shapeval, 905 scale = scaleval, log = TRUE)) 906 } 907 try.this <- 908 grid.search2(scale.grid, shape.grid, 909 objfun = perks.Loglikfun2, 910 y = yvec, w = wvec, 911 ret.objfun = TRUE) # Last value is the loglik 912 if (!length( .iscale )) 913 matC[, spp.] <- try.this["Value1"] 914 if (!length( .ishape )) 915 matH[, spp.] <- try.this["Value2"] 916 } # spp. 917 918 etastart <- 919 cbind(theta2eta(matC, .lscale , .escale ), 920 theta2eta(matH, .lshape , .eshape ))[, 921 interleave.VGAM(M, M1 = M1)] 922 } # End of !length(etastart) 923 }), list( .lscale = lscale, .lshape = lshape, 924 .eshape = eshape, .escale = escale, 925 .gshape = gshape, .gscale = gscale, 926 .ishape = ishape, .iscale = iscale 927 ))), 928 929 linkinv = eval(substitute(function(eta, extra = NULL) { 930 Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale ) 931 Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape ) 932 933 qperks(p = 0.5, shape = Shape, scale = Scale) 934 }, list( .lscale = lscale, .lshape = lshape, 935 .escale = escale, .eshape = eshape ))), 936 last = eval(substitute(expression({ 937 938 misc$link <- 939 c(rep_len( .lscale , ncoly), 940 rep_len( .lshape , ncoly))[interleave.VGAM(M, M1 = M1)] 941 temp.names <- c(mynames1, mynames2)[interleave.VGAM(M, M1 = M1)] 942 names(misc$link) <- temp.names 943 944 misc$earg <- vector("list", M) 945 names(misc$earg) <- temp.names 946 for (ii in 1:ncoly) { 947 misc$earg[[M1*ii-1]] <- .escale 948 misc$earg[[M1*ii ]] <- .eshape 949 } 950 951 952 misc$M1 <- M1 953 misc$expected <- TRUE 954 misc$multipleResponses <- TRUE 955 misc$nsimEIM <- .nsimEIM 956 }), list( .lscale = lscale, .lshape = lshape, 957 .escale = escale, .eshape = eshape, 958 .nsimEIM = nsimEIM ))), 959 loglikelihood = eval(substitute( 960 function(mu, y, w, residuals = FALSE, eta, 961 extra = NULL, 962 summation = TRUE) { 963 Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale ) 964 Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape ) 965 if (residuals) { 966 stop("loglikelihood residuals not implemented yet") 967 } else { 968 ll.elts <- c(w) * dperks(x = y, shape = Shape, 969 scale = Scale, log = TRUE) 970 if (summation) { 971 sum(ll.elts) 972 } else { 973 ll.elts 974 } 975 } 976 }, list( .lscale = lscale, .lshape = lshape, 977 .escale = escale, .eshape = eshape ))), 978 vfamily = c("perks"), 979 validparams = eval(substitute(function(eta, y, extra = NULL) { 980 Scale <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE], 981 .lscale , .escale ) 982 Shape <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE], 983 .lshape , .eshape ) 984 okay1 <- all(is.finite(Scale)) && all(0 < Scale) && 985 all(is.finite(Shape)) && all(0 < Shape) 986 okay1 987 }, list( .lscale = lscale, .lshape = lshape, 988 .escale = escale, .eshape = eshape) )), 989 990 991 992 993 994 simslot = eval(substitute( 995 function(object, nsim) { 996 997 pwts <- if (length(pwts <- object@prior.weights) > 0) 998 pwts else weights(object, type = "prior") 999 if (any(pwts != 1)) 1000 warning("ignoring prior weights") 1001 eta <- predict(object) 1002 Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale ) 1003 Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape ) 1004 rperks(nsim * length(Scale), shape = Shape, scale = Scale) 1005 }, list( .lscale = lscale, .lshape = lshape, 1006 .escale = escale, .eshape = eshape ))), 1007 1008 1009 1010 1011 1012 1013 1014 deriv = eval(substitute(expression({ 1015 M1 <- 2 1016 scale <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE], 1017 .lscale , .escale ) 1018 shape <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE], 1019 .lshape , .eshape ) 1020 1021 1022 temp2 <- exp(y * scale) 1023 temp3 <- 1 + shape * temp2 1024 dl.dshape <- 1 / shape + 1 / (scale * (1 + shape)) - 1025 (1 + 1 / scale) * temp2 / temp3 1026 dl.dscale <- y - log1p(shape) / scale^2 + 1027 log1p(shape * temp2) / scale^2 - 1028 (1 + 1 / scale) * shape * y * temp2 / temp3 1029 1030 dshape.deta <- dtheta.deta(shape, .lshape , .eshape ) 1031 dscale.deta <- dtheta.deta(scale, .lscale , .escale ) 1032 1033 dthetas.detas <- cbind(dscale.deta, dshape.deta) 1034 myderiv <- c(w) * cbind(dl.dscale, dl.dshape) * dthetas.detas 1035 myderiv[, interleave.VGAM(M, M1 = M1)] 1036 }), list( .lscale = lscale, .lshape = lshape, 1037 .escale = escale, .eshape = eshape ))), 1038 1039 1040 weight = eval(substitute(expression({ 1041 1042 NOS <- M / M1 1043 dThetas.detas <- dthetas.detas[, interleave.VGAM(M, M1 = M1)] 1044 1045 wz <- matrix(0.0, n, M + M - 1) # wz is 'tridiagonal' 1046 1047 ind1 <- iam(NA, NA, M = M1, both = TRUE, diag = TRUE) 1048 1049 1050 for (spp. in 1:NOS) { 1051 run.varcov <- 0 1052 Scale <- scale[, spp.] 1053 Shape <- shape[, spp.] 1054 1055 1056 1057 1058 if (FALSE && intercept.only && .oim.mean ) { 1059 1060 stop("this is wrong") 1061 temp8 <- (1 + Shape * exp(Scale * y[, spp.]))^2 1062 nd2l.dadb <- 2 * y[, spp.] * exp(Scale * y[, spp.]) / temp8 1063 1064 nd2l.dada <- 1 / Shape^2 + 1 / (1 + Shape)^2 - 1065 2 * exp(2 * Scale * y[, spp.]) / temp8 1066 1067 nd2l.dbdb <- 2 * Shape * y[, spp.]^2 * 1068 exp(Scale * y[, spp.]) / temp8 1069 1070 ave.oim11 <- weighted.mean(nd2l.dada, w[, spp.]) 1071 ave.oim12 <- weighted.mean(nd2l.dadb, w[, spp.]) 1072 ave.oim22 <- weighted.mean(nd2l.dbdb, w[, spp.]) 1073 run.varcov <- cbind(ave.oim11, ave.oim22, ave.oim12) 1074 } else { 1075 1076 for (ii in 1:( .nsimEIM )) { 1077 ysim <- rperks(n = n, shape = Shape, scale = Scale) 1078if (ii < 3) { 1079} 1080 1081 temp2 <- exp(ysim * Scale) 1082 temp3 <- 1 + Shape * temp2 1083 dl.dshape <- 1 / Shape + 1 / (Scale * (1 + Shape)) - 1084 (1 + 1 / Scale) * temp2 / temp3 1085 dl.dscale <- ysim - log1p(Shape) / Scale^2 + 1086 log1p(Shape * temp2) / Scale^2 - 1087 (1 + 1 / Scale) * Shape * ysim * temp2 / temp3 1088 1089 1090 temp7 <- cbind(dl.dscale, dl.dshape) 1091if (ii < 3) { 1092} 1093 run.varcov <- run.varcov + 1094 temp7[, ind1$row.index] * 1095 temp7[, ind1$col.index] 1096 } 1097 run.varcov <- cbind(run.varcov / .nsimEIM ) 1098 1099 } 1100 1101 1102 1103 wz1 <- if (intercept.only) 1104 matrix(colMeans(run.varcov), 1105 nrow = n, ncol = ncol(run.varcov), byrow = TRUE) else 1106 run.varcov 1107 1108 wz1 <- wz1 * dThetas.detas[, M1 * (spp. - 1) + ind1$row] * 1109 dThetas.detas[, M1 * (spp. - 1) + ind1$col] 1110 1111 1112 for (jay in 1:M1) 1113 for (kay in jay:M1) { 1114 cptr <- iam((spp. - 1) * M1 + jay, 1115 (spp. - 1) * M1 + kay, 1116 M = M) 1117 wz[, cptr] <- wz1[, iam(jay, kay, M = M1)] 1118 } 1119 } # End of for (spp.) loop 1120 1121 1122 1123 w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = M / M1) 1124 }), list( .lscale = lscale, 1125 .escale = escale, 1126 .nsimEIM = nsimEIM, .oim.mean = oim.mean )))) 1127} # perks() 1128 1129 1130 1131 1132 1133 1134 1135 1136dmakeham <- function(x, scale = 1, shape, epsilon = 0, log = FALSE) { 1137 1138 if (!is.logical(log.arg <- log) || length(log) != 1) 1139 stop("bad input for argument 'log'") 1140 rm(log) 1141 1142 LLL <- max(length(x), length(shape), length(scale), length(epsilon)) 1143 if (length(x) != LLL) x <- rep_len(x, LLL) 1144 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 1145 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 1146 if (length(epsilon) != LLL) epsilon <- rep_len(epsilon, LLL) 1147 1148 index0 <- (x < 0) 1149 ans <- log(epsilon * exp(-x * scale) + shape) + 1150 x * (scale - epsilon) - 1151 (shape / scale) * expm1(x * scale) 1152 ans[index0] <- log(0) 1153 ans[x == Inf] <- log(0) 1154 if (log.arg) { 1155 } else { 1156 ans <- exp(ans) 1157 ans[index0] <- 0 1158 ans[x == Inf] <- 0 1159 } 1160 ans[shape <= 0 | scale <= 0 | epsilon < 0] <- NaN 1161 ans 1162} 1163 1164 1165 1166pmakeham <- function(q, scale = 1, shape, epsilon = 0, 1167 lower.tail = TRUE, log.p = FALSE) { 1168 1169 1170 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 1171 stop("bad input for argument 'lower.tail'") 1172 1173 if (!is.logical(log.p) || length(log.p) != 1) 1174 stop("bad input for argument 'log.p'") 1175 1176 LLL <- max(length(q), length(shape), length(scale), length(epsilon)) 1177 if (length(q) != LLL) q <- rep_len(q, LLL) 1178 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 1179 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 1180 if (length(epsilon) != LLL) epsilon <- rep_len(epsilon, LLL) 1181 1182 if (lower.tail) { 1183 if (log.p) { 1184 ans <- log(-expm1(-q * epsilon - (shape / scale) * 1185 expm1(scale * q))) 1186 ans[q <= 0 ] <- -Inf 1187 ans[q == Inf] <- 0 1188 } else { 1189 ans <- -expm1(-q * epsilon - (shape / scale) * expm1(scale * q)) 1190 ans[q <= 0] <- 0 1191 ans[q == Inf] <- 1 1192 } 1193 } else { 1194 if (log.p) { 1195 ans <- -q * epsilon - (shape / scale) * expm1(scale * q) 1196 ans[q <= 0] <- 0 1197 ans[q == Inf] <- -Inf 1198 } else { 1199 ans <- exp(-q * epsilon - (shape / scale) * expm1(scale * q)) 1200 ans[q <= 0] <- 1 1201 ans[q == Inf] <- 0 1202 } 1203 } 1204 1205 ans[shape <= 0 | scale <= 0 | epsilon < 0] <- NaN 1206 ans 1207} 1208 1209 1210 1211qmakeham <- function(p, scale = 1, shape, epsilon = 0, 1212 lower.tail = TRUE, log.p = FALSE) { 1213 1214 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 1215 stop("bad input for argument 'lower.tail'") 1216 1217 if (!is.logical(log.p) || length(log.p) != 1) 1218 stop("bad input for argument 'log.p'") 1219 1220 LLL <- max(length(p), length(shape), length(scale), length(epsilon)) 1221 if (length(p) != LLL) p <- rep_len(p, LLL) 1222 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 1223 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 1224 if (length(epsilon) != LLL) epsilon <- rep_len(epsilon, LLL) 1225 1226 1227 if (lower.tail) { 1228 if (log.p) { 1229 ln.p <- p 1230 ans <- shape / (scale * epsilon) - log(-expm1(ln.p)) / epsilon - 1231 lambertW((shape / epsilon) * exp(shape / epsilon) * 1232 exp(log(-expm1(ln.p)) * (-scale / epsilon))) / scale 1233 ans[ln.p == 0] <- Inf 1234 ans[ln.p > 0] <- NaN 1235 } else { 1236 ans <- shape / (scale * epsilon) - log1p(-p) / epsilon - 1237 lambertW((shape / epsilon) * exp(shape / epsilon) * 1238 exp( (-scale / epsilon) * log1p(-p) )) / scale 1239 ans[p < 0] <- NaN 1240 ans[p == 0] <- 0 1241 ans[p == 1] <- Inf 1242 ans[p > 1] <- NaN 1243 } 1244 } else { 1245 if (log.p) { 1246 ln.p <- p 1247 ans <- shape / (scale * epsilon) - ln.p / epsilon - 1248 lambertW((shape / epsilon) * exp(shape / epsilon) * 1249 exp(ln.p * (-scale / epsilon))) / scale 1250 ans[ln.p == -Inf] <- Inf 1251 ans[ln.p > 0] <- NaN 1252 } else { 1253 ans <- shape / (scale * epsilon) - log(p) / epsilon - 1254 lambertW((shape / epsilon) * exp(shape / epsilon) * 1255 p^(-scale / epsilon)) / scale 1256 ans[p < 0] <- NaN 1257 ans[p == 0] <- Inf 1258 ans[p == 1] <- 0 1259 ans[p > 1] <- NaN 1260 } 1261 } 1262 1263 ans[epsilon == 0] <- 1264 qgompertz(p = p[epsilon == 0], 1265 shape = shape[epsilon == 0], 1266 scale = scale[epsilon == 0], 1267 lower.tail = lower.tail, 1268 log.p = log.p) 1269 1270 ans[shape <= 0 | scale <= 0 | epsilon < 0] <- NaN 1271 ans 1272} 1273 1274 1275 1276rmakeham <- function(n, scale = 1, shape, epsilon = 0) { 1277 qmakeham(runif(n), scale = scale, shape = shape, epsilon = epsilon) 1278} 1279 1280 1281 1282 1283makeham.control <- function(save.weights = TRUE, ...) { 1284 list(save.weights = save.weights) 1285} 1286 1287 1288 makeham <- 1289 function(lscale = "loglink", lshape = "loglink", lepsilon = "loglink", 1290 iscale = NULL, ishape = NULL, iepsilon = NULL, # 0.3, 1291 gscale = exp(-5:5), 1292 gshape = exp(-5:5), 1293 gepsilon = exp(-4:1), 1294 nsimEIM = 500, 1295 oim.mean = TRUE, 1296 zero = NULL, nowarning = FALSE) { 1297 1298 1299 1300 1301 1302 1303 1304 lepsil <- lepsilon 1305 iepsil <- iepsilon 1306 1307 1308 lshape <- as.list(substitute(lshape)) 1309 eshape <- link2list(lshape) 1310 lshape <- attr(eshape, "function.name") 1311 1312 lscale <- as.list(substitute(lscale)) 1313 escale <- link2list(lscale) 1314 lscale <- attr(escale, "function.name") 1315 1316 lepsil <- as.list(substitute(lepsil)) 1317 eepsil <- link2list(lepsil) 1318 lepsil <- attr(eepsil, "function.name") 1319 1320 if (!is.Numeric(nsimEIM, length.arg = 1, 1321 integer.valued = TRUE)) 1322 stop("bad input for argument 'nsimEIM'") 1323 if (nsimEIM <= 50) 1324 warning("argument 'nsimEIM' should be an integer ", 1325 "greater than 50, say") 1326 1327 1328 if (length(ishape)) 1329 if (!is.Numeric(ishape, positive = TRUE)) 1330 stop("argument 'ishape' values must be positive") 1331 if (length(iscale)) 1332 if (!is.Numeric(iscale, positive = TRUE)) 1333 stop("argument 'iscale' values must be positive") 1334 if (length(iepsil)) 1335 if (!is.Numeric(iepsil, positive = TRUE)) 1336 stop("argument 'iepsil' values must be positive") 1337 1338 1339 1340 1341 1342 if (!is.logical(oim.mean) || length(oim.mean) != 1) 1343 stop("bad input for argument 'oim.mean'") 1344 1345 1346 1347 1348 new("vglmff", 1349 blurb = c("Makeham distribution\n\n", 1350 "Links: ", 1351 namesof("scale", lscale, escale), ", ", 1352 namesof("shape", lshape, eshape), ", ", 1353 namesof("epsilon", lepsil, eepsil), "\n", 1354 "Median: qmakeham(p = 0.5, scale, shape, epsilon)"), 1355 1356 constraints = eval(substitute(expression({ 1357 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 1358 predictors.names = predictors.names, 1359 M1 = 3) 1360 }), list( .zero = zero ))), 1361 1362 infos = eval(substitute(function(...) { 1363 list(M1 = 3, 1364 Q1 = 1, 1365 nsimEIM = .nsimEIM, 1366 parameters.names = c("scale", "shape"), 1367 zero = .zero ) 1368 }, list( .zero = zero, 1369 .nsimEIM = nsimEIM ))), 1370 initialize = eval(substitute(expression({ 1371 1372 1373 temp5 <- 1374 w.y.check(w = w, y = y, 1375 Is.positive.y = TRUE, 1376 ncol.w.max = Inf, 1377 ncol.y.max = Inf, 1378 out.wy = TRUE, 1379 colsyperw = 1, 1380 maximize = TRUE) 1381 w <- temp5$w 1382 y <- temp5$y 1383 1384 1385 ncoly <- ncol(y) 1386 1387 M1 <- 3 1388 extra$ncoly <- ncoly 1389 extra$M1 <- M1 1390 M <- M1 * ncoly 1391 1392 1393 mynames1 <- param.names("scale", ncoly, skip1 = TRUE) 1394 mynames2 <- param.names("shape", ncoly, skip1 = TRUE) 1395 mynames3 <- param.names("epsilon", ncoly, skip1 = TRUE) 1396 predictors.names <- 1397 c(namesof(mynames1, .lscale , .escale , tag = FALSE), 1398 namesof(mynames2, .lshape , .eshape , tag = FALSE), 1399 namesof(mynames3, .lepsil , .eepsil , tag = FALSE))[ 1400 interleave.VGAM(M, M1 = M1)] 1401 1402 1403 if (!length(etastart)) { 1404 1405 matC <- matrix(if (length( .iscale )) .iscale else 0 + NA, 1406 n, ncoly, byrow = TRUE) 1407 matH <- matrix(if (length( .ishape )) .ishape else 0 + NA, 1408 n, ncoly, byrow = TRUE) 1409 1410 matE <- matrix(if (length( .iepsil )) .iepsil else 0.3, 1411 n, ncoly, byrow = TRUE) 1412 1413 1414 shape.grid <- unique(sort(c( .gshape ))) 1415 scale.grid <- unique(sort(c( .gscale ))) 1416 1417 1418 1419 for (spp. in 1:ncoly) { 1420 yvec <- y[, spp.] 1421 wvec <- w[, spp.] 1422 1423 1424 makeham.Loglikfun2 <- function(scaleval, shapeval, 1425 y, x, w, extraargs) { 1426 sum(c(w) * dmakeham(x = y, shape = shapeval, 1427 epsilon = extraargs$Epsil, 1428 scale = scaleval, log = TRUE)) 1429 } 1430 try.this <- 1431 grid.search2(scale.grid, shape.grid, 1432 objfun = makeham.Loglikfun2, 1433 y = yvec, w = wvec, 1434 extraargs = list(Epsilon = matE[1, spp.]), 1435 ret.objfun = TRUE) # Last value is the loglik 1436 if (!length( .iscale )) 1437 matC[, spp.] <- try.this["Value1"] 1438 if (!length( .ishape )) 1439 matH[, spp.] <- try.this["Value2"] 1440 } # spp. 1441 1442 1443 1444 1445 1446 epsil.grid <- c( .gepsil ) 1447 for (spp. in 1:ncoly) { 1448 yvec <- y[, spp.] 1449 wvec <- w[, spp.] 1450 1451 makeham.Loglikfun2 <- function(epsilval, y, x, w, extraargs) { 1452 ans <- 1453 sum(c(w) * dmakeham(x = y, shape = extraargs$Shape, 1454 epsilon = epsilval, 1455 scale = extraargs$Scale, log = TRUE)) 1456 ans 1457 } 1458 Init.epsil <- 1459 grid.search(epsil.grid, 1460 objfun = makeham.Loglikfun2, 1461 y = yvec, x = x, w = wvec, 1462 extraargs = list(Shape = matH[1, spp.], 1463 Scale = matC[1, spp.])) 1464 1465 matE[, spp.] <- Init.epsil 1466 } # spp. 1467 1468 1469 etastart <- cbind(theta2eta(matC, .lscale , .escale ), 1470 theta2eta(matH, .lshape , .eshape ), 1471 theta2eta(matE, .lepsil , .eepsil ))[, 1472 interleave.VGAM(M, M1 = M1)] 1473 } # End of !length(etastart) 1474 }), list( 1475 .lshape = lshape, .lscale = lscale, .lepsil = lepsil, 1476 .eshape = eshape, .escale = escale, .eepsil = eepsil, 1477 .gshape = gshape, .gscale = gscale, .gepsil = gepsilon, 1478 .ishape = ishape, .iscale = iscale, .iepsil = iepsil 1479 ))), 1480 1481 linkinv = eval(substitute(function(eta, extra = NULL) { 1482 scale <- eta2theta(eta[, c(TRUE, FALSE, FALSE)], .lscale , .escale ) 1483 shape <- eta2theta(eta[, c(FALSE, TRUE, FALSE)], .lshape , .eshape ) 1484 epsil <- eta2theta(eta[, c(FALSE, FALSE, TRUE)], .lepsil , .eepsil ) 1485 qmakeham(p = 0.5, scale = scale, shape = shape, epsil = epsil) 1486 }, list( 1487 .lshape = lshape, .lscale = lscale, .lepsil = lepsil, 1488 .eshape = eshape, .escale = escale, .eepsil = eepsil 1489 ))), 1490 last = eval(substitute(expression({ 1491 M1 <- extra$M1 1492 misc$link <- 1493 c(rep_len( .lscale , ncoly), 1494 rep_len( .lshape , ncoly), 1495 rep_len( .lepsil , ncoly))[interleave.VGAM(M, M1 = M1)] 1496 temp.names <- c(mynames1, mynames2, mynames3)[ 1497 interleave.VGAM(M, M1 = M1)] 1498 names(misc$link) <- temp.names 1499 1500 misc$earg <- vector("list", M) 1501 names(misc$earg) <- temp.names 1502 for (ii in 1:ncoly) { 1503 misc$earg[[M1*ii-2]] <- .escale 1504 misc$earg[[M1*ii-1]] <- .eshape 1505 misc$earg[[M1*ii ]] <- .eepsil 1506 } 1507 1508 misc$M1 <- M1 1509 misc$expected <- TRUE 1510 misc$multipleResponses <- TRUE 1511 misc$nsimEIM <- .nsimEIM 1512 }), list( 1513 .lshape = lshape, .lscale = lscale, .lepsil = lepsil, 1514 .eshape = eshape, .escale = escale, .eepsil = eepsil, 1515 .nsimEIM = nsimEIM ))), 1516 loglikelihood = eval(substitute( 1517 function(mu, y, w, residuals = FALSE, eta, 1518 extra = NULL, 1519 summation = TRUE) { 1520 scale <- eta2theta(eta[, c(TRUE, FALSE, FALSE)], .lscale , .escale ) 1521 shape <- eta2theta(eta[, c(FALSE, TRUE, FALSE)], .lshape , .eshape ) 1522 epsil <- eta2theta(eta[, c(FALSE, FALSE, TRUE)], .lepsil , .eepsil ) 1523 if (residuals) { 1524 stop("loglikelihood residuals not implemented yet") 1525 } else { 1526 ll.elts <- c(w) * dmakeham(x = y, scale = scale, shape = shape, 1527 epsil = epsil, log = TRUE) 1528 if (summation) { 1529 sum(ll.elts) 1530 } else { 1531 ll.elts 1532 } 1533 } 1534 }, list( .lshape = lshape, .lscale = lscale, .lepsil = lepsil, 1535 .eshape = eshape, .escale = escale, .eepsil = eepsil ))), 1536 vfamily = c("makeham"), 1537 # ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 1538 validparams = eval(substitute(function(eta, y, extra = NULL) { 1539 Scale <- eta2theta(eta[, c(TRUE, FALSE, FALSE), drop = FALSE], 1540 .lscale , .escale ) 1541 shape <- eta2theta(eta[, c(FALSE, TRUE, FALSE), drop = FALSE], 1542 .lshape , .eshape ) 1543 epsil <- eta2theta(eta[, c(FALSE, FALSE, TRUE), drop = FALSE], 1544 .lepsil , .eepsil ) 1545 okay1 <- all(is.finite(Scale)) && all(0 < Scale) && 1546 all(is.finite(shape)) && all(0 < shape) && 1547 all(is.finite(epsil)) && all(0 < epsil) 1548 okay1 1549 }, list( .lshape = lshape, .lscale = lscale, .lepsil = lepsil, 1550 .eshape = eshape, .escale = escale, .eepsil = eepsil ))), 1551 1552 1553 1554 1555 simslot = eval(substitute( 1556 function(object, nsim) { 1557 1558 pwts <- if (length(pwts <- object@prior.weights) > 0) 1559 pwts else weights(object, type = "prior") 1560 if (any(pwts != 1)) 1561 warning("ignoring prior weights") 1562 eta <- predict(object) 1563 Scale <- eta2theta(eta[, c(TRUE, FALSE, FALSE), drop = FALSE], 1564 .lscale , .escale ) 1565 shape <- eta2theta(eta[, c(FALSE, TRUE, FALSE), drop = FALSE], 1566 .lshape , .eshape ) 1567 epsil <- eta2theta(eta[, c(FALSE, FALSE, TRUE), drop = FALSE], 1568 .lepsil , .eepsil ) 1569 rmakeham(nsim * length(Scale), 1570 scale = c(Scale), shape = c(shape), epsilon = c(epsil)) 1571 }, list( .lshape = lshape, .lscale = lscale, .lepsil = lepsil, 1572 .eshape = eshape, .escale = escale, .eepsil = eepsil ))), 1573 1574 1575 1576 1577 deriv = eval(substitute(expression({ 1578 scale <- eta2theta(eta[, c(TRUE, FALSE, FALSE), drop = FALSE], 1579 .lscale , .escale ) 1580 shape <- eta2theta(eta[, c(FALSE, TRUE, FALSE), drop = FALSE], 1581 .lshape , .eshape ) 1582 epsil <- eta2theta(eta[, c(FALSE, FALSE, TRUE), drop = FALSE], 1583 .lepsil , .eepsil ) 1584 1585 temp2 <- exp(y * scale) 1586 temp3 <- epsil + shape * temp2 1587 dl.dshape <- temp2 / temp3 - expm1(y * scale) / scale 1588 dl.dscale <- shape * y * temp2 / temp3 + 1589 shape * expm1(y * scale) / scale^2 - 1590 shape * y * temp2 / scale 1591 dl.depsil <- 1 / temp3 - y 1592 1593 dshape.deta <- dtheta.deta(shape, .lshape , .eshape ) 1594 dscale.deta <- dtheta.deta(scale, .lscale , .escale ) 1595 depsil.deta <- dtheta.deta(epsil, .lepsil , .eepsil ) 1596 1597 dthetas.detas <- cbind(dscale.deta, dshape.deta, depsil.deta) 1598 myderiv <- c(w) * cbind(dl.dscale, 1599 dl.dshape, 1600 dl.depsil) * dthetas.detas 1601 myderiv[, interleave.VGAM(M, M1 = M1)] 1602 }), list( .lshape = lshape, .lscale = lscale, .lepsil = lepsil, 1603 .eshape = eshape, .escale = escale, .eepsil = eepsil ))), 1604 1605 weight = eval(substitute(expression({ 1606 NOS <- M / M1 1607 dThetas.detas <- dthetas.detas[, interleave.VGAM(M, M1 = M1)] 1608 wz <- matrix(0.0, n, M + M - 1 + M - 2) # wz has half-bandwidth 3 1609 1610 ind1 <- iam(NA, NA, M = M1, both = TRUE, diag = TRUE) # Use SFS 1611 1612 for (spp. in 1:NOS) { 1613 run.varcov <- 0 1614 Shape <- shape[, spp.] 1615 Scale <- scale[, spp.] 1616 Epsil <- epsil[, spp.] 1617 1618 for (ii in 1:( .nsimEIM )) { 1619 ysim <- rmakeham(n, scale = Scale, shape = Shape, epsil = Epsil) 1620 temp2 <- exp(ysim * Scale) 1621 temp3 <- Epsil + Shape * temp2 1622 dl.dshape <- temp2 / temp3 - expm1(ysim * Scale) / Scale 1623 dl.dscale <- Shape * ysim * temp2 / temp3 + 1624 Shape * expm1(ysim * Scale) / Scale^2 - 1625 Shape * ysim * temp2 / Scale 1626 dl.depsil <- 1 / temp3 - ysim 1627 1628 temp7 <- cbind(dl.dscale, dl.dshape, dl.depsil) 1629 run.varcov <- run.varcov + temp7[, ind1$row.index] * 1630 temp7[, ind1$col.index] 1631 } 1632 run.varcov <- cbind(run.varcov / .nsimEIM ) 1633 1634 1635 wz1 <- if (intercept.only) 1636 matrix(colMeans(run.varcov, na.rm = TRUE), 1637 nrow = n, ncol = ncol(run.varcov), byrow = TRUE) else 1638 run.varcov 1639 1640 wz1 <- wz1 * dThetas.detas[, M1 * (spp. - 1) + ind1$row] * 1641 dThetas.detas[, M1 * (spp. - 1) + ind1$col] 1642 for (jay in 1:M1) 1643 for (kay in jay:M1) { # Now copy wz1 into wz 1644 cptr <- iam((spp. - 1) * M1 + jay, 1645 (spp. - 1) * M1 + kay, M = M) 1646 wz[, cptr] <- wz1[, iam(jay, kay, M = M1)] 1647 } 1648 } # End of for (spp.) loop 1649 w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = M / M1) 1650 }), list( .lshape = lshape, .lscale = lscale, .lepsil = lepsil, 1651 .eshape = eshape, .escale = escale, .eepsil = eepsil, 1652 .nsimEIM = nsimEIM, .oim.mean = oim.mean )))) 1653} # makeham() 1654 1655 1656 1657 1658 1659 1660 1661 1662dgompertz <- function(x, scale = 1, shape, log = FALSE) { 1663 1664 if (!is.logical(log.arg <- log) || length(log) != 1) 1665 stop("bad input for argument 'log'") 1666 rm(log) 1667 1668 LLL <- max(length(x), length(shape), length(scale)) 1669 if (length(x) != LLL) x <- rep_len(x, LLL) 1670 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 1671 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 1672 1673 1674 index0 <- (x < 0) 1675 index1 <- abs(x * scale) < 0.1 & is.finite(x * scale) 1676 ans <- log(shape) + x * scale - (shape / scale) * (exp(x * scale) - 1) 1677 ans[index1] <- log(shape[index1]) + x[index1] * scale[index1] - 1678 (shape[index1] / scale[index1]) * 1679 expm1(x[index1] * scale[index1]) 1680 ans[index0] <- log(0) 1681 ans[x == Inf] <- log(0) 1682 if (log.arg) { 1683 } else { 1684 ans <- exp(ans) 1685 ans[index0] <- 0 1686 ans[x == Inf] <- 0 1687 } 1688 ans[shape <= 0 | scale <= 0] <- NaN 1689 ans 1690} 1691 1692 1693 1694pgompertz <- function(q, scale = 1, shape, 1695 lower.tail = TRUE, log.p = FALSE) { 1696 1697 1698 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 1699 stop("bad input for argument 'lower.tail'") 1700 1701 if (!is.logical(log.p) || length(log.p) != 1) 1702 stop("bad input for argument 'log.p'") 1703 1704 LLL <- max(length(q), length(shape), length(scale)) 1705 if (length(q) != LLL) q <- rep_len(q, LLL) 1706 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 1707 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 1708 1709 1710 if (lower.tail) { 1711 if (log.p) { 1712 ans <- log1p(-exp((-shape / scale) * expm1(scale * q))) 1713 ans[q <= 0 ] <- -Inf 1714 ans[q == Inf] <- 0 1715 } else { 1716 ans <- -expm1((-shape / scale) * expm1(scale * q)) 1717 ans[q <= 0] <- 0 1718 ans[q == Inf] <- 1 1719 } 1720 } else { 1721 if (log.p) { 1722 ans <- (-shape / scale) * expm1(scale * q) 1723 ans[q <= 0] <- 0 1724 ans[q == Inf] <- -Inf 1725 } else { 1726 ans <- exp((-shape / scale) * expm1(scale * q)) 1727 ans[q <= 0] <- 1 1728 ans[q == Inf] <- 0 1729 } 1730 } 1731 ans[shape <= 0 | scale <= 0] <- NaN 1732 ans 1733} 1734 1735 1736 1737qgompertz <- function(p, scale = 1, shape, 1738 lower.tail = TRUE, log.p = FALSE) { 1739 1740 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 1741 stop("bad input for argument 'lower.tail'") 1742 if (!is.logical(log.p) || length(log.p) != 1) 1743 stop("bad input for argument 'log.p'") 1744 1745 LLL <- max(length(p), length(shape), length(scale)) 1746 if (length(p) != LLL) p <- rep_len(p, LLL) 1747 if (length(shape) != LLL) shape <- rep_len(shape, LLL) 1748 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 1749 1750 if (lower.tail) { 1751 if (log.p) { 1752 ln.p <- p 1753 ans <- log1p((-scale / shape) * log(-expm1(ln.p))) / scale 1754 ans[ln.p > 0] <- NaN 1755 } else { 1756 ans <- log1p((-scale / shape) * log1p(-p)) / scale 1757 ans[p < 0] <- NaN 1758 ans[p == 0] <- 0 1759 ans[p == 1] <- Inf 1760 ans[p > 1] <- NaN 1761 } 1762 } else { 1763 if (log.p) { 1764 ln.p <- p 1765 ans <- log1p((-scale / shape) * ln.p) / scale 1766 ans[ln.p > 0] <- NaN 1767 } else { 1768 ans <- log1p((-scale / shape) * log(p)) / scale 1769 ans[p < 0] <- NaN 1770 ans[p == 0] <- Inf 1771 ans[p == 1] <- 0 1772 ans[p > 1] <- NaN 1773 } 1774 } 1775 ans[shape <= 0 | scale <= 0] <- NaN 1776 ans 1777} 1778 1779 1780 1781 1782 1783rgompertz <- function(n, scale = 1, shape) { 1784 qgompertz(runif(n), scale = scale, shape = shape) 1785} 1786 1787 1788 1789 1790 1791 1792 1793gompertz.control <- function(save.weights = TRUE, ...) { 1794 list(save.weights = save.weights) 1795} 1796 1797 1798 gompertz <- 1799 function(lscale = "loglink", lshape = "loglink", 1800 iscale = NULL, ishape = NULL, 1801 nsimEIM = 500, 1802 zero = NULL, nowarning = FALSE) { 1803 1804 1805 1806 1807 1808 lshape <- as.list(substitute(lshape)) 1809 eshape <- link2list(lshape) 1810 lshape <- attr(eshape, "function.name") 1811 1812 lscale <- as.list(substitute(lscale)) 1813 escale <- link2list(lscale) 1814 lscale <- attr(escale, "function.name") 1815 1816 1817 1818 if (!is.Numeric(nsimEIM, length.arg = 1, 1819 integer.valued = TRUE)) 1820 stop("bad input for argument 'nsimEIM'") 1821 if (nsimEIM <= 50) 1822 warning("argument 'nsimEIM' should be an integer ", 1823 "greater than 50, say") 1824 1825 1826 if (length(ishape)) 1827 if (!is.Numeric(ishape, positive = TRUE)) 1828 stop("argument 'ishape' values must be positive") 1829 if (length(iscale)) 1830 if (!is.Numeric(iscale, positive = TRUE)) 1831 stop("argument 'iscale' values must be positive") 1832 1833 1834 1835 1836 1837 new("vglmff", 1838 blurb = c("Gompertz distribution\n\n", 1839 "Links: ", 1840 namesof("scale", lscale, escale ), ", ", 1841 namesof("shape", lshape, eshape ), "\n", 1842 "Median: scale * log(2 - 1 / shape)"), 1843 1844 constraints = eval(substitute(expression({ 1845 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 1846 predictors.names = predictors.names, 1847 M1 = 2) 1848 }), list( .zero = zero ))), 1849 1850 infos = eval(substitute(function(...) { 1851 list(M1 = 2, 1852 Q1 = 1, 1853 nsimEIM = .nsimEIM, 1854 parameters.names = c("scale", "shape"), 1855 zero = .zero ) 1856 }, list( .zero = zero, 1857 .nsimEIM = nsimEIM ))), 1858 initialize = eval(substitute(expression({ 1859 1860 temp5 <- 1861 w.y.check(w = w, y = y, 1862 Is.positive.y = TRUE, 1863 ncol.w.max = Inf, 1864 ncol.y.max = Inf, 1865 out.wy = TRUE, 1866 colsyperw = 1, 1867 maximize = TRUE) 1868 w <- temp5$w 1869 y <- temp5$y 1870 1871 1872 1873 ncoly <- ncol(y) 1874 M1 <- 2 1875 extra$ncoly <- ncoly 1876 extra$M1 <- M1 1877 M <- M1 * ncoly 1878 1879 1880 mynames1 <- param.names("scale", ncoly, skip1 = TRUE) 1881 mynames2 <- param.names("shape", ncoly, skip1 = TRUE) 1882 predictors.names <- 1883 c(namesof(mynames1, .lscale , .escale , tag = FALSE), 1884 namesof(mynames2, .lshape , .eshape , tag = FALSE))[ 1885 interleave.VGAM(M, M1 = M1)] 1886 1887 1888 1889 if (!length(etastart)) { 1890 1891 matH <- matrix(if (length( .ishape )) .ishape else 0 + NA, 1892 n, ncoly, byrow = TRUE) 1893 matC <- matrix(if (length( .iscale )) .iscale else 0 + NA, 1894 n, ncoly, byrow = TRUE) 1895 1896 shape.grid <- c(exp(-seq(4, 0.1, len = 07)), 1, 1897 exp( seq(0.1, 4, len = 07))) 1898 scale.grid <- c(exp(-seq(4, 0.1, len = 07)), 1, 1899 exp( seq(0.1, 4, len = 07))) 1900 1901 for (spp. in 1:ncoly) { 1902 yvec <- y[, spp.] 1903 wvec <- w[, spp.] 1904 1905 1906 gompertz.Loglikfun <- function(scaleval, y, x, w, extraargs) { 1907 ans <- 1908 sum(c(w) * dgompertz(x = y, shape = extraargs$Shape, 1909 scale = scaleval, log = TRUE)) 1910 ans 1911 } 1912 1913 mymat <- matrix(-1, length(shape.grid), 2) 1914 for (jlocal in seq_along(shape.grid)) { 1915 mymat[jlocal, ] <- 1916 grid.search(scale.grid, 1917 objfun = gompertz.Loglikfun, 1918 y = yvec, x = x, w = wvec, 1919 ret.objfun = TRUE, 1920 extraargs = list(Shape = shape.grid[jlocal])) 1921 } 1922 index.shape <- which(mymat[, 2] == max(mymat[, 2]))[1] 1923 1924 if (!length( .ishape )) 1925 matH[, spp.] <- shape.grid[index.shape] 1926 if (!length( .iscale )) 1927 matC[, spp.] <- mymat[index.shape, 1] 1928 } # spp. 1929 1930 etastart <- cbind(theta2eta(matC, .lscale , .escale ), 1931 theta2eta(matH, .lshape , .eshape ))[, 1932 interleave.VGAM(M, M1 = M1)] 1933 } # End of !length(etastart) 1934 }), list( .lshape = lshape, .lscale = lscale, 1935 .eshape = eshape, .escale = escale, 1936 .ishape = ishape, .iscale = iscale 1937 ))), 1938 1939 linkinv = eval(substitute(function(eta, extra = NULL) { 1940 scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale ) 1941 shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape ) 1942 log1p((scale / shape) * log(2)) / scale 1943 }, list( .lshape = lshape, .lscale = lscale, 1944 .eshape = eshape, .escale = escale ))), 1945 last = eval(substitute(expression({ 1946 M1 <- extra$M1 1947 misc$link <- 1948 c(rep_len( .lscale , ncoly), 1949 rep_len( .lshape , ncoly))[interleave.VGAM(M, M1 = M1)] 1950 temp.names <- c(mynames1, mynames2)[interleave.VGAM(M, M1 = M1)] 1951 names(misc$link) <- temp.names 1952 1953 misc$earg <- vector("list", M) 1954 names(misc$earg) <- temp.names 1955 for (ii in 1:ncoly) { 1956 misc$earg[[M1*ii-1]] <- .escale 1957 misc$earg[[M1*ii ]] <- .eshape 1958 } 1959 1960 misc$M1 <- M1 1961 misc$expected <- TRUE 1962 misc$multipleResponses <- TRUE 1963 misc$nsimEIM <- .nsimEIM 1964 }), list( .lshape = lshape, .lscale = lscale, 1965 .eshape = eshape, .escale = escale, 1966 .nsimEIM = nsimEIM ))), 1967 loglikelihood = eval(substitute( 1968 function(mu, y, w, residuals = FALSE, eta, 1969 extra = NULL, 1970 summation = TRUE) { 1971 scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale ) 1972 shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape ) 1973 if (residuals) { 1974 stop("loglikelihood residuals not implemented yet") 1975 } else { 1976 ll.elts <- c(w) * dgompertz(x = y, scale = scale, 1977 shape = shape, log = TRUE) 1978 if (summation) { 1979 sum(ll.elts) 1980 } else { 1981 ll.elts 1982 } 1983 } 1984 }, list( .lshape = lshape, .lscale = lscale, 1985 .eshape = eshape, .escale = escale ))), 1986 vfamily = c("gompertz"), 1987 # ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 1988 validparams = eval(substitute(function(eta, y, extra = NULL) { 1989 Scale <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE], .lscale , 1990 .escale ) 1991 shape <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE], .lshape , 1992 .eshape ) 1993 okay1 <- all(is.finite(Scale)) && all(0 < Scale) && 1994 all(is.finite(shape)) && all(0 < shape) 1995 okay1 1996 }, list( .lshape = lshape, .lscale = lscale, 1997 .eshape = eshape, .escale = escale ))), 1998 1999 2000 2001 2002 2003 simslot = eval(substitute( 2004 function(object, nsim) { 2005 2006 pwts <- if (length(pwts <- object@prior.weights) > 0) 2007 pwts else weights(object, type = "prior") 2008 if (any(pwts != 1)) 2009 warning("ignoring prior weights") 2010 eta <- predict(object) 2011 Scale <- eta2theta(eta[, c(TRUE, FALSE)], .lscale , .escale ) 2012 Shape <- eta2theta(eta[, c(FALSE, TRUE)], .lshape , .eshape ) 2013 rgompertz(nsim * length(Scale), 2014 shape = c(Shape), scale = c(Scale)) 2015 }, list( .lshape = lshape, .lscale = lscale, 2016 .eshape = eshape, .escale = escale ))), 2017 2018 2019 2020 deriv = eval(substitute(expression({ 2021 M1 <- 2 2022 scale <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE], .lscale , 2023 .escale ) 2024 shape <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE], .lshape , 2025 .eshape ) 2026 2027 2028 temp2 <- exp(y * scale) 2029 temp4 <- -expm1(y * scale) 2030 dl.dshape <- 1 / shape + temp4 / scale 2031 dl.dscale <- y * (1 - shape * temp2 / scale) - 2032 shape * temp4 / scale^2 2033 2034 dscale.deta <- dtheta.deta(scale, .lscale , .escale ) 2035 dshape.deta <- dtheta.deta(shape, .lshape , .eshape ) 2036 2037 dthetas.detas <- cbind(dscale.deta, dshape.deta) 2038 myderiv <- c(w) * cbind(dl.dscale, dl.dshape) * dthetas.detas 2039 myderiv[, interleave.VGAM(M, M1 = M1)] 2040 }), list( .lshape = lshape, .lscale = lscale, 2041 .eshape = eshape, .escale = escale ))), 2042 2043 2044 weight = eval(substitute(expression({ 2045 2046 NOS <- M / M1 2047 dThetas.detas <- dthetas.detas[, interleave.VGAM(M, M1 = M1)] 2048 2049 wz <- matrix(0.0, n, M + M - 1) # wz is 'tridiagonal' 2050 2051 ind1 <- iam(NA, NA, M = M1, both = TRUE, diag = TRUE) 2052 2053 2054 for (spp. in 1:NOS) { 2055 run.varcov <- 0 2056 Shape <- shape[, spp.] 2057 Scale <- scale[, spp.] 2058 2059 for (ii in 1:( .nsimEIM )) { 2060 ysim <- rgompertz(n = n, shape = Shape, scale = Scale) 2061if (ii < 3) { 2062} 2063 2064 temp2 <- exp(ysim * scale) 2065 temp4 <- -expm1(ysim * scale) 2066 dl.dshape <- 1 / shape + temp4 / scale 2067 dl.dscale <- ysim * (1 - shape * temp2 / scale) - 2068 shape * temp4 / scale^2 2069 2070 2071 temp7 <- cbind(dl.dscale, dl.dshape) 2072 run.varcov <- run.varcov + 2073 temp7[, ind1$row.index] * 2074 temp7[, ind1$col.index] 2075 } 2076 run.varcov <- cbind(run.varcov / .nsimEIM ) 2077 2078 wz1 <- if (intercept.only) 2079 matrix(colMeans(run.varcov), 2080 nrow = n, ncol = ncol(run.varcov), byrow = TRUE) else 2081 run.varcov 2082 2083 wz1 <- wz1 * dThetas.detas[, M1 * (spp. - 1) + ind1$row] * 2084 dThetas.detas[, M1 * (spp. - 1) + ind1$col] 2085 2086 2087 for (jay in 1:M1) 2088 for (kay in jay:M1) { 2089 cptr <- iam((spp. - 1) * M1 + jay, 2090 (spp. - 1) * M1 + kay, 2091 M = M) 2092 wz[, cptr] <- wz1[, iam(jay, kay, M = M1)] 2093 } 2094 } # End of for (spp.) loop 2095 2096 2097 2098 w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = M / M1) 2099 }), list( .lscale = lscale, 2100 .escale = escale, 2101 .nsimEIM = nsimEIM )))) 2102} # gompertz() 2103 2104 2105 2106 2107 2108 2109 dmoe <- function (x, alpha = 1, lambda = 1, log = FALSE) { 2110 if (!is.logical(log.arg <- log) || length(log) != 1) 2111 stop("bad input for argument 'log'") 2112 rm(log) 2113 2114 LLL <- max(length(x), 2115 length(alpha), 2116 length(lambda)) 2117 if (length(x) != LLL) x <- rep_len(x, LLL) 2118 if (length(alpha) != LLL) alpha <- rep_len(alpha, LLL) 2119 if (length(lambda) != LLL) lambda <- rep_len(lambda, LLL) 2120 2121 index0 <- (x < 0) 2122 if (log.arg) { 2123 ans <- log(lambda) + (lambda * x) - 2124 2 * log(expm1(lambda * x) + alpha) 2125 ans[index0] <- log(0) 2126 } else { 2127 ans <- lambda * exp(lambda * x) / (expm1(lambda * x) + alpha)^2 2128 ans[index0] <- 0 2129 } 2130 ans[alpha <= 0 | lambda <= 0] <- NaN 2131 ans 2132} 2133 2134 2135 2136 pmoe <- function (q, alpha = 1, lambda = 1) { 2137 ret <- ifelse(alpha <= 0 | lambda <= 0, NaN, 2138 1 - 1 / (expm1(lambda * q) + alpha)) 2139 ret[q < log(2 - alpha) / lambda] <- 0 2140 ret 2141} 2142 2143 2144 2145qmoe <- function (p, alpha = 1, lambda = 1) { 2146 ifelse(p < 0 | p > 1 | alpha <= 0 | lambda <= 0, NaN, 2147 log1p(-alpha + 1 / (1 - p)) / lambda) 2148} 2149 2150 2151 2152rmoe <- function (n, alpha = 1, lambda = 1) { 2153 qmoe(p = runif(n), alpha = alpha, lambda = lambda) 2154} 2155 2156 2157 2158 2159exponential.mo.control <- function(save.weights = TRUE, ...) { 2160 list(save.weights = save.weights) 2161} 2162 2163 2164 2165 2166 exponential.mo <- 2167 function(lalpha = "loglink", llambda = "loglink", 2168 ealpha = list(), elambda = list(), 2169 ialpha = 1, ilambda = NULL, 2170 imethod = 1, 2171 nsimEIM = 200, 2172 zero = NULL) { 2173 2174 stop("fundamentally unable to estimate the parameters as ", 2175 "the support of the density depends on the parameters") 2176 2177 2178 lalpha <- as.list(substitute(lalpha)) 2179 ealpha <- link2list(lalpha) 2180 lalpha <- attr(ealpha, "function.name") 2181 2182 llambda <- as.list(substitute(llambda)) 2183 elambda <- link2list(llambda) 2184 llambda <- attr(elambda, "function.name") 2185 2186 lalpha0 <- lalpha 2187 ealpha0 <- ealpha 2188 ialpha0 <- ialpha 2189 2190 2191 2192 if (!is.Numeric(nsimEIM, length.arg = 1, 2193 integer.valued = TRUE)) 2194 stop("bad input for argument 'nsimEIM'") 2195 if (nsimEIM <= 50) 2196 warning("argument 'nsimEIM' should be an integer ", 2197 "greater than 50, say") 2198 2199 if (length(ialpha0)) 2200 if (!is.Numeric(ialpha0, positive = TRUE)) 2201 stop("argument 'ialpha' values must be positive") 2202 if (length(ilambda)) 2203 if (!is.Numeric(ilambda, positive = TRUE)) 2204 stop("argument 'ilambda' values must be positive") 2205 2206 2207 if (!is.Numeric(imethod, length.arg = 1, 2208 integer.valued = TRUE, positive = TRUE) || 2209 imethod > 2) 2210 stop("argument 'imethod' must be 1 or 2") 2211 2212 2213 2214 new("vglmff", 2215 blurb = c("Marshall-Olkin exponential distribution\n\n", 2216 "Links: ", 2217 namesof("alpha", lalpha0, ealpha0 ), ", ", 2218 namesof("lambda", llambda, elambda ), "\n", 2219 "Median: log(3 - alpha) / lambda"), 2220 2221 constraints = eval(substitute(expression({ 2222 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 2223 predictors.names = predictors.names, 2224 M1 = 2) 2225 }), list( .zero = zero ))), 2226 2227 infos = eval(substitute(function(...) { 2228 list(M1 = 2, 2229 Q1 = 1, 2230 nsimEIM = .nsimEIM, 2231 parameters.names = c("alpha", "lambda"), 2232 zero = .zero ) 2233 }, list( .zero = zero, 2234 .nsimEIM = nsimEIM ))), 2235 initialize = eval(substitute(expression({ 2236 2237 temp5 <- 2238 w.y.check(w = w, y = y, 2239 Is.positive.y = TRUE, 2240 ncol.w.max = Inf, 2241 ncol.y.max = Inf, 2242 out.wy = TRUE, 2243 colsyperw = 1, 2244 maximize = TRUE) 2245 w <- temp5$w 2246 y <- temp5$y 2247 2248 2249 2250 ncoly <- ncol(y) 2251 2252 M1 <- 2 2253 extra$ncoly <- ncoly 2254 extra$M1 <- M1 2255 M <- M1 * ncoly 2256 2257 2258 mynames1 <- param.names("alpha", ncoly, skip1 = TRUE) 2259 mynames2 <- param.names("lambda", ncoly, skip1 = TRUE) 2260 predictors.names <- 2261 c(namesof(mynames1, .lalpha0 , .ealpha0 , tag = FALSE), 2262 namesof(mynames2, .llambda , .elambda , tag = FALSE))[ 2263 interleave.VGAM(M, M1 = M1)] 2264 2265 2266 2267 if (!length(etastart)) { 2268 2269 matL <- matrix(if (length( .ilambda )) .ilambda else 0, 2270 n, ncoly, byrow = TRUE) 2271 matA <- matrix(if (length( .ialpha0 )) .ialpha0 else 0, 2272 n, ncoly, byrow = TRUE) 2273 2274 2275 for (spp. in 1:ncoly) { 2276 yvec <- y[, spp.] 2277 2278 moexpon.Loglikfun <- function(lambdaval, y, x, w, extraargs) { 2279 ans <- 2280 sum(c(w) * log(dmoe(x = y, alpha = extraargs$alpha, 2281 lambda = lambdaval))) 2282 ans 2283 } 2284 Alpha.init <- .ialpha0 2285 lambda.grid <- seq(0.1, 10.0, len = 21) 2286 Lambda.init <- grid.search(lambda.grid, 2287 objfun = moexpon.Loglikfun, 2288 y = y, x = x, w = w, 2289 extraargs = list(alpha = Alpha.init)) 2290 2291 if (length(mustart)) { 2292 Lambda.init <- Lambda.init / (1 - Phimat.init) 2293 } 2294 2295 if (!length( .ialpha0 )) 2296 matA[, spp.] <- Alpha0.init 2297 if (!length( .ilambda )) 2298 matL[, spp.] <- Lambda.init 2299 } # spp. 2300 2301 etastart <- cbind(theta2eta(matA, .lalpha0, .ealpha0 ), 2302 theta2eta(matL, .llambda, .elambda ))[, 2303 interleave.VGAM(M, M1 = M1)] 2304 mustart <- NULL # Since etastart has been computed. 2305 } # End of !length(etastart) 2306 }), list( .lalpha0 = lalpha0, .llambda = llambda, 2307 .ealpha0 = ealpha0, .elambda = elambda, 2308 .ialpha0 = ialpha0, .ilambda = ilambda, 2309 .imethod = imethod 2310 ))), 2311 2312 linkinv = eval(substitute(function(eta, extra = NULL) { 2313 alpha0 <- eta2theta(eta[, c(TRUE, FALSE)], .lalpha0 , .ealpha0 ) 2314 lambda <- eta2theta(eta[, c(FALSE, TRUE)], .llambda , .elambda ) 2315 log(3 - alpha0) / lambda 2316 }, list( .lalpha0 = lalpha0, .llambda = llambda, 2317 .ealpha0 = ealpha0, .elambda = elambda ))), 2318 last = eval(substitute(expression({ 2319 M1 <- extra$M1 2320 misc$link <- 2321 c(rep_len( .lalpha0 , ncoly), 2322 rep_len( .llambda , ncoly))[interleave.VGAM(M, M1 = M1)] 2323 temp.names <- c(mynames1, mynames2)[interleave.VGAM(M, M1 = M1)] 2324 names(misc$link) <- temp.names 2325 2326 misc$earg <- vector("list", M) 2327 names(misc$earg) <- temp.names 2328 for (ii in 1:ncoly) { 2329 misc$earg[[M1*ii-1]] <- .ealpha0 2330 misc$earg[[M1*ii ]] <- .elambda 2331 } 2332 2333 misc$M1 <- M1 2334 misc$imethod <- .imethod 2335 misc$expected <- TRUE 2336 misc$multipleResponses <- TRUE 2337 misc$nsimEIM <- .nsimEIM 2338 }), list( .lalpha0 = lalpha0, .llambda = llambda, 2339 .ealpha0 = ealpha0, .elambda = elambda, 2340 .nsimEIM = nsimEIM, 2341 .imethod = imethod ))), 2342 loglikelihood = eval(substitute( 2343 function(mu, y, w, residuals = FALSE, eta, 2344 extra = NULL, 2345 summation = TRUE) { 2346 alpha0 <- eta2theta(eta[, c(TRUE, FALSE)], .lalpha0 , .ealpha0 ) 2347 lambda <- eta2theta(eta[, c(FALSE, TRUE)], .llambda , .elambda ) 2348 if (residuals) { 2349 stop("loglikelihood residuals not implemented yet") 2350 } else { 2351 ll.elts <- c(w) * log(dmoe(x = y, alpha = alpha0, 2352 lambda = lambda)) 2353 if (summation) { 2354 sum(ll.elts) 2355 } else { 2356 ll.elts 2357 } 2358 } 2359 }, list( .lalpha0 = lalpha0, .llambda = llambda, 2360 .ealpha0 = ealpha0, .elambda = elambda ))), 2361 vfamily = c("exponential.mo"), 2362 validparams = eval(substitute(function(eta, y, extra = NULL) { 2363 alpha0 <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE], .lalpha0 , 2364 .ealpha0 ) 2365 lambda <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE], .llambda , 2366 .elambda ) 2367 okay1 <- all(is.finite(alpha0)) && all(0 < alpha0) && 2368 all(is.finite(lambda)) && all(0 < lambda) 2369 okay1 2370 }, list( .lalpha0 = lalpha0, .llambda = llambda, 2371 .ealpha0 = ealpha0, .elambda = elambda ))), 2372 2373 2374 deriv = eval(substitute(expression({ 2375 M1 <- 2 2376 alpha0 <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE], .lalpha0 , 2377 .ealpha0 ) 2378 lambda <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE], .llambda , 2379 .elambda ) 2380 2381 temp2 <- (expm1(lambda * y) + alpha0) 2382 dl.dalpha0 <- -2 / temp2 2383 dl.dlambda <- 1 / lambda + y - 2 * y * exp(lambda * y) / temp2 2384 2385 dalpha0.deta <- dtheta.deta(alpha0, .lalpha0 , .ealpha0 ) 2386 dlambda.deta <- dtheta.deta(lambda, .llambda , .elambda ) 2387 2388 dthetas.detas <- cbind(dalpha0.deta, 2389 dlambda.deta) 2390 myderiv <- c(w) * cbind(dl.dalpha0, dl.dlambda) * dthetas.detas 2391 myderiv[, interleave.VGAM(M, M1 = M1)] 2392 }), list( .lalpha0 = lalpha0, .llambda = llambda, 2393 .ealpha0 = ealpha0, .elambda = elambda ))), 2394 2395 2396 weight = eval(substitute(expression({ 2397 2398 NOS <- M / M1 2399 dThetas.detas <- dthetas.detas[, interleave.VGAM(M, M1 = M1)] 2400 2401 wz <- matrix(0.0, n, M + M - 1) # wz is 'tridiagonal' 2402 2403 ind1 <- iam(NA, NA, M = M1, both = TRUE, diag = TRUE) 2404 2405 2406 for (spp. in 1:NOS) { 2407 run.varcov <- 0 2408 Alph <- alpha0[, spp.] 2409 Lamb <- lambda[, spp.] 2410 2411 for (ii in 1:( .nsimEIM )) { 2412 ysim <- rmoe(n = n, alpha = Alph, lambda = Lamb) 2413if (ii < 3) { 2414} 2415 2416 temp2 <- (expm1(lambda * ysim) + alpha0) 2417 dl.dalpha0 <- -2 / temp2 2418 dl.dlambda <- 1 / lambda + ysim - 2419 2 * ysim * exp(lambda * ysim) / temp2 2420 2421 2422 temp3 <- cbind(dl.dalpha0, dl.dlambda) 2423 run.varcov <- run.varcov + 2424 temp3[, ind1$row.index] * 2425 temp3[, ind1$col.index] 2426 } 2427 run.varcov <- cbind(run.varcov / .nsimEIM) 2428 2429 wz1 <- if (intercept.only) 2430 matrix(colMeans(run.varcov), 2431 nrow = n, ncol = ncol(run.varcov), byrow = TRUE) else 2432 run.varcov 2433 2434 wz1 <- wz1 * dThetas.detas[, M1 * (spp. - 1) + ind1$row] * 2435 dThetas.detas[, M1 * (spp. - 1) + ind1$col] 2436 2437 2438 for (jay in 1:M1) 2439 for (kay in jay:M1) { 2440 cptr <- iam((spp. - 1) * M1 + jay, 2441 (spp. - 1) * M1 + kay, 2442 M = M) 2443 wz[, cptr] <- wz1[, iam(jay, kay, M = M1)] 2444 } 2445 } # End of for (spp.) loop 2446 2447 2448 2449 2450 w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = M / M1) 2451 }), list( .llambda = llambda, 2452 .elambda = elambda, 2453 .nsimEIM = nsimEIM )))) 2454} # exponential.mo() 2455 2456 2457 2458 2459 2460 2461 2462 2463 2464 genbetaII.Loglikfun4 <- 2465 function(scaleval, shape1.a, shape2.p, shape3.q, 2466 y, x, w, extraargs) { 2467 sum(c(w) * dgenbetaII(x = y, scale = scaleval, 2468 shape1.a = shape1.a, 2469 shape2.p = shape2.p, 2470 shape3.q = shape3.q, log = TRUE)) 2471 } 2472 2473 2474 2475 genbetaII <- function(lscale = "loglink", 2476 lshape1.a = "loglink", 2477 lshape2.p = "loglink", 2478 lshape3.q = "loglink", 2479 iscale = NULL, 2480 ishape1.a = NULL, 2481 ishape2.p = NULL, 2482 ishape3.q = NULL, 2483 lss = TRUE, 2484 gscale = exp(-5:5), 2485 gshape1.a = exp(-5:5), 2486 gshape2.p = exp(-5:5), 2487 gshape3.q = exp(-5:5), 2488 zero = "shape") { 2489 2490 2491 2492 2493 2494 if (length(lss) != 1 && !is.logical(lss)) 2495 stop("Argument 'lss' not specified correctly") 2496 2497 2498 if (length(iscale) && !is.Numeric(iscale, positive = TRUE)) 2499 stop("Bad input for argument 'iscale'") 2500 2501 if (length(ishape1.a) && !is.Numeric(ishape1.a, positive = TRUE)) 2502 stop("Bad input for argument 'ishape1.a'") 2503 2504 if (length(ishape2.p) && !is.Numeric(ishape2.p, positive = TRUE)) 2505 stop("Bad input for argument 'ishape2.p'") 2506 2507 if (length(ishape3.q) && !is.Numeric(ishape3.q, positive = TRUE)) 2508 stop("Bad input for argument 'ishape3.q'") 2509 2510 2511 2512 lscale <- as.list(substitute(lscale)) 2513 escale <- link2list(lscale) 2514 lscale <- attr(escale, "function.name") 2515 2516 lshape1.a <- as.list(substitute(lshape1.a)) 2517 eshape1.a <- link2list(lshape1.a) 2518 lshape1.a <- attr(eshape1.a, "function.name") 2519 2520 lshape2.p <- as.list(substitute(lshape2.p)) 2521 eshape2.p <- link2list(lshape2.p) 2522 lshape2.p <- attr(eshape2.p, "function.name") 2523 2524 lshape3.q <- as.list(substitute(lshape3.q)) 2525 eshape3.q <- link2list(lshape3.q) 2526 lshape3.q <- attr(eshape3.q, "function.name") 2527 2528 2529 new("vglmff", 2530 blurb = 2531 c("Generalized Beta II distribution \n\n", 2532 "Links: ", 2533 ifelse (lss, 2534 namesof("scale" , lscale , earg = escale), 2535 namesof("shape1.a", lshape1.a, earg = eshape1.a)), ", ", 2536 ifelse (lss, 2537 namesof("shape1.a", lshape1.a, earg = eshape1.a), 2538 namesof("scale" , lscale , earg = escale)), ", ", 2539 namesof("shape2.p" , lshape2.p, earg = eshape2.p), ", ", 2540 namesof("shape3.q" , lshape3.q, earg = eshape3.q), "\n", 2541 "Mean: scale * gamma(shape2.p + 1/shape1.a) * ", 2542 "gamma(shape3.q - 1/shape1.a) / ", 2543 "(gamma(shape2.p) * gamma(shape3.q))"), 2544 constraints = eval(substitute(expression({ 2545 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 2546 predictors.names = predictors.names, 2547 M1 = 4) 2548 }), list( .zero = zero ))), 2549 infos = eval(substitute(function(...) { 2550 list(M1 = 4, 2551 Q1 = 1, 2552 expected = TRUE, 2553 zero = .zero , 2554 multipleResponses = TRUE, 2555 parameters.names = if ( .lss ) 2556 c("scale", "shape1.a", "shape2.p", "shape3.q") else 2557 c("shape1.a", "scale", "shape2.p", "shape3.q"), 2558 lscale = .lscale , lshape1.a = .lshape1.a , 2559 escale = .escale , eshape1.a = .eshape1.a , 2560 lshape2.p = .lshape2.p , lshape3.q = .lshape3.q , 2561 eshape2.p = .eshape2.p , eshape3.q = .eshape3.q ) 2562 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 2563 .escale = escale , .eshape1.a = eshape1.a, 2564 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 2565 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q, 2566 .lss = lss , 2567 .zero = zero ))), 2568 initialize = eval(substitute(expression({ 2569 temp5 <- w.y.check(w = w, y = y, 2570 Is.positive.y = TRUE, 2571 ncol.w.max = Inf, 2572 ncol.y.max = Inf, 2573 out.wy = TRUE, 2574 colsyperw = 1, 2575 maximize = TRUE) 2576 y <- temp5$y 2577 w <- temp5$w 2578 M1 <- 4 # Number of parameters for one response 2579 NOS <- ncoly <- ncol(y) 2580 M <- M1*ncol(y) 2581 2582 scaL.names <- param.names("scale", NOS, skip1 = TRUE) 2583 sha1.names <- param.names("shape1.a", NOS, skip1 = TRUE) 2584 sha2.names <- param.names("shape2.p", NOS, skip1 = TRUE) 2585 sha3.names <- param.names("shape3.q", NOS, skip1 = TRUE) 2586 2587 predictors.names <- c( 2588 if ( .lss ) { 2589 c(namesof(scaL.names , .lscale , earg = .escale , 2590 tag = FALSE), 2591 namesof(sha1.names , .lshape1.a , earg = .eshape1.a , 2592 tag = FALSE)) 2593 } else { 2594 c(namesof(sha1.names , .lshape1.a , earg = .eshape1.a , 2595 tag = FALSE), 2596 namesof(scaL.names , .lscale , earg = .escale , 2597 tag = FALSE)) 2598 }, 2599 namesof(sha2.names , .lshape2.p , earg = .eshape2.p , tag = FALSE), 2600 namesof(sha3.names , .lshape3.q , earg = .eshape3.q , tag = FALSE)) 2601 predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)] 2602 2603 if (!length(etastart)) { 2604 sc.init <- 2605 aa.init <- 2606 pp.init <- 2607 qq.init <- matrix(NA_real_, n, NOS) 2608 2609 for (spp. in 1:NOS) { # For each response 'y_spp.'... do: 2610 yvec <- y[, spp.] 2611 wvec <- w[, spp.] 2612 2613 gscale <- .gscale 2614 gshape1.a <- .gshape1.a 2615 gshape2.p <- .gshape2.p 2616 gshape3.q <- .gshape3.q 2617 if (length( .iscale )) gscale <- rep_len( .iscale , NOS) 2618 if (length( .ishape1.a )) gshape1.a <- rep_len( .ishape1.a , NOS) 2619 if (length( .ishape2.p )) gshape2.p <- rep_len( .ishape2.p , NOS) 2620 if (length( .ishape3.q )) gshape3.q <- rep_len( .ishape3.q , NOS) 2621 2622 2623 try.this <- 2624 grid.search4(gscale, gshape1.a, gshape2.p, gshape3.q, 2625 objfun = genbetaII.Loglikfun4, # x = x, 2626 y = yvec, w = wvec, # extraargs = NULL, 2627 ret.objfun = TRUE) # Last value is the loglik 2628 sc.init[, spp.] <- try.this["Value1"] 2629 aa.init[, spp.] <- try.this["Value2"] 2630 pp.init[, spp.] <- try.this["Value3"] 2631 qq.init[, spp.] <- try.this["Value4"] 2632 } # End of for (spp. ...) 2633 2634 2635 finite.mean <- 1 < aa.init * qq.init 2636 COP.use <- 1.15 2637 while (FALSE && any(!finite.mean)) { 2638 qq.init[!finite.mean] <- 0.1 + qq.init[!finite.mean] * COP.use 2639 aa.init[!finite.mean] <- 0.1 + aa.init[!finite.mean] * COP.use 2640 finite.mean <- 1 < aa.init * qq.init 2641 } 2642 2643 etastart <- 2644 cbind(if ( .lss ) 2645 cbind(theta2eta(sc.init, .lscale , earg = .escale ), 2646 theta2eta(aa.init, .lshape1.a , earg = .eshape1.a )) else 2647 cbind(theta2eta(aa.init, .lshape1.a , earg = .eshape1.a ), 2648 theta2eta(sc.init, .lscale , earg = .escale )), 2649 theta2eta(pp.init , .lshape2.p , earg = .eshape2.p ), 2650 theta2eta(qq.init , .lshape3.q , earg = .eshape3.q )) 2651 etastart <- etastart[, interleave.VGAM(M, M1 = M1)] 2652 } # End of etastart. 2653 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 2654 .escale = escale , .eshape1.a = eshape1.a, 2655 .iscale = iscale , .ishape1.a = ishape1.a, 2656 .gscale = gscale , .gshape1.a = gshape1.a, 2657 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 2658 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q, 2659 .ishape2.p = ishape2.p, .ishape3.q = ishape3.q, 2660 .gshape2.p = gshape2.p, .gshape3.q = gshape3.q, 2661 .lss = lss ))), 2662 linkinv = eval(substitute(function(eta, extra = NULL) { 2663 M1 <- 4 2664 NOS <- ncol(eta) / M1 2665 if ( .lss ) { 2666 Scale <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE], 2667 .lscale , earg = .escale ) 2668 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 2669 .lshape1.a , earg = .eshape1.a ) 2670 } else { 2671 aa <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE], 2672 .lshape1.a , earg = .eshape1.a ) 2673 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 2674 .lscale , earg = .escale ) 2675 } 2676 parg <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 2677 .lshape2.p , earg = .eshape2.p ) 2678 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 2679 .lshape3.q , earg = .eshape3.q ) 2680 2681 ans <- cbind(Scale * exp(lgamma(parg + 1/aa) + 2682 lgamma(qq - 1/aa) - 2683 lgamma(parg) - lgamma(qq))) 2684 ans 2685 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 2686 .escale = escale , .eshape1.a = eshape1.a, 2687 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 2688 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q, 2689 .lss = lss ))), 2690 last = eval(substitute(expression({ 2691 M1 <- 4 2692 2693 misc$link <- c(rep_len( if ( .lss ) .lscale else .lshape1.a , ncoly), 2694 rep_len( if ( .lss ) .lshape1.a else .lscale , ncoly), 2695 rep_len( .lshape2.p , ncoly), 2696 rep_len( .lshape3.q , ncoly))[ 2697 interleave.VGAM(M, M1 = M1)] 2698 temp.names <- if ( .lss ) { 2699 c(scaL.names, sha1.names, sha2.names, sha3.names) 2700 } else { 2701 c(sha1.names, scaL.names, sha2.names, sha3.names) 2702 } 2703 names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)] 2704 2705 misc$earg <- vector("list", M) 2706 names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)] 2707 for (ii in 1:ncoly) { 2708 if ( .lss ) { 2709 misc$earg[[M1*ii-3]] <- .escale 2710 misc$earg[[M1*ii-2]] <- .eshape1.a 2711 } else { 2712 misc$earg[[M1*ii-3]] <- .eshape1.a 2713 misc$earg[[M1*ii-2]] <- .escale 2714 } 2715 misc$earg[[M1*ii-1]] <- .eshape2.p 2716 misc$earg[[M1*ii ]] <- .eshape3.q 2717 } 2718 2719 misc$expected <- TRUE 2720 misc$multipleResponses <- TRUE 2721 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 2722 .escale = escale , .eshape1.a = eshape1.a, 2723 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 2724 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q, 2725 .lss = lss ))), 2726 loglikelihood = eval(substitute( 2727 function(mu, y, w, residuals = FALSE, 2728 eta, extra = NULL, summation = TRUE) { 2729 M1 <- 4 2730 NOS <- ncol(eta)/M1 2731 if ( .lss ) { 2732 Scale <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE], 2733 .lscale , earg = .escale ) 2734 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 2735 .lshape1.a , earg = .eshape1.a ) 2736 } else { 2737 aa <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE], 2738 .lshape1.a , earg = .eshape1.a ) 2739 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 2740 .lscale , earg = .escale ) 2741 } 2742 parg <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 2743 .lshape2.p , earg = .eshape2.p ) 2744 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 2745 .lshape3.q , earg = .eshape3.q ) 2746 if (residuals) { 2747 stop("loglikelihood residuals not implemented yet") 2748 } else { 2749 ll.elts <- 2750 c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = aa, 2751 shape2.p = parg, shape3.q = qq, log = TRUE) 2752 if (summation) { 2753 sum(ll.elts) 2754 } else { 2755 ll.elts 2756 } 2757 } 2758 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 2759 .escale = escale , .eshape1.a = eshape1.a, 2760 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 2761 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q, 2762 .lss = lss ))), 2763 vfamily = c("genbetaII"), 2764 2765 2766 2767 2768 validparams = eval(substitute(function(eta, y, extra = NULL) { 2769 M1 <- 4 2770 NOS <- ncol(eta) / M1 2771 if ( .lss ) { 2772 Scale <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE], 2773 .lscale , earg = .escale ) 2774 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 2775 .lshape1.a , earg = .eshape1.a) 2776 } else { 2777 aa <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE], 2778 .lshape1.a , earg = .eshape1.a) 2779 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 2780 .lscale , earg = .escale ) 2781 } 2782 parg <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 2783 .lshape2.p , earg = .eshape2.p) 2784 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 2785 .lshape3.q , earg = .eshape3.q) 2786 2787 2788 okay1 <- all(is.finite(Scale)) && all(Scale > 0) && 2789 all(is.finite(aa )) && all(aa > 0) && 2790 all(is.finite(parg )) && all(parg > 0) && 2791 all(is.finite(qq )) && all(qq > 0) 2792 okay.support <- if (okay1) all(-aa * parg < 1 & 1 < aa * qq) else TRUE 2793 if (!okay.support) 2794 warning("parameter constraint -a*p < 1 < a*q has been violated; ", 2795 "solution may be at the boundary of the ", 2796 "parameter space.") 2797 okay1 && okay.support 2798 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 2799 .escale = escale , .eshape1.a = eshape1.a, 2800 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 2801 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q, 2802 .lss = lss ))), 2803 2804 2805 2806 2807 deriv = eval(substitute(expression({ 2808 M1 <- 4 2809 NOS <- ncol(eta)/M1 # Needed for summary() 2810 if ( .lss ) { 2811 Scale <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE], 2812 .lscale , earg = .escale ) 2813 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 2814 .lshape1.a , earg = .eshape1.a) 2815 } else { 2816 aa <- eta2theta(eta[, M1*(1:NOS) - 3, drop = FALSE], 2817 .lshape1.a , earg = .eshape1.a) 2818 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 2819 .lscale , earg = .escale ) 2820 } 2821 parg <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 2822 .lshape2.p , earg = .eshape2.p) 2823 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 2824 .lshape3.q , earg = .eshape3.q) 2825 temp1 <- log(y/Scale) 2826 temp2 <- (y/Scale)^aa 2827 temp3 <- digamma(parg + qq) 2828 temp3a <- digamma(parg) 2829 temp3b <- digamma(qq) 2830 temp4 <- log1p(temp2) 2831 2832 dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2)) 2833 dl.da <- 1/aa + parg * temp1 - (parg+qq) * temp1 / (1+1/temp2) 2834 dl.dp <- aa * temp1 + temp3 - temp3a - temp4 2835 dl.dq <- temp3 - temp3b - temp4 2836 2837 dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale ) 2838 da.deta <- dtheta.deta(aa, .lshape1.a , earg = .eshape1.a ) 2839 dp.deta <- dtheta.deta(parg, .lshape2.p , earg = .eshape2.p ) 2840 dq.deta <- dtheta.deta(qq, .lshape3.q , earg = .eshape3.q ) 2841 2842 myderiv <- if ( .lss ) { 2843 c(w) * cbind(dl.dscale * dscale.deta, 2844 dl.da * da.deta, 2845 dl.dp * dp.deta, 2846 dl.dq * dq.deta) 2847 } else { 2848 c(w) * cbind(dl.da * da.deta, 2849 dl.dscale * dscale.deta, 2850 dl.dp * dp.deta, 2851 dl.dq * dq.deta) 2852 } 2853 myderiv[, interleave.VGAM(M, M1 = M1)] 2854 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 2855 .escale = escale , .eshape1.a = eshape1.a, 2856 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 2857 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q, 2858 .lss = lss ))), 2859 weight = eval(substitute(expression({ 2860 temp5 <- trigamma(parg + qq) 2861 temp5a <- trigamma(parg) 2862 temp5b <- trigamma(qq) 2863 2864 ned2l.da <- (1 + parg + qq + parg * qq * (temp5a + temp5b + 2865 (temp3b - temp3a + (parg-qq)/(parg*qq))^2 - 2866 (parg^2 + qq^2) / (parg*qq)^2)) / (aa^2 * (1+parg+qq)) 2867 ned2l.dscale <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2) 2868 ned2l.dp <- temp5a - temp5 2869 ned2l.dq <- temp5b - temp5 2870 ned2l.dascale <- (parg - qq - parg * qq * 2871 (temp3a -temp3b)) / (Scale*(1 + parg+qq)) 2872 ned2l.dap <- -(qq * (temp3a -temp3b) -1) / (aa*(parg+qq)) 2873 ned2l.daq <- -(parg * (temp3b -temp3a) -1) / (aa*(parg+qq)) 2874 ned2l.dscalep <- aa * qq / (Scale*(parg+qq)) 2875 ned2l.dscaleq <- -aa * parg / (Scale*(parg+qq)) 2876 ned2l.dpq <- -temp5 2877 wz <- if ( .lss ) { 2878 array(c(c(w) * ned2l.dscale * dscale.deta^2, 2879 c(w) * ned2l.da * da.deta^2, 2880 c(w) * ned2l.dp * dp.deta^2, 2881 c(w) * ned2l.dq * dq.deta^2, 2882 c(w) * ned2l.dascale * da.deta * dscale.deta, 2883 c(w) * ned2l.dap * da.deta * dp.deta, 2884 c(w) * ned2l.dpq * dp.deta * dq.deta, 2885 c(w) * ned2l.dscalep * dscale.deta * dp.deta, 2886 c(w) * ned2l.daq * da.deta * dq.deta, 2887 c(w) * ned2l.dscaleq * dscale.deta * dq.deta), 2888 dim = c(n, M/M1, M1*(M1+1)/2)) 2889 } else { 2890 array(c(c(w) * ned2l.da * da.deta^2, 2891 c(w) * ned2l.dscale * dscale.deta^2, 2892 c(w) * ned2l.dp * dp.deta^2, 2893 c(w) * ned2l.dq * dq.deta^2, 2894 c(w) * ned2l.dascale * da.deta * dscale.deta, 2895 c(w) * ned2l.dscalep * dscale.deta * dp.deta, 2896 c(w) * ned2l.dpq * dp.deta * dq.deta, 2897 c(w) * ned2l.dap * da.deta * dp.deta, 2898 c(w) * ned2l.dscaleq * dscale.deta * dq.deta, 2899 c(w) * ned2l.daq * da.deta * dq.deta), 2900 dim = c(n, M/M1, M1*(M1+1)/2)) 2901 } 2902 wz <- arwz2wz(wz, M = M, M1 = M1) 2903 wz 2904 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 2905 .escale = escale , .eshape1.a = eshape1.a, 2906 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 2907 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q, 2908 .lss = lss )))) 2909} 2910 2911 2912 2913 2914 2915 2916 2917 2918 2919 2920dgenbetaII <- function(x, scale = 1, shape1.a, shape2.p, shape3.q, 2921 log = FALSE) { 2922 2923 if (!is.logical(log.arg <- log) || length(log) != 1) 2924 stop("Bad input for argument 'log'") 2925 rm(log) 2926 2927 2928 logden <- log(shape1.a) + (shape1.a * shape2.p - 1) * log(abs(x)) - 2929 shape1.a * shape2.p * log(scale) - 2930 lbeta(shape2.p, shape3.q) - 2931 (shape2.p + shape3.q) * log1p((abs(x)/scale)^shape1.a) 2932 2933 2934 if (any(x <= 0) || any(is.infinite(x))) { 2935 LLL <- max(length(x), length(scale), 2936 length(shape1.a), length(shape2.p), length(shape3.q)) 2937 if (length(x) != LLL) x <- rep_len(x, LLL) 2938 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 2939 if (length(shape1.a) != LLL) shape1.a <- rep_len(shape1.a, LLL) 2940 if (length(shape2.p) != LLL) shape2.p <- rep_len(shape2.p, LLL) 2941 if (length(shape3.q) != LLL) shape3.q <- rep_len(shape3.q, LLL) 2942 2943 logden[is.infinite(x)] <- log(0) 2944 logden[x < 0] <- log(0) 2945 x.eq.0 <- !is.na(x) & (x == 0) 2946 if (any(x.eq.0)) { 2947 axp <- shape1.a[x.eq.0] * shape2.p[x.eq.0] 2948 logden[x.eq.0 & axp < 1] <- log(Inf) 2949 ind5 <- x.eq.0 & axp == 1 2950 logden[ind5] <- log(shape1.a[ind5]) - 2951 shape1.a[ind5] * shape2.p[ind5] * log(scale[ind5]) - 2952 lbeta(shape2.p[ind5], shape3.q[ind5]) - 2953 (shape2.p[ind5] + shape3.q[ind5]) * 2954 log1p((0/scale[ind5])^shape1.a[ind5]) 2955 logden[x.eq.0 & axp > 1] <- log(0) 2956 } 2957 } 2958 2959 if (log.arg) logden else exp(logden) 2960} 2961 2962 2963 2964 2965 2966rsinmad <- function(n, scale = 1, shape1.a, shape3.q) 2967 qsinmad(runif(n), shape1.a = shape1.a, scale = scale, 2968 shape3.q = shape3.q) 2969 2970 2971rlomax <- function(n, scale = 1, shape3.q) 2972 rsinmad(n, scale = scale, shape1.a = 1, shape3.q = shape3.q) 2973 2974 2975rfisk <- function(n, scale = 1, shape1.a) 2976 rsinmad(n, scale = scale, shape1.a = shape1.a, shape3.q = 1) 2977 2978 2979rparalogistic <- function(n, scale = 1, shape1.a) 2980 rsinmad(n, scale = scale, shape1.a = shape1.a, shape3.q = shape1.a) 2981 2982 2983rdagum <- function(n, scale = 1, shape1.a, shape2.p) 2984 qdagum(runif(n), scale = scale, shape1.a = shape1.a, 2985 shape2.p = shape2.p) 2986 2987 2988rinv.lomax <- function(n, scale = 1, shape2.p) 2989 rdagum(n, scale = scale, shape1.a = 1, shape2.p = shape2.p) 2990 2991 2992rinv.paralogistic <- function(n, scale = 1, shape1.a) 2993 rdagum(n, scale = scale, shape1.a = shape1.a, shape2.p = shape1.a) 2994 2995 2996 2997 2998qsinmad <- function(p, scale = 1, shape1.a, shape3.q, 2999 lower.tail = TRUE, 3000 log.p = FALSE) { 3001 3002 3003 3004 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 3005 stop("bad input for argument 'lower.tail'") 3006 3007 if (!is.logical(log.p) || length(log.p) != 1) 3008 stop("bad input for argument 'log.p'") 3009 3010 3011 3012 LLL <- max(length(p), length(shape1.a), length(scale), 3013 length(shape3.q)) 3014 if (length(p) != LLL) p <- rep_len(p, LLL) 3015 if (length(shape1.a) != LLL) shape1.a <- rep_len(shape1.a, LLL) 3016 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 3017 if (length(shape3.q) != LLL) shape3.q <- rep_len(shape3.q, LLL) 3018 3019 3020 if (lower.tail) { 3021 if (log.p) { 3022 ln.p <- p 3023 ans <- scale * expm1((-1/shape3.q) * log(-expm1(ln.p)))^(1/shape1.a) 3024 } else { 3025 ans <- scale * expm1((-1/shape3.q) * log1p(-p))^(1/shape1.a) 3026 ans[p == 0] <- 0 3027 ans[p == 1] <- Inf 3028 } 3029 } else { 3030 if (log.p) { 3031 ln.p <- p 3032 ans <- scale * expm1(-ln.p / shape3.q)^(1/shape1.a) 3033 } else { 3034 ans <- scale * expm1(-log(p) / shape3.q)^(1/shape1.a) 3035 ans[p == 0] <- Inf 3036 ans[p == 1] <- 0 3037 } 3038 } 3039 3040 ans[scale <= 0 | shape1.a <= 0 | shape3.q <= 0] <- NaN 3041 ans 3042} 3043 3044 3045 3046 3047qlomax <- function(p, scale = 1, shape3.q, 3048 lower.tail = TRUE, log.p = FALSE) 3049 qsinmad(p, shape1.a = 1, scale = scale, shape3.q = shape3.q, 3050 lower.tail = lower.tail, log.p = log.p) 3051 3052qfisk <- function(p, scale = 1, shape1.a, 3053 lower.tail = TRUE, log.p = FALSE) 3054 qsinmad(p, shape1.a = shape1.a, scale = scale, shape3.q = 1, 3055 lower.tail = lower.tail, log.p = log.p) 3056 3057qparalogistic <- function(p, scale = 1, shape1.a, 3058 lower.tail = TRUE, log.p = FALSE) 3059 qsinmad(p, shape1.a = shape1.a, scale = scale, 3060 shape3.q = shape1.a, ## 20150121 KaiH; add shape3.q = shape1.a 3061 lower.tail = lower.tail, log.p = log.p) 3062 3063 3064 3065 3066 3067qdagum <- function(p, scale = 1, shape1.a, shape2.p, 3068 lower.tail = TRUE, log.p = FALSE) { 3069 3070 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 3071 stop("bad input for argument 'lower.tail'") 3072 if (!is.logical(log.p) || length(log.p) != 1) 3073 stop("bad input for argument 'log.p'") 3074 3075 3076 LLL <- max(length(p), length(shape1.a), length(scale), 3077 length(shape2.p)) 3078 if (length(p) != LLL) p <- rep_len(p, LLL) 3079 if (length(shape1.a) != LLL) shape1.a <- rep_len(shape1.a, LLL) 3080 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 3081 if (length(shape2.p) != LLL) shape2.p <- rep_len(shape2.p, LLL) 3082 3083 if (lower.tail) { 3084 if (log.p) { 3085 ln.p <- p 3086 ans <- scale * (expm1(-ln.p/shape2.p))^(-1/shape1.a) 3087 ans[ln.p > 0] <- NaN 3088 } else { 3089 ans <- scale * (expm1(-log(p)/shape2.p))^(-1/shape1.a) 3090 ans[p < 0] <- NaN 3091 ans[p == 0] <- 0 3092 ans[p == 1] <- Inf 3093 ans[p > 1] <- NaN 3094 } 3095 } else { 3096 if (log.p) { 3097 ln.p <- p 3098 ans <- scale * (expm1(-log(-expm1(ln.p))/shape2.p))^(-1/shape1.a) 3099 ans[ln.p > 0] <- NaN 3100 } else { 3101 ans <- scale * (expm1(-log1p(-p)/shape2.p))^(-1/shape1.a) 3102 ans[p < 0] <- NaN 3103 ans[p == 0] <- Inf 3104 ans[p == 1] <- 0 3105 ans[p > 1] <- NaN 3106 } 3107 } 3108 3109 ans[scale <= 0 | shape1.a <= 0 | shape2.p <= 0] <- NaN 3110 ans 3111} 3112 3113 3114 3115 3116qinv.lomax <- function(p, scale = 1, shape2.p, 3117 lower.tail = TRUE, log.p = FALSE) 3118 qdagum(p, scale = scale, shape1.a = 1, shape2.p = shape2.p, 3119 lower.tail = lower.tail, log.p = log.p) 3120 3121 3122qinv.paralogistic <- function(p, scale = 1, shape1.a, 3123 lower.tail = TRUE, log.p = FALSE) 3124 qdagum(p, scale = scale, shape1.a = shape1.a, 3125 shape2.p = shape1.a, ## 20150121 Kai; add shape2.p = shape1.a 3126 lower.tail = lower.tail, log.p = log.p) 3127 3128 3129 3130 3131psinmad <- function(q, scale = 1, shape1.a, shape3.q, 3132 lower.tail = TRUE, log.p = FALSE) { 3133 3134 3135 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 3136 stop("bad input for argument 'lower.tail'") 3137 3138 if (!is.logical(log.p) || length(log.p) != 1) 3139 stop("bad input for argument 'log.p'") 3140 3141 3142 LLL <- max(length(q), length(shape1.a), length(scale), 3143 length(shape3.q)) 3144 if (length(q) != LLL) q <- rep_len(q, LLL) 3145 if (length(shape1.a) != LLL) shape1.a <- rep_len(shape1.a, LLL) 3146 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 3147 if (length(shape3.q) != LLL) shape3.q <- rep_len(shape3.q, LLL) 3148 3149 # 20150121 KaiH 3150 if (lower.tail) { 3151 if (log.p) { 3152 ans <- log1p(-(1 + (q / scale)^shape1.a)^(-shape3.q)) 3153 ans[q <= 0 ] <- -Inf 3154 ans[q == Inf] <- 0 3155 } else { 3156 ans <- exp(log1p(-(1 + (q / scale)^shape1.a)^(-shape3.q))) 3157 ans[q <= 0] <- 0 3158 ans[q == Inf] <- 1 3159 } 3160 } else { 3161 if (log.p) { 3162 ans <- (-shape3.q) * log1p((q / scale)^shape1.a) 3163 ans[q <= 0] <- 0 3164 ans[q == Inf] <- -Inf 3165 } else { 3166 ans <- (1 + (q / scale)^shape1.a)^(-shape3.q) 3167 ans[q <= 0] <- 1 3168 ans[q == Inf] <- 0 3169 } 3170 } 3171 ans[scale <= 0 | shape1.a <= 0 | shape3.q <= 0] <- NaN 3172 ans 3173} 3174 3175 3176 3177 3178 3179 3180plomax <- function(q, scale = 1, shape3.q, # Change the order 3181 lower.tail = TRUE, log.p = FALSE) 3182 psinmad(q, shape1.a = 1, scale = scale, shape3.q = shape3.q, 3183 lower.tail = lower.tail, log.p = log.p) 3184 3185 3186pfisk <- function(q, scale = 1, shape1.a, 3187 lower.tail = TRUE, log.p = FALSE) 3188 psinmad(q, shape1.a = shape1.a, scale = scale, shape3.q = 1, 3189 lower.tail = lower.tail, log.p = log.p) 3190 3191 3192pparalogistic <- function(q, scale = 1, shape1.a, # Change the order 3193 lower.tail = TRUE, log.p = FALSE) 3194 psinmad(q, shape1.a = shape1.a, scale = scale, 3195 shape3.q = shape1.a, # Add shape3.q = shape1.a 3196 lower.tail = lower.tail, log.p = log.p) 3197 3198 3199 3200 3201 3202 3203pdagum <- function(q, scale = 1, shape1.a, shape2.p, 3204 lower.tail = TRUE, log.p = FALSE) { 3205 3206 3207 3208 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 3209 stop("bad input for argument 'lower.tail'") 3210 3211 if (!is.logical(log.p) || length(log.p) != 1) 3212 stop("bad input for argument 'log.p'") 3213 3214 3215 LLL <- max(length(q), length(shape1.a), length(scale), 3216 length(shape2.p)) 3217 if (length(q) != LLL) q <- rep_len(q, LLL) 3218 if (length(shape1.a) != LLL) shape1.a <- rep_len(shape1.a, LLL) 3219 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 3220 if (length(shape2.p) != LLL) shape2.p <- rep_len(shape2.p, LLL) 3221 3222 3223 if (lower.tail) { 3224 if (log.p) { 3225 ans <- (-shape2.p) * log1p((q/scale)^(-shape1.a)) 3226 ans[q <= 0 ] <- -Inf 3227 ans[q == Inf] <- 0 3228 } else { 3229 ans <- exp( (-shape2.p) * log1p((q/scale)^(-shape1.a)) ) 3230 ans[q <= 0] <- 0 3231 ans[q == Inf] <- 1 3232 } 3233 } else { 3234 if (log.p) { 3235 ans <- log1p(-(1 + (q/scale)^(-shape1.a))^(-shape2.p)) 3236 ans[q <= 0] <- 0 3237 ans[q == Inf] <- -Inf 3238 } else { 3239 stop("unfinished") 3240 ans[q <= 0] <- 1 3241 ans[q == Inf] <- 0 3242 } 3243 } 3244 3245 ans[shape1.a <= 0 | scale <= 0 | shape2.p <= 0] <- NaN 3246 ans 3247} 3248 3249 3250 3251 3252 3253pinv.lomax <- function(q, scale = 1, shape2.p, 3254 lower.tail = TRUE, log.p = FALSE) 3255 pdagum(q, scale = scale, shape1.a = 1, shape2.p = shape2.p, 3256 lower.tail = lower.tail, log.p = log.p) 3257 3258 3259pinv.paralogistic <- function(q, scale = 1, shape1.a, 3260 lower.tail = TRUE, log.p = FALSE) 3261 pdagum(q, scale = scale, shape1.a = shape1.a, 3262 shape2.p = shape1.a, 3263 lower.tail = lower.tail, log.p = log.p) 3264 3265 3266 3267 3268 3269dbetaII <- function(x, scale = 1, shape2.p, shape3.q, log = FALSE) 3270 dgenbetaII(x = x, scale = scale, shape1.a = 1, 3271 shape2.p = shape2.p, shape3.q = shape3.q, log = log) 3272 3273 3274 3275 3276dsinmad <- function(x, scale = 1, shape1.a, shape3.q, log = FALSE) { 3277 3278 if (!is.logical(log.arg <- log) || length(log) != 1) 3279 stop("bad input for argument 'log'") 3280 rm(log) 3281 3282 LLL <- max(length(x), length(shape1.a), 3283 length(scale), length(shape3.q)) 3284 if (length(x) != LLL) x <- rep_len(x, LLL) 3285 if (length(shape1.a) != LLL) shape1.a <- rep_len(shape1.a, LLL) 3286 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 3287 if (length(shape3.q) != LLL) shape3.q <- rep_len(shape3.q, LLL) 3288 3289 Loglik <- rep_len(log(0), LLL) 3290 xok <- (x > 0) & !is.na(x) # Avoids log(x) if x<0, and handles NAs 3291 Loglik[xok] <- log(shape1.a[xok]) + log(shape3.q[xok]) + 3292 (shape1.a[xok]-1) * log(x[xok]) - 3293 shape1.a[xok] * log(scale[xok]) - 3294 (1 + shape3.q[xok]) * log1p((x[xok]/scale[xok])^shape1.a[xok]) 3295 x.eq.0 <- (x == 0) & !is.na(x) 3296 Loglik[x.eq.0] <- log(shape1.a[x.eq.0]) + log(shape3.q[x.eq.0]) - 3297 shape1.a[x.eq.0] * log(scale[x.eq.0]) 3298 Loglik[is.na(x)] <- NA 3299 Loglik[is.nan(x)] <- NaN 3300 Loglik[x == Inf] <- log(0) 3301 3302 if (log.arg) Loglik else exp(Loglik) 3303} 3304 3305 3306 3307dlomax <- function(x, scale = 1, shape3.q, log = FALSE) 3308 dsinmad(x, scale = scale, shape1.a = 1, shape3.q = shape3.q, log = log) 3309 3310 3311 3312dfisk <- function(x, scale = 1, shape1.a, log = FALSE) 3313 dsinmad(x, scale = scale, shape1.a = shape1.a, shape3.q = 1, log = log) 3314 3315 3316 3317dparalogistic <- function(x, scale = 1, shape1.a, log = FALSE) 3318 dsinmad(x, scale = scale, shape1.a = shape1.a, shape3.q = shape1.a, 3319 log = log) 3320 3321 3322 3323ddagum <- function(x, scale = 1, shape1.a, shape2.p, log = FALSE) { 3324 if (!is.logical(log.arg <- log) || length(log) != 1) 3325 stop("bad input for argument 'log'") 3326 rm(log) 3327 3328 LLL <- max(length(x), 3329 length(shape1.a), 3330 length(scale), 3331 length(shape2.p)) 3332 if (length(x) != LLL) x <- rep_len(x, LLL) 3333 if (length(shape1.a) != LLL) shape1.a <- rep_len(shape1.a, LLL) 3334 if (length(scale) != LLL) scale <- rep_len(scale, LLL) 3335 if (length(shape2.p) != LLL) shape2.p <- rep_len(shape2.p, LLL) 3336 3337 Loglik <- rep_len(log(0), LLL) 3338 xok <- (x > 0) & !is.na(x) # Avoids log(x) if x<0, and handles NAs 3339 Loglik[xok] <- log(shape1.a[xok]) + 3340 log(shape2.p[xok]) + 3341 (shape1.a[xok] * shape2.p[xok]-1) * log( x[xok]) - 3342 shape1.a[xok] * shape2.p[xok] * log(scale[xok]) - 3343 (1 + shape2.p[xok]) * 3344 log1p((x[xok]/scale[xok])^shape1.a[xok]) 3345 Loglik[shape2.p <= 0] <- NaN 3346 3347 x.eq.0 <- (x == 0) & !is.na(x) 3348 Loglik[x.eq.0] <- log(shape1.a[x.eq.0]) + 3349 log(shape2.p[x.eq.0]) - 3350 shape1.a[x.eq.0] * shape2.p[x.eq.0] * 3351 log(scale[x.eq.0]) 3352 Loglik[is.na(x)] <- NA 3353 Loglik[is.nan(x)] <- NaN 3354 Loglik[x == Inf] <- log(0) 3355 3356 if (log.arg) Loglik else exp(Loglik) 3357} 3358 3359 3360 3361dinv.lomax <- function(x, scale = 1, shape2.p, log = FALSE) 3362 ddagum(x, scale = scale, shape1.a = 1, shape2.p = shape2.p, 3363 log = log) 3364 3365 3366dinv.paralogistic <- function(x, scale = 1, shape1.a, log = FALSE) 3367 ddagum(x, scale = scale, shape1.a = shape1.a, shape2.p = shape1.a, 3368 log = log) 3369 3370 3371 3372 3373 3374 3375 3376 3377 sinmad <- function(lscale = "loglink", 3378 lshape1.a = "loglink", 3379 lshape3.q = "loglink", 3380 iscale = NULL, 3381 ishape1.a = NULL, 3382 ishape3.q = NULL, 3383 imethod = 1, 3384 lss = TRUE, 3385 gscale = exp(-5:5), 3386 gshape1.a = exp(-5:5), 3387 gshape3.q = exp(-5:5), 3388 probs.y = c(0.25, 0.50, 0.75), 3389 zero = "shape") { 3390 3391 3392 3393 3394 if (length(lss) != 1 && !is.logical(lss)) 3395 stop("Argument 'lss' not specified correctly") 3396 3397 if (!is.Numeric(imethod, length.arg = 1, 3398 integer.valued = TRUE, 3399 positive = TRUE) || imethod > 2) 3400 stop("Bad input for argument 'imethod'") 3401 3402 if (length(iscale) && !is.Numeric(iscale, positive = TRUE)) 3403 stop("Bad input for argument 'iscale'") 3404 3405 if (length(ishape1.a) && !is.Numeric(ishape1.a, positive = TRUE)) 3406 stop("Bad input for argument 'ishape1.a'") 3407 3408 if (length(ishape3.q) && !is.Numeric(ishape3.q, positive = TRUE)) 3409 stop("Bad input for argument 'ishape3.q'") 3410 3411 if (length(probs.y) < 2 || max(probs.y) > 1 || 3412 !is.Numeric(probs.y, positive = TRUE)) 3413 stop("Bad input for argument 'probs.y'") 3414 3415 lscale <- as.list(substitute(lscale)) 3416 escale <- link2list(lscale) 3417 lscale <- attr(escale, "function.name") 3418 3419 lshape1.a <- as.list(substitute(lshape1.a)) 3420 eshape1.a <- link2list(lshape1.a) 3421 lshape1.a <- attr(eshape1.a, "function.name") 3422 3423 lshape3.q <- as.list(substitute(lshape3.q)) 3424 eshape3.q <- link2list(lshape3.q) 3425 lshape3.q <- attr(eshape3.q, "function.name") 3426 3427 3428 new("vglmff", 3429 blurb = 3430 c("Singh-Maddala distribution \n\n", 3431 "Links: ", 3432 ifelse (lss, 3433 namesof("scale" , lscale , earg = escale), 3434 namesof("shape1.a", lshape1.a, earg = eshape1.a)), ", ", 3435 ifelse (lss, 3436 namesof("shape1.a", lshape1.a, earg = eshape1.a), 3437 namesof("scale" , lscale , earg = escale)), ", ", 3438 namesof("shape3.q" , lshape3.q, earg = eshape3.q), "\n", 3439 "Mean: scale * gamma(shape2.p + 1/shape1.a) * ", 3440 "gamma(shape3.q - 1/shape1.a) / ", 3441 "gamma(shape3.q)"), 3442 constraints = eval(substitute(expression({ 3443 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 3444 predictors.names = predictors.names, 3445 M1 = 3) 3446 }), list( .zero = zero ))), 3447 infos = eval(substitute(function(...) { 3448 list(M1 = 3, 3449 Q1 = 1, 3450 expected = TRUE, 3451 zero = .zero , 3452 multipleResponses = TRUE, 3453 parameters.names = if ( .lss ) 3454 c("scale", "shape1.a", "shape3.q") else 3455 c("shape1.a", "scale", "shape3.q"), 3456 lscale = .lscale , lshape1.a = .lshape1.a , 3457 escale = .escale , eshape1.a = .eshape1.a , 3458 lshape3.q = .lshape3.q , 3459 eshape3.q = .eshape3.q ) 3460 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 3461 .escale = escale , .eshape1.a = eshape1.a, 3462 .lshape3.q = lshape3.q, 3463 .eshape3.q = eshape3.q, 3464 .lss = lss , 3465 .zero = zero ))), 3466 initialize = eval(substitute(expression({ 3467 temp5 <- w.y.check(w = w, y = y, 3468 Is.positive.y = TRUE, 3469 ncol.w.max = Inf, 3470 ncol.y.max = Inf, 3471 out.wy = TRUE, 3472 colsyperw = 1, 3473 maximize = TRUE) 3474 y <- temp5$y 3475 w <- temp5$w 3476 M1 <- 3 # Number of parameters for one response 3477 NOS <- ncoly <- ncol(y) 3478 M <- M1*ncol(y) 3479 3480 scaL.names <- param.names("scale", NOS, skip1 = TRUE) 3481 sha1.names <- param.names("shape1.a", NOS, skip1 = TRUE) 3482 sha3.names <- param.names("shape3.q", NOS, skip1 = TRUE) 3483 3484 predictors.names <- c( 3485 if ( .lss ) { 3486 c(namesof(scaL.names , .lscale , earg = .escale , 3487 tag = FALSE), 3488 namesof(sha1.names , .lshape1.a , earg = .eshape1.a , 3489 tag = FALSE)) 3490 } else { 3491 c(namesof(sha1.names , .lshape1.a , earg = .eshape1.a , 3492 tag = FALSE), 3493 namesof(scaL.names , .lscale , earg = .escale , 3494 tag = FALSE)) 3495 }, 3496 namesof(sha3.names , .lshape3.q , earg = .eshape3.q , tag = FALSE)) 3497 predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)] 3498 3499 if (!length(etastart)) { 3500 sc.init <- 3501 aa.init <- 3502 qq.init <- matrix(NA_real_, n, NOS) 3503 3504 for (spp. in 1:NOS) { # For each response 'y_spp.'... do: 3505 yvec <- y[, spp.] 3506 wvec <- w[, spp.] 3507 3508 if ( .imethod == 1 ) { 3509 gscale <- .gscale 3510 gshape1.a <- .gshape1.a 3511 gshape3.q <- .gshape3.q 3512 if (length( .iscale )) gscale <- rep_len( .iscale , NOS) 3513 if (length( .ishape1.a )) gshape1.a <- rep_len( .ishape1.a , NOS) 3514 if (length( .ishape3.q )) gshape3.q <- rep_len( .ishape3.q , NOS) 3515 3516 3517 3518 try.this <- 3519 grid.search4(gscale, gshape1.a, vov3 = 1, gshape3.q, 3520 objfun = genbetaII.Loglikfun4, # x = x, 3521 y = yvec, w = wvec, # extraargs = NULL, 3522 ret.objfun = TRUE) # Last value is the loglik 3523 sc.init[, spp.] <- try.this["Value1"] 3524 aa.init[, spp.] <- try.this["Value2"] 3525 qq.init[, spp.] <- try.this["Value4"] 3526 } else { # .imethod == 2 3527 qvec <- .probs.y 3528 ishape3.q <- if (length( .ishape3.q )) .ishape3.q else 1 3529 xvec <- log( (1-qvec)^(-1/ ishape3.q ) - 1 ) 3530 fit0 <- lsfit(x = xvec, y = log(quantile(yvec, qvec))) 3531 sc.init[, spp.] <- if (length( .iscale )) .iscale else 3532 exp(fit0$coef[1]) 3533 aa.init[, spp.] <- if (length( .ishape1.a )) .ishape1.a else 3534 1/fit0$coef[2] 3535 qq.init[, spp.] <- ishape3.q 3536 } 3537 } # End of for (spp. ...) 3538 3539 3540 3541 3542 finite.mean <- 1 < aa.init * qq.init 3543 COP.use <- 1.15 3544 while (FALSE && any(!finite.mean)) { 3545 qq.init[!finite.mean] <- 0.1 + qq.init[!finite.mean] * COP.use 3546 aa.init[!finite.mean] <- 0.1 + aa.init[!finite.mean] * COP.use 3547 finite.mean <- 1 < aa.init * qq.init 3548 } 3549 3550 etastart <- 3551 cbind(if ( .lss ) 3552 cbind(theta2eta(sc.init, .lscale , earg = .escale ), 3553 theta2eta(aa.init, .lshape1.a , earg = .eshape1.a )) else 3554 cbind(theta2eta(aa.init, .lshape1.a , earg = .eshape1.a ), 3555 theta2eta(sc.init, .lscale , earg = .escale )), 3556 theta2eta(qq.init , .lshape3.q , earg = .eshape3.q )) 3557 etastart <- etastart[, interleave.VGAM(M, M1 = M1)] 3558 } # End of etastart. 3559 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 3560 .escale = escale , .eshape1.a = eshape1.a, 3561 .iscale = iscale , .ishape1.a = ishape1.a, 3562 .gscale = gscale , .gshape1.a = gshape1.a, 3563 .lshape3.q = lshape3.q, 3564 .eshape3.q = eshape3.q, 3565 .ishape3.q = ishape3.q, 3566 .gshape3.q = gshape3.q, 3567 .imethod = imethod , .probs.y = probs.y, 3568 .lss = lss ))), 3569 linkinv = eval(substitute(function(eta, extra = NULL) { 3570 M1 <- 3 3571 NOS <- ncol(eta)/M1 3572 if ( .lss ) { 3573 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 3574 .lscale , earg = .escale ) 3575 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 3576 .lshape1.a , earg = .eshape1.a ) 3577 } else { 3578 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 3579 .lshape1.a , earg = .eshape1.a ) 3580 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 3581 .lscale , earg = .escale ) 3582 } 3583 parg <- 1 3584 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 3585 .lshape3.q , earg = .eshape3.q ) 3586 3587 ans <- Scale * exp(lgamma(parg + 1/aa) + 3588 lgamma(qq - 1/aa) - lgamma(parg) - lgamma(qq)) 3589 ans[parg + 1/aa <= 0] <- NA 3590 ans[qq - 1/aa <= 0] <- NA 3591 ans[aa <= 0] <- NA 3592 ans[Scale <= 0] <- NA 3593 ans[qq <= 0] <- NA 3594 ans 3595 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 3596 .escale = escale , .eshape1.a = eshape1.a, 3597 .lshape3.q = lshape3.q, 3598 .eshape3.q = eshape3.q, 3599 .lss = lss ))), 3600 last = eval(substitute(expression({ 3601 M1 <- 3 3602 3603 misc$link <- c(rep_len( if ( .lss ) .lscale else .lshape1.a , ncoly), 3604 rep_len( if ( .lss ) .lshape1.a else .lscale , ncoly), 3605 rep_len( .lshape3.q , ncoly))[ 3606 interleave.VGAM(M, M1 = M1)] 3607 temp.names <- if ( .lss ) { 3608 c(scaL.names, sha1.names, sha3.names) 3609 } else { 3610 c(sha1.names, scaL.names, sha3.names) 3611 } 3612 names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)] 3613 3614 misc$earg <- vector("list", M) 3615 names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)] 3616 for (ii in 1:ncoly) { 3617 if ( .lss ) { 3618 misc$earg[[M1*ii-2]] <- .escale 3619 misc$earg[[M1*ii-1]] <- .eshape1.a 3620 } else { 3621 misc$earg[[M1*ii-2]] <- .eshape1.a 3622 misc$earg[[M1*ii-1]] <- .escale 3623 } 3624 misc$earg[[M1*ii ]] <- .eshape3.q 3625 } 3626 3627 misc$expected <- TRUE 3628 misc$multipleResponses <- TRUE 3629 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 3630 .escale = escale , .eshape1.a = eshape1.a, 3631 .lshape3.q = lshape3.q, 3632 .eshape3.q = eshape3.q, 3633 .lss = lss ))), 3634 loglikelihood = eval(substitute( 3635 function(mu, y, w, residuals = FALSE, 3636 eta, extra = NULL, summation = TRUE) { 3637 M1 <- 3 3638 NOS <- ncol(eta)/M1 3639 if ( .lss ) { 3640 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 3641 .lscale , earg = .escale ) 3642 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 3643 .lshape1.a , earg = .eshape1.a ) 3644 } else { 3645 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 3646 .lshape1.a , earg = .eshape1.a ) 3647 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 3648 .lscale , earg = .escale ) 3649 } 3650 parg <- 1 3651 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 3652 .lshape3.q , earg = .eshape3.q ) 3653 if (residuals) { 3654 stop("loglikelihood residuals not implemented yet") 3655 } else { 3656 ll.elts <- 3657 c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = aa, 3658 shape2.p = parg, shape3.q = qq, log = TRUE) 3659 if (summation) { 3660 sum(ll.elts) 3661 } else { 3662 ll.elts 3663 } 3664 } 3665 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 3666 .escale = escale , .eshape1.a = eshape1.a, 3667 .lshape3.q = lshape3.q, 3668 .eshape3.q = eshape3.q, 3669 .lss = lss ))), 3670 vfamily = c("sinmad"), 3671 simslot = eval(substitute( 3672 function(object, nsim) { 3673 3674 pwts <- if (length(pwts <- object@prior.weights) > 0) 3675 pwts else weights(object, type = "prior") 3676 if (any(pwts != 1)) 3677 warning("ignoring prior weights") 3678 3679 eta <- predict(object) 3680 M1 <- 3 3681 NOS <- ncol(eta)/M1 3682 if ( .lss ) { 3683 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 3684 .lscale , earg = .escale ) 3685 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 3686 .lshape1.a , earg = .eshape1.a ) 3687 } else { 3688 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 3689 .lshape1.a , earg = .eshape1.a ) 3690 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 3691 .lscale , earg = .escale ) 3692 } 3693 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 3694 .lshape3.q , earg = .eshape3.q ) 3695 3696 rsinmad(nsim * length(qq), shape1.a = aa, scale = Scale, 3697 shape3.q = qq) 3698 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 3699 .escale = escale , .eshape1.a = eshape1.a, 3700 .lshape3.q = lshape3.q, 3701 .eshape3.q = eshape3.q 3702 ))), 3703 3704 validparams = eval(substitute(function(eta, y, extra = NULL) { 3705 M1 <- 3 3706 NOS <- ncol(eta) / M1 3707 if ( .lss ) { 3708 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 3709 .lscale , earg = .escale ) 3710 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 3711 .lshape1.a , earg = .eshape1.a) 3712 } else { 3713 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 3714 .lshape1.a , earg = .eshape1.a) 3715 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 3716 .lscale , earg = .escale ) 3717 } 3718 parg <- 1 3719 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 3720 .lshape3.q , earg = .eshape3.q) 3721 3722 okay1 <- all(is.finite(Scale)) && all(Scale > 0) && 3723 all(is.finite(aa )) && all(aa > 0) && 3724 all(is.finite(parg )) && all(parg > 0) && 3725 all(is.finite(qq )) && all(qq > 0) 3726 okay.support <- if (okay1) all(-aa * parg < 1 & 1 < aa * qq) else TRUE 3727 if (!okay.support) 3728 warning("parameter constraint -a < 1 < a*q has been violated; ", 3729 "solution may be at the boundary of the ", 3730 "parameter space.") 3731 okay1 && okay.support 3732 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 3733 .escale = escale , .eshape1.a = eshape1.a, 3734 .lshape3.q = lshape3.q, 3735 .eshape3.q = eshape3.q, 3736 .lss = lss ))), 3737 3738 3739 3740 3741 deriv = eval(substitute(expression({ 3742 M1 <- 3 3743 NOS <- ncol(eta)/M1 # Needed for summary() 3744 if ( .lss ) { 3745 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 3746 .lscale , earg = .escale ) 3747 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 3748 .lshape1.a , earg = .eshape1.a) 3749 } else { 3750 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 3751 .lshape1.a , earg = .eshape1.a) 3752 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 3753 .lscale , earg = .escale ) 3754 } 3755 parg <- 1 3756 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 3757 .lshape3.q , earg = .eshape3.q) 3758 temp1 <- log(y/Scale) 3759 temp2 <- (y/Scale)^aa 3760 temp3 <- digamma(parg + qq) 3761 temp3a <- digamma(parg) 3762 temp3b <- digamma(qq) 3763 temp4 <- log1p(temp2) 3764 3765 dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2)) 3766 dl.da <- 1/aa + parg * temp1 - (parg+qq) * temp1 / (1+1/temp2) 3767 dl.dq <- temp3 - temp3b - temp4 3768 3769 dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale ) 3770 da.deta <- dtheta.deta(aa, .lshape1.a , earg = .eshape1.a ) 3771 dq.deta <- dtheta.deta(qq, .lshape3.q , earg = .eshape3.q ) 3772 3773 myderiv <- if ( .lss ) { 3774 c(w) * cbind(dl.dscale * dscale.deta, 3775 dl.da * da.deta, 3776 dl.dq * dq.deta) 3777 } else { 3778 c(w) * cbind(dl.da * da.deta, 3779 dl.dscale * dscale.deta, 3780 dl.dq * dq.deta) 3781 } 3782 myderiv[, interleave.VGAM(M, M1 = M1)] 3783 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 3784 .escale = escale , .eshape1.a = eshape1.a, 3785 .lshape3.q = lshape3.q, 3786 .eshape3.q = eshape3.q, 3787 .lss = lss ))), 3788 weight = eval(substitute(expression({ 3789 temp5 <- trigamma(parg + qq) 3790 temp5a <- trigamma(parg) 3791 temp5b <- trigamma(qq) 3792 3793 ned2l.da <- (1 + parg + qq + parg * qq * (temp5a + temp5b + 3794 (temp3b - temp3a + (parg-qq)/(parg*qq))^2 - 3795 (parg^2 + qq^2) / (parg*qq)^2)) / (aa^2 * (1+parg+qq)) 3796 ned2l.dscale <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2) 3797 ned2l.dq <- temp5b - temp5 3798 ned2l.dascale <- (parg - qq - parg * qq * 3799 (temp3a -temp3b)) / (Scale*(1 + parg+qq)) 3800 ned2l.daq <- -(parg * (temp3b -temp3a) -1) / (aa*(parg+qq)) 3801 ned2l.dscaleq <- -aa * parg / (Scale*(parg+qq)) 3802 wz <- if ( .lss ) { 3803 array(c(c(w) * ned2l.dscale * dscale.deta^2, 3804 c(w) * ned2l.da * da.deta^2, 3805 c(w) * ned2l.dq * dq.deta^2, 3806 c(w) * ned2l.dascale * da.deta * dscale.deta, 3807 c(w) * ned2l.daq * da.deta * dq.deta, 3808 c(w) * ned2l.dscaleq * dscale.deta * dq.deta), 3809 dim = c(n, M/M1, M1*(M1+1)/2)) 3810 } else { 3811 array(c(c(w) * ned2l.da * da.deta^2, 3812 c(w) * ned2l.dscale * dscale.deta^2, 3813 c(w) * ned2l.dq * dq.deta^2, 3814 c(w) * ned2l.dascale * da.deta * dscale.deta, 3815 c(w) * ned2l.dscaleq * dscale.deta * dq.deta, 3816 c(w) * ned2l.daq * da.deta * dq.deta), 3817 dim = c(n, M/M1, M1*(M1+1)/2)) 3818 } 3819 wz <- arwz2wz(wz, M = M, M1 = M1) 3820 wz 3821 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 3822 .escale = escale , .eshape1.a = eshape1.a, 3823 .lshape3.q = lshape3.q, 3824 .eshape3.q = eshape3.q, 3825 .lss = lss )))) 3826} 3827 3828 3829 3830 3831 3832 3833 3834 3835 3836 3837 3838 3839 dagum <- 3840 function(lscale = "loglink", 3841 lshape1.a = "loglink", 3842 lshape2.p = "loglink", 3843 iscale = NULL, 3844 ishape1.a = NULL, 3845 ishape2.p = NULL, 3846 imethod = 1, 3847 lss = TRUE, 3848 gscale = exp(-5:5), 3849 gshape1.a = seq(0.75, 4, by = 0.25), # exp(-5:5), 3850 gshape2.p = exp(-5:5), 3851 probs.y = c(0.25, 0.50, 0.75), 3852 zero = "shape") { 3853 3854 3855 3856 3857 3858 if (length(lss) != 1 && !is.logical(lss)) 3859 stop("Argument 'lss' not specified correctly") 3860 3861 if (!is.Numeric(imethod, length.arg = 1, 3862 integer.valued = TRUE, 3863 positive = TRUE) || imethod > 2) 3864 stop("Bad input for argument 'imethod'") 3865 3866 if (length(iscale) && !is.Numeric(iscale, positive = TRUE)) 3867 stop("Bad input for argument 'iscale'") 3868 3869 if (length(ishape1.a) && !is.Numeric(ishape1.a, positive = TRUE)) 3870 stop("Bad input for argument 'ishape1.a'") 3871 3872 if (length(ishape2.p) && !is.Numeric(ishape2.p, positive = TRUE)) 3873 stop("Bad input for argument 'ishape2.p'") 3874 3875 if (length(probs.y) < 2 || max(probs.y) > 1 || 3876 !is.Numeric(probs.y, positive = TRUE)) 3877 stop("Bad input for argument 'probs.y'") 3878 3879 3880 lscale <- as.list(substitute(lscale)) 3881 escale <- link2list(lscale) 3882 lscale <- attr(escale, "function.name") 3883 3884 lshape1.a <- as.list(substitute(lshape1.a)) 3885 eshape1.a <- link2list(lshape1.a) 3886 lshape1.a <- attr(eshape1.a, "function.name") 3887 3888 lshape2.p <- as.list(substitute(lshape2.p)) 3889 eshape2.p <- link2list(lshape2.p) 3890 lshape2.p <- attr(eshape2.p, "function.name") 3891 3892 3893 new("vglmff", 3894 blurb = 3895 c("Dagum distribution \n\n", 3896 "Links: ", 3897 ifelse (lss, 3898 namesof("scale" , lscale , earg = escale), 3899 namesof("shape1.a", lshape1.a, earg = eshape1.a)), ", ", 3900 ifelse (lss, 3901 namesof("shape1.a", lshape1.a, earg = eshape1.a), 3902 namesof("scale" , lscale , earg = escale)), ", ", 3903 namesof("shape2.p" , lshape2.p, earg = eshape2.p), "\n", 3904 "Mean: scale * gamma(shape2.p + 1/shape1.a) * ", 3905 "gamma(1 - 1/shape1.a) / ", 3906 "gamma(shape2.p)"), 3907 constraints = eval(substitute(expression({ 3908 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 3909 predictors.names = predictors.names, 3910 M1 = 3) 3911 }), list( .zero = zero ))), 3912 infos = eval(substitute(function(...) { 3913 list(M1 = 3, 3914 Q1 = 1, 3915 expected = TRUE, 3916 zero = .zero , 3917 multipleResponses = TRUE, 3918 parameters.names = 3919 if ( .lss ) c("scale", "shape1.a", "shape2.p") else 3920 c("shape1.a", "scale", "shape2.p"), 3921 lscale = .lscale , lshape1.a = .lshape1.a , 3922 escale = .escale , eshape1.a = .eshape1.a , 3923 lshape2.p = .lshape2.p , 3924 eshape2.p = .eshape2.p ) 3925 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 3926 .escale = escale , .eshape1.a = eshape1.a, 3927 .lshape2.p = lshape2.p, 3928 .eshape2.p = eshape2.p, 3929 .lss = lss , 3930 .zero = zero ))), 3931 initialize = eval(substitute(expression({ 3932 temp5 <- w.y.check(w = w, y = y, 3933 Is.positive.y = TRUE, 3934 ncol.w.max = Inf, 3935 ncol.y.max = Inf, 3936 out.wy = TRUE, 3937 colsyperw = 1, 3938 maximize = TRUE) 3939 y <- temp5$y 3940 w <- temp5$w 3941 M1 <- 3 # Number of parameters for one response 3942 NOS <- ncoly <- ncol(y) 3943 M <- M1*ncol(y) 3944 3945 scaL.names <- param.names("scale", NOS, skip1 = TRUE) 3946 sha1.names <- param.names("shape1.a", NOS, skip1 = TRUE) 3947 sha2.names <- param.names("shape2.p", NOS, skip1 = TRUE) 3948 3949 predictors.names <- c( 3950 if ( .lss ) { 3951 c(namesof(scaL.names , .lscale , earg = .escale , tag = FALSE), 3952 namesof(sha1.names , .lshape1.a , earg = .eshape1.a , tag = FALSE)) 3953 } else { 3954 c(namesof(sha1.names , .lshape1.a , earg = .eshape1.a , tag = FALSE), 3955 namesof(scaL.names , .lscale , earg = .escale , tag = FALSE)) 3956 }, 3957 namesof(sha2.names , .lshape2.p , earg = .eshape2.p , tag = FALSE)) 3958 predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)] 3959 3960 if (!length(etastart)) { 3961 sc.init <- 3962 aa.init <- 3963 pp.init <- matrix(NA_real_, n, NOS) 3964 3965 for (spp. in 1:NOS) { # For each response 'y_spp.'... do: 3966 yvec <- y[, spp.] 3967 wvec <- w[, spp.] 3968 3969 if ( .imethod == 1 ) { 3970 gscale <- .gscale 3971 gshape1.a <- .gshape1.a 3972 gshape2.p <- .gshape2.p 3973 if (length( .iscale )) gscale <- rep_len( .iscale , NOS) 3974 if (length( .ishape1.a )) gshape1.a <- rep_len( .ishape1.a , NOS) 3975 if (length( .ishape2.p )) gshape2.p <- rep_len( .ishape2.p , NOS) 3976 3977 3978 3979 try.this <- 3980 grid.search4(gscale, gshape1.a, gshape2.p, vov4 = 1, 3981 objfun = genbetaII.Loglikfun4, # x = x, 3982 y = yvec, w = wvec, # extraargs = NULL, 3983 ret.objfun = TRUE) # Last value is the loglik 3984 sc.init[, spp.] <- try.this["Value1"] 3985 aa.init[, spp.] <- try.this["Value2"] 3986 pp.init[, spp.] <- try.this["Value3"] 3987 } else { # .imethod == 2 3988 qvec <- .probs.y 3989 ishape2.p <- if (length( .ishape2.p )) .ishape2.p else 1 3990 xvec <- log( qvec^(-1/ ishape2.p) - 1 ) 3991 fit0 <- lsfit(x = xvec, y = log(quantile(yvec, qvec))) 3992 sc.init[, spp.] <- if (length( .iscale )) .iscale else 3993 exp(fit0$coef[1]) 3994 aa.init[, spp.] <- if (length( .ishape1.a )) .ishape1.a else 3995 -1/fit0$coef[2] 3996 pp.init[, spp.] <- ishape2.p 3997 } 3998 } # End of for (spp. ...) 3999 4000 4001 4002 finite.mean <- 1 < aa.init 4003 COP.use <- 1.15 4004 while (FALSE && any(!finite.mean)) { 4005 aa.init[!finite.mean] <- 0.1 + aa.init[!finite.mean] * COP.use 4006 finite.mean <- 1 < aa.init 4007 } 4008 4009 etastart <- 4010 cbind(if ( .lss ) 4011 cbind(theta2eta(sc.init, .lscale , earg = .escale ), 4012 theta2eta(aa.init, .lshape1.a , earg = .eshape1.a )) else 4013 cbind(theta2eta(aa.init, .lshape1.a , earg = .eshape1.a ), 4014 theta2eta(sc.init, .lscale , earg = .escale )), 4015 theta2eta(pp.init , .lshape2.p , earg = .eshape2.p )) 4016 etastart <- etastart[, interleave.VGAM(M, M1 = M1)] 4017 } # End of etastart. 4018 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 4019 .escale = escale , .eshape1.a = eshape1.a, 4020 .iscale = iscale , .ishape1.a = ishape1.a, 4021 .gscale = gscale , .gshape1.a = gshape1.a, 4022 .lshape2.p = lshape2.p, 4023 .eshape2.p = eshape2.p, 4024 .ishape2.p = ishape2.p, 4025 .gshape2.p = gshape2.p, 4026 .imethod = imethod , .probs.y = probs.y, 4027 .lss = lss ))), 4028 linkinv = eval(substitute(function(eta, extra = NULL) { 4029 M1 <- 3 4030 NOS <- ncol(eta)/M1 4031 if ( .lss ) { 4032 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4033 .lscale , earg = .escale ) 4034 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4035 .lshape1.a , earg = .eshape1.a ) 4036 } else { 4037 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4038 .lshape1.a , earg = .eshape1.a ) 4039 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4040 .lscale , earg = .escale ) 4041 } 4042 parg <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4043 .lshape2.p , earg = .eshape2.p ) 4044 4045 qq <- 1 4046 4047 ans <- Scale * exp(lgamma(parg + 1/aa) + 4048 lgamma(qq - 1/aa) - lgamma(parg) - lgamma(qq)) 4049 ans[parg + 1/aa <= 0] <- NA 4050 ans[qq - 1/aa <= 0] <- NA 4051 ans[aa <= 0] <- NA 4052 ans[Scale <= 0] <- NA 4053 ans[parg <= 0] <- NA 4054 ans 4055 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 4056 .escale = escale , .eshape1.a = eshape1.a, 4057 .lshape2.p = lshape2.p, 4058 .eshape2.p = eshape2.p, 4059 .lss = lss ))), 4060 last = eval(substitute(expression({ 4061 M1 <- 3 4062 4063 misc$link <- c(rep_len( if ( .lss ) .lscale else .lshape1.a , ncoly), 4064 rep_len( if ( .lss ) .lshape1.a else .lscale , ncoly), 4065 rep_len( .lshape2.p , ncoly))[ 4066 interleave.VGAM(M, M1 = M1)] 4067 temp.names <- if ( .lss ) { 4068 c(scaL.names, sha1.names, sha2.names) 4069 } else { 4070 c(sha1.names, scaL.names, sha2.names) 4071 } 4072 names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)] 4073 4074 misc$earg <- vector("list", M) 4075 names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)] 4076 for (ii in 1:ncoly) { 4077 if ( .lss ) { 4078 misc$earg[[M1*ii-2]] <- .escale 4079 misc$earg[[M1*ii-1]] <- .eshape1.a 4080 } else { 4081 misc$earg[[M1*ii-2]] <- .eshape1.a 4082 misc$earg[[M1*ii-1]] <- .escale 4083 } 4084 misc$earg[[M1*ii ]] <- .eshape2.p 4085 } 4086 4087 misc$expected <- TRUE 4088 misc$multipleResponses <- TRUE 4089 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 4090 .escale = escale , .eshape1.a = eshape1.a, 4091 .lshape2.p = lshape2.p, 4092 .eshape2.p = eshape2.p, 4093 .lss = lss ))), 4094 loglikelihood = eval(substitute( 4095 function(mu, y, w, residuals = FALSE, 4096 eta, extra = NULL, summation = TRUE) { 4097 M1 <- 3 4098 NOS <- ncol(eta)/M1 4099 if ( .lss ) { 4100 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4101 .lscale , earg = .escale ) 4102 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4103 .lshape1.a , earg = .eshape1.a ) 4104 } else { 4105 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4106 .lshape1.a , earg = .eshape1.a ) 4107 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4108 .lscale , earg = .escale ) 4109 } 4110 parg <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4111 .lshape2.p , earg = .eshape2.p ) 4112 if (residuals) { 4113 stop("loglikelihood residuals not implemented yet") 4114 } else { 4115 ll.elts <- 4116 c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = aa, 4117 shape2.p = parg, shape3.q = 1, log = TRUE) 4118 if (summation) { 4119 sum(ll.elts) 4120 } else { 4121 ll.elts 4122 } 4123 } 4124 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 4125 .escale = escale , .eshape1.a = eshape1.a, 4126 .lshape2.p = lshape2.p, 4127 .eshape2.p = eshape2.p, 4128 .lss = lss ))), 4129 vfamily = c("dagum"), 4130 4131 simslot = eval(substitute( 4132 function(object, nsim) { 4133 4134 pwts <- if (length(pwts <- object@prior.weights) > 0) 4135 pwts else weights(object, type = "prior") 4136 if (any(pwts != 1)) 4137 warning("ignoring prior weights") 4138 eta <- predict(object) 4139 if ( .lss ) { 4140 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4141 .lscale , earg = .escale ) 4142 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4143 .lshape1.a , earg = .eshape1.a) 4144 } else { 4145 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4146 .lshape1.a , earg = .eshape1.a) 4147 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4148 .lscale , earg = .escale ) 4149 } 4150 parg <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4151 .lshape2.p , earg = .eshape2.p) 4152 qq <- 1 4153 rdagum(nsim * length(parg), shape1.a = aa, scale = Scale, 4154 shape2.p = parg) 4155 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 4156 .escale = escale , .eshape1.a = eshape1.a, 4157 .lshape2.p = lshape2.p, 4158 .eshape2.p = eshape2.p, 4159 .lss = lss ))), 4160 4161 4162 4163 validparams = eval(substitute(function(eta, y, extra = NULL) { 4164 M1 <- 3 4165 NOS <- ncol(eta)/M1 # Needed for summary() 4166 if ( .lss ) { 4167 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4168 .lscale , earg = .escale ) 4169 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4170 .lshape1.a , earg = .eshape1.a) 4171 } else { 4172 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4173 .lshape1.a , earg = .eshape1.a) 4174 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4175 .lscale , earg = .escale ) 4176 } 4177 parg <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4178 .lshape2.p , earg = .eshape2.p) 4179 qq <- 1 4180 4181 4182 okay1 <- all(is.finite(Scale)) && all(Scale > 0) && 4183 all(is.finite(aa )) && all(aa > 0) && 4184 all(is.finite(parg )) && all(parg > 0) && 4185 all(is.finite(qq )) && all(qq > 0) 4186 okay.support <- if (okay1) all(-aa * parg < 1 & 1 < aa * qq) else TRUE 4187 if (!okay.support) 4188 warning("parameter constraint -a*p < 1 < a has been violated; ", 4189 "solution may be at the boundary of the ", 4190 "parameter space.") 4191 okay1 && okay.support 4192 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 4193 .escale = escale , .eshape1.a = eshape1.a, 4194 .lshape2.p = lshape2.p, 4195 .eshape2.p = eshape2.p, 4196 .lss = lss ))), 4197 4198 4199 4200 4201 deriv = eval(substitute(expression({ 4202 M1 <- 3 4203 NOS <- ncol(eta)/M1 # Needed for summary() 4204 if ( .lss ) { 4205 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4206 .lscale , earg = .escale ) 4207 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4208 .lshape1.a , earg = .eshape1.a) 4209 } else { 4210 aa <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4211 .lshape1.a , earg = .eshape1.a) 4212 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4213 .lscale , earg = .escale ) 4214 } 4215 parg <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4216 .lshape2.p , earg = .eshape2.p) 4217 qq <- 1 4218 temp1 <- log(y/Scale) 4219 temp2 <- (y/Scale)^aa 4220 temp3 <- digamma(parg + qq) 4221 temp3a <- digamma(parg) 4222 temp3b <- digamma(qq) 4223 temp4 <- log1p(temp2) 4224 4225 dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2)) 4226 dl.da <- 1/aa + parg * temp1 - (parg+qq) * temp1 / (1+1/temp2) 4227 dl.dp <- aa * temp1 + temp3 - temp3a - temp4 4228 4229 dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale ) 4230 da.deta <- dtheta.deta(aa, .lshape1.a , earg = .eshape1.a ) 4231 dp.deta <- dtheta.deta(parg, .lshape2.p , earg = .eshape2.p ) 4232 4233 myderiv <- if ( .lss ) { 4234 c(w) * cbind(dl.dscale * dscale.deta, 4235 dl.da * da.deta, 4236 dl.dp * dp.deta) 4237 } else { 4238 c(w) * cbind(dl.da * da.deta, 4239 dl.dscale * dscale.deta, 4240 dl.dp * dp.deta) 4241 } 4242 myderiv[, interleave.VGAM(M, M1 = M1)] 4243 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 4244 .escale = escale , .eshape1.a = eshape1.a, 4245 .lshape2.p = lshape2.p, 4246 .eshape2.p = eshape2.p, 4247 .lss = lss ))), 4248 weight = eval(substitute(expression({ 4249 temp5 <- trigamma(parg + qq) 4250 temp5a <- trigamma(parg) 4251 temp5b <- trigamma(qq) 4252 4253 ned2l.da <- (1 + parg + qq + parg * qq * (temp5a + temp5b + 4254 (temp3b - temp3a + (parg-qq)/(parg*qq))^2 - 4255 (parg^2 + qq^2) / (parg*qq)^2)) / (aa^2 * (1+parg+qq)) 4256 ned2l.dscale <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2) 4257 ned2l.dp <- temp5a - temp5 4258 ned2l.dascale <- (parg - qq - parg * qq * 4259 (temp3a -temp3b)) / (Scale*(1 + parg+qq)) 4260 ned2l.dap <- -(qq * (temp3a -temp3b) -1) / (aa*(parg+qq)) 4261 ned2l.dscalep <- aa * qq / (Scale*(parg+qq)) 4262 wz <- if ( .lss ) { 4263 array(c(c(w) * ned2l.dscale * dscale.deta^2, 4264 c(w) * ned2l.da * da.deta^2, 4265 c(w) * ned2l.dp * dp.deta^2, 4266 c(w) * ned2l.dascale * da.deta * dscale.deta, 4267 c(w) * ned2l.dap * da.deta * dp.deta, 4268 c(w) * ned2l.dscalep * dscale.deta * dp.deta), 4269 dim = c(n, M/M1, M1*(M1+1)/2)) 4270 } else { 4271 array(c(c(w) * ned2l.da * da.deta^2, 4272 c(w) * ned2l.dscale * dscale.deta^2, 4273 c(w) * ned2l.dp * dp.deta^2, 4274 c(w) * ned2l.dascale * da.deta * dscale.deta, 4275 c(w) * ned2l.dscalep * dscale.deta * dp.deta, 4276 c(w) * ned2l.dap * da.deta * dp.deta), 4277 dim = c(n, M/M1, M1*(M1+1)/2)) 4278 } 4279 wz <- arwz2wz(wz, M = M, M1 = M1) 4280 wz 4281 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 4282 .escale = escale , .eshape1.a = eshape1.a, 4283 .lshape2.p = lshape2.p, 4284 .eshape2.p = eshape2.p, 4285 .lss = lss )))) 4286} 4287 4288 4289 4290 4291 4292 4293 betaII <- 4294 function(lscale = "loglink", 4295 lshape2.p = "loglink", 4296 lshape3.q = "loglink", 4297 iscale = NULL, 4298 ishape2.p = NULL, 4299 ishape3.q = NULL, 4300 imethod = 1, 4301 gscale = exp(-5:5), 4302 gshape2.p = exp(-5:5), 4303 gshape3.q = seq(0.75, 4, by = 0.25), # exp(-5:5), 4304 probs.y = c(0.25, 0.50, 0.75), 4305 zero = "shape") { 4306 4307 4308 4309 4310 if (!is.Numeric(imethod, length.arg = 1, 4311 integer.valued = TRUE, 4312 positive = TRUE) || imethod > 2) 4313 stop("Bad input for argument 'imethod'") 4314 4315 if (length(iscale ) && !is.Numeric(iscale , positive = TRUE)) 4316 stop("Bad input for argument 'iscale'") 4317 4318 if (length(ishape2.p) && !is.Numeric(ishape2.p, positive = TRUE)) 4319 stop("Bad input for argument 'ishape2.p'") 4320 4321 if (length(ishape3.q) && !is.Numeric(ishape3.q, positive = TRUE)) 4322 stop("Bad input for argument 'ishape3.q'") 4323 4324 if (length(probs.y) < 2 || max(probs.y) > 1 || 4325 !is.Numeric(probs.y, positive = TRUE)) 4326 stop("Bad input for argument 'probs.y'") 4327 4328 4329 lscale <- as.list(substitute(lscale)) 4330 escale <- link2list(lscale) 4331 lscale <- attr(escale, "function.name") 4332 4333 lshape2.p <- as.list(substitute(lshape2.p)) 4334 eshape2.p <- link2list(lshape2.p) 4335 lshape2.p <- attr(eshape2.p, "function.name") 4336 4337 lshape3.q <- as.list(substitute(lshape3.q)) 4338 eshape3.q <- link2list(lshape3.q) 4339 lshape3.q <- attr(eshape3.q, "function.name") 4340 4341 4342 new("vglmff", 4343 blurb = 4344 c("Beta II distribution \n\n", 4345 "Links: ", 4346 namesof("scale" , lscale , earg = escale ), ", ", 4347 namesof("shape2.p" , lshape2.p, earg = eshape2.p), ", ", 4348 namesof("shape3.q" , lshape3.q, earg = eshape3.q), "\n", 4349 "Mean: scale * gamma(shape2.p + 1) * ", 4350 "gamma(shape3.q - 1) / ", 4351 "(gamma(shape2.p) * gamma(shape3.q))"), 4352 constraints = eval(substitute(expression({ 4353 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 4354 predictors.names = predictors.names, 4355 M1 = 3) 4356 }), list( .zero = zero ))), 4357 infos = eval(substitute(function(...) { 4358 list(M1 = 3, 4359 Q1 = 1, 4360 expected = TRUE, 4361 zero = .zero , 4362 multipleResponses = TRUE, 4363 parameters.names = c("scale", "shape2.p", "shape3.q"), 4364 lscale = .lscale , 4365 escale = .escale , 4366 lshape2.p = .lshape2.p , lshape3.q = .lshape3.q , 4367 eshape2.p = .eshape2.p , eshape3.q = .eshape3.q ) 4368 }, list( .lscale = lscale , 4369 .escale = escale , 4370 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 4371 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q, 4372 .zero = zero ))), 4373 initialize = eval(substitute(expression({ 4374 temp5 <- w.y.check(w = w, y = y, 4375 Is.positive.y = TRUE, 4376 ncol.w.max = Inf, 4377 ncol.y.max = Inf, 4378 out.wy = TRUE, 4379 colsyperw = 1, 4380 maximize = TRUE) 4381 y <- temp5$y 4382 w <- temp5$w 4383 M1 <- 3 # Number of parameters for one response 4384 NOS <- ncoly <- ncol(y) 4385 M <- M1*ncol(y) 4386 4387 scaL.names <- param.names("scale", NOS, skip1 = TRUE) 4388 sha2.names <- param.names("shape2.p", NOS, skip1 = TRUE) 4389 sha3.names <- param.names("shape3.q", NOS, skip1 = TRUE) 4390 4391 predictors.names <- 4392 c(namesof(scaL.names , .lscale , earg = .escale , tag = FALSE), 4393 namesof(sha2.names , .lshape2.p , earg = .eshape2.p , tag = FALSE), 4394 namesof(sha3.names , .lshape3.q , earg = .eshape3.q , tag = FALSE)) 4395 predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)] 4396 4397 if (!length(etastart)) { 4398 sc.init <- 4399 pp.init <- 4400 qq.init <- matrix(NA_real_, n, NOS) 4401 4402 for (spp. in 1:NOS) { # For each response 'y_spp.'... do: 4403 yvec <- y[, spp.] 4404 wvec <- w[, spp.] 4405 4406 if ( .imethod == 1 ) { 4407 gscale <- .gscale 4408 gshape2.p <- .gshape2.p 4409 gshape3.q <- .gshape3.q 4410 if (length( .iscale )) gscale <- rep_len( .iscale , NOS) 4411 if (length( .ishape2.p )) gshape2.p <- rep_len( .ishape2.p , NOS) 4412 if (length( .ishape3.q )) gshape3.q <- rep_len( .ishape3.q , NOS) 4413 4414 4415 4416 4417 try.this <- 4418 grid.search4(gscale, vov2 = 1, gshape2.p, gshape3.q, 4419 objfun = genbetaII.Loglikfun4, # x = x, 4420 y = yvec, w = wvec, # extraargs = NULL, 4421 ret.objfun = TRUE) # Last value is the loglik 4422 sc.init[, spp.] <- try.this["Value1"] 4423 pp.init[, spp.] <- try.this["Value3"] 4424 qq.init[, spp.] <- try.this["Value4"] 4425 } else { # .imethod == 2 4426 4427 sc.init[, spp.] <- if (length( .iscale )) .iscale else { 4428 qvec <- .probs.y 4429 ishape3.q <- if (length( .ishape3.q )) .ishape3.q else 1 4430 xvec <- log( (1-qvec)^(-1/ ishape3.q ) - 1 ) 4431 fit0 <- lsfit(x = xvec, y = log(quantile(yvec, qvec ))) 4432 exp(fit0$coef[1]) 4433 } 4434 pp.init[, spp.] <- if (length( .ishape2.p )) .ishape2.p else 1.0 4435 qq.init[, spp.] <- if (length( .ishape3.q )) .ishape3.q else 1.0 4436 } 4437 } # End of for (spp. ...) 4438 4439 4440 4441 4442 finite.mean <- 1 < qq.init 4443 COP.use <- 1.15 4444 while (any(!finite.mean)) { 4445 qq.init[!finite.mean] <- 0.1 + qq.init[!finite.mean] * COP.use 4446 finite.mean <- 1 < qq.init 4447 } 4448 4449 etastart <- 4450 cbind(theta2eta(sc.init , .lscale , earg = .escale ), 4451 theta2eta(pp.init , .lshape2.p , earg = .eshape2.p ), 4452 theta2eta(qq.init , .lshape3.q , earg = .eshape3.q )) 4453 etastart <- etastart[, interleave.VGAM(M, M1 = M1)] 4454 } # End of etastart. 4455 }), list( .lscale = lscale , 4456 .escale = escale , 4457 .iscale = iscale , 4458 .gscale = gscale , 4459 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 4460 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q, 4461 .ishape2.p = ishape2.p, .ishape3.q = ishape3.q, 4462 .gshape2.p = gshape2.p, .gshape3.q = gshape3.q, 4463 .imethod = imethod , .probs.y = probs.y 4464 ))), 4465 linkinv = eval(substitute(function(eta, extra = NULL) { 4466 M1 <- 3 4467 NOS <- ncol(eta)/M1 4468 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4469 .lscale , earg = .escale ) 4470 aa <- 1 4471 parg <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4472 .lshape2.p , earg = .eshape2.p ) 4473 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4474 .lshape3.q , earg = .eshape3.q ) 4475 4476 ans <- Scale * exp(lgamma(parg + 1/aa) + 4477 lgamma(qq - 1/aa) - lgamma(parg) - lgamma(qq)) 4478 ans[parg + 1/aa <= 0] <- NA 4479 ans[qq - 1/aa <= 0] <- NA 4480 ans[Scale <= 0] <- NA 4481 ans[parg <= 0] <- NA 4482 ans[qq <= 0] <- NA 4483 ans 4484 }, list( .lscale = lscale , 4485 .escale = escale , 4486 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 4487 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q 4488 ))), 4489 last = eval(substitute(expression({ 4490 M1 <- 3 4491 4492 misc$link <- 4493 c(rep_len( .lscale , ncoly), 4494 rep_len( .lshape2.p , ncoly), 4495 rep_len( .lshape3.q , ncoly))[interleave.VGAM(M, M1 = M1)] 4496 temp.names <- c(scaL.names, sha2.names, sha3.names) 4497 names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)] 4498 4499 misc$earg <- vector("list", M) 4500 names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)] 4501 for (ii in 1:ncoly) { 4502 misc$earg[[M1*ii-2]] <- .escale 4503 misc$earg[[M1*ii-1]] <- .eshape2.p 4504 misc$earg[[M1*ii ]] <- .eshape3.q 4505 } 4506 4507 misc$expected <- TRUE 4508 misc$multipleResponses <- TRUE 4509 }), list( .lscale = lscale , 4510 .escale = escale , 4511 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 4512 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q 4513 ))), 4514 loglikelihood = eval(substitute( 4515 function(mu, y, w, residuals = FALSE, 4516 eta, extra = NULL, summation = TRUE) { 4517 M1 <- 3 4518 NOS <- ncol(eta)/M1 4519 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4520 .lscale , earg = .escale ) 4521 aa <- 1 4522 parg <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4523 .lshape2.p , earg = .eshape2.p ) 4524 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4525 .lshape3.q , earg = .eshape3.q ) 4526 if (residuals) { 4527 stop("loglikelihood residuals not implemented yet") 4528 } else { 4529 ll.elts <- 4530 c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = aa, 4531 shape2.p = parg, shape3.q = qq, log = TRUE) 4532 if (summation) { 4533 sum(ll.elts) 4534 } else { 4535 ll.elts 4536 } 4537 } 4538 }, list( .lscale = lscale , 4539 .escale = escale , 4540 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 4541 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q 4542 ))), 4543 vfamily = c("betaII"), 4544 4545 validparams = eval(substitute(function(eta, y, extra = NULL) { 4546 M1 <- 3 4547 NOS <- ncol(eta)/M1 4548 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4549 .lscale , earg = .escale ) 4550 aa <- 1 4551 parg <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4552 .lshape2.p , earg = .eshape2.p ) 4553 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4554 .lshape3.q , earg = .eshape3.q ) 4555 4556 4557 okay1 <- all(is.finite(Scale)) && all(Scale > 0) && 4558 all(is.finite(aa )) && all(aa > 0) && 4559 all(is.finite(parg )) && all(parg > 0) && 4560 all(is.finite(qq )) && all(qq > 0) 4561 okay.support <- if (okay1) all(-aa * parg < 1 & 1 < aa * qq) else TRUE 4562 if (!okay.support) 4563 warning("parameter constraint -p < 1 < q has been violated; ", 4564 "solution may be at the boundary of the ", 4565 "parameter space.") 4566 okay1 && okay.support 4567 }, list( .lscale = lscale , 4568 .escale = escale , 4569 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 4570 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q ))), 4571 4572 4573 4574 4575 deriv = eval(substitute(expression({ 4576 M1 <- 3 4577 NOS <- ncol(eta)/M1 # Needed for summary() 4578 Scale <- eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 4579 .lscale , earg = .escale ) 4580 aa <- 1 4581 parg <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4582 .lshape2.p , earg = .eshape2.p ) 4583 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4584 .lshape3.q , earg = .eshape3.q ) 4585 temp1 <- log(y/Scale) 4586 temp2 <- (y/Scale)^aa 4587 temp3 <- digamma(parg + qq) 4588 temp3a <- digamma(parg) 4589 temp3b <- digamma(qq) 4590 temp4 <- log1p(temp2) 4591 4592 dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2)) 4593 dl.dp <- aa * temp1 + temp3 - temp3a - temp4 4594 dl.dq <- temp3 - temp3b - temp4 4595 4596 dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale ) 4597 dp.deta <- dtheta.deta(parg, .lshape2.p , earg = .eshape2.p ) 4598 dq.deta <- dtheta.deta(qq, .lshape3.q , earg = .eshape3.q ) 4599 4600 myderiv <- 4601 c(w) * cbind(dl.dscale * dscale.deta, 4602 dl.dp * dp.deta, 4603 dl.dq * dq.deta) 4604 myderiv[, interleave.VGAM(M, M1 = M1)] 4605 }), list( .lscale = lscale , 4606 .escale = escale , 4607 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 4608 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q 4609 ))), 4610 weight = eval(substitute(expression({ 4611 temp5 <- trigamma(parg + qq) 4612 temp5a <- trigamma(parg) 4613 temp5b <- trigamma(qq) 4614 4615 ned2l.dscale <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2) 4616 ned2l.dp <- temp5a - temp5 4617 ned2l.dq <- temp5b - temp5 4618 ned2l.dscalep <- aa * qq / (Scale*(parg+qq)) 4619 ned2l.dscaleq <- -aa * parg / (Scale*(parg+qq)) 4620 ned2l.dpq <- -temp5 4621 wz <- 4622 array(c(c(w) * ned2l.dscale * dscale.deta^2, 4623 c(w) * ned2l.dp * dp.deta^2, 4624 c(w) * ned2l.dq * dq.deta^2, 4625 c(w) * ned2l.dscalep * dscale.deta * dp.deta, 4626 c(w) * ned2l.dpq * dp.deta * dq.deta, # Switched!! 4627 c(w) * ned2l.dscaleq * dscale.deta * dq.deta), 4628 dim = c(n, M/M1, M1*(M1+1)/2)) 4629 wz <- arwz2wz(wz, M = M, M1 = M1) 4630 wz 4631 }), list( .lscale = lscale , 4632 .escale = escale , 4633 .lshape2.p = lshape2.p, .lshape3.q = lshape3.q, 4634 .eshape2.p = eshape2.p, .eshape3.q = eshape3.q 4635 )))) 4636} 4637 4638 4639 4640 4641 4642 4643 lomax <- 4644 function(lscale = "loglink", 4645 lshape3.q = "loglink", 4646 iscale = NULL, 4647 ishape3.q = NULL, 4648 imethod = 1, 4649 gscale = exp(-5:5), 4650 gshape3.q = seq(0.75, 4, by = 0.25), # exp(-5:5), 4651 probs.y = c(0.25, 0.50, 0.75), 4652 zero = "shape") { 4653 4654 4655 4656 4657 if (!is.Numeric(imethod, length.arg = 1, 4658 integer.valued = TRUE, 4659 positive = TRUE) || imethod > 2) 4660 stop("Bad input for argument 'imethod'") 4661 4662 if (length(iscale) && !is.Numeric(iscale, positive = TRUE)) 4663 stop("Bad input for argument 'iscale'") 4664 4665 if (length(ishape3.q) && !is.Numeric(ishape3.q, positive = TRUE)) 4666 stop("Bad input for argument 'ishape3.q'") 4667 4668 if (length(probs.y) < 2 || max(probs.y) > 1 || 4669 !is.Numeric(probs.y, positive = TRUE)) 4670 stop("Bad input for argument 'probs.y'") 4671 4672 4673 lscale <- as.list(substitute(lscale)) 4674 escale <- link2list(lscale) 4675 lscale <- attr(escale, "function.name") 4676 4677 lshape3.q <- as.list(substitute(lshape3.q)) 4678 eshape3.q <- link2list(lshape3.q) 4679 lshape3.q <- attr(eshape3.q, "function.name") 4680 4681 4682 new("vglmff", 4683 blurb = 4684 c("Lomax distribution \n\n", 4685 "Links: ", 4686 namesof("scale" , lscale , earg = escale ), ", ", 4687 namesof("shape3.q" , lshape3.q, earg = eshape3.q), "\n", 4688 "Mean: scale / (shape3.q - 1)"), 4689 constraints = eval(substitute(expression({ 4690 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 4691 predictors.names = predictors.names, 4692 M1 = 2) 4693 }), list( .zero = zero ))), 4694 infos = eval(substitute(function(...) { 4695 list(M1 = 2, 4696 Q1 = 1, 4697 expected = TRUE, 4698 zero = .zero , 4699 multipleResponses = TRUE, 4700 parameters.names = c("scale", "shape3.q"), 4701 lscale = .lscale , 4702 escale = .escale , 4703 lshape3.q = .lshape3.q , 4704 eshape3.q = .eshape3.q ) 4705 }, list( .lscale = lscale , 4706 .escale = escale , 4707 .lshape3.q = lshape3.q, 4708 .eshape3.q = eshape3.q, 4709 .zero = zero ))), 4710 initialize = eval(substitute(expression({ 4711 temp5 <- w.y.check(w = w, y = y, 4712 Is.positive.y = TRUE, 4713 ncol.w.max = Inf, 4714 ncol.y.max = Inf, 4715 out.wy = TRUE, 4716 colsyperw = 1, 4717 maximize = TRUE) 4718 y <- temp5$y 4719 w <- temp5$w 4720 M1 <- 2 # Number of parameters for one response 4721 NOS <- ncoly <- ncol(y) 4722 M <- M1*ncol(y) 4723 4724 scaL.names <- param.names("scale", NOS, skip1 = TRUE) 4725 sha3.names <- param.names("shape3.q", NOS, skip1 = TRUE) 4726 4727 predictors.names <- 4728 c(namesof(scaL.names , .lscale , earg = .escale , tag = FALSE), 4729 namesof(sha3.names , .lshape3.q , earg = .eshape3.q , tag = FALSE)) 4730 predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)] 4731 4732 if (!length(etastart)) { 4733 sc.init <- 4734 qq.init <- matrix(NA_real_, n, NOS) 4735 4736 for (spp. in 1:NOS) { # For each response 'y_spp.'... do: 4737 yvec <- y[, spp.] 4738 wvec <- w[, spp.] 4739 4740 if ( .imethod == 1) { 4741 gscale <- .gscale 4742 gshape3.q <- .gshape3.q 4743 if (length( .iscale )) gscale <- rep_len( .iscale , NOS) 4744 if (length( .ishape3.q )) gshape3.q <- rep_len( .ishape3.q , NOS) 4745 4746 4747 4748 try.this <- 4749 grid.search4(gscale, vov2 = 1, vov3 = 1, gshape3.q, 4750 objfun = genbetaII.Loglikfun4, # x = x, 4751 y = yvec, w = wvec, # extraargs = NULL, 4752 ret.objfun = TRUE) # Last value is the loglik 4753 sc.init[, spp.] <- try.this["Value1"] 4754 qq.init[, spp.] <- try.this["Value4"] 4755 } else { # .imethod == 2 4756 qvec <- .probs.y 4757 iscale <- if (length( .iscale )) .iscale else 1 4758 xvec <- log1p( quantile(yvec / iscale, probs = qvec) ) 4759 fit0 <- lsfit(x = xvec, y = -log1p(-qvec), intercept = FALSE) 4760 sc.init[, spp.] <- iscale 4761 qq.init[, spp.] <- if (length( .ishape3.q )) .ishape3.q else 4762 fit0$coef 4763 } 4764 } # End of for (spp. ...) 4765 4766 finite.mean <- 1 < qq.init 4767 COP.use <- 1.15 4768 while (FALSE && any(!finite.mean)) { 4769 qq.init[!finite.mean] <- 0.1 + qq.init[!finite.mean] * COP.use 4770 finite.mean <- 1 < qq.init 4771 } 4772 4773 etastart <- 4774 cbind(theta2eta(sc.init, .lscale , earg = .escale ), 4775 theta2eta(qq.init, .lshape3.q , earg = .eshape3.q )) 4776 etastart <- etastart[, interleave.VGAM(M, M1 = M1)] 4777 } # End of etastart. 4778 }), list( .lscale = lscale , 4779 .escale = escale , 4780 .iscale = iscale , 4781 .gscale = gscale , 4782 .lshape3.q = lshape3.q, 4783 .eshape3.q = eshape3.q, 4784 .ishape3.q = ishape3.q, 4785 .gshape3.q = gshape3.q, 4786 .imethod = imethod , .probs.y = probs.y 4787 ))), 4788 linkinv = eval(substitute(function(eta, extra = NULL) { 4789 M1 <- 2 4790 NOS <- ncol(eta)/M1 4791 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4792 .lscale , earg = .escale ) 4793 aa <- 1 4794 parg <- 1 4795 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4796 .lshape3.q , earg = .eshape3.q ) 4797 4798 ans <- Scale * exp(lgamma(parg + 1/aa) + 4799 lgamma(qq - 1/aa) - lgamma(parg) - lgamma(qq)) 4800 ans[parg + 1/aa <= 0] <- NA 4801 ans[qq - 1/aa <= 0] <- NA 4802 ans[Scale <= 0] <- NA 4803 ans[qq <= 0] <- NA 4804 ans 4805 }, list( .lscale = lscale , 4806 .escale = escale , 4807 .lshape3.q = lshape3.q, 4808 .eshape3.q = eshape3.q 4809 ))), 4810 last = eval(substitute(expression({ 4811 M1 <- 2 4812 4813 misc$link <- 4814 c(rep_len( .lscale , ncoly), 4815 rep_len( .lshape3.q , ncoly))[interleave.VGAM(M, M1 = M1)] 4816 temp.names <- 4817 c(scaL.names, sha3.names) 4818 names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)] 4819 4820 misc$earg <- vector("list", M) 4821 names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)] 4822 for (ii in 1:ncoly) { 4823 misc$earg[[M1*ii-1]] <- .escale 4824 misc$earg[[M1*ii ]] <- .eshape3.q 4825 } 4826 4827 misc$expected <- TRUE 4828 misc$multipleResponses <- TRUE 4829 }), list( .lscale = lscale , 4830 .escale = escale , 4831 .lshape3.q = lshape3.q, 4832 .eshape3.q = eshape3.q 4833 ))), 4834 loglikelihood = eval(substitute( 4835 function(mu, y, w, residuals = FALSE, 4836 eta, extra = NULL, summation = TRUE) { 4837 M1 <- 2 4838 NOS <- ncol(eta)/M1 4839 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4840 .lscale , earg = .escale ) 4841 aa <- 1 4842 parg <- 1 4843 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4844 .lshape3.q , earg = .eshape3.q ) 4845 if (residuals) { 4846 stop("loglikelihood residuals not implemented yet") 4847 } else { 4848 ll.elts <- 4849 c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = aa, 4850 shape2.p = parg, shape3.q = qq, log = TRUE) 4851 if (summation) { 4852 sum(ll.elts) 4853 } else { 4854 ll.elts 4855 } 4856 } 4857 }, list( .lscale = lscale , 4858 .escale = escale , 4859 .lshape3.q = lshape3.q, 4860 .eshape3.q = eshape3.q 4861 ))), 4862 vfamily = c("lomax"), 4863 simslot = eval(substitute( 4864 function(object, nsim) { 4865 4866 pwts <- if (length(pwts <- object@prior.weights) > 0) 4867 pwts else weights(object, type = "prior") 4868 if (any(pwts != 1)) 4869 warning("ignoring prior weights") 4870 4871 eta <- predict(object) 4872 M1 <- 2 4873 NOS <- ncol(eta)/M1 4874 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4875 .lscale , earg = .escale ) 4876 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4877 .lshape3.q , earg = .eshape3.q ) 4878 4879 rlomax(nsim * length(qq), scale = Scale, shape3.q = qq) 4880 }, list( .lscale = lscale , 4881 .escale = escale , 4882 .lshape3.q = lshape3.q, 4883 .eshape3.q = eshape3.q 4884 ))), 4885 4886 validparams = eval(substitute(function(eta, y, extra = NULL) { 4887 M1 <- 2 4888 NOS <- ncol(eta) / M1 4889 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4890 .lscale , earg = .escale ) 4891 aa <- 1 4892 parg <- 1 4893 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4894 .lshape3.q , earg = .eshape3.q) 4895 okay1 <- all(is.finite(Scale)) && all(Scale > 0) && 4896 all(is.finite(qq )) && all(qq > 0) 4897 okay.support <- if (okay1) all(-aa * parg < 1 & 1 < aa * qq) else TRUE 4898 if (!okay.support) 4899 warning("parameter constraint 1 < q has been violated; ", 4900 "solution may be at the boundary of the ", 4901 "parameter space.") 4902 okay1 && okay.support 4903 }, list( .lscale = lscale , 4904 .escale = escale , 4905 .lshape3.q = lshape3.q, 4906 .eshape3.q = eshape3.q ))), 4907 4908 4909 4910 4911 deriv = eval(substitute(expression({ 4912 M1 <- 2 4913 NOS <- ncol(eta)/M1 # Needed for summary() 4914 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 4915 .lscale , earg = .escale ) 4916 aa <- 1 4917 parg <- 1 4918 qq <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 4919 .lshape3.q , earg = .eshape3.q) 4920 temp1 <- log(y/Scale) 4921 temp2 <- (y/Scale)^aa 4922 temp3 <- digamma(parg + qq) 4923 temp3a <- digamma(parg) 4924 temp3b <- digamma(qq) 4925 temp4 <- log1p(temp2) 4926 4927 dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2)) 4928 dl.dq <- temp3 - temp3b - temp4 4929 4930 dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale ) 4931 dq.deta <- dtheta.deta(qq, .lshape3.q , earg = .eshape3.q ) 4932 4933 myderiv <- 4934 c(w) * cbind(dl.dscale * dscale.deta, 4935 dl.dq * dq.deta) 4936 myderiv[, interleave.VGAM(M, M1 = M1)] 4937 }), list( .lscale = lscale , 4938 .escale = escale , 4939 .lshape3.q = lshape3.q, 4940 .eshape3.q = eshape3.q 4941 ))), 4942 weight = eval(substitute(expression({ 4943 temp5 <- trigamma(parg + qq) 4944 temp5a <- trigamma(parg) 4945 temp5b <- trigamma(qq) 4946 4947 ned2l.dscale <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2) 4948 ned2l.dq <- temp5b - temp5 4949 ned2l.dscaleq <- -aa * parg / (Scale*(parg+qq)) 4950 wz <- 4951 array(c(c(w) * ned2l.dscale * dscale.deta^2, 4952 c(w) * ned2l.dq * dq.deta^2, 4953 c(w) * ned2l.dscaleq * dscale.deta * dq.deta), 4954 dim = c(n, M/M1, M1*(M1+1)/2)) 4955 wz <- arwz2wz(wz, M = M, M1 = M1) 4956 wz 4957 }), list( .lscale = lscale , 4958 .escale = escale , 4959 .lshape3.q = lshape3.q, 4960 .eshape3.q = eshape3.q 4961 )))) 4962} 4963 4964 4965 4966 4967 4968 4969 4970 4971 4972 4973 fisk <- 4974 function(lscale = "loglink", 4975 lshape1.a = "loglink", 4976 iscale = NULL, 4977 ishape1.a = NULL, 4978 imethod = 1, 4979 lss = TRUE, 4980 gscale = exp(-5:5), 4981 gshape1.a = seq(0.75, 4, by = 0.25), # exp(-5:5), 4982 probs.y = c(0.25, 0.50, 0.75), 4983 zero = "shape") { 4984 4985 4986 4987 4988 4989 if (length(lss) != 1 && !is.logical(lss)) 4990 stop("Argument 'lss' not specified correctly") 4991 4992 if (!is.Numeric(imethod, length.arg = 1, 4993 integer.valued = TRUE, 4994 positive = TRUE) || imethod > 2) 4995 stop("Bad input for argument 'imethod'") 4996 4997 if (length(iscale) && !is.Numeric(iscale, positive = TRUE)) 4998 stop("Bad input for argument 'iscale'") 4999 5000 if (length(ishape1.a) && !is.Numeric(ishape1.a, positive = TRUE)) 5001 stop("Bad input for argument 'ishape1.a'") 5002 5003 if (length(probs.y) < 2 || max(probs.y) > 1 || 5004 !is.Numeric(probs.y, positive = TRUE)) 5005 stop("Bad input for argument 'probs.y'") 5006 5007 5008 lscale <- as.list(substitute(lscale)) 5009 escale <- link2list(lscale) 5010 lscale <- attr(escale, "function.name") 5011 5012 lshape1.a <- as.list(substitute(lshape1.a)) 5013 eshape1.a <- link2list(lshape1.a) 5014 lshape1.a <- attr(eshape1.a, "function.name") 5015 5016 5017 new("vglmff", 5018 blurb = 5019 c("Fisk distribution \n\n", 5020 "Links: ", 5021 ifelse (lss, 5022 namesof("scale" , lscale , earg = escale), 5023 namesof("shape1.a", lshape1.a, earg = eshape1.a)), ", ", 5024 ifelse (lss, 5025 namesof("shape1.a", lshape1.a, earg = eshape1.a), 5026 namesof("scale" , lscale , earg = escale)), "\n", 5027 "Mean: scale * gamma(1 + 1/shape1.a) * ", 5028 "gamma(1 - 1/shape1.a)"), 5029 constraints = eval(substitute(expression({ 5030 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 5031 predictors.names = predictors.names, 5032 M1 = 2) 5033 }), list( .zero = zero ))), 5034 infos = eval(substitute(function(...) { 5035 list(M1 = 2, 5036 Q1 = 1, 5037 expected = TRUE, 5038 zero = .zero , 5039 multipleResponses = TRUE, 5040 parameters.names = if ( .lss ) 5041 c("scale", "shape1.a") else 5042 c("shape1.a", "scale"), 5043 lscale = .lscale , lshape1.a = .lshape1.a , 5044 escale = .escale , eshape1.a = .eshape1.a ) 5045 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 5046 .escale = escale , .eshape1.a = eshape1.a, 5047 .lss = lss , 5048 .zero = zero ))), 5049 initialize = eval(substitute(expression({ 5050 temp5 <- w.y.check(w = w, y = y, 5051 Is.positive.y = TRUE, 5052 ncol.w.max = Inf, 5053 ncol.y.max = Inf, 5054 out.wy = TRUE, 5055 colsyperw = 1, 5056 maximize = TRUE) 5057 y <- temp5$y 5058 w <- temp5$w 5059 M1 <- 2 # Number of parameters for one response 5060 NOS <- ncoly <- ncol(y) 5061 M <- M1*ncol(y) 5062 5063 scaL.names <- param.names("scale", NOS, skip1 = TRUE) 5064 sha1.names <- param.names("shape1.a", NOS, skip1 = TRUE) 5065 5066 predictors.names <- if ( .lss ) { 5067 c(namesof(scaL.names , .lscale , earg = .escale , tag = FALSE), 5068 namesof(sha1.names , .lshape1.a , earg = .eshape1.a , tag = FALSE)) 5069 } else { 5070 c(namesof(sha1.names , .lshape1.a , earg = .eshape1.a , tag = FALSE), 5071 namesof(scaL.names , .lscale , earg = .escale , tag = FALSE)) 5072 } 5073 predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)] 5074 5075 if (!length(etastart)) { 5076 sc.init <- 5077 aa.init <- matrix(NA_real_, n, NOS) 5078 5079 for (spp. in 1:NOS) { # For each response 'y_spp.'... do: 5080 yvec <- y[, spp.] 5081 wvec <- w[, spp.] 5082 5083 if ( .imethod == 1 ) { 5084 gscale <- .gscale 5085 gshape1.a <- .gshape1.a 5086 if (length( .iscale )) gscale <- rep_len( .iscale , NOS) 5087 if (length( .ishape1.a )) gshape1.a <- rep_len( .ishape1.a , NOS) 5088 5089 5090 5091 5092 try.this <- 5093 grid.search4(gscale, gshape1.a, vov3 = 1, vov4 = 1, 5094 objfun = genbetaII.Loglikfun4, # x = x, 5095 y = yvec, w = wvec, # extraargs = NULL, 5096 ret.objfun = TRUE) # Last value is the loglik 5097 sc.init[, spp.] <- try.this["Value1"] 5098 aa.init[, spp.] <- try.this["Value2"] 5099 } else { # .imethod == 2 5100 qvec <- .probs.y 5101 iscale <- if (length( .iscale )) .iscale else 1 5102 xvec <- log( quantile(yvec / iscale, probs = qvec) ) 5103 fit0 <- lsfit(x = xvec, y = logitlink(qvec), intercept = FALSE) 5104 sc.init[, spp.] <- iscale 5105 aa.init[, spp.] <- if (length( .ishape1.a )) .ishape1.a else 5106 fit0$coef 5107 } 5108 } # End of for (spp. ...) 5109 5110 5111 finite.mean <- 1 < aa.init 5112 COP.use <- 1.15 5113 while (FALSE && any(!finite.mean)) { 5114 aa.init[!finite.mean] <- 0.1 + aa.init[!finite.mean] * COP.use 5115 finite.mean <- 1 < aa.init 5116 } 5117 5118 etastart <- 5119 if ( .lss ) 5120 cbind(theta2eta(sc.init, .lscale , earg = .escale ), 5121 theta2eta(aa.init, .lshape1.a , earg = .eshape1.a )) else 5122 cbind(theta2eta(aa.init, .lshape1.a , earg = .eshape1.a ), 5123 theta2eta(sc.init, .lscale , earg = .escale )) 5124 etastart <- etastart[, interleave.VGAM(M, M1 = M1)] 5125 } # End of etastart. 5126 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 5127 .escale = escale , .eshape1.a = eshape1.a, 5128 .iscale = iscale , .ishape1.a = ishape1.a, 5129 .gscale = gscale , .gshape1.a = gshape1.a, 5130 .imethod = imethod , .probs.y = probs.y, 5131 .lss = lss ))), 5132 linkinv = eval(substitute(function(eta, extra = NULL) { 5133 M1 <- 2 5134 NOS <- ncol(eta)/M1 5135 if ( .lss ) { 5136 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5137 .lscale , earg = .escale ) 5138 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5139 .lshape1.a , earg = .eshape1.a ) 5140 } else { 5141 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5142 .lshape1.a , earg = .eshape1.a ) 5143 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5144 .lscale , earg = .escale ) 5145 } 5146 parg <- 1 5147 qq <- 1 5148 5149 ans <- Scale * exp(lgamma(parg + 1/aa) + 5150 lgamma(qq - 1/aa) - lgamma(parg) - lgamma(qq)) 5151 ans[parg + 1/aa <= 0] <- NA 5152 ans[qq - 1/aa <= 0] <- NA 5153 ans[Scale <= 0] <- NA 5154 ans[aa <= 0] <- NA 5155 ans 5156 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 5157 .escale = escale , .eshape1.a = eshape1.a, 5158 .lss = lss ))), 5159 last = eval(substitute(expression({ 5160 M1 <- 2 5161 5162 misc$link <- c(rep_len( if ( .lss ) .lscale else .lshape1.a , ncoly), 5163 rep_len( if ( .lss ) .lshape1.a else .lscale , ncoly))[ 5164 interleave.VGAM(M, M1 = M1)] 5165 temp.names <- if ( .lss ) { 5166 c(scaL.names, sha1.names) 5167 } else { 5168 c(sha1.names, scaL.names) 5169 } 5170 names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)] 5171 5172 misc$earg <- vector("list", M) 5173 names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)] 5174 for (ii in 1:ncoly) 5175 if ( .lss ) { 5176 misc$earg[[M1*ii-1]] <- .escale 5177 misc$earg[[M1*ii ]] <- .eshape1.a 5178 } else { 5179 misc$earg[[M1*ii-1]] <- .eshape1.a 5180 misc$earg[[M1*ii ]] <- .escale 5181 } 5182 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 5183 .escale = escale , .eshape1.a = eshape1.a, 5184 .lss = lss ))), 5185 loglikelihood = eval(substitute( 5186 function(mu, y, w, residuals = FALSE, 5187 eta, extra = NULL, summation = TRUE) { 5188 M1 <- 2 5189 NOS <- ncol(eta)/M1 5190 if ( .lss ) { 5191 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5192 .lscale , earg = .escale ) 5193 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5194 .lshape1.a , earg = .eshape1.a ) 5195 } else { 5196 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5197 .lshape1.a , earg = .eshape1.a ) 5198 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5199 .lscale , earg = .escale ) 5200 } 5201 if (residuals) { 5202 stop("loglikelihood residuals not implemented yet") 5203 } else { 5204 ll.elts <- 5205 c(w) * dfisk(x = y, scale = Scale, shape1.a = aa, log = TRUE) 5206 if (summation) { 5207 sum(ll.elts) 5208 } else { 5209 ll.elts 5210 } 5211 } 5212 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 5213 .escale = escale , .eshape1.a = eshape1.a, 5214 .lss = lss ))), 5215 vfamily = c("fisk"), 5216 simslot = eval(substitute( 5217 function(object, nsim) { 5218 5219 pwts <- if (length(pwts <- object@prior.weights) > 0) 5220 pwts else weights(object, type = "prior") 5221 if (any(pwts != 1)) 5222 warning("ignoring prior weights") 5223 5224 eta <- predict(object) 5225 M1 <- 2 5226 NOS <- ncol(eta)/M1 5227 if ( .lss ) { 5228 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5229 .lscale , earg = .escale ) 5230 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5231 .lshape1.a , earg = .eshape1.a ) 5232 } else { 5233 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5234 .lshape1.a , earg = .eshape1.a ) 5235 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5236 .lscale , earg = .escale ) 5237 } 5238 rfisk(nsim * length(aa), shape1.a = aa, scale = Scale) 5239 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 5240 .escale = escale , .eshape1.a = eshape1.a, 5241 .lss = lss ))), 5242 5243 validparams = eval(substitute(function(eta, y, extra = NULL) { 5244 M1 <- 2 5245 NOS <- ncol(eta) / M1 5246 if ( .lss ) { 5247 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5248 .lscale , earg = .escale ) 5249 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5250 .lshape1.a , earg = .eshape1.a) 5251 } else { 5252 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5253 .lshape1.a , earg = .eshape1.a) 5254 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5255 .lscale , earg = .escale ) 5256 } 5257 parg <- 1 5258 qq <- 1 5259 5260 okay1 <- all(is.finite(Scale)) && all(Scale > 0) && 5261 all(is.finite(aa )) && all(aa > 0) && 5262 all(is.finite(parg )) && all(parg > 0) && 5263 all(is.finite(qq )) && all(qq > 0) 5264 okay.support <- if (okay1) all(-aa * parg < 1 & 1 < aa * qq) else TRUE 5265 if (!okay.support) 5266 warning("parameter constraint -a < 1 < a has been violated; ", 5267 "solution may be at the boundary of the ", 5268 "parameter space.") 5269 okay1 && okay.support 5270 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 5271 .escale = escale , .eshape1.a = eshape1.a, 5272 .lss = lss ))), 5273 5274 5275 5276 5277 deriv = eval(substitute(expression({ 5278 M1 <- 2 5279 NOS <- ncol(eta)/M1 # Needed for summary() 5280 if ( .lss ) { 5281 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5282 .lscale , earg = .escale ) 5283 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5284 .lshape1.a , earg = .eshape1.a) 5285 } else { 5286 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5287 .lshape1.a , earg = .eshape1.a) 5288 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5289 .lscale , earg = .escale ) 5290 } 5291 parg <- 1 5292 qq <- 1 5293 temp1 <- log(y/Scale) 5294 temp2 <- (y/Scale)^aa 5295 5296 dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2)) 5297 dl.da <- 1/aa + parg * temp1 - (parg+qq) * temp1 / (1+1/temp2) 5298 5299 dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale ) 5300 da.deta <- dtheta.deta(aa, .lshape1.a , earg = .eshape1.a ) 5301 5302 myderiv <- if ( .lss ) { 5303 c(w) * cbind(dl.dscale * dscale.deta, 5304 dl.da * da.deta) 5305 } else { 5306 c(w) * cbind(dl.da * da.deta, 5307 dl.dscale * dscale.deta) 5308 } 5309 myderiv[, interleave.VGAM(M, M1 = M1)] 5310 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 5311 .escale = escale , .eshape1.a = eshape1.a, 5312 .lss = lss ))), 5313 weight = eval(substitute(expression({ 5314 ned2l.da <- (1 + (-6 + pi^2) / 9) / aa^2 5315 ned2l.dscale <- ((aa / Scale)^2) / 3 5316 5317 5318 wz <- matrix(0, n, M) # Diagonal 5319 if ( .lss ) { 5320 wz[, c(TRUE, FALSE)] <- ned2l.dscale * dscale.deta^2 5321 wz[, c(FALSE, TRUE)] <- ned2l.da * da.deta^2 5322 } else { 5323 wz[, c(TRUE, FALSE)] <- ned2l.da * da.deta^2 5324 wz[, c(FALSE, TRUE)] <- ned2l.dscale * dscale.deta^2 5325 } 5326 c(w) * wz 5327 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 5328 .escale = escale , .eshape1.a = eshape1.a, 5329 .lss = lss )))) 5330} 5331 5332 5333 5334 5335 5336 5337 5338 inv.lomax <- function(lscale = "loglink", 5339 lshape2.p = "loglink", 5340 iscale = NULL, 5341 ishape2.p = NULL, 5342 imethod = 1, 5343 gscale = exp(-5:5), 5344 gshape2.p = exp(-5:5), 5345 probs.y = c(0.25, 0.50, 0.75), 5346 zero = "shape2.p") { 5347 5348 5349 5350 5351 5352 5353 if (!is.Numeric(imethod, length.arg = 1, 5354 integer.valued = TRUE, 5355 positive = TRUE) || imethod > 2) 5356 stop("Bad input for argument 'imethod'") 5357 5358 if (length(iscale) && !is.Numeric(iscale, positive = TRUE)) 5359 stop("Bad input for argument 'iscale'") 5360 5361 if (length(ishape2.p) && !is.Numeric(ishape2.p, positive = TRUE)) 5362 stop("Bad input for argument 'ishape2.p'") 5363 5364 if (length(probs.y) < 2 || max(probs.y) > 1 || 5365 !is.Numeric(probs.y, positive = TRUE)) 5366 stop("Bad input for argument 'probs.y'") 5367 5368 5369 lscale <- as.list(substitute(lscale)) 5370 escale <- link2list(lscale) 5371 lscale <- attr(escale, "function.name") 5372 5373 lshape2.p <- as.list(substitute(lshape2.p)) 5374 eshape2.p <- link2list(lshape2.p) 5375 lshape2.p <- attr(eshape2.p, "function.name") 5376 5377 5378 new("vglmff", 5379 blurb = 5380 c("Inverse Lomax distribution \n\n", 5381 "Links: ", 5382 namesof("scale" , lscale , earg = escale), ", ", 5383 namesof("shape2.p" , lshape2.p, earg = eshape2.p), "\n", 5384 "Mean: does not exist"), 5385 constraints = eval(substitute(expression({ 5386 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 5387 predictors.names = predictors.names, 5388 M1 = 2) 5389 }), list( .zero = zero ))), 5390 infos = eval(substitute(function(...) { 5391 list(M1 = 2, 5392 Q1 = 1, 5393 expected = TRUE, 5394 zero = .zero , 5395 multipleResponses = TRUE, 5396 parameters.names = c("scale", "shape2.p"), 5397 lscale = .lscale , 5398 escale = .escale , 5399 lshape2.p = .lshape2.p , 5400 eshape2.p = .eshape2.p ) 5401 }, list( .lscale = lscale , 5402 .escale = escale , 5403 .lshape2.p = lshape2.p, 5404 .eshape2.p = eshape2.p, 5405 .zero = zero ))), 5406 initialize = eval(substitute(expression({ 5407 temp5 <- w.y.check(w = w, y = y, 5408 Is.positive.y = TRUE, 5409 ncol.w.max = Inf, 5410 ncol.y.max = Inf, 5411 out.wy = TRUE, 5412 colsyperw = 1, 5413 maximize = TRUE) 5414 y <- temp5$y 5415 w <- temp5$w 5416 M1 <- 2 # Number of parameters for one response 5417 NOS <- ncoly <- ncol(y) 5418 M <- M1*ncol(y) 5419 5420 scaL.names <- param.names("scale", NOS, skip1 = TRUE) 5421 sha2.names <- param.names("shape2.p", NOS, skip1 = TRUE) 5422 5423 predictors.names <- 5424 c(namesof(scaL.names , .lscale , earg = .escale , tag = FALSE), 5425 namesof(sha2.names , .lshape2.p , earg = .eshape2.p , tag = FALSE)) 5426 predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)] 5427 5428 if (!length(etastart)) { 5429 sc.init <- 5430 pp.init <- matrix(NA_real_, n, NOS) 5431 5432 for (spp. in 1:NOS) { # For each response 'y_spp.'... do: 5433 yvec <- y[, spp.] 5434 wvec <- w[, spp.] 5435 5436 if ( .imethod == 1 ) { 5437 gscale <- .gscale 5438 gshape2.p <- .gshape2.p 5439 if (length( .iscale )) gscale <- rep_len( .iscale , NOS) 5440 if (length( .ishape2.p )) gshape2.p <- rep_len( .ishape2.p , NOS) 5441 5442 5443 5444 5445 5446 try.this <- 5447 grid.search4(gscale, vov2 = 1, gshape2.p, vov4 = 1, 5448 objfun = genbetaII.Loglikfun4, # x = x, 5449 y = yvec, w = wvec, # extraargs = NULL, 5450 ret.objfun = TRUE) # Last value is the loglik 5451 sc.init[, spp.] <- try.this["Value1"] 5452 pp.init[, spp.] <- try.this["Value3"] 5453 } else { # .imethod == 2 5454 qvec <- .probs.y 5455 ishape2.p <- if (length( .ishape2.p )) .ishape2.p else 1 5456 xvec <- log( qvec^(-1/ ishape2.p) - 1 ) 5457 fit0 <- lsfit(x = xvec, y = log(quantile(yvec, qvec))) 5458 sc.init[, spp.] <- if (length( .iscale )) .iscale else 5459 exp(fit0$coef[1]) 5460 pp.init[, spp.] <- ishape2.p 5461 } 5462 } # End of for (spp. ...) 5463 5464 5465 etastart <- 5466 cbind(theta2eta(sc.init, .lscale , earg = .escale ), 5467 theta2eta(pp.init, .lshape2.p , earg = .eshape2.p )) 5468 etastart <- etastart[, interleave.VGAM(M, M1 = M1)] 5469 } # End of etastart. 5470 }), list( .lscale = lscale , 5471 .escale = escale , 5472 .iscale = iscale , 5473 .gscale = gscale , 5474 .lshape2.p = lshape2.p, 5475 .eshape2.p = eshape2.p, 5476 .ishape2.p = ishape2.p, 5477 .gshape2.p = gshape2.p, 5478 .imethod = imethod , .probs.y = probs.y 5479 ))), 5480 linkinv = eval(substitute(function(eta, extra = NULL) { 5481 5482 M1 <- 2 5483 NOS <- ncol(eta)/M1 5484 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5485 .lscale , earg = .escale ) 5486 parg <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5487 .lshape2.p , earg = .eshape2.p ) 5488 5489 qinv.lomax(p = 0.5, scale = Scale, shape2.p = parg) 5490 }, list( .lscale = lscale , 5491 .escale = escale , 5492 .lshape2.p = lshape2.p, 5493 .eshape2.p = eshape2.p 5494 ))), 5495 last = eval(substitute(expression({ 5496 M1 <- 2 5497 5498 misc$link <- 5499 c(rep_len( .lscale , ncoly), 5500 rep_len( .lshape2.p , ncoly))[interleave.VGAM(M, M1 = M1)] 5501 temp.names <- c(scaL.names, sha2.names) 5502 names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)] 5503 5504 misc$earg <- vector("list", M) 5505 names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)] 5506 for (ii in 1:ncoly) { 5507 misc$earg[[M1*ii-1]] <- .escale 5508 misc$earg[[M1*ii ]] <- .eshape2.p 5509 } 5510 5511 misc$expected <- TRUE 5512 misc$multipleResponses <- TRUE 5513 }), list( .lscale = lscale , 5514 .escale = escale , 5515 .lshape2.p = lshape2.p, 5516 .eshape2.p = eshape2.p 5517 ))), 5518 loglikelihood = eval(substitute( 5519 function(mu, y, w, residuals = FALSE, 5520 eta, extra = NULL, summation = TRUE) { 5521 M1 <- 2 5522 NOS <- ncol(eta)/M1 5523 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5524 .lscale , earg = .escale ) 5525 parg <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5526 .lshape2.p , earg = .eshape2.p ) 5527 if (residuals) { 5528 stop("loglikelihood residuals not implemented yet") 5529 } else { 5530 ll.elts <- 5531 c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = 1, 5532 shape2.p = parg, shape3.q = 1, log = TRUE) 5533 if (summation) { 5534 sum(ll.elts) 5535 } else { 5536 ll.elts 5537 } 5538 } 5539 }, list( .lscale = lscale , 5540 .escale = escale , 5541 .lshape2.p = lshape2.p, 5542 .eshape2.p = eshape2.p 5543 ))), 5544 vfamily = c("inv.lomax"), 5545 5546 simslot = eval(substitute( 5547 function(object, nsim) { 5548 5549 pwts <- if (length(pwts <- object@prior.weights) > 0) 5550 pwts else weights(object, type = "prior") 5551 if (any(pwts != 1)) 5552 warning("ignoring prior weights") 5553 eta <- predict(object) 5554 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5555 .lscale , earg = .escale ) 5556 parg <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5557 .lshape2.p , earg = .eshape2.p ) 5558 aa <- 1 5559 qq <- 1 5560 rinv.lomax(nsim * length(Scale), scale = Scale, shape2.p = parg) 5561 }, list( .lscale = lscale , 5562 .escale = escale , 5563 .lshape2.p = lshape2.p, 5564 .eshape2.p = eshape2.p 5565 ))), 5566 5567 5568 5569 validparams = eval(substitute(function(eta, y, extra = NULL) { 5570 M1 <- 2 5571 NOS <- ncol(eta) / M1 5572 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5573 .lscale , earg = .escale ) 5574 parg <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5575 .lshape2.p , earg = .eshape2.p ) 5576 qq <- 1 5577 aa <- 1 5578 okay1 <- all(is.finite(Scale)) && all(Scale > 0) && 5579 all(is.finite(aa )) && all(aa > 0) && 5580 all(is.finite(parg )) && all(parg > 0) && 5581 all(is.finite(qq )) && all(qq > 0) 5582 okay.support <- if (okay1) all(-aa * parg < 1 ) else TRUE 5583 if (!okay.support) 5584 warning("parameter constraint -a*p < 1 has been violated; ", 5585 "solution may be at the boundary of the ", 5586 "parameter space.") 5587 okay1 && okay.support 5588 }, list( .lscale = lscale , 5589 .escale = escale , 5590 .lshape2.p = lshape2.p, 5591 .eshape2.p = eshape2.p 5592 ))), 5593 5594 5595 5596 5597 deriv = eval(substitute(expression({ 5598 M1 <- 2 5599 NOS <- ncol(eta)/M1 # Needed for summary() 5600 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5601 .lscale , earg = .escale ) 5602 parg <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5603 .lshape2.p , earg = .eshape2.p ) 5604 qq <- 1 5605 aa <- 1 5606 temp1 <- log(y/Scale) 5607 temp2 <- (y/Scale)^aa 5608 temp3 <- digamma(parg + qq) 5609 temp3a <- digamma(parg) 5610 temp3b <- digamma(qq) 5611 temp4 <- log1p(temp2) 5612 5613 dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2)) 5614 dl.dp <- aa * temp1 + temp3 - temp3a - temp4 5615 5616 dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale ) 5617 dp.deta <- dtheta.deta(parg, .lshape2.p , earg = .eshape2.p ) 5618 5619 myderiv <- 5620 c(w) * cbind(dl.dscale * dscale.deta, 5621 dl.dp * dp.deta) 5622 myderiv[, interleave.VGAM(M, M1 = M1)] 5623 }), list( .lscale = lscale , 5624 .escale = escale , 5625 .lshape2.p = lshape2.p, 5626 .eshape2.p = eshape2.p 5627 ))), 5628 weight = eval(substitute(expression({ 5629 temp5 <- trigamma(parg + qq) 5630 temp5a <- trigamma(parg) 5631 temp5b <- trigamma(qq) 5632 5633 ned2l.dscale <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2) 5634 ned2l.dp <- temp5a - temp5 5635 ned2l.dscalep <- aa * qq / (Scale*(parg+qq)) 5636 wz <- 5637 array(c(c(w) * ned2l.dscale * dscale.deta^2, 5638 c(w) * ned2l.dp * dp.deta^2, 5639 c(w) * ned2l.dscalep * dscale.deta * dp.deta), 5640 dim = c(n, M/M1, M1*(M1+1)/2)) 5641 wz <- arwz2wz(wz, M = M, M1 = M1) 5642 wz 5643 }), list( .lscale = lscale , 5644 .escale = escale , 5645 .lshape2.p = lshape2.p, 5646 .eshape2.p = eshape2.p 5647 )))) 5648} 5649 5650 5651 5652 5653 5654 5655 paralogistic <- 5656 function(lscale = "loglink", 5657 lshape1.a = "loglink", 5658 iscale = NULL, 5659 ishape1.a = NULL, 5660 imethod = 1, 5661 lss = TRUE, 5662 gscale = exp(-5:5), 5663 gshape1.a = seq(0.75, 4, by = 0.25), # exp(-5:5), 5664 probs.y = c(0.25, 0.50, 0.75), 5665 zero = "shape") { 5666 5667 5668 5669 5670 5671 if (length(lss) != 1 && !is.logical(lss)) 5672 stop("Argument 'lss' not specified correctly") 5673 5674 if (!is.Numeric(imethod, length.arg = 1, 5675 integer.valued = TRUE, 5676 positive = TRUE) || imethod > 2) 5677 stop("Bad input for argument 'imethod'") 5678 5679 if (length(iscale) && !is.Numeric(iscale, positive = TRUE)) 5680 stop("Bad input for argument 'iscale'") 5681 5682 if (length(ishape1.a) && !is.Numeric(ishape1.a, positive = TRUE)) 5683 stop("Bad input for argument 'ishape1.a'") 5684 5685 if (length(probs.y) < 2 || max(probs.y) > 1 || 5686 !is.Numeric(probs.y, positive = TRUE)) 5687 stop("Bad input for argument 'probs.y'") 5688 5689 5690 lscale <- as.list(substitute(lscale)) 5691 escale <- link2list(lscale) 5692 lscale <- attr(escale, "function.name") 5693 5694 lshape1.a <- as.list(substitute(lshape1.a)) 5695 eshape1.a <- link2list(lshape1.a) 5696 lshape1.a <- attr(eshape1.a, "function.name") 5697 5698 5699 new("vglmff", 5700 blurb = 5701 c("Paralogistic distribution \n\n", 5702 "Links: ", 5703 ifelse (lss, 5704 namesof("scale" , lscale , earg = escale), 5705 namesof("shape1.a", lshape1.a, earg = eshape1.a)), ", ", 5706 ifelse (lss, 5707 namesof("shape1.a", lshape1.a, earg = eshape1.a), 5708 namesof("scale" , lscale , earg = escale)), "\n", 5709 "Mean: scale * gamma(1 + 1/shape1.a) * ", 5710 "gamma(shape1.a - 1/shape1.a) / ", 5711 "gamma(shape1.a)"), 5712 constraints = eval(substitute(expression({ 5713 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 5714 predictors.names = predictors.names, 5715 M1 = 2) 5716 }), list( .zero = zero ))), 5717 infos = eval(substitute(function(...) { 5718 list(M1 = 2, 5719 Q1 = 1, 5720 expected = TRUE, 5721 zero = .zero , 5722 multipleResponses = TRUE, 5723 parameters.names = if ( .lss ) 5724 c("scale", "shape1.a") else 5725 c("shape1.a", "scale"), 5726 lscale = .lscale , lshape1.a = .lshape1.a , 5727 escale = .escale , eshape1.a = .eshape1.a ) 5728 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 5729 .escale = escale , .eshape1.a = eshape1.a, 5730 .lss = lss , 5731 .zero = zero ))), 5732 initialize = eval(substitute(expression({ 5733 temp5 <- w.y.check(w = w, y = y, 5734 Is.positive.y = TRUE, 5735 ncol.w.max = Inf, 5736 ncol.y.max = Inf, 5737 out.wy = TRUE, 5738 colsyperw = 1, 5739 maximize = TRUE) 5740 y <- temp5$y 5741 w <- temp5$w 5742 M1 <- 2 # Number of parameters for one response 5743 NOS <- ncoly <- ncol(y) 5744 M <- M1*ncol(y) 5745 5746 scaL.names <- param.names("scale", NOS, skip1 = TRUE) 5747 sha1.names <- param.names("shape1.a", NOS, skip1 = TRUE) 5748 5749 predictors.names <- 5750 if ( .lss ) { 5751 c(namesof(scaL.names , .lscale , earg = .escale , 5752 tag = FALSE), 5753 namesof(sha1.names , .lshape1.a , earg = .eshape1.a , 5754 tag = FALSE)) 5755 } else { 5756 c(namesof(sha1.names , .lshape1.a , earg = .eshape1.a , 5757 tag = FALSE), 5758 namesof(scaL.names , .lscale , earg = .escale , 5759 tag = FALSE)) 5760 } 5761 predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)] 5762 5763 if (!length(etastart)) { 5764 sc.init <- 5765 aa.init <- matrix(NA_real_, n, NOS) 5766 5767 for (spp. in 1:NOS) { # For each response 'y_spp.'... do: 5768 yvec <- y[, spp.] 5769 wvec <- w[, spp.] 5770 5771 if ( .imethod == 1 ) { 5772 gscale <- .gscale 5773 gshape1.a <- .gshape1.a 5774 if (length( .iscale )) gscale <- rep_len( .iscale , NOS) 5775 if (length( .ishape1.a )) gshape1.a <- rep_len( .ishape1.a , NOS) 5776 5777 5778 paralogistic.Loglikfun2 <- 5779 function(scaleval, shape1.a, 5780 y, x, w, extraargs) { 5781 sum(c(w) * dgenbetaII(x = y, scale = scaleval, 5782 shape1.a = shape1.a, 5783 shape2.p = 1, 5784 shape3.q = shape1.a, log = TRUE)) 5785 } 5786 5787 try.this <- 5788 grid.search2(gscale, gshape1.a, # vov3 = 1, vov4 = gshape1.a, 5789 objfun = paralogistic.Loglikfun2, # x = x, 5790 y = yvec, w = wvec, # extraargs = NULL, 5791 ret.objfun = TRUE) # Last value is the loglik 5792 sc.init[, spp.] <- try.this["Value1"] 5793 aa.init[, spp.] <- try.this["Value2"] 5794 } else { # .imethod == 2 5795 qvec <- .probs.y 5796 ishape3.q <- if (length( .ishape1.a )) .ishape1.a else 1 5797 xvec <- log( (1-qvec)^(-1/ ishape3.q) - 1 ) 5798 fit0 <- lsfit(x = xvec, y = log(quantile(yvec, qvec))) 5799 sc.init[, spp.] <- if (length( .iscale )) .iscale else 5800 exp(fit0$coef[1]) 5801 aa.init[, spp.] <- if (length( .ishape1.a )) .ishape1.a else 5802 1/fit0$coef[2] 5803 } 5804 } # End of for (spp. ...) 5805 5806 finite.mean <- 1 < aa.init * aa.init 5807 COP.use <- 1.15 5808 while (FALSE && any(!finite.mean)) { 5809 aa.init[!finite.mean] <- 0.1 + aa.init[!finite.mean] * COP.use 5810 finite.mean <- 1 < aa.init * aa.init 5811 } 5812 5813 etastart <- if ( .lss ) 5814 cbind(theta2eta(sc.init, .lscale , earg = .escale ), 5815 theta2eta(aa.init, .lshape1.a , earg = .eshape1.a )) else 5816 cbind(theta2eta(aa.init, .lshape1.a , earg = .eshape1.a ), 5817 theta2eta(sc.init, .lscale , earg = .escale )) 5818 etastart <- etastart[, interleave.VGAM(M, M1 = M1)] 5819 } # End of etastart. 5820 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 5821 .escale = escale , .eshape1.a = eshape1.a, 5822 .iscale = iscale , .ishape1.a = ishape1.a, 5823 .gscale = gscale , .gshape1.a = gshape1.a, 5824 .imethod = imethod , .probs.y = probs.y, 5825 .lss = lss ))), 5826 linkinv = eval(substitute(function(eta, extra = NULL) { 5827 M1 <- 2 5828 NOS <- ncol(eta)/M1 5829 if ( .lss ) { 5830 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5831 .lscale , earg = .escale ) 5832 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5833 .lshape1.a , earg = .eshape1.a ) 5834 } else { 5835 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5836 .lshape1.a , earg = .eshape1.a ) 5837 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5838 .lscale , earg = .escale ) 5839 } 5840 parg <- 1 5841 qq <- aa 5842 5843 ans <- Scale * exp(lgamma(parg + 1/aa) + 5844 lgamma(qq - 1/aa) - lgamma(parg) - lgamma(qq)) 5845 ans[parg + 1/aa <= 0] <- NA 5846 ans[qq - 1/aa <= 0] <- NA 5847 ans[aa <= 0] <- NA 5848 ans[Scale <= 0] <- NA 5849 ans 5850 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 5851 .escale = escale , .eshape1.a = eshape1.a, 5852 .lss = lss ))), 5853 last = eval(substitute(expression({ 5854 M1 <- 2 5855 5856 misc$link <- c(rep_len(if ( .lss ) .lscale else .lshape1.a , ncoly), 5857 rep_len(if ( .lss ) .lshape1.a else .lscale , ncoly))[ 5858 interleave.VGAM(M, M1 = M1)] 5859 temp.names <- if ( .lss ) { 5860 c(scaL.names, sha1.names) 5861 } else { 5862 c(sha1.names, scaL.names) 5863 } 5864 names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)] 5865 5866 misc$earg <- vector("list", M) 5867 names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)] 5868 for (ii in 1:ncoly) 5869 if ( .lss ) { 5870 misc$earg[[M1*ii-1]] <- .escale 5871 misc$earg[[M1*ii ]] <- .eshape1.a 5872 } else { 5873 misc$earg[[M1*ii-1]] <- .eshape1.a 5874 misc$earg[[M1*ii ]] <- .escale 5875 } 5876 5877 misc$expected <- TRUE 5878 misc$multipleResponses <- TRUE 5879 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 5880 .escale = escale , .eshape1.a = eshape1.a, 5881 .lss = lss ))), 5882 loglikelihood = eval(substitute( 5883 function(mu, y, w, residuals = FALSE, 5884 eta, extra = NULL, summation = TRUE) { 5885 M1 <- 2 5886 NOS <- ncol(eta)/M1 5887 if ( .lss ) { 5888 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5889 .lscale , earg = .escale ) 5890 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5891 .lshape1.a , earg = .eshape1.a ) 5892 } else { 5893 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5894 .lshape1.a , earg = .eshape1.a ) 5895 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5896 .lscale , earg = .escale ) 5897 } 5898 parg <- 1 5899 qq <- aa 5900 if (residuals) { 5901 stop("loglikelihood residuals not implemented yet") 5902 } else { 5903 ll.elts <- 5904 c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = aa, 5905 shape2.p = parg, shape3.q = aa, log = TRUE) 5906 if (summation) { 5907 sum(ll.elts) 5908 } else { 5909 ll.elts 5910 } 5911 } 5912 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 5913 .escale = escale , .eshape1.a = eshape1.a, 5914 .lss = lss ))), 5915 vfamily = c("paralogistic"), 5916 simslot = eval(substitute( 5917 function(object, nsim) { 5918 5919 pwts <- if (length(pwts <- object@prior.weights) > 0) 5920 pwts else weights(object, type = "prior") 5921 if (any(pwts != 1)) 5922 warning("ignoring prior weights") 5923 5924 eta <- predict(object) 5925 M1 <- 2 5926 NOS <- ncol(eta)/M1 5927 if ( .lss ) { 5928 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5929 .lscale , earg = .escale ) 5930 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5931 .lshape1.a , earg = .eshape1.a ) 5932 } else { 5933 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5934 .lshape1.a , earg = .eshape1.a ) 5935 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5936 .lscale , earg = .escale ) 5937 } 5938 rparalogistic(nsim * length(Scale), shape1.a = aa, scale = Scale) 5939 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 5940 .escale = escale , .eshape1.a = eshape1.a, 5941 .lss = lss ))), 5942 5943 validparams = eval(substitute(function(eta, y, extra = NULL) { 5944 M1 <- 2 5945 NOS <- ncol(eta) / M1 5946 if ( .lss ) { 5947 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5948 .lscale , earg = .escale ) 5949 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5950 .lshape1.a , earg = .eshape1.a) 5951 } else { 5952 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5953 .lshape1.a , earg = .eshape1.a) 5954 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5955 .lscale , earg = .escale ) 5956 } 5957 parg <- 1 5958 qq <- aa 5959 okay1 <- all(is.finite(Scale)) && all(Scale > 0) && 5960 all(is.finite(aa )) && all(aa > 0) && 5961 all(is.finite(parg )) && all(parg > 0) && 5962 all(is.finite(qq )) && all(qq > 0) 5963 okay.support <- if (okay1) all(-aa * parg < 1 & 1 < aa * qq) else TRUE 5964 if (!okay.support) 5965 warning("parameter constraint -a < 1 < a*a has been violated; ", 5966 "solution may be at the boundary of the ", 5967 "parameter space.") 5968 okay1 && okay.support 5969 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 5970 .escale = escale , .eshape1.a = eshape1.a, 5971 .lss = lss ))), 5972 5973 5974 5975 5976 deriv = eval(substitute(expression({ 5977 M1 <- 2 5978 NOS <- ncol(eta)/M1 # Needed for summary() 5979 if ( .lss ) { 5980 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5981 .lscale , earg = .escale ) 5982 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5983 .lshape1.a , earg = .eshape1.a) 5984 } else { 5985 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 5986 .lshape1.a , earg = .eshape1.a) 5987 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 5988 .lscale , earg = .escale ) 5989 } 5990 parg <- 1 5991 qq <- aa 5992 5993 temp1 <- log(y/Scale) 5994 temp2 <- (y/Scale)^aa 5995 temp3 <- digamma(parg + qq) 5996 temp3a <- digamma(parg) 5997 temp3b <- digamma(qq) 5998 temp4 <- log1p(temp2) 5999 6000 dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2)) 6001 dl.da <- 1/aa + parg * temp1 - (parg+qq) * temp1 / (1+1/temp2) 6002 6003 dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale ) 6004 da.deta <- dtheta.deta(aa, .lshape1.a , earg = .eshape1.a ) 6005 6006 myderiv <- if ( .lss ) { 6007 c(w) * cbind(dl.dscale * dscale.deta, 6008 dl.da * da.deta) 6009 } else { 6010 c(w) * cbind(dl.da * da.deta, 6011 dl.dscale * dscale.deta) 6012 } 6013 myderiv[, interleave.VGAM(M, M1 = M1)] 6014 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 6015 .escale = escale , .eshape1.a = eshape1.a, 6016 .lss = lss ))), 6017 weight = eval(substitute(expression({ 6018 temp5 <- trigamma(parg + qq) 6019 temp5a <- trigamma(parg) 6020 temp5b <- trigamma(qq) 6021 6022 ned2l.da <- (1 + parg + qq + parg * qq * (temp5a + temp5b + 6023 (temp3b - temp3a + (parg-qq)/(parg*qq))^2 - 6024 (parg^2 + qq^2) / (parg*qq)^2)) / (aa^2 * (1+parg+qq)) 6025 ned2l.dscale <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2) 6026 ned2l.dascale <- (parg - qq - parg * qq * 6027 (temp3a -temp3b)) / (Scale*(1 + parg+qq)) 6028 wz <- if ( .lss ) { 6029 array(c(c(w) * ned2l.dscale * dscale.deta^2, 6030 c(w) * ned2l.da * da.deta^2, 6031 c(w) * ned2l.dascale * da.deta * dscale.deta), 6032 dim = c(n, M/M1, M1*(M1+1)/2)) 6033 } else { 6034 array(c(c(w) * ned2l.da * da.deta^2, 6035 c(w) * ned2l.dscale * dscale.deta^2, 6036 c(w) * ned2l.dascale * da.deta * dscale.deta), 6037 dim = c(n, M/M1, M1*(M1+1)/2)) 6038 } 6039 wz <- arwz2wz(wz, M = M, M1 = M1) 6040 wz 6041 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 6042 .escale = escale , .eshape1.a = eshape1.a, 6043 .lss = lss )))) 6044} 6045 6046 6047 6048 6049 6050 6051 6052 6053 6054 6055 6056 inv.paralogistic <- 6057 function(lscale = "loglink", 6058 lshape1.a = "loglink", 6059 iscale = NULL, 6060 ishape1.a = NULL, 6061 imethod = 1, 6062 lss = TRUE, 6063 gscale = exp(-5:5), 6064 gshape1.a = seq(0.75, 4, by = 0.25), # exp(-5:5), 6065 probs.y = c(0.25, 0.50, 0.75), 6066 zero = "shape") { 6067 6068 6069 6070 6071 if (length(lss) != 1 && !is.logical(lss)) 6072 stop("Argument 'lss' not specified correctly") 6073 6074 if (!is.Numeric(imethod, length.arg = 1, 6075 integer.valued = TRUE, 6076 positive = TRUE) || imethod > 2) 6077 stop("Bad input for argument 'imethod'") 6078 6079 if (length(iscale) && !is.Numeric(iscale, positive = TRUE)) 6080 stop("Bad input for argument 'iscale'") 6081 6082 if (length(ishape1.a) && !is.Numeric(ishape1.a, positive = TRUE)) 6083 stop("Bad input for argument 'ishape1.a'") 6084 6085 if (length(probs.y) < 2 || max(probs.y) > 1 || 6086 !is.Numeric(probs.y, positive = TRUE)) 6087 stop("Bad input for argument 'probs.y'") 6088 6089 6090 lscale <- as.list(substitute(lscale)) 6091 escale <- link2list(lscale) 6092 lscale <- attr(escale, "function.name") 6093 6094 lshape1.a <- as.list(substitute(lshape1.a)) 6095 eshape1.a <- link2list(lshape1.a) 6096 lshape1.a <- attr(eshape1.a, "function.name") 6097 6098 6099 new("vglmff", 6100 blurb = 6101 c("Inverse paralogistic distribution \n\n", 6102 "Links: ", 6103 ifelse (lss, 6104 namesof("scale" , lscale , earg = escale), 6105 namesof("shape1.a", lshape1.a, earg = eshape1.a)), ", ", 6106 ifelse (lss, 6107 namesof("shape1.a", lshape1.a, earg = eshape1.a), 6108 namesof("scale" , lscale , earg = escale)), "\n", 6109 "Mean: scale * gamma(shape1.a + 1/shape1.a) * ", 6110 "gamma(1 - 1/shape1.a) / ", 6111 "gamma(shape1.a)"), 6112 constraints = eval(substitute(expression({ 6113 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 6114 predictors.names = predictors.names, 6115 M1 = 2) 6116 }), list( .zero = zero ))), 6117 infos = eval(substitute(function(...) { 6118 list(M1 = 2, 6119 Q1 = 1, 6120 expected = TRUE, 6121 zero = .zero , 6122 multipleResponses = TRUE, 6123 parameters.names = if ( .lss ) 6124 c("scale", "shape1.a") else 6125 c("shape1.a", "scale"), 6126 lscale = .lscale , lshape1.a = .lshape1.a , 6127 escale = .escale , eshape1.a = .eshape1.a ) 6128 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 6129 .escale = escale , .eshape1.a = eshape1.a, 6130 .lss = lss , 6131 .zero = zero ))), 6132 initialize = eval(substitute(expression({ 6133 temp5 <- w.y.check(w = w, y = y, 6134 Is.positive.y = TRUE, 6135 ncol.w.max = Inf, 6136 ncol.y.max = Inf, 6137 out.wy = TRUE, 6138 colsyperw = 1, 6139 maximize = TRUE) 6140 y <- temp5$y 6141 w <- temp5$w 6142 M1 <- 2 # Number of parameters for one response 6143 NOS <- ncoly <- ncol(y) 6144 M <- M1*ncol(y) 6145 6146 scaL.names <- param.names("scale", NOS, skip1 = TRUE) 6147 sha1.names <- param.names("shape1.a", NOS, skip1 = TRUE) 6148 6149 predictors.names <- if ( .lss ) { 6150 c(namesof(scaL.names , .lscale , earg = .escale , tag = FALSE), 6151 namesof(sha1.names , .lshape1.a , earg = .eshape1.a , tag = FALSE)) 6152 } else { 6153 c(namesof(sha1.names , .lshape1.a , earg = .eshape1.a , tag = FALSE), 6154 namesof(scaL.names , .lscale , earg = .escale , tag = FALSE)) 6155 } 6156 predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)] 6157 6158 if (!length(etastart)) { 6159 sc.init <- 6160 aa.init <- matrix(NA_real_, n, NOS) 6161 6162 for (spp. in 1:NOS) { # For each response 'y_spp.'... do: 6163 yvec <- y[, spp.] 6164 wvec <- w[, spp.] 6165 6166 if ( .imethod == 1 ) { 6167 gscale <- .gscale 6168 gshape1.a <- .gshape1.a 6169 if (length( .iscale )) gscale <- rep_len( .iscale , NOS) 6170 if (length( .ishape1.a )) gshape1.a <- rep_len( .ishape1.a , NOS) 6171 6172 6173 inv.paralogistic.Loglikfun2 <- 6174 function(scaleval, shape1.a, 6175 y, x, w, extraargs) { 6176 sum(c(w) * dgenbetaII(x = y, scale = scaleval, 6177 shape1.a = shape1.a, 6178 shape2.p = shape1.a, 6179 shape3.q = 1, log = TRUE)) 6180 } 6181 6182 6183 try.this <- 6184 grid.search2(gscale, gshape1.a, # vov3 = 1, vov4 = gshape1.a, 6185 objfun = inv.paralogistic.Loglikfun2, # x = x, 6186 y = yvec, w = wvec, # extraargs = NULL, 6187 ret.objfun = TRUE) # Last value is the loglik 6188 sc.init[, spp.] <- try.this["Value1"] 6189 aa.init[, spp.] <- try.this["Value2"] 6190 } else { # .imethod == 2 6191 qvec <- .probs.y 6192 ishape2.p <- if (length( .ishape1.a )) .ishape1.a else 1 6193 xvec <- log( qvec^(-1/ ishape2.p) - 1 ) 6194 fit0 <- lsfit(x = xvec, y = log(quantile(yvec, qvec))) 6195 sc.init[, spp.] <- if (length( .iscale )) .iscale else 6196 exp(fit0$coef[1]) 6197 aa.init[, spp.] <- if (length( .ishape1.a )) .ishape1.a else 6198 -1/fit0$coef[2] 6199 } 6200 } # End of for (spp. ...) 6201 6202 6203 6204 finite.mean <- 1 < aa.init 6205 COP.use <- 1.15 6206 while (FALSE && any(!finite.mean)) { 6207 aa.init[!finite.mean] <- 0.1 + aa.init[!finite.mean] * COP.use 6208 finite.mean <- 1 < aa.init 6209 } 6210 6211 etastart <- if ( .lss ) 6212 cbind(theta2eta(sc.init, .lscale , earg = .escale ), 6213 theta2eta(aa.init, .lshape1.a , earg = .eshape1.a )) else 6214 cbind(theta2eta(aa.init, .lshape1.a , earg = .eshape1.a ), 6215 theta2eta(sc.init, .lscale , earg = .escale )) 6216 etastart <- etastart[, interleave.VGAM(M, M1 = M1)] 6217 } # End of etastart. 6218 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 6219 .escale = escale , .eshape1.a = eshape1.a, 6220 .iscale = iscale , .ishape1.a = ishape1.a, 6221 .gscale = gscale , .gshape1.a = gshape1.a, 6222 .imethod = imethod , .probs.y = probs.y, 6223 .lss = lss ))), 6224 linkinv = eval(substitute(function(eta, extra = NULL) { 6225 M1 <- 2 6226 NOS <- ncol(eta)/M1 6227 if ( .lss ) { 6228 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 6229 .lscale , earg = .escale ) 6230 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 6231 .lshape1.a , earg = .eshape1.a ) 6232 } else { 6233 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 6234 .lshape1.a , earg = .eshape1.a ) 6235 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 6236 .lscale , earg = .escale ) 6237 } 6238 parg <- aa 6239 qq <- 1 6240 6241 ans <- Scale * exp(lgamma(parg + 1/aa) + 6242 lgamma(qq - 1/aa) - lgamma(parg) - lgamma(qq)) 6243 ans[parg + 1/aa <= 0] <- NA 6244 ans[qq - 1/aa <= 0] <- NA 6245 ans[aa <= 0] <- NA 6246 ans[Scale <= 0] <- NA 6247 ans 6248 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 6249 .escale = escale , .eshape1.a = eshape1.a, 6250 .lss = lss ))), 6251 last = eval(substitute(expression({ 6252 M1 <- 2 6253 6254 misc$link <- c(rep_len(if ( .lss ) .lscale else .lshape1.a , ncoly), 6255 rep_len(if ( .lss ) .lshape1.a else .lscale , ncoly))[ 6256 interleave.VGAM(M, M1 = M1)] 6257 temp.names <- if ( .lss ) { 6258 c(scaL.names, sha1.names) 6259 } else { 6260 c(sha1.names, scaL.names) 6261 } 6262 names(misc$link) <- temp.names[interleave.VGAM(M, M1 = M1)] 6263 6264 misc$earg <- vector("list", M) 6265 names(misc$earg) <- temp.names[interleave.VGAM(M, M1 = M1)] 6266 for (ii in 1:ncoly) 6267 if ( .lss ) { 6268 misc$earg[[M1*ii-1]] <- .escale 6269 misc$earg[[M1*ii ]] <- .eshape1.a 6270 } else { 6271 misc$earg[[M1*ii-1]] <- .eshape1.a 6272 misc$earg[[M1*ii ]] <- .escale 6273 } 6274 6275 misc$expected <- TRUE 6276 misc$multipleResponses <- TRUE 6277 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 6278 .escale = escale , .eshape1.a = eshape1.a, 6279 .lss = lss ))), 6280 loglikelihood = eval(substitute( 6281 function(mu, y, w, residuals = FALSE, 6282 eta, extra = NULL, summation = TRUE) { 6283 M1 <- 2 6284 NOS <- ncol(eta)/M1 6285 if ( .lss ) { 6286 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 6287 .lscale , earg = .escale ) 6288 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 6289 .lshape1.a , earg = .eshape1.a ) 6290 } else { 6291 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 6292 .lshape1.a , earg = .eshape1.a ) 6293 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 6294 .lscale , earg = .escale ) 6295 } 6296 parg <- aa 6297 if (residuals) { 6298 stop("loglikelihood residuals not implemented yet") 6299 } else { 6300 ll.elts <- 6301 c(w) * dgenbetaII(x = y, scale = Scale, shape1.a = aa, 6302 shape2.p = aa, shape3.q = 1, log = TRUE) 6303 if (summation) { 6304 sum(ll.elts) 6305 } else { 6306 ll.elts 6307 } 6308 } 6309 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 6310 .escale = escale , .eshape1.a = eshape1.a, 6311 .lss = lss ))), 6312 vfamily = c("inv.paralogistic"), 6313 6314 simslot = eval(substitute( 6315 function(object, nsim) { 6316 6317 pwts <- if (length(pwts <- object@prior.weights) > 0) 6318 pwts else weights(object, type = "prior") 6319 if (any(pwts != 1)) 6320 warning("ignoring prior weights") 6321 eta <- predict(object) 6322 if ( .lss ) { 6323 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 6324 .lscale , earg = .escale ) 6325 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 6326 .lshape1.a , earg = .eshape1.a) 6327 } else { 6328 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 6329 .lshape1.a , earg = .eshape1.a) 6330 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 6331 .lscale , earg = .escale ) 6332 } 6333 parg <- aa 6334 qq <- 1 6335 rinv.paralogistic(nsim * length(Scale), shape1.a = aa, scale = Scale) 6336 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 6337 .escale = escale , .eshape1.a = eshape1.a, 6338 .lss = lss ))), 6339 6340 6341 6342 validparams = eval(substitute(function(eta, y, extra = NULL) { 6343 M1 <- 2 6344 NOS <- ncol(eta) / M1 6345 if ( .lss ) { 6346 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 6347 .lscale , earg = .escale ) 6348 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 6349 .lshape1.a , earg = .eshape1.a) 6350 } else { 6351 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 6352 .lshape1.a , earg = .eshape1.a) 6353 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 6354 .lscale , earg = .escale ) 6355 } 6356 parg <- aa 6357 qq <- 1 6358 okay1 <- all(is.finite(Scale)) && all(Scale > 0) && 6359 all(is.finite(aa )) && all(aa > 0) && 6360 all(is.finite(parg )) && all(parg > 0) && 6361 all(is.finite(qq )) && all(qq > 0) 6362 okay.support <- if (okay1) all(-aa * parg < 1 & 1 < aa * qq) else TRUE 6363 if (!okay.support) 6364 warning("parameter constraint -a*a < 1 < a has been violated; ", 6365 "solution may be at the boundary of the ", 6366 "parameter space.") 6367 okay1 && okay.support 6368 }, list( .lscale = lscale , .lshape1.a = lshape1.a, 6369 .escale = escale , .eshape1.a = eshape1.a, 6370 .lss = lss ))), 6371 6372 6373 6374 6375 deriv = eval(substitute(expression({ 6376 M1 <- 2 6377 NOS <- ncol(eta)/M1 # Needed for summary() 6378 if ( .lss ) { 6379 Scale <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 6380 .lscale , earg = .escale ) 6381 aa <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 6382 .lshape1.a , earg = .eshape1.a) 6383 } else { 6384 aa <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 6385 .lshape1.a , earg = .eshape1.a) 6386 Scale <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 6387 .lscale , earg = .escale ) 6388 } 6389 parg <- aa 6390 qq <- 1 6391 temp1 <- log(y/Scale) 6392 temp2 <- (y/Scale)^aa 6393 temp3 <- digamma(parg + qq) 6394 temp3a <- digamma(parg) 6395 temp3b <- digamma(qq) 6396 temp4 <- log1p(temp2) 6397 6398 dl.dscale <- (aa/Scale) * (-parg + (parg+qq) / (1+1/temp2)) 6399 dl.da <- 1/aa + parg * temp1 - (parg+qq) * temp1 / (1+1/temp2) 6400 6401 dscale.deta <- dtheta.deta(Scale, .lscale , earg = .escale ) 6402 da.deta <- dtheta.deta(aa, .lshape1.a , earg = .eshape1.a ) 6403 6404 myderiv <- if ( .lss ) { 6405 c(w) * cbind(dl.dscale * dscale.deta, 6406 dl.da * da.deta) 6407 } else { 6408 c(w) * cbind(dl.da * da.deta, 6409 dl.dscale * dscale.deta) 6410 } 6411 myderiv[, interleave.VGAM(M, M1 = M1)] 6412 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 6413 .escale = escale , .eshape1.a = eshape1.a, 6414 .lss = lss ))), 6415 weight = eval(substitute(expression({ 6416 temp5 <- trigamma(parg + qq) 6417 temp5a <- trigamma(parg) 6418 temp5b <- trigamma(qq) 6419 6420 ned2l.da <- (1 + parg + qq + parg * qq * (temp5a + temp5b + 6421 (temp3b - temp3a + (parg-qq)/(parg*qq))^2 - 6422 (parg^2 + qq^2) / (parg*qq)^2)) / (aa^2 * (1+parg+qq)) 6423 ned2l.dscale <- (aa^2) * parg * qq / ((1+parg+qq) * Scale^2) 6424 ned2l.dascale <- (parg - qq - parg * qq * 6425 (temp3a -temp3b)) / (Scale*(1 + parg+qq)) 6426 wz <- if ( .lss ) { 6427 array(c(c(w) * ned2l.dscale * dscale.deta^2, 6428 c(w) * ned2l.da * da.deta^2, 6429 c(w) * ned2l.dascale * da.deta * dscale.deta), 6430 dim = c(n, M/M1, M1*(M1+1)/2)) 6431 } else { 6432 array(c(c(w) * ned2l.da * da.deta^2, 6433 c(w) * ned2l.dscale * dscale.deta^2, 6434 c(w) * ned2l.dascale * da.deta * dscale.deta), 6435 dim = c(n, M/M1, M1*(M1+1)/2)) 6436 } 6437 wz <- arwz2wz(wz, M = M, M1 = M1) 6438 wz 6439 }), list( .lscale = lscale , .lshape1.a = lshape1.a, 6440 .escale = escale , .eshape1.a = eshape1.a, 6441 .lss = lss )))) 6442} 6443 6444 6445 6446 6447 6448 6449 6450 6451 6452 6453 6454 6455 6456 6457 6458 6459 6460 6461