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 dtrinorm <- 17 function(x1, x2, x3, mean1 = 0, mean2 = 0, mean3 = 0, 18 var1 = 1, var2 = 1, var3 = 1, 19 cov12 = 0, cov23 = 0, cov13 = 0, 20 log = FALSE) { 21 22 if (!is.logical(log.arg <- log) || length(log) != 1) 23 stop("bad input for argument 'log'") 24 rm(log) 25 26 M <- 3 27 n <- max(length(x1), length(x2), length(x3), 28 length(mean1), length(mean2), length(mean3), 29 length(var1 ), length(var2 ), length(var3 ), 30 length(cov12), length(cov23), length(cov13)) 31 32 sd1 <- sqrt(var1) 33 sd2 <- sqrt(var2) 34 sd3 <- sqrt(var3) 35 rho12 <- cov12 / (sd1 * sd2) 36 rho13 <- cov13 / (sd1 * sd3) 37 rho23 <- cov23 / (sd2 * sd3) 38 39 bbb <- 1 - rho12^2 - rho13^2 - rho23^2 + 2 * rho12 * rho13 * rho23 40 logdet <- 2 * (log(sd1) + log(sd2) + log(sd3)) + log(bbb) 41 42 Sigmainv <- matrix(0, n, dimm(3)) # sum(3:1) 43 Sigmainv[, iam(1, 1, M = M)] <- (1 - rho23^2) / (bbb * sd1^2) 44 Sigmainv[, iam(2, 2, M = M)] <- (1 - rho13^2) / (bbb * sd2^2) 45 Sigmainv[, iam(3, 3, M = M)] <- (1 - rho12^2) / (bbb * sd3^2) 46 Sigmainv[, iam(1, 2, M = M)] <- (rho13 * rho23 - rho12) / ( 47 sd1 * sd2 * bbb) 48 Sigmainv[, iam(2, 3, M = M)] <- (rho12 * rho13 - rho23) / ( 49 sd2 * sd3 * bbb) 50 Sigmainv[, iam(1, 3, M = M)] <- (rho12 * rho23 - rho13) / ( 51 sd1 * sd3 * bbb) 52 53 ymatt <- rbind(x1 - mean1, x2 - mean2, x3 - mean3) 54 dim(ymatt) <- c(nrow(ymatt), 1, ncol(ymatt)) # For mux5() 55 qform <- mux5(x = ymatt, cc = Sigmainv, M = 3, matrix.arg = TRUE) 56 57 logpdf <- -1.5 * log(2 * pi) - 0.5 * logdet - 0.5 * c(qform) 58 logpdf[is.infinite(x1) | is.infinite(x2) | is.infinite(x3)] <- log(0) 59 60 if (log.arg) logpdf else exp(logpdf) 61} 62 63 64 65rtrinorm <- function(n, mean1 = 0, mean2 = 0, mean3 = 0, 66 var1 = 1, var2 = 1, var3 = 1, 67 cov12 = 0, cov23 = 0, cov13 = 0) { 68 69 Y1 <- rnorm(n, mean1, sqrt(var1)) 70 ans2 <- rbinorm(n, 71 mean1 = mean2 + cov12 * (Y1 - mean1) / var1, 72 mean2 = mean3 + cov13 * (Y1 - mean1) / var1, 73 var1 = var2 - cov12 * cov12 / var1, 74 var2 = var3 - cov13 * cov13 / var1, 75 cov12 = cov23 - cov12 * cov13 / var1) 76 77 ans <- cbind(Y1, ans2) 78 colnames(ans) <- paste("X", 1:3, sep = "") 79 ans 80} 81 82 83 84 85 86trinormal.control <- 87 function(summary.HDEtest = FALSE, # Overwrites the summary() default. 88 ...) { 89 list(summary.HDEtest = summary.HDEtest) 90} 91 92 93 94 trinormal <- 95 function( 96 zero = c("sd", "rho"), 97 eq.mean = FALSE, 98 eq.sd = FALSE, 99 eq.cor = FALSE, 100 lmean1 = "identitylink", 101 lmean2 = "identitylink", 102 lmean3 = "identitylink", 103 lsd1 = "loglink", 104 lsd2 = "loglink", 105 lsd3 = "loglink", 106 lrho12 = "rhobitlink", 107 lrho23 = "rhobitlink", 108 lrho13 = "rhobitlink", 109 imean1 = NULL, imean2 = NULL, imean3 = NULL, 110 isd1 = NULL, isd2 = NULL, isd3 = NULL, 111 irho12 = NULL, irho23 = NULL, irho13 = NULL, 112 imethod = 1) { 113 114 115 lmean1 <- as.list(substitute(lmean1)) 116 emean1 <- link2list(lmean1) 117 lmean1 <- attr(emean1, "function.name") 118 119 lmean2 <- as.list(substitute(lmean2)) 120 emean2 <- link2list(lmean2) 121 lmean2 <- attr(emean2, "function.name") 122 123 lmean3 <- as.list(substitute(lmean3)) 124 emean3 <- link2list(lmean3) 125 lmean3 <- attr(emean3, "function.name") 126 127 lsd1 <- as.list(substitute(lsd1)) 128 esd1 <- link2list(lsd1) 129 lsd1 <- attr(esd1, "function.name") 130 131 lsd2 <- as.list(substitute(lsd2)) 132 esd2 <- link2list(lsd2) 133 lsd2 <- attr(esd2, "function.name") 134 135 lsd3 <- as.list(substitute(lsd3)) 136 esd3 <- link2list(lsd3) 137 lsd3 <- attr(esd3, "function.name") 138 139 lrho12 <- as.list(substitute(lrho12)) 140 erho12 <- link2list(lrho12) 141 lrho12 <- attr(erho12, "function.name") 142 143 lrho23 <- as.list(substitute(lrho23)) 144 erho23 <- link2list(lrho23) 145 lrho23 <- attr(erho23, "function.name") 146 147 lrho13 <- as.list(substitute(lrho13)) 148 erho13 <- link2list(lrho13) 149 lrho13 <- attr(erho13, "function.name") 150 151 152 if (!is.logical(eq.mean) || length(eq.mean) != 1) 153 stop("argument 'eq.mean' must be a single logical") 154 if (!is.logical(eq.sd) || length(eq.sd) != 1) 155 stop("argument 'eq.sd' must be a single logical") 156 if (!is.logical(eq.cor) || length(eq.cor) != 1) 157 stop("argument 'eq.cor' must be a single logical") 158 159 160 if (!is.Numeric(imethod, length.arg = 1, 161 integer.valued = TRUE, positive = TRUE) || 162 imethod > 2) 163 stop("argument 'imethod' must be 1 or 2") 164 165 new("vglmff", 166 blurb = c("Trivariate normal distribution\n", 167 "Links: ", 168 namesof("mean1", lmean1, earg = emean1 ), ", ", 169 namesof("mean2", lmean2, earg = emean2 ), ", ", 170 namesof("mean3", lmean3, earg = emean3 ), ", ", 171 namesof("sd1", lsd1, earg = esd1 ), ", ", 172 namesof("sd2", lsd2, earg = esd2 ), ", ", 173 namesof("sd3", lsd3, earg = esd3 ), ",\n", 174 " ", 175 namesof("rho12", lrho12, earg = erho12 ), ", ", 176 namesof("rho23", lrho23, earg = erho23 ), ", ", 177 namesof("rho13", lrho13, earg = erho13 )), 178 constraints = eval(substitute(expression({ 179 180 constraints.orig <- constraints 181 M1 <- 9 182 NOS <- M / M1 183 184 cm1.m <- 185 cmk.m <- kronecker(diag(NOS), rbind(diag(3), matrix(0, 6, 3))) 186 con.m <- cm.VGAM(kronecker(diag(NOS), eijfun(1:3, 9)), 187 x = x, 188 bool = .eq.mean , # 189 constraints = constraints.orig, 190 apply.int = TRUE, 191 cm.default = cmk.m, 192 cm.intercept.default = cm1.m) 193 194 195 cm1.s <- 196 cmk.s <- kronecker(diag(NOS), 197 rbind(matrix(0, 3, 3), diag(3), matrix(0, 3, 3))) 198 con.s <- cm.VGAM(kronecker(diag(NOS), eijfun(4:6, 9)), 199 x = x, 200 bool = .eq.sd , # 201 constraints = constraints.orig, 202 apply.int = TRUE, 203 cm.default = cmk.s, 204 cm.intercept.default = cm1.s) 205 206 207 208 cm1.r <- 209 cmk.r <- kronecker(diag(NOS), 210 rbind(matrix(0, 3, 3), matrix(0, 3, 3), diag(3))) 211 con.r <- cm.VGAM(kronecker(diag(NOS), eijfun(7:9, 9)), 212 x = x, 213 bool = .eq.cor , # 214 constraints = constraints.orig, 215 apply.int = TRUE, 216 cm.default = cmk.r, 217 cm.intercept.default = cm1.r) 218 219 220 221 con.use <- con.m 222 for (klocal in seq_along(con.m)) { 223 con.use[[klocal]] <- 224 cbind(con.m[[klocal]], 225 con.s[[klocal]], 226 con.r[[klocal]]) 227 } 228 229 230 231 232 constraints <- con.use 233 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 234 predictors.names = predictors.names, 235 M1 = M1) 236 }), list( .zero = zero, 237 .eq.sd = eq.sd, 238 .eq.mean = eq.mean, 239 .eq.cor = eq.cor ))), 240 241 infos = eval(substitute(function(...) { 242 list(M1 = 9, 243 Q1 = 3, 244 expected = TRUE, 245 multipleResponses = FALSE, 246 parameters.names = c("mean1", "mean2", "mean3", 247 "sd1", "sd2", "sd3", 248 "rho12", "rho13", "rho23"), 249 eq.cor = .eq.cor , 250 eq.mean = .eq.mean , 251 eq.sd = .eq.sd , 252 zero = .zero ) 253 }, list( .zero = zero, 254 .eq.cor = eq.cor, 255 .eq.mean = eq.mean, 256 .eq.sd = eq.sd ))), 257 258 initialize = eval(substitute(expression({ 259 Q1 <- 3 260 temp5 <- 261 w.y.check(w = w, y = y, 262 ncol.y.max = Q1, 263 ncol.w.max = 1, 264 ncol.y.min = Q1, 265 out.wy = TRUE, 266 colsyperw = Q1, 267 maximize = TRUE) 268 w <- temp5$w 269 y <- temp5$y 270 271 272 273 predictors.names <- c( 274 namesof("mean1", .lmean1 , earg = .emean1 , short = TRUE), 275 namesof("mean2", .lmean2 , earg = .emean2 , short = TRUE), 276 namesof("mean3", .lmean3 , earg = .emean3 , short = TRUE), 277 namesof("sd1", .lsd1 , earg = .esd1 , short = TRUE), 278 namesof("sd2", .lsd2 , earg = .esd2 , short = TRUE), 279 namesof("sd3", .lsd3 , earg = .esd3 , short = TRUE), 280 namesof("rho12", .lrho12 , earg = .erho12 , short = TRUE), 281 namesof("rho23", .lrho23 , earg = .erho23 , short = TRUE), 282 namesof("rho13", .lrho13 , earg = .erho13 , short = TRUE)) 283 284 extra$colnames.y <- colnames(y) 285 286 if (!length(etastart)) { 287 imean1 <- rep_len(if (length( .imean1 )) .imean1 else 288 weighted.mean(y[, 1], w = w), n) 289 imean2 <- rep_len(if (length( .imean2 )) .imean2 else 290 weighted.mean(y[, 2], w = w), n) 291 imean3 <- rep_len(if (length( .imean3 )) .imean3 else 292 weighted.mean(y[, 3], w = w), n) 293 isd1 <- rep_len(if (length( .isd1 )) .isd1 else sd(y[, 1]), n) 294 isd2 <- rep_len(if (length( .isd2 )) .isd2 else sd(y[, 2]), n) 295 isd3 <- rep_len(if (length( .isd3 )) .isd3 else sd(y[, 3]), n) 296 irho12 <- rep_len(if (length( .irho12 )) .irho12 else 297 cor(y[, 1], y[, 2]), n) 298 irho23 <- rep_len(if (length( .irho23 )) .irho23 else 299 cor(y[, 2], y[, 3]), n) 300 irho13 <- rep_len(if (length( .irho13 )) .irho13 else 301 cor(y[, 1], y[, 3]), n) 302 303 if ( .imethod == 2) { 304 imean1 <- abs(imean1) + 0.01 305 imean2 <- abs(imean2) + 0.01 306 imean3 <- abs(imean3) + 0.01 307 } 308 etastart <- 309 cbind(theta2eta(imean1, .lmean1 , earg = .emean1 ), 310 theta2eta(imean2, .lmean2 , earg = .emean2 ), 311 theta2eta(imean3, .lmean3 , earg = .emean3 ), 312 theta2eta(isd1, .lsd1 , earg = .esd1 ), 313 theta2eta(isd2, .lsd2 , earg = .esd2 ), 314 theta2eta(isd3, .lsd3 , earg = .esd3 ), 315 theta2eta(irho12, .lrho12 , earg = .erho12 ), 316 theta2eta(irho23, .lrho23 , earg = .erho23 ), 317 theta2eta(irho13, .lrho13 , earg = .erho13 )) 318 } 319 }), list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3, 320 .emean1 = emean1, .emean2 = emean2, .emean3 = emean3, 321 .imean1 = imean1, .imean2 = imean2, .imean3 = imean3, 322 .lsd1 = lsd1 , .lsd2 = lsd2 , .lsd3 = lsd3 , 323 .esd1 = esd1 , .esd2 = esd2 , .esd3 = esd3 , 324 .isd1 = isd1, .isd2 = isd2, .isd3 = isd3, 325 .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23, 326 .erho12 = erho12, .erho13 = erho13, .erho23 = erho23, 327 .irho12 = irho12, .irho13 = irho13, .irho23 = irho23, 328 .imethod = imethod ))), 329 linkinv = eval(substitute(function(eta, extra = NULL) { 330 NOS <- ncol(eta) / c(M1 = 9) 331 mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 ) 332 mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 ) 333 mean3 <- eta2theta(eta[, 3], .lmean3 , earg = .emean3 ) 334 fv.mat <- cbind(mean1, mean2, mean3) 335 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 336 } , list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3, 337 .emean1 = emean1, .emean2 = emean2, .emean3 = emean3, 338 .lsd1 = lsd1 , .lsd2 = lsd2 , .lsd3 = lsd3 , 339 .esd1 = esd1 , .esd2 = esd2 , .esd3 = esd3 , 340 .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 ))), 341 342 last = eval(substitute(expression({ 343 misc$link <- c("mean1" = .lmean1 , 344 "mean2" = .lmean2 , 345 "mean3" = .lmean3 , 346 "sd1" = .lsd1 , 347 "sd2" = .lsd2 , 348 "sd3" = .lsd3 , 349 "rho12" = .lrho12 , 350 "rho23" = .lrho23 , 351 "rho13" = .lrho13 ) 352 353 misc$earg <- list("mean1" = .emean1 , 354 "mean2" = .emean2 , 355 "mean3" = .emean3 , 356 "sd1" = .esd1 , 357 "sd2" = .esd2 , 358 "sd3" = .esd3 , 359 "rho12" = .erho12 , 360 "rho23" = .erho23 , 361 "rho13" = .erho13 ) 362 }), list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3, 363 .emean1 = emean1, .emean2 = emean2, .emean3 = emean3, 364 .lsd1 = lsd1 , .lsd2 = lsd2 , .lsd3 = lsd3 , 365 .esd1 = esd1 , .esd2 = esd2 , .esd3 = esd3 , 366 .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23, 367 .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 ))), 368 loglikelihood = eval(substitute( 369 function(mu, y, w, residuals = FALSE, eta, 370 extra = NULL, 371 summation = TRUE) { 372 mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 ) 373 mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 ) 374 mean3 <- eta2theta(eta[, 3], .lmean3 , earg = .emean3 ) 375 sd1 <- eta2theta(eta[, 4], .lsd1 , earg = .esd1 ) 376 sd2 <- eta2theta(eta[, 5], .lsd2 , earg = .esd2 ) 377 sd3 <- eta2theta(eta[, 6], .lsd3 , earg = .esd3 ) 378 Rho12 <- eta2theta(eta[, 7], .lrho12 , earg = .erho12 ) 379 Rho23 <- eta2theta(eta[, 8], .lrho23 , earg = .erho23 ) 380 Rho13 <- eta2theta(eta[, 9], .lrho13 , earg = .erho13 ) 381 382 if (residuals) { 383 stop("loglikelihood residuals not implemented yet") 384 } else { 385 ll.elts <- 386 c(w) * dtrinorm(x1 = y[, 1], x2 = y[, 2], x3 = y[, 3], 387 mean1 = mean1, mean2 = mean2, mean3 = mean3, 388 var1 = sd1^2, var2 = sd2^2, var3 = sd3^2, 389 cov12 = Rho12 * sd1 * sd2, 390 cov23 = Rho23 * sd2 * sd3, 391 cov13 = Rho13 * sd1 * sd3, 392 log = TRUE) 393 if (summation) { 394 sum(ll.elts) 395 } else { 396 ll.elts 397 } 398 } 399 } , list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3, 400 .emean1 = emean1, .emean2 = emean2, .emean3 = emean3, 401 .lsd1 = lsd1 , .lsd2 = lsd2 , .lsd3 = lsd3 , 402 .esd1 = esd1 , .esd2 = esd2 , .esd3 = esd3 , 403 .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23, 404 .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 ))), 405 vfamily = c("trinormal"), 406 validparams = eval(substitute(function(eta, y, extra = NULL) { 407 mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 ) 408 mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 ) 409 mean3 <- eta2theta(eta[, 3], .lmean3 , earg = .emean3 ) 410 sd1 <- eta2theta(eta[, 4], .lsd1 , earg = .esd1 ) 411 sd2 <- eta2theta(eta[, 5], .lsd2 , earg = .esd2 ) 412 sd3 <- eta2theta(eta[, 6], .lsd3 , earg = .esd3 ) 413 Rho12 <- eta2theta(eta[, 7], .lrho12 , earg = .erho12 ) 414 Rho23 <- eta2theta(eta[, 8], .lrho23 , earg = .erho23 ) 415 Rho13 <- eta2theta(eta[, 9], .lrho13 , earg = .erho13 ) 416 okay1 <- all(is.finite(mean1)) && 417 all(is.finite(mean2)) && 418 all(is.finite(mean3)) && 419 all(is.finite(sd1 )) && all(0 < sd1) && 420 all(is.finite(sd2 )) && all(0 < sd2) && 421 all(is.finite(sd3 )) && all(0 < sd3) && 422 all(is.finite(Rho12)) && all(abs(Rho12) < 1) && 423 all(is.finite(Rho23)) && all(abs(Rho23) < 1) && 424 all(is.finite(Rho13)) && all(abs(Rho13) < 1) && 425 all(0 < 1 - Rho12^2 - Rho13^2 - Rho23^2 + 426 2 * Rho12 * Rho13 * Rho23) 427 okay1 428 } , list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3, 429 .emean1 = emean1, .emean2 = emean2, .emean3 = emean3, 430 .lsd1 = lsd1 , .lsd2 = lsd2 , .lsd3 = lsd3 , 431 .esd1 = esd1 , .esd2 = esd2 , .esd3 = esd3 , 432 .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23, 433 .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 ))), 434 435 436 437 simslot = eval(substitute( 438 function(object, nsim) { 439 440 pwts <- if (length(pwts <- object@prior.weights) > 0) 441 pwts else weights(object, type = "prior") 442 if (any(pwts != 1)) 443 warning("ignoring prior weights") 444 eta <- predict(object) 445 mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 ) 446 mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 ) 447 mean3 <- eta2theta(eta[, 3], .lmean3 , earg = .emean3 ) 448 sd1 <- eta2theta(eta[, 4], .lsd1 , earg = .esd1 ) 449 sd2 <- eta2theta(eta[, 5], .lsd2 , earg = .esd2 ) 450 sd3 <- eta2theta(eta[, 6], .lsd3 , earg = .esd3 ) 451 Rho12 <- eta2theta(eta[, 7], .lrho12 , earg = .erho12 ) 452 Rho23 <- eta2theta(eta[, 8], .lrho23 , earg = .erho23 ) 453 Rho13 <- eta2theta(eta[, 9], .lrho13 , earg = .erho13 ) 454 rtrinorm(nsim * length(sd1), 455 mean1 = mean1, mean2 = mean2, mean3 = mean3, 456 var1 = sd1^2, var2 = sd2^2, var3 = sd3^2, 457 cov12 = Rho12 * sd1 * sd2, 458 cov23 = Rho23 * sd2 * sd3, 459 cov13 = Rho13 * sd1 * sd3) 460 } , list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3, 461 .emean1 = emean1, .emean2 = emean2, .emean3 = emean3, 462 .lsd1 = lsd1 , .lsd2 = lsd2 , .lsd3 = lsd3 , 463 .esd1 = esd1 , .esd2 = esd2 , .esd3 = esd3 , 464 .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23, 465 .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 ))), 466 467 468 469 470 deriv = eval(substitute(expression({ 471 mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 ) 472 mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 ) 473 mean3 <- eta2theta(eta[, 3], .lmean3 , earg = .emean3 ) 474 sd1 <- eta2theta(eta[, 4], .lsd1 , earg = .esd1 ) 475 sd2 <- eta2theta(eta[, 5], .lsd2 , earg = .esd2 ) 476 sd3 <- eta2theta(eta[, 6], .lsd3 , earg = .esd3 ) 477 rho12 <- eta2theta(eta[, 7], .lrho12 , earg = .erho12 ) 478 rho23 <- eta2theta(eta[, 8], .lrho23 , earg = .erho23 ) 479 rho13 <- eta2theta(eta[, 9], .lrho13 , earg = .erho13 ) 480 bbb <- 1 - rho12^2 - rho13^2 - rho23^2 + 2 * rho12 * rho13 * rho23 481 482 483 484 Sigmainv <- matrix(0, n, dimm(3)) # sum(3:1) 485 Sigmainv[, iam(1, 1, M = 3)] <- (1 - rho23^2) / (bbb * sd1^2) 486 Sigmainv[, iam(2, 2, M = 3)] <- (1 - rho13^2) / (bbb * sd2^2) 487 Sigmainv[, iam(3, 3, M = 3)] <- (1 - rho12^2) / (bbb * sd3^2) 488 Sigmainv[, iam(1, 2, M = 3)] <- (rho13 * rho23 - rho12) / ( 489 sd1 * sd2 * bbb) 490 Sigmainv[, iam(2, 3, M = 3)] <- (rho12 * rho13 - rho23) / ( 491 sd2 * sd3 * bbb) 492 Sigmainv[, iam(1, 3, M = 3)] <- (rho12 * rho23 - rho13) / ( 493 sd1 * sd3 * bbb) 494 495 496 dem <- bbb * (sd1 * sd2 * sd3)^2 497 ymat.cen <- y - cbind(mean1, mean2, mean3) # Usual dimensions n x 3 498 ymatt.cen <- t(ymat.cen) 499 dim(ymatt.cen) <- c(nrow(ymatt.cen), 1, ncol(ymatt.cen)) # 4 mux5() 500 dl.dmeans <- mux22(t(Sigmainv), ymat.cen, M = 3, as.matrix = TRUE) 501 502 503 504 SI.sd1 <- Sigmainv * 0 505 SI.sd1[, iam(1, 1, M = 3)] <- -2 * Sigmainv[, iam(1, 1, M = 3)] / sd1 506 SI.sd1[, iam(2, 2, M = 3)] <- 0 507 SI.sd1[, iam(3, 3, M = 3)] <- 0 508 SI.sd1[, iam(1, 2, M = 3)] <- -1 * Sigmainv[, iam(1, 2, M = 3)] / sd1 509 SI.sd1[, iam(2, 3, M = 3)] <- 0 510 SI.sd1[, iam(1, 3, M = 3)] <- -1 * Sigmainv[, iam(1, 3, M = 3)] / sd1 511 512 SI.sd2 <- Sigmainv * 0 513 SI.sd2[, iam(2, 2, M = 3)] <- -2 * Sigmainv[, iam(2, 2, M = 3)] / sd2 514 SI.sd2[, iam(1, 1, M = 3)] <- 0 515 SI.sd2[, iam(3, 3, M = 3)] <- 0 516 SI.sd2[, iam(1, 2, M = 3)] <- -1 * Sigmainv[, iam(1, 2, M = 3)] / sd2 517 SI.sd2[, iam(1, 3, M = 3)] <- 0 518 SI.sd2[, iam(2, 3, M = 3)] <- -1 * Sigmainv[, iam(2, 3, M = 3)] / sd2 519 520 SI.sd3 <- Sigmainv * 0 521 SI.sd3[, iam(3, 3, M = 3)] <- -2 * Sigmainv[, iam(3, 3, M = 3)] / sd3 522 SI.sd3[, iam(2, 2, M = 3)] <- 0 523 SI.sd3[, iam(1, 1, M = 3)] <- 0 524 SI.sd3[, iam(1, 3, M = 3)] <- -1 * Sigmainv[, iam(1, 3, M = 3)] / sd3 525 SI.sd3[, iam(1, 2, M = 3)] <- 0 526 SI.sd3[, iam(2, 3, M = 3)] <- -1 * Sigmainv[, iam(2, 3, M = 3)] / sd3 527 528 529 dl.dsd1 <- -1 / sd1 - 0.5 * 530 c(mux5(x = ymatt.cen, cc = SI.sd1, M = 3, matrix.arg = TRUE)) 531 dl.dsd2 <- -1 / sd2 - 0.5 * 532 c(mux5(x = ymatt.cen, cc = SI.sd2, M = 3, matrix.arg = TRUE)) 533 dl.dsd3 <- -1 / sd3 - 0.5 * 534 c(mux5(x = ymatt.cen, cc = SI.sd3, M = 3, matrix.arg = TRUE)) 535 536 537 dbbb.drho12 <- 2 * (rho13 * rho23 - rho12) 538 dbbb.drho23 <- 2 * (rho12 * rho13 - rho23) 539 dbbb.drho13 <- 2 * (rho12 * rho23 - rho13) 540 SI.rho12 <- Sigmainv * 0 541 SI.rho12[, iam(1, 1, M = 3)] <- 542 -1 * Sigmainv[, iam(1, 1, M = 3)] * dbbb.drho12 / bbb 543 SI.rho12[, iam(2, 2, M = 3)] <- 544 -1 * Sigmainv[, iam(2, 2, M = 3)] * dbbb.drho12 / bbb 545 SI.rho12[, iam(3, 3, M = 3)] <- 546 (-2 * rho12 - (1 - rho12^2) * dbbb.drho12 / bbb) / (bbb * sd3^2) 547 SI.rho12[, iam(1, 2, M = 3)] <- 548 (-1 - (rho13 * rho23 - rho12) * dbbb.drho12 / bbb) / ( 549 bbb * sd1 * sd2) 550 SI.rho12[, iam(2, 3, M = 3)] <- 551 (rho13 - (rho12 * rho13 - rho23) * dbbb.drho12 / bbb) / ( 552 bbb * sd2 * sd3) 553 SI.rho12[, iam(1, 3, M = 3)] <- 554 (rho23 - (rho12 * rho23 - rho13) * dbbb.drho12 / bbb) / ( 555 bbb * sd1 * sd3) 556 557 558 SI.rho23 <- Sigmainv * 0 559 SI.rho23[, iam(1, 1, M = 3)] <- 560 (-2 * rho23 - (1 - rho23^2) * dbbb.drho23 / bbb) / (bbb * sd1^2) 561 SI.rho23[, iam(2, 2, M = 3)] <- 562 -1 * Sigmainv[, iam(2, 2, M = 3)] * dbbb.drho23 / bbb 563 SI.rho23[, iam(3, 3, M = 3)] <- 564 -1 * Sigmainv[, iam(3, 3, M = 3)] * dbbb.drho23 / bbb 565 SI.rho23[, iam(1, 2, M = 3)] <- 566 (rho13 - (rho13 * rho23 - rho12) * dbbb.drho23 / bbb) / ( 567 bbb * sd1 * sd2) 568 SI.rho23[, iam(2, 3, M = 3)] <- 569 (-1 - (rho12 * rho13 - rho23) * dbbb.drho23 / bbb) / ( 570 bbb * sd2 * sd3) 571 SI.rho23[, iam(1, 3, M = 3)] <- 572 (rho12 - (rho12 * rho23 - rho13) * dbbb.drho23 / bbb) / ( 573 bbb * sd1 * sd3) 574 575 SI.rho13 <- Sigmainv * 0 576 SI.rho13[, iam(1, 1, M = 3)] <- 577 -1 * Sigmainv[, iam(1, 1, M = 3)] * dbbb.drho13 / bbb 578 SI.rho13[, iam(2, 2, M = 3)] <- 579 (-2 * rho13 - (1 - rho13^2) * dbbb.drho13 / bbb) / (bbb * sd2^2) 580 SI.rho13[, iam(3, 3, M = 3)] <- 581 -1 * Sigmainv[, iam(3, 3, M = 3)] * dbbb.drho13 / bbb 582 SI.rho13[, iam(1, 2, M = 3)] <- 583 (rho23 - (rho13 * rho23 - rho12) * dbbb.drho13 / bbb) / ( 584 bbb * sd1 * sd2) 585 SI.rho13[, iam(2, 3, M = 3)] <- 586 (rho12 - (rho12 * rho13 - rho23) * dbbb.drho13 / bbb) / ( 587 bbb * sd2 * sd3) 588 SI.rho13[, iam(1, 3, M = 3)] <- 589 (-1 - (rho12 * rho23 - rho13) * dbbb.drho13 / bbb) / ( 590 bbb * sd1 * sd3) 591 592 dl.drho12 <- -0.5 * dbbb.drho12 / bbb - 0.5 * 593 c(mux5(x = ymatt.cen, cc = SI.rho12, M = 3, matrix.arg = TRUE)) 594 dl.drho23 <- -0.5 * dbbb.drho23 / bbb - 0.5 * 595 c(mux5(x = ymatt.cen, cc = SI.rho23, M = 3, matrix.arg = TRUE)) 596 dl.drho13 <- -0.5 * dbbb.drho13 / bbb - 0.5 * 597 c(mux5(x = ymatt.cen, cc = SI.rho13, M = 3, matrix.arg = TRUE)) 598 599 600 dmean1.deta <- dtheta.deta(mean1, .lmean1 ) 601 dmean2.deta <- dtheta.deta(mean2, .lmean2 ) 602 dmean3.deta <- dtheta.deta(mean3, .lmean3 ) 603 dsd1.deta <- dtheta.deta(sd1 , .lsd1 ) 604 dsd2.deta <- dtheta.deta(sd2 , .lsd2 ) 605 dsd3.deta <- dtheta.deta(sd3 , .lsd3 ) 606 drho12.deta <- dtheta.deta(rho12, .lrho12 ) 607 drho23.deta <- dtheta.deta(rho23, .lrho23 ) 608 drho13.deta <- dtheta.deta(rho13, .lrho13 ) 609 dThetas.detas <- cbind(dmean1.deta, 610 dmean2.deta, 611 dmean3.deta, 612 dsd1.deta, 613 dsd2.deta, 614 dsd3.deta, 615 drho12.deta, 616 drho23.deta, 617 drho13.deta) 618 c(w) * cbind(dl.dmeans, # dl.dmeans[, 1:3], 619 dl.dsd1, 620 dl.dsd2, 621 dl.dsd3, 622 dl.drho12, 623 dl.drho23, 624 dl.drho13) * dThetas.detas 625 }), list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3, 626 .emean1 = emean1, .emean2 = emean2, .emean3 = emean3, 627 .lsd1 = lsd1 , .lsd2 = lsd2 , .lsd3 = lsd3 , 628 .esd1 = esd1 , .esd2 = esd2 , .esd3 = esd3 , 629 .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23, 630 .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 ))), 631 632 weight = eval(substitute(expression({ 633 wz <- matrix(0, n, dimm(M)) 634 wz[, iam(1, 1, M)] <- Sigmainv[, iam(1, 1, M = 3)] 635 wz[, iam(2, 2, M)] <- Sigmainv[, iam(2, 2, M = 3)] 636 wz[, iam(3, 3, M)] <- Sigmainv[, iam(3, 3, M = 3)] 637 wz[, iam(1, 2, M)] <- Sigmainv[, iam(1, 2, M = 3)] 638 wz[, iam(2, 3, M)] <- Sigmainv[, iam(2, 3, M = 3)] 639 wz[, iam(1, 3, M)] <- Sigmainv[, iam(1, 3, M = 3)] 640 641 642if (FALSE) { 643 wz[, iam(4, 4, M)] <- -1 / sd1^2 + 644 (1 - rho23^2 + 2 * bbb) / (bbb * sd1^2) 645 wz[, iam(5, 5, M)] <- -1 / sd2^2 + 646 (1 - rho13^2 + 2 * bbb) / (bbb * sd2^2) 647 wz[, iam(6, 6, M)] <- -1 / sd3^2 + 648 (1 - rho12^2 + 2 * bbb) / (bbb * sd3^2) 649 wz[, iam(4, 5, M)] <- 0 - 650 rho12 * (rho13 * rho23 - rho12) / (sd1 * sd2 * bbb) 651 wz[, iam(5, 6, M)] <- 0 - 652 rho23 * (rho12 * rho13 - rho23) / (sd2 * sd3 * bbb) 653 wz[, iam(4, 6, M)] <- 0 - 654 rho13 * (rho12 * rho23 - rho13) / (sd1 * sd3 * bbb) 655} 656 657if (FALSE) { 658 d2bbb.drho12.12 <- -2 659 d2bbb.drho23.23 <- -2 660 d2bbb.drho13.13 <- -2 661 d2bbb.drho12.13 <- 2 * rho23 662 d2bbb.drho12.23 <- 2 * rho13 663 d2bbb.drho13.23 <- 2 * rho12 664 wz[, iam(7, 7, M)] <- 665 0.5 * (d2bbb.drho12.12 - dbbb.drho12 * dbbb.drho12 / bbb) / bbb 666 wz[, iam(8, 8, M)] <- 667 0.5 * (d2bbb.drho23.23 - dbbb.drho23 * dbbb.drho23 / bbb) / bbb 668 wz[, iam(9, 9, M)] <- 669 0.5 * (d2bbb.drho13.13 - dbbb.drho13 * dbbb.drho13 / bbb) / bbb 670 wz[, iam(7, 8, M)] <- 671 0.5 * (d2bbb.drho12.23 - dbbb.drho12 * dbbb.drho23 / bbb) / bbb 672 wz[, iam(7, 9, M)] <- 673 0.5 * (d2bbb.drho12.13 - dbbb.drho12 * dbbb.drho13 / bbb) / bbb 674 wz[, iam(8, 9, M)] <- 675 0.5 * (d2bbb.drho13.23 - dbbb.drho13 * dbbb.drho23 / bbb) / bbb 676} 677 678 679 680 mux43mat <- function(A, B, C, D, aa, bb) { 681 682 s <- rep(0, length(A[, 1])) 683 for (i1 in 1:3) 684 for (i2 in 1:3) 685 for (i3 in 1:3) 686 s <- s + A[, iam(aa, i1, M = 3)] * 687 B[, iam(i1, i2, M = 3)] * 688 C[, iam(i2, i3, M = 3)] * 689 D[, iam(i3, bb, M = 3)] 690 s 691 } # mux43mat 692 693 694 695 Sigma <- matrix(0, n, dimm(3)) # sum(3:1) 696 Sigma[, iam(1, 1, M = 3)] <- sd1^2 697 Sigma[, iam(2, 2, M = 3)] <- sd2^2 698 Sigma[, iam(3, 3, M = 3)] <- sd3^2 699 Sigma[, iam(1, 2, M = 3)] <- rho12 * sd1 * sd2 700 Sigma[, iam(2, 3, M = 3)] <- rho23 * sd2 * sd3 701 Sigma[, iam(1, 3, M = 3)] <- rho13 * sd1 * sd3 702 703 704 705 for (ii in 1:3) 706 wz[, iam(4, 4, M)] <- 707 wz[, iam(4, 4, M)] + 708 0.5 * mux43mat(Sigma, SI.sd1, Sigma, SI.sd1, ii, ii) 709 710 for (ii in 1:3) 711 wz[, iam(5, 5, M)] <- 712 wz[, iam(5, 5, M)] + 713 0.5 * mux43mat(Sigma, SI.sd2, Sigma, SI.sd2, ii, ii) 714 715 for (ii in 1:3) 716 wz[, iam(6, 6, M)] <- 717 wz[, iam(6, 6, M)] + 718 0.5 * mux43mat(Sigma, SI.sd3, Sigma, SI.sd3, ii, ii) 719 720 for (ii in 1:3) 721 wz[, iam(4, 5, M)] <- 722 wz[, iam(4, 5, M)] + 723 0.5 * mux43mat(Sigma, SI.sd1, Sigma, SI.sd2, ii, ii) 724 725 for (ii in 1:3) 726 wz[, iam(5, 6, M)] <- 727 wz[, iam(5, 6, M)] + 728 0.5 * mux43mat(Sigma, SI.sd2, Sigma, SI.sd3, ii, ii) 729 730 for (ii in 1:3) 731 wz[, iam(4, 6, M)] <- 732 wz[, iam(4, 6, M)] + 733 0.5 * mux43mat(Sigma, SI.sd1, Sigma, SI.sd3, ii, ii) 734 735 736 737 738 739 740 for (ii in 1:3) 741 wz[, iam(4, 7, M)] <- 742 wz[, iam(4, 7, M)] + 743 0.5 * mux43mat(Sigma, SI.sd1, Sigma, SI.rho12, ii, ii) 744 745 for (ii in 1:3) 746 wz[, iam(4, 8, M)] <- 747 wz[, iam(4, 8, M)] + 748 0.5 * mux43mat(Sigma, SI.sd1, Sigma, SI.rho23, ii, ii) 749 750 for (ii in 1:3) 751 wz[, iam(4, 9, M)] <- 752 wz[, iam(4, 9, M)] + 753 0.5 * mux43mat(Sigma, SI.sd1, Sigma, SI.rho13, ii, ii) 754 755 for (ii in 1:3) 756 wz[, iam(5, 7, M)] <- 757 wz[, iam(5, 7, M)] + 758 0.5 * mux43mat(Sigma, SI.sd2, Sigma, SI.rho12, ii, ii) 759 760 for (ii in 1:3) 761 wz[, iam(5, 8, M)] <- 762 wz[, iam(5, 8, M)] + 763 0.5 * mux43mat(Sigma, SI.sd2, Sigma, SI.rho23, ii, ii) 764 765 for (ii in 1:3) 766 wz[, iam(5, 9, M)] <- 767 wz[, iam(5, 9, M)] + 768 0.5 * mux43mat(Sigma, SI.sd2, Sigma, SI.rho13, ii, ii) 769 770 for (ii in 1:3) 771 wz[, iam(6, 7, M)] <- 772 wz[, iam(6, 7, M)] + 773 0.5 * mux43mat(Sigma, SI.sd3, Sigma, SI.rho12, ii, ii) 774 775 for (ii in 1:3) 776 wz[, iam(6, 8, M)] <- 777 wz[, iam(6, 8, M)] + 778 0.5 * mux43mat(Sigma, SI.sd3, Sigma, SI.rho23, ii, ii) 779 780 for (ii in 1:3) 781 wz[, iam(6, 9, M)] <- 782 wz[, iam(6, 9, M)] + 783 0.5 * mux43mat(Sigma, SI.sd3, Sigma, SI.rho13, ii, ii) 784 785 786 787 788 for (ii in 1:3) 789 wz[, iam(7, 7, M)] <- 790 wz[, iam(7, 7, M)] + 791 0.5 * mux43mat(Sigma, SI.rho12, Sigma, SI.rho12, ii, ii) 792 793 for (ii in 1:3) 794 wz[, iam(8, 8, M)] <- 795 wz[, iam(8, 8, M)] + 796 0.5 * mux43mat(Sigma, SI.rho23, Sigma, SI.rho23, ii, ii) 797 798 for (ii in 1:3) 799 wz[, iam(9, 9, M)] <- 800 wz[, iam(9, 9, M)] + 801 0.5 * mux43mat(Sigma, SI.rho13, Sigma, SI.rho13, ii, ii) 802 803 for (ii in 1:3) 804 wz[, iam(7, 8, M)] <- 805 wz[, iam(7, 8, M)] + 806 0.5 * mux43mat(Sigma, SI.rho12, Sigma, SI.rho23, ii, ii) 807 808 for (ii in 1:3) 809 wz[, iam(8, 9, M)] <- 810 wz[, iam(8, 9, M)] + 811 0.5 * mux43mat(Sigma, SI.rho23, Sigma, SI.rho13, ii, ii) 812 813 for (ii in 1:3) 814 wz[, iam(7, 9, M)] <- 815 wz[, iam(7, 9, M)] + 816 0.5 * mux43mat(Sigma, SI.rho12, Sigma, SI.rho13, ii, ii) 817 818 819 820 ind5 <- iam(NA, NA, M = M, both = TRUE, diag = TRUE) 821 wz <- wz * dThetas.detas[, ind5$row.index] * 822 dThetas.detas[, ind5$col.index] 823 c(w) * wz 824 }), list( .lmean1 = lmean1, .lmean2 = lmean2, .lmean3 = lmean3, 825 .emean1 = emean1, .emean2 = emean2, .emean3 = emean3, 826 .lsd1 = lsd1 , .lsd2 = lsd2 , .lsd3 = lsd3 , 827 .esd1 = esd1 , .esd2 = esd2 , .esd3 = esd3 , 828 .lrho12 = lrho12, .lrho13 = lrho13, .lrho23 = lrho23, 829 .erho12 = erho12, .erho13 = erho13, .erho23 = erho23 )))) 830} # trinormal 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846dbiclaytoncop <- function(x1, x2, apar = 0, log = FALSE) { 847 if (!is.logical(log.arg <- log) || length(log) != 1) 848 stop("bad input for argument 'log'") 849 rm(log) 850 851 A <- x1^(-apar) + x2^(-apar) - 1 852 logdensity <- log1p(apar) - 853 (1 + apar) * (log(x1) + log(x2)) - 854 (2 + 1 / apar) * log(abs(A)) # Avoid warning 855 856 out.square <- (x1 < 0) | (x1 > 1) | (x2 < 0) | (x2 > 1) 857 logdensity[out.square] <- log(0.0) 858 859 860 index0 <- (rep_len(apar, length(A)) < sqrt(.Machine$double.eps)) 861 if (any(index0)) 862 logdensity[index0] <- log(1.0) 863 864 865 index1 <- (rep_len(apar, length(A)) < 0.0) | (A < 0.0) 866 if (any(index1)) 867 logdensity[index1] <- NaN 868 869 870 871 872 873 874 if (log.arg) logdensity else exp(logdensity) 875} 876 877 878 879rbiclaytoncop <- function(n, apar = 0) { 880 if (any(apar < 0)) 881 stop("argument 'apar' must be greater or equal to 0") 882 883 u1 <- runif(n = n) 884 v2 <- runif(n = n) 885 886 u2 <- (u1^(-apar) * 887 (v2^(-apar / (1 + apar)) - 1) + 1)^(-1 / apar) 888 889 890 index0 <- (rep_len(apar, length(u1)) < sqrt(.Machine$double.eps)) 891 if (any(index0)) 892 u2[index0] <- runif(sum(index0)) 893 894 cbind(u1, u2) 895} # dbiclaytoncop 896 897 898 899 biclaytoncop <- function(lapar = "loglink", 900 iapar = NULL, 901 imethod = 1, 902 parallel = FALSE, 903 zero = NULL) { 904 905 apply.parint <- TRUE 906 907 908 lapar <- as.list(substitute(lapar)) 909 eapar <- link2list(lapar) 910 lapar <- attr(eapar, "function.name") 911 912 913 if (length(iapar) && any(iapar <= 0)) 914 stop("argument 'iapar' must have values in (0, Inf)") 915 916 917 918 if (!is.Numeric(imethod, length.arg = 1, 919 integer.valued = TRUE, positive = TRUE) || 920 imethod > 3) 921 stop("argument 'imethod' must be 1 or 2 or 3") 922 923 924 925 new("vglmff", 926 blurb = c("Bivariate Clayton copula distribution)\n", 927 "Links: ", namesof("apar", lapar, earg = eapar)), 928 929 constraints = eval(substitute(expression({ 930 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 931 bool = .parallel , 932 constraints = constraints, 933 apply.int = .apply.parint ) 934 935 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 936 predictors.names = predictors.names, 937 M1 = 1) 938 }), list( .zero = zero, 939 .apply.parint = apply.parint, 940 .parallel = parallel ))), 941 942 infos = eval(substitute(function(...) { 943 list(M1 = 1, 944 Q1 = 2, 945 apply.parint = .apply.parint , 946 parameters.names = c("apar"), 947 lapar = .lapar , 948 parallel = .parallel , 949 zero = .zero ) 950 }, list( .zero = zero, 951 .apply.parint = apply.parint, 952 .lapar = lapar, 953 .parallel = parallel ))), 954 955 initialize = eval(substitute(expression({ 956 M1 <- 1 957 Q1 <- 2 958 959 temp5 <- 960 w.y.check(w = w, y = y, 961 Is.positive.y = TRUE, 962 ncol.w.max = Inf, 963 ncol.y.max = Inf, 964 ncol.y.min = Q1, 965 out.wy = TRUE, 966 colsyperw = Q1, 967 maximize = TRUE) 968 969 w <- temp5$w 970 y <- temp5$y 971 972 973 ncoly <- ncol(y) 974 extra$ncoly <- ncoly 975 extra$M1 <- M1 976 extra$Q1 <- Q1 977 M <- M1 * (ncoly / Q1) 978 mynames1 <- param.names("apar", M / M1, skip1 = TRUE) 979 predictors.names <- 980 namesof(mynames1, .lapar , earg = .eapar , short = TRUE) 981 982 extra$colnames.y <- colnames(y) 983 984 if (!length(etastart)) { 985 986 apar.init <- matrix(if (length( .iapar )) .iapar else NA_real_, 987 n, M / M1, byrow = TRUE) 988 989 if (!length( .iapar )) 990 for (spp. in 1:(M / M1)) { 991 ymatj <- y[, (Q1 * spp. - 1):(Q1 * spp.)] 992 993 994 apar.init0 <- if ( .imethod == 1) { 995 k.tau <- kendall.tau(ymatj[, 1], ymatj[, 2], exact = FALSE, 996 max.n = 500) 997 998 max(0.1, 2 * k.tau / (1 - k.tau)) # Must be positive 999 } else if ( .imethod == 2) { 1000 spearman.rho <- max(0.05, cor(ymatj[, 1], 1001 ymatj[, 2], meth = "spearman")) 1002 rhobitlink(spearman.rho) 1003 } else { 1004 pearson.rho <- max(0.05, cor(ymatj[, 1], ymatj[, 2])) 1005 rhobitlink(pearson.rho) 1006 } 1007 1008 if (anyNA(apar.init[, spp.])) 1009 apar.init[, spp.] <- apar.init0 1010 } 1011 1012 etastart <- theta2eta(apar.init, .lapar , earg = .eapar ) 1013 } 1014 }), list( .imethod = imethod, 1015 .lapar = lapar, 1016 .eapar = eapar, 1017 .iapar = iapar ))), 1018 linkinv = eval(substitute(function(eta, extra = NULL) { 1019 NOS <- NCOL(eta) / c(M1 = 1) 1020 Q1 <- 2 1021 fv.mat <- matrix(0.5, NROW(eta), NOS * Q1) 1022 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 1023 } , list( .lapar = lapar, 1024 .eapar = eapar ))), 1025 1026 last = eval(substitute(expression({ 1027 M1 <- extra$M1 1028 Q1 <- extra$Q1 1029 misc$link <- rep_len( .lapar , M) 1030 temp.names <- mynames1 1031 names(misc$link) <- temp.names 1032 1033 misc$earg <- vector("list", M) 1034 names(misc$earg) <- temp.names 1035 for (ii in 1:M) { 1036 misc$earg[[ii]] <- .eapar 1037 } 1038 1039 misc$M1 <- M1 1040 misc$Q1 <- Q1 1041 misc$imethod <- .imethod 1042 misc$expected <- TRUE 1043 misc$parallel <- .parallel 1044 misc$apply.parint <- .apply.parint 1045 misc$multipleResponses <- TRUE 1046 }) , list( .imethod = imethod, 1047 .parallel = parallel, .apply.parint = apply.parint, 1048 .lapar = lapar, 1049 .eapar = eapar ))), 1050 loglikelihood = eval(substitute( 1051 function(mu, y, w, residuals = FALSE, eta, 1052 extra = NULL, 1053 summation = TRUE) { 1054 Alpha <- eta2theta(eta, .lapar , earg = .eapar ) 1055 1056 if (residuals) { 1057 stop("loglikelihood residuals not implemented yet") 1058 } else { 1059 1060 ll.elts <- 1061 c(w) * dbiclaytoncop(x1 = c(y[, c(TRUE, FALSE)]), 1062 x2 = c(y[, c(FALSE, TRUE)]), 1063 apar = c(Alpha), log = TRUE) 1064 if (summation) { 1065 sum(ll.elts) 1066 } else { 1067 ll.elts 1068 } 1069 } 1070 } , list( .lapar = lapar, .eapar = eapar, .imethod = imethod ))), 1071 vfamily = c("biclaytoncop"), 1072 validparams = eval(substitute(function(eta, y, extra = NULL) { 1073 Alpha <- eta2theta(eta, .lapar , earg = .eapar ) 1074 okay1 <- all(is.finite(Alpha)) && all(0 < Alpha) 1075 okay1 1076 } , list( .lapar = lapar, .eapar = eapar, .imethod = imethod ))), 1077 1078 simslot = eval(substitute( 1079 function(object, nsim) { 1080 pwts <- if (length(pwts <- object@prior.weights) > 0) 1081 pwts else weights(object, type = "prior") 1082 if (any(pwts != 1)) 1083 warning("ignoring prior weights") 1084 eta <- predict(object) 1085 Alpha <- eta2theta(eta, .lapar , earg = .eapar ) 1086 rbiclaytoncop(nsim * length(Alpha), 1087 apar = c(Alpha)) 1088 } , list( .lapar = lapar, 1089 .eapar = eapar ))), 1090 1091 1092 deriv = eval(substitute(expression({ 1093 Alpha <- eta2theta(eta, .lapar , earg = .eapar ) 1094 Yindex1 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) - 1 1095 Yindex2 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) 1096 1097 1098 1099 1100 1101 AA <- y[, Yindex1]^(-Alpha) + y[, Yindex2]^(-Alpha) - 1 1102 dAA.dapar <- -y[, Yindex1]^(-Alpha) * log(y[, Yindex1]) - 1103 y[, Yindex2]^(-Alpha) * log(y[, Yindex2]) 1104 dl.dapar <- 1 / (1 + Alpha) - log(y[, Yindex1] * y[, Yindex2]) - 1105 dAA.dapar / AA * (2 + 1 / Alpha ) + log(AA) / Alpha^2 1106 1107 1108 1109 dapar.deta <- dtheta.deta(Alpha, .lapar , earg = .eapar ) 1110 1111 dl.deta <- c(w) * cbind(dl.dapar) * dapar.deta 1112 dl.deta 1113 }), list( .lapar = lapar, 1114 .eapar = eapar, 1115 .imethod = imethod ))), 1116 1117 weight = eval(substitute(expression({ 1118 par <- Alpha + 1 1119 denom1 <- (3 * par -2) * (2 * par - 1) 1120 denom2 <- 2 * (par - 1) 1121 v1 <- trigamma(1 / denom2) 1122 v2 <- trigamma(par / denom2) 1123 v3 <- trigamma((2 * par - 1) / denom2) 1124 Rho. <- (1 + par * (v1 - v2) / denom2 + 1125 (v2 - v3) / denom2) / denom1 1126 1127 ned2l.dapar <- 1 / par^2 + 2 / (par * (par - 1) * (2 * par - 1)) + 1128 4 * par / (3 * par - 2) - 1129 2 * (2 * par - 1) * Rho. / (par - 1) 1130 1131 wz <- ned2l.dapar * dapar.deta^2 1132 c(w) * wz 1133 }), list( .lapar = lapar, 1134 .eapar = eapar, 1135 .imethod = imethod )))) 1136} # biclaytoncop 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154dbistudentt <- function(x1, x2, df, rho = 0, log = FALSE) { 1155 1156 1157 1158 1159 if (!is.logical(log.arg <- log) || length(log) != 1) 1160 stop("bad input for argument 'log'") 1161 rm(log) 1162 1163 logdensity <- 1164 -(df/2 + 1) * log1p( 1165 (x1^2 + x2^2 - 2 * rho * x1 * x2) / (df * (1 - rho^2))) - 1166 log(2 * pi) - 0.5 * log1p(-rho^2) # - 1167 1168 logdensity[df <= 0] <- NaN # Not picked up by dt(). 1169 1170 logdensity[is.infinite(x1) | is.infinite(x2)] <- log(0) 1171 1172 if (log.arg) logdensity else exp(logdensity) 1173} 1174 1175 1176 1177 1178if (FALSE) 1179bistudent.deriv.dof <- function(u, v, nu, rho) { 1180 1181 1182 t1 <- qt(u, nu, 1, 0) 1183 t2 <- qt(v, nu, 1, 0) 1184 t3 <- -(nu + 2.0) / 2.0 1185 t10 <- nu * (1.0 - rho * rho) 1186 t4 <- -2.0 * t1 * t2 / t10 1187 t11 <- (t1 * t1 + t2 * t2 - 2.0 * rho * t1 * t2) 1188 t5 <- 2.0 * t11 * rho / t10 / (1.0 - rho * rho) 1189 t6 <- 1.0 + (t11 / t10) 1190 t7 <- rho / (1.0 - rho * rho) 1191 out <- (t3 * (t4 + t5) / t6 + t7) 1192} 1193 1194 1195 1196 1197 1198 1199 1200 bistudentt <- 1201 function(ldf = "logloglink", 1202 lrho = "rhobitlink", 1203 idf = NULL, 1204 irho = NULL, 1205 imethod = 1, 1206 parallel = FALSE, 1207 zero = "rho") { 1208 1209 1210 1211 1212 apply.parint <- TRUE 1213 1214 ldof <- as.list(substitute(ldf)) 1215 edof <- link2list(ldof) 1216 ldof <- attr(edof, "function.name") 1217 1218 lrho <- as.list(substitute(lrho)) 1219 erho <- link2list(lrho) 1220 lrho <- attr(erho, "function.name") 1221 1222 1223 idof <- idf 1224 if (length(idof) && 1225 any(idof <= 1)) 1226 stop("argument 'idf' must have values in (1,Inf)") 1227 1228 1229 if (length(irho) && 1230 any(abs(irho) >= 1)) 1231 stop("argument 'irho' must have values in (-1,1)") 1232 1233 1234 1235 if (!is.Numeric(imethod, length.arg = 1, 1236 integer.valued = TRUE, positive = TRUE) || 1237 imethod > 2) 1238 stop("argument 'imethod' must be 1 or 2") 1239 1240 new("vglmff", 1241 blurb = c("Bivariate student-t distribution\n", 1242 "Links: ", 1243 namesof("df", ldof, earg = edof), ", ", 1244 namesof("rho", lrho, earg = erho)), 1245 1246 constraints = eval(substitute(expression({ 1247 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 1248 bool = .parallel , 1249 constraints = constraints, 1250 apply.int = .apply.parint ) 1251 1252 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 1253 predictors.names = predictors.names, 1254 M1 = 2) 1255 }), list( .zero = zero, 1256 .apply.parint = apply.parint, 1257 .parallel = parallel ))), 1258 1259 infos = eval(substitute(function(...) { 1260 list(M1 = 2, 1261 Q1 = 2, 1262 parameters.names = c("df", "rho"), 1263 apply.parint = .apply.parint , 1264 parallel = .parallel , 1265 zero = .zero ) 1266 }, list( .zero = zero, 1267 .apply.parint = apply.parint, 1268 .parallel = parallel ))), 1269 1270 initialize = eval(substitute(expression({ 1271 M1 <- 2 1272 Q1 <- 2 1273 1274 temp5 <- 1275 w.y.check(w = w, y = y, 1276 ncol.w.max = Inf, 1277 ncol.y.max = Inf, 1278 ncol.y.min = Q1, 1279 out.wy = TRUE, 1280 colsyperw = Q1, 1281 maximize = TRUE) 1282 1283 w <- temp5$w 1284 y <- temp5$y 1285 1286 1287 ncoly <- ncol(y) 1288 extra$ncoly <- ncoly 1289 extra$M1 <- M1 1290 extra$Q1 <- Q1 1291 M <- M1 * (ncoly / Q1) 1292 mynames1 <- param.names("df", M / M1, skip1 = TRUE) 1293 mynames2 <- param.names("rho", M / M1, skip1 = TRUE) 1294 predictors.names <- c( 1295 namesof(mynames1, .ldof , earg = .edof , short = TRUE), 1296 namesof(mynames2, .lrho , earg = .erho , short = TRUE))[ 1297 interleave.VGAM(M, M1 = M1)] 1298 1299 1300 extra$colnames.y <- colnames(y) 1301 1302 if (!length(etastart)) { 1303 1304 dof.init <- matrix(if (length( .idof )) .idof else 0 + NA, 1305 n, M / M1, byrow = TRUE) 1306 rho.init <- matrix(if (length( .irho )) .irho else 0 + NA, 1307 n, M / M1, byrow = TRUE) 1308 1309 if (!length( .idof ) || !length( .irho )) 1310 for (spp. in 1:(M / M1)) { 1311 ymatj <- y[, (M1 * spp. - 1):(M1 * spp.)] 1312 1313 1314 dof.init0 <- if ( .imethod == 1) { 1315 1316 1317 2 + rexp(n = 1, rate = 0.1) 1318 } else { 1319 10 1320 } 1321 1322 if (anyNA(dof.init[, spp.])) 1323 dof.init[, spp.] <- dof.init0 1324 1325 1326 rho.init0 <- if ( .imethod == 2) { 1327 runif(n, min = -1 + 0.1, max = 1 - 0.1) 1328 } else { 1329 cor(ymatj[, 1], ymatj[, 2]) 1330 } 1331 1332 if (anyNA(rho.init[, spp.])) 1333 rho.init[, spp.] <- rho.init0 1334 1335 } 1336 1337 etastart <- 1338 cbind(theta2eta(dof.init, .ldof , earg = .edof ), 1339 theta2eta(rho.init, .lrho , earg = .erho )) 1340 1341 etastart <- etastart[, interleave.VGAM(M, M1 = M1)] 1342 1343 } 1344 }), list( .imethod = imethod, 1345 .lrho = lrho, .ldof = ldof, 1346 .erho = erho, .edof = edof, 1347 .idof = idof, .irho = irho ))), 1348 linkinv = eval(substitute(function(eta, extra = NULL) { 1349 NOS <- ncol(eta) / c(M1 = 2) 1350 Q1 <- 2 1351 fv.mat <- matrix(0, nrow(eta), Q1 * NOS) 1352 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 1353 } , list( .lrho = lrho, .ldof = ldof, 1354 .erho = erho, .edof = edof ))), 1355 1356 last = eval(substitute(expression({ 1357 1358 M1 <- extra$M1 1359 Q1 <- extra$Q1 1360 misc$link <- 1361 c(rep_len( .ldof , M / M1), 1362 rep_len( .lrho , M / M1))[interleave.VGAM(M, M1 = M1)] 1363 temp.names <- c(mynames1, mynames2)[interleave.VGAM(M, M1 = M1)] 1364 names(misc$link) <- temp.names 1365 1366 misc$earg <- vector("list", M) 1367 names(misc$earg) <- temp.names 1368 for (ii in 1:(M / M1)) { 1369 misc$earg[[M1*ii-1]] <- .edof 1370 misc$earg[[M1*ii ]] <- .erho 1371 } 1372 1373 misc$M1 <- M1 1374 misc$Q1 <- Q1 1375 misc$imethod <- .imethod 1376 misc$expected <- TRUE 1377 misc$parallel <- .parallel 1378 misc$apply.parint <- .apply.parint 1379 misc$multipleResponses <- TRUE 1380 1381 }) , list( .imethod = imethod, 1382 .parallel = parallel, 1383 .apply.parint = apply.parint, 1384 .lrho = lrho, .ldof = ldof, 1385 .erho = erho, .edof = edof ))), 1386 loglikelihood = eval(substitute( 1387 function(mu, y, w, residuals = FALSE, eta, 1388 extra = NULL, 1389 summation = TRUE) { 1390 Dof <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE], 1391 .ldof , earg = .edof ) 1392 Rho <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE], 1393 .lrho , earg = .erho ) 1394 1395 if (residuals) { 1396 stop("loglikelihood residuals not implemented yet") 1397 } else { 1398 Yindex1 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) - 1 1399 Yindex2 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) 1400 ll.elts <- 1401 c(w) * dbistudentt(x1 = y[, Yindex1, drop = FALSE], 1402 x2 = y[, Yindex2, drop = FALSE], 1403 df = Dof, 1404 rho = Rho, log = TRUE) 1405 if (summation) { 1406 sum(ll.elts) 1407 } else { 1408 ll.elts 1409 } 1410 } 1411 }, list( .lrho = lrho, .ldof = ldof, 1412 .erho = erho, .edof = edof, 1413 .imethod = imethod ))), 1414 vfamily = c("bistudentt"), 1415 validparams = eval(substitute(function(eta, y, extra = NULL) { 1416 Dof <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE], 1417 .ldof , earg = .edof ) 1418 Rho <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE], 1419 .lrho , earg = .erho ) 1420 okay1 <- all(is.finite(Dof)) && all(0 < Dof) && 1421 all(is.finite(Rho)) && all(abs(Rho) < 1) 1422 okay1 1423 }, list( .lrho = lrho, .ldof = ldof, 1424 .erho = erho, .edof = edof, 1425 .imethod = imethod ))), 1426 deriv = eval(substitute(expression({ 1427 M1 <- Q1 <- 2 1428 Dof <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE], 1429 .ldof , earg = .edof ) 1430 Rho <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE], 1431 .lrho , earg = .erho ) 1432 Yindex1 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) - 1 1433 Yindex2 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) 1434 1435 1436 x1 <- c(y[, Yindex1]) # Convert into a vector 1437 x2 <- c(y[, Yindex2]) 1438 1439 dee3 <- deriv3( ~ 1440 -(Dof/2 + 1) * log(1 + 1441 (x1^2 + x2^2 - 2 * Rho * x1 * x2) / (Dof * (1 - Rho^2))) - 1442 log(2 * pi) - 0.5 * log(1 - Rho^2), 1443 namevec = c("Dof", "Rho"), hessian = FALSE) 1444 eval.d3 <- eval(dee3) 1445 1446 dl.dthetas <- attr(eval.d3, "gradient") 1447 1448 dl.ddof <- matrix(dl.dthetas[, "Dof"], n, length(Yindex1)) 1449 dl.drho <- matrix(dl.dthetas[, "Rho"], n, length(Yindex2)) 1450 1451 1452 if (FALSE) { 1453 dd <- cbind(y, Rho, Dof) 1454 pp <- apply(dd, 1, function(x) 1455 BiCopPDF(x[1], x[2], family = 2, x[3], x[4])) 1456 alt.dl.ddof <- apply(dd, 1, function(x) 1457 BiCopDeriv(x[1], x[2], family = 2, 1458 x[3], x[4], "par2")) / pp 1459 alt.dl.drho <- apply(dd, 1, function(x) 1460 BiCopDeriv(x[1], x[2], family = 2, 1461 x[3], x[4], "par")) / pp 1462 1463 1464 1465 } 1466 1467 1468 1469 1470 1471 ddof.deta <- dtheta.deta(Dof, .ldof , earg = .edof ) 1472 drho.deta <- dtheta.deta(Rho, .lrho , earg = .erho ) 1473 1474 ans <- c(w) * cbind(dl.ddof * ddof.deta, 1475 dl.drho * drho.deta) 1476 ans <- ans[, interleave.VGAM(M, M1 = M1)] 1477 ans 1478 }), list( .lrho = lrho, .ldof = ldof, 1479 .erho = erho, .edof = edof, 1480 .imethod = imethod ))), 1481 1482 weight = eval(substitute(expression({ 1483 wz11 <- beta(2, Dof / 2) / Dof - 1484 beta(3, Dof / 2) * (Dof + 2) / (4 * Dof) 1485 wz12 <- -Rho / (2 * (1 - Rho^2)) * (beta(2, Dof / 2) - 1486 beta(3, Dof / 2) * (Dof + 2) / 2) 1487 wz22 <- (1 + Rho^2) / (1 - Rho^2)^2 + 1488 (Dof^2 + 2 * Dof) * Rho^2 * 1489 beta(3, Dof / 2) / (4 * (1 - Rho^2)^2) 1490 wz22 <- wz22 + (Dof^2 + 2 * Dof) * (2 - 3 * Rho^2 + Rho^6) * 1491 beta(3, Dof / 2) / (16 * (1 - Rho^2)^4) 1492 wz22 <- wz22 + (Dof^2 + 2 * Dof) * (1 + Rho^2) * # Replace - by + ??? 1493 beta(2, Dof / 2) / (4 * (1 - Rho^2)^2) # denom == 4 or 2 ??? 1494 ned2l.ddof2 <- wz11 1495 ned2l.ddofrho <- wz12 1496 ned2l.drho2 <- wz22 1497 1498 wz <- array(c(c(w) * ned2l.ddof2 * ddof.deta^2, 1499 c(w) * ned2l.drho2 * drho.deta^2, 1500 c(w) * ned2l.ddofrho * ddof.deta * drho.deta), 1501 dim = c(n, M / M1, 3)) 1502 wz <- arwz2wz(wz, M = M, M1 = M1) 1503 wz 1504 }), list( .lrho = lrho, .ldof = ldof, 1505 .erho = erho, .edof = edof, 1506 .imethod = imethod )))) 1507} 1508 1509 1510 1511 1512 1513 1514 1515dbinormcop <- function(x1, x2, rho = 0, log = FALSE) { 1516 if (!is.logical(log.arg <- log) || length(log) != 1) 1517 stop("bad input for argument 'log'") 1518 rm(log) 1519 1520 x1 <- qnorm(x1) 1521 x2 <- qnorm(x2) 1522 1523 logdensity <- (2 * rho * x1 * x2 - 1524 rho^2 * (x1^2 + x2^2)) / (2 * (1 - rho^2)) - 1525 0.5 * log1p(-rho^2) 1526 1527 if (log.arg) logdensity else exp(logdensity) 1528} 1529 1530 1531 1532 1533pbinormcop <- function(q1, q2, rho = 0) { 1534 1535 if (!is.Numeric(q1, positive = TRUE) || 1536 any(q1 >= 1)) 1537 stop("bad input for argument 'q1'") 1538 if (!is.Numeric(q2, positive = TRUE) || 1539 any(q2 >= 1)) 1540 stop("bad input for argument 'q2'") 1541 if (!is.Numeric(rho) || 1542 any(abs(rho) >= 1)) 1543 stop("bad input for argument 'rho'") 1544 1545 pbinorm(q1 = qnorm(q1), q2 = qnorm(q2), cov12 = rho) 1546} 1547 1548 1549rbinormcop <- function(n, rho = 0 #, inverse = FALSE 1550 ) { 1551 1552 inverse <- FALSE 1553 ymat <- rbinorm(n = n, cov12 = rho) 1554 if (inverse) { 1555 ymat 1556 } else { 1557 cbind(y1 = pnorm(ymat[, 1]), 1558 y2 = pnorm(ymat[, 2])) 1559 } 1560} 1561 1562 1563 1564 1565 1566 binormalcop <- function(lrho = "rhobitlink", 1567 irho = NULL, 1568 imethod = 1, 1569 parallel = FALSE, 1570 zero = NULL) { 1571 1572 1573 1574 apply.parint <- TRUE 1575 1576 1577 lrho <- as.list(substitute(lrho)) 1578 erho <- link2list(lrho) 1579 lrho <- attr(erho, "function.name") 1580 1581 1582 if (length(irho) && 1583 any(abs(irho) >= 1)) 1584 stop("argument 'irho' must have values in (-1,1)") 1585 1586 1587 1588 if (!is.Numeric(imethod, length.arg = 1, 1589 integer.valued = TRUE, positive = TRUE) || 1590 imethod > 3) 1591 stop("argument 'imethod' must be 1 or 2 or 3") 1592 1593 new("vglmff", 1594 blurb = c("Gaussian copula ", 1595 "(based on the bivariate normal distribution)\n", 1596 "Links: ", 1597 namesof("rho", lrho, earg = erho)), 1598 1599 constraints = eval(substitute(expression({ 1600 constraints <- cm.VGAM(matrix(1, M, 1), x = x, 1601 bool = .parallel , 1602 constraints = constraints, 1603 apply.int = .apply.parint ) 1604 1605 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 1606 predictors.names = predictors.names, 1607 M1 = 1) 1608 }), list( .zero = zero, 1609 .apply.parint = apply.parint, 1610 .parallel = parallel ))), 1611 1612 infos = eval(substitute(function(...) { 1613 list(M1 = 1, 1614 Q1 = 2, 1615 parameters.names = c("rho"), 1616 apply.parint = .apply.parint , 1617 parallel = .parallel , 1618 zero = .zero ) 1619 }, list( .zero = zero, 1620 .apply.parint = apply.parint, 1621 .parallel = parallel ))), 1622 1623 initialize = eval(substitute(expression({ 1624 M1 <- 1 1625 Q1 <- 2 1626 1627 temp5 <- 1628 w.y.check(w = w, y = y, 1629 Is.positive.y = TRUE, 1630 ncol.w.max = Inf, 1631 ncol.y.max = Inf, 1632 ncol.y.min = Q1, 1633 out.wy = TRUE, 1634 colsyperw = Q1, 1635 maximize = TRUE) 1636 1637 w <- temp5$w 1638 y <- temp5$y 1639 1640 1641 ncoly <- ncol(y) 1642 extra$ncoly <- ncoly 1643 extra$M1 <- M1 1644 extra$Q1 <- Q1 1645 M <- M1 * (ncoly / Q1) 1646 mynames1 <- param.names("rho", M / M1, skip1 = TRUE) 1647 predictors.names <- c( 1648 namesof(mynames1, .lrho , earg = .erho , short = TRUE)) 1649 1650 1651 extra$colnames.y <- colnames(y) 1652 1653 if (!length(etastart)) { 1654 1655 rho.init <- matrix(if (length( .irho )) .irho else 0 + NA, 1656 n, M / M1, byrow = TRUE) 1657 1658 if (!length( .irho )) 1659 for (spp. in 1:(M / M1)) { 1660 ymatj <- y[, (Q1 * spp. - 1):(Q1 * spp.)] 1661 1662 1663 rho.init0 <- if ( .imethod == 1) { 1664 sin(kendall.tau(ymatj[, 1], ymatj[, 2], 1665 exact = FALSE, 1666 max.n = 200) * pi / 2) 1667 } else if ( .imethod == 2) { 1668 sin(cor(ymatj[, 1], ymatj[, 2], 1669 method = "spearman") * pi / 6) * 2 1670 } else { 1671 cor(ymatj[, 1], ymatj[, 2]) 1672 } 1673 1674 1675 1676 1677 1678 if (anyNA(rho.init[, spp.])) 1679 rho.init[, spp.] <- rho.init0 1680 } 1681 1682 etastart <- theta2eta(rho.init, .lrho , earg = .erho ) 1683 } 1684 }), list( .imethod = imethod, 1685 .lrho = lrho, 1686 .erho = erho, 1687 .irho = irho ))), 1688 linkinv = eval(substitute(function(eta, extra = NULL) { 1689 NOS <- NCOL(eta) / c(M1 = 1) 1690 Q1 <- 2 1691 fv.mat <- matrix(0.5, NROW(eta), NOS * Q1) 1692 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 1693 } , list( .lrho = lrho, 1694 .erho = erho ))), 1695 1696 last = eval(substitute(expression({ 1697 M1 <- extra$M1 1698 Q1 <- extra$Q1 1699 misc$link <- rep_len( .lrho , M) 1700 temp.names <- mynames1 1701 names(misc$link) <- temp.names 1702 1703 misc$earg <- vector("list", M) 1704 names(misc$earg) <- temp.names 1705 for (ii in 1:M) { 1706 misc$earg[[ii]] <- .erho 1707 } 1708 1709 misc$M1 <- M1 1710 misc$Q1 <- Q1 1711 misc$imethod <- .imethod 1712 misc$expected <- TRUE 1713 misc$parallel <- .parallel 1714 misc$apply.parint <- .apply.parint 1715 misc$multipleResponses <- TRUE 1716 1717 }) , list( .imethod = imethod, 1718 .parallel = parallel, 1719 .apply.parint = apply.parint, 1720 .lrho = lrho, 1721 .erho = erho ))), 1722 loglikelihood = eval(substitute( 1723 function(mu, y, w, residuals = FALSE, eta, 1724 extra = NULL, 1725 summation = TRUE) { 1726 Rho <- eta2theta(eta, .lrho , earg = .erho ) 1727 1728 if (residuals) { 1729 stop("loglikelihood residuals not implemented yet") 1730 } else { 1731 Yindex1 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) - 1 1732 Yindex2 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) 1733 ll.elts <- 1734 c(w) * dbinormcop(x1 = y[, Yindex1, drop = FALSE], 1735 x2 = y[, Yindex2, drop = FALSE], 1736 rho = Rho, log = TRUE) 1737 if (summation) { 1738 sum(ll.elts) 1739 } else { 1740 ll.elts 1741 } 1742 } 1743 } , list( .lrho = lrho, 1744 .erho = erho, 1745 .imethod = imethod ))), 1746 vfamily = c("binormalcop"), 1747 validparams = eval(substitute(function(eta, y, extra = NULL) { 1748 Rho <- eta2theta(eta, .lrho , earg = .erho ) 1749 okay1 <- all(is.finite(Rho)) && all(abs(Rho) < 1) 1750 okay1 1751 }, list( .lrho = lrho, .erho = erho, .imethod = imethod ))), 1752 1753 1754 1755 simslot = eval(substitute( 1756 function(object, nsim) { 1757 1758 pwts <- if (length(pwts <- object@prior.weights) > 0) 1759 pwts else weights(object, type = "prior") 1760 if (any(pwts != 1)) 1761 warning("ignoring prior weights") 1762 eta <- predict(object) 1763 Rho <- eta2theta(eta, .lrho , earg = .erho ) 1764 rbinormcop(nsim * length(Rho), 1765 rho = c(Rho)) 1766 } , list( .lrho = lrho, 1767 .erho = erho ))), 1768 1769 1770 1771 deriv = eval(substitute(expression({ 1772 Rho <- eta2theta(eta, .lrho , earg = .erho ) 1773 Yindex1 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) - 1 1774 Yindex2 <- extra$Q1 * (1:(extra$ncoly/extra$Q1)) 1775 1776 temp7 <- 1 - Rho^2 1777 q.y <- qnorm(y) 1778 1779 dl.drho <- ((1 + Rho^2) * q.y[, Yindex1] * q.y[, Yindex2] - 1780 Rho * (q.y[, Yindex1]^2 + q.y[, Yindex2]^2)) / temp7^2 + 1781 Rho / temp7 1782 1783 drho.deta <- dtheta.deta(Rho, .lrho , earg = .erho ) 1784 1785 c(w) * cbind(dl.drho) * drho.deta 1786 }), list( .lrho = lrho, 1787 .erho = erho, 1788 .imethod = imethod ))), 1789 1790 weight = eval(substitute(expression({ 1791 ned2l.drho <- (1 + Rho^2) / temp7^2 1792 wz <- ned2l.drho * drho.deta^2 1793 c(w) * wz 1794 }), list( .lrho = lrho, 1795 .erho = erho, 1796 .imethod = imethod )))) 1797} 1798 1799 1800 1801 1802 1803 1804 1805 1806 1807 1808 1809bilogistic.control <- function(save.weights = TRUE, ...) { 1810 list(save.weights = save.weights) 1811} 1812 1813 1814 bilogistic <- function(llocation = "identitylink", 1815 lscale = "loglink", 1816 iloc1 = NULL, iscale1 = NULL, 1817 iloc2 = NULL, iscale2 = NULL, 1818 imethod = 1, 1819 nsimEIM = 250, 1820 zero = NULL) { 1821 1822 llocat <- as.list(substitute(llocation)) 1823 elocat <- link2list(llocat) 1824 llocat <- attr(elocat, "function.name") 1825 1826 lscale <- as.list(substitute(lscale)) 1827 escale <- link2list(lscale) 1828 lscale <- attr(escale, "function.name") 1829 1830 1831 1832 if (!is.Numeric(nsimEIM, length.arg = 1, 1833 integer.valued = TRUE) || 1834 nsimEIM <= 50) 1835 warning("'nsimEIM' should be an integer greater than 50") 1836 1837 1838 1839 if (!is.Numeric(imethod, length.arg = 1, 1840 integer.valued = TRUE, positive = TRUE) || 1841 imethod > 2) stop("argument 'imethod' must be 1 or 2") 1842 1843 new("vglmff", 1844 blurb = c("Bivariate logistic distribution\n\n", 1845 "Link: ", 1846 namesof("location1", llocat, elocat), ", ", 1847 namesof("scale1", lscale, escale), ", ", 1848 namesof("location2", llocat, elocat), ", ", 1849 namesof("scale2", lscale, escale), "\n", "\n", 1850 "Means: location1, location2"), 1851 constraints = eval(substitute(expression({ 1852 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 1853 predictors.names = predictors.names, 1854 M1 = 4) 1855 }), list( .zero = zero))), 1856 1857 1858 infos = eval(substitute(function(...) { 1859 list(M1 = 4, 1860 Q1 = 2, 1861 expected = FALSE, 1862 parameters.names = 1863 c("location1", "scale1", "location2", "scale2"), 1864 multipleResponses = FALSE, 1865 zero = .zero ) 1866 }, list( .zero = zero 1867 ))), 1868 1869 1870 initialize = eval(substitute(expression({ 1871 1872 temp5 <- 1873 w.y.check(w = w, y = y, 1874 ncol.w.max = 1, 1875 ncol.y.max = 2, 1876 ncol.y.min = 2, 1877 out.wy = TRUE, 1878 colsyperw = 2, 1879 maximize = TRUE) 1880 w <- temp5$w 1881 y <- temp5$y 1882 1883 extra$colnames.y <- colnames(y) 1884 1885 1886 predictors.names <- 1887 c(namesof("location1", .llocat, .elocat , tag = FALSE), 1888 namesof("scale1", .lscale, .escale , tag = FALSE), 1889 namesof("location2", .llocat, .elocat , tag = FALSE), 1890 namesof("scale2", .lscale, .escale , tag = FALSE)) 1891 1892 if (!length(etastart)) { 1893 if ( .imethod == 1) { 1894 locat.init1 <- y[, 1] 1895 scale.init1 <- sqrt(3) * sd(y[, 1]) / pi 1896 locat.init2 <- y[, 2] 1897 scale.init2 <- sqrt(3) * sd(y[, 2]) / pi 1898 } else { 1899 locat.init1 <- median(rep(y[, 1], w)) 1900 locat.init2 <- median(rep(y[, 2], w)) 1901 const4 <- sqrt(3) / (sum(w) * pi) 1902 scale.init1 <- const4 * sum(c(w) *(y[, 1] - locat.init1)^2) 1903 scale.init2 <- const4 * sum(c(w) *(y[, 2] - locat.init2)^2) 1904 } 1905 loc1.init <- if (length( .iloc1 )) rep_len( .iloc1 , n) else 1906 rep_len(locat.init1, n) 1907 loc2.init <- if (length( .iloc2 )) rep_len( .iloc2 , n) else 1908 rep_len(locat.init2, n) 1909 scale1.init <- if (length( .iscale1 )) rep_len( .iscale1 , n) else 1910 rep_len(1, n) 1911 scale2.init <- if (length( .iscale2 )) rep_len( .iscale2 , n) else 1912 rep_len(1, n) 1913 1914 if ( .llocat == "loglink") 1915 locat.init1 <- abs(locat.init1) + 0.001 1916 if ( .llocat == "loglink") 1917 locat.init2 <- abs(locat.init2) + 0.001 1918 1919 etastart <- 1920 cbind(theta2eta(locat.init1, .llocat , .elocat ), 1921 theta2eta(scale1.init, .lscale , .escale ), 1922 theta2eta(locat.init2, .llocat , .elocat ), 1923 theta2eta(scale2.init, .lscale , .escale )) 1924 } 1925 }), list(.imethod = imethod, 1926 .iloc1 = iloc1, .iloc2 = iloc2, 1927 .llocat = llocat, .lscale = lscale, 1928 .elocat = elocat, .escale = escale, 1929 .iscale1 = iscale1, .iscale2 = iscale2))), 1930 linkinv = function(eta, extra = NULL) { 1931 NOS <- NCOL(eta) / c(M1 = 4) 1932 fv.mat <- eta[, 1:2] 1933 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 1934 }, 1935 last = eval(substitute(expression({ 1936 misc$link <- c(location1 = .llocat , scale1 = .lscale , 1937 location2 = .llocat , scale2 = .lscale ) 1938 1939 misc$earg <- list(location1 = .elocat , scale1 = .escale , 1940 location2 = .elocat , scale2 = .escale ) 1941 1942 }), list( .llocat = llocat, .lscale = lscale, 1943 .elocat = elocat, .escale = escale ))), 1944 loglikelihood = eval(substitute( 1945 function(mu, y, w, residuals = FALSE, eta, 1946 extra = NULL, 1947 summation = TRUE) { 1948 locat1 <- eta2theta(eta[, 1], .llocat , .elocat ) 1949 Scale1 <- eta2theta(eta[, 2], .lscale , .escale ) 1950 locat2 <- eta2theta(eta[, 3], .llocat , .elocat ) 1951 Scale2 <- eta2theta(eta[, 4], .lscale , .escale ) 1952 1953 zedd1 <- (y[, 1]-locat1) / Scale1 1954 zedd2 <- (y[, 2]-locat2) / Scale2 1955 1956 if (residuals) { 1957 stop("loglikelihood residuals not implemented yet") 1958 } else { 1959 ll.elts <- 1960 c(w) * (-zedd1 - zedd2 - 1961 3 * log1p(exp(-zedd1) + exp(-zedd2)) - 1962 log(Scale1) - log(Scale2)) 1963 if (summation) { 1964 sum(ll.elts) 1965 } else { 1966 ll.elts 1967 } 1968 } 1969 }, list( .llocat = llocat, .lscale = lscale, 1970 .elocat = elocat, .escale = escale ))), 1971 vfamily = c("bilogistic"), 1972 validparams = eval(substitute(function(eta, y, extra = NULL) { 1973 locat1 <- eta2theta(eta[, 1], .llocat , .elocat ) 1974 Scale1 <- eta2theta(eta[, 2], .lscale , .escale ) 1975 locat2 <- eta2theta(eta[, 3], .llocat , .elocat ) 1976 Scale2 <- eta2theta(eta[, 4], .lscale , .escale ) 1977 okay1 <- all(is.finite(locat1)) && 1978 all(is.finite(Scale1)) && all(0 < Scale1) && 1979 all(is.finite(locat2)) && 1980 all(is.finite(Scale2)) && all(0 < Scale2) 1981 okay1 1982 }, list( .llocat = llocat, .lscale = lscale, 1983 .elocat = elocat, .escale = escale ))), 1984 1985 1986 simslot = eval(substitute( 1987 function(object, nsim) { 1988 1989 pwts <- if (length(pwts <- object@prior.weights) > 0) 1990 pwts else weights(object, type = "prior") 1991 if (any(pwts != 1)) 1992 warning("ignoring prior weights") 1993 eta <- predict(object) 1994 locat1 <- eta2theta(eta[, 1], .llocat , .elocat ) 1995 Scale1 <- eta2theta(eta[, 2], .lscale , .escale ) 1996 locat2 <- eta2theta(eta[, 3], .llocat , .elocat ) 1997 Scale2 <- eta2theta(eta[, 4], .lscale , .escale ) 1998 rbilogis(nsim * length(locat1), 1999 loc1 = locat1, scale1 = Scale1, 2000 loc2 = locat2, scale2 = Scale2) 2001 }, list( .llocat = llocat, .lscale = lscale, 2002 .elocat = elocat, .escale = escale ))), 2003 2004 2005 2006 2007 deriv = eval(substitute(expression({ 2008 locat1 <- eta2theta(eta[, 1], .llocat , .elocat ) 2009 Scale1 <- eta2theta(eta[, 2], .lscale , .escale ) 2010 locat2 <- eta2theta(eta[, 3], .llocat , .elocat ) 2011 Scale2 <- eta2theta(eta[, 4], .lscale , .escale ) 2012 2013 zedd1 <- (y[, 1] - locat1) / Scale1 2014 zedd2 <- (y[, 2] - locat2) / Scale2 2015 ezedd1 <- exp(-zedd1) 2016 ezedd2 <- exp(-zedd2) 2017 denom <- 1 + ezedd1 + ezedd2 2018 2019 dl.dlocat1 <- (1 - 3 * ezedd1 / denom) / Scale1 2020 dl.dlocat2 <- (1 - 3 * ezedd2 / denom) / Scale2 2021 dl.dscale1 <- (zedd1 - 1 - 3 * ezedd1 * zedd1 / denom) / Scale1 2022 dl.dscale2 <- (zedd2 - 1 - 3 * ezedd2 * zedd2 / denom) / Scale2 2023 2024 dlocat1.deta <- dtheta.deta(locat1, .llocat , .elocat ) 2025 dlocat2.deta <- dtheta.deta(locat2, .llocat , .elocat ) 2026 dscale1.deta <- dtheta.deta(Scale1, .lscale , .escale ) 2027 dscale2.deta <- dtheta.deta(Scale2, .lscale , .escale ) 2028 2029 derivnew <- c(w) * cbind(dl.dlocat1 * dlocat1.deta, 2030 dl.dscale1 * dscale1.deta, 2031 dl.dlocat2 * dlocat2.deta, 2032 dl.dscale2 * dscale2.deta) 2033 derivnew 2034 }), list( .llocat = llocat, .lscale = lscale, 2035 .elocat = elocat, .escale = escale ))), 2036 weight = eval(substitute(expression({ 2037 run.varcov <- 0 2038 ind1 <- iam(NA_real_, NA_real_, M = M, both = TRUE, diag = TRUE) 2039 for (ii in 1:( .nsimEIM )) { 2040 ysim <- rbilogis( .nsimEIM * length(locat1), 2041 loc1 = locat1, scale1 = Scale1, 2042 loc2 = locat2, scale2 = Scale2) 2043 2044 zedd1 <- (ysim[, 1] - locat1) / Scale1 2045 zedd2 <- (ysim[, 2] - locat2) / Scale2 2046 ezedd1 <- exp(-zedd1) 2047 ezedd2 <- exp(-zedd2) 2048 denom <- 1 + ezedd1 + ezedd2 2049 2050 dl.dlocat1 <- (1 - 3 * ezedd1 / denom) / Scale1 2051 dl.dlocat2 <- (1 - 3 * ezedd2 / denom) / Scale2 2052 dl.dscale1 <- (zedd1 - 1 - 3 * ezedd1 * zedd1 / denom) / Scale1 2053 dl.dscale2 <- (zedd2 - 1 - 3 * ezedd2 * zedd2 / denom) / Scale2 2054 2055 2056 rm(ysim) 2057 temp3 <- cbind(dl.dlocat1, 2058 dl.dscale1, 2059 dl.dlocat2, 2060 dl.dscale2) 2061 run.varcov <- run.varcov + temp3[, ind1$row] * temp3[, ind1$col] 2062 } # ii 2063 run.varcov <- run.varcov / .nsimEIM 2064 wz <- if (intercept.only) 2065 matrix(colMeans(run.varcov, na.rm = FALSE), 2066 n, ncol(run.varcov), byrow = TRUE) else run.varcov 2067 dthetas.detas <- cbind(dlocat1.deta, 2068 dscale1.deta, 2069 dlocat2.deta, 2070 dscale2.deta) 2071 wz <- wz * dthetas.detas[, ind1$row] * dthetas.detas[, ind1$col] 2072 c(w) * wz 2073 }), list( .lscale = lscale, 2074 .escale = escale, 2075 .llocat = llocat, 2076 .nsimEIM = nsimEIM )))) 2077} 2078 2079 2080 2081 2082 2083 2084dbilogis <- function(x1, x2, loc1 = 0, scale1 = 1, 2085 loc2 = 0, scale2 = 1, log = FALSE) { 2086 if (!is.logical(log.arg <- log) || length(log) != 1) 2087 stop("bad input for argument 'log'") 2088 rm(log) 2089 2090 2091 2092 2093 L <- max(length(x1), length(x2), 2094 length(loc1), length(loc2), 2095 length(scale1), length(scale2)) 2096 if (length(x1 ) != L) x1 <- rep_len(x1, L) 2097 if (length(x2 ) != L) x2 <- rep_len(x2, L) 2098 if (length(loc1 ) != L) loc1 <- rep_len(loc1, L) 2099 if (length(loc2 ) != L) loc2 <- rep_len(loc2, L) 2100 if (length(scale1) != L) scale1 <- rep_len(scale1, L) 2101 if (length(scale2) != L) scale2 <- rep_len(scale2, L) 2102 zedd1 <- (x1 - loc1) / scale1 2103 zedd2 <- (x2 - loc2) / scale2 2104 2105 2106 2107 2108 logdensity <- log(2) - zedd1 - zedd2 - log(scale1) - 2109 log(scale1) - 3 * log1p(exp(-zedd1) + exp(-zedd2)) 2110 2111 2112 logdensity[x1 == -Inf | x2 == -Inf] <- log(0) # 20141216 KaiH 2113 2114 2115 if (log.arg) logdensity else exp(logdensity) 2116} 2117 2118 2119 2120pbilogis <- 2121 function(q1, q2, loc1 = 0, scale1 = 1, loc2 = 0, scale2 = 1) { 2122 2123 ans <- 1 / (1 + exp(-(q1-loc1)/scale1) + exp(-(q2-loc2)/scale2)) 2124 ans[scale1 <= 0] <- NA 2125 ans[scale2 <= 0] <- NA 2126 ans 2127} 2128 2129 2130 2131rbilogis <- function(n, loc1 = 0, scale1 = 1, loc2 = 0, scale2 = 1) { 2132 2133 2134 y1 <- rlogis(n = n, location = loc1, scale = scale1) 2135 ezedd1 <- exp(-(y1-loc1)/scale1) 2136 y2 <- loc2 - scale2 * 2137 log(1/sqrt(runif(n) / (1 + ezedd1)^2) - 1 - ezedd1) 2138 ans <- cbind(y1, y2) 2139 ans[scale2 <= 0, ] <- NA 2140 ans 2141} 2142 2143 2144 2145 2146 freund61 <- 2147 function(la = "loglink", 2148 lap = "loglink", 2149 lb = "loglink", 2150 lbp = "loglink", 2151 ia = NULL, iap = NULL, ib = NULL, ibp = NULL, 2152 independent = FALSE, 2153 zero = NULL) { 2154 la <- as.list(substitute(la)) 2155 ea <- link2list(la) 2156 la <- attr(ea, "function.name") 2157 2158 lap <- as.list(substitute(lap)) 2159 eap <- link2list(lap) 2160 lap <- attr(eap, "function.name") 2161 2162 lb <- as.list(substitute(lb)) 2163 eb <- link2list(lb) 2164 lb <- attr(eb, "function.name") 2165 2166 2167 lbp <- as.list(substitute(lbp)) 2168 ebp <- link2list(lbp) 2169 lbp <- attr(ebp, "function.name") 2170 2171 2172 2173 new("vglmff", 2174 blurb = c("Freund (1961) bivariate exponential distribution\n", 2175 "Links: ", 2176 namesof("a", la, earg = ea ), ", ", 2177 namesof("ap", lap, earg = eap), ", ", 2178 namesof("b", lb, earg = eb ), ", ", 2179 namesof("bp", lbp, earg = ebp)), 2180 constraints = eval(substitute(expression({ 2181 M1 <- 4 2182 Q1 <- 2 2183 constraints <- cm.VGAM(matrix(c(1, 1,0,0, 0,0, 1, 1), M, 2), x = x, 2184 bool = .independent , 2185 constraints = constraints, 2186 apply.int = TRUE) 2187 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 2188 predictors.names = predictors.names, 2189 M1 = 4) 2190 }), list( .independent = independent, .zero = zero))), 2191 2192 2193 2194 infos = eval(substitute(function(...) { 2195 list(M1 = 4, 2196 Q1 = 2, 2197 expected = TRUE, 2198 multipleResponses = FALSE, 2199 parameters.names = c("a", "ap", "b", "bp"), 2200 la = .la , 2201 lap = .lap , 2202 lb = .lb , 2203 lbp = .lbp , 2204 independent = .independent , 2205 zero = .zero ) 2206 }, list( .zero = zero, 2207 .la = la , 2208 .lap = lap , 2209 .lb = lb , 2210 .lbp = lbp , 2211 .independent = independent ))), 2212 2213 2214 initialize = eval(substitute(expression({ 2215 2216 temp5 <- 2217 w.y.check(w = w, y = y, 2218 ncol.w.max = 1, 2219 ncol.y.max = 2, 2220 ncol.y.min = 2, 2221 out.wy = TRUE, 2222 colsyperw = 2, 2223 maximize = TRUE) 2224 w <- temp5$w 2225 y <- temp5$y 2226 2227 2228 predictors.names <- 2229 c(namesof("a", .la , earg = .ea , short = TRUE), 2230 namesof("ap", .lap , earg = .eap , short = TRUE), 2231 namesof("b", .lb , earg = .eb , short = TRUE), 2232 namesof("bp", .lbp , earg = .ebp , short = TRUE)) 2233 extra$y1.lt.y2 = y[, 1] < y[, 2] 2234 2235 if (!(arr <- sum(extra$y1.lt.y2)) || arr == n) 2236 stop("identifiability problem: either all y1<y2 or y2<y1") 2237 2238 if (!length(etastart)) { 2239 sumx <- sum(y[ extra$y1.lt.y2, 1]); 2240 sumxp <- sum(y[!extra$y1.lt.y2, 1]) 2241 sumy <- sum(y[ extra$y1.lt.y2, 2]); 2242 sumyp <- sum(y[!extra$y1.lt.y2, 2]) 2243 2244 if (FALSE) { # Noise: 2245 arr <- min(arr + n/10, n*0.95) 2246 sumx <- sumx * 1.1; sumxp <- sumxp * 1.2; 2247 sumy <- sumy * 1.2; sumyp <- sumyp * 1.3; 2248 } 2249 ainit <- if (length( .ia )) rep_len( .ia , n) else 2250 arr / (sumx + sumyp) 2251 apinit <- if (length( .iap )) rep_len( .iap , n) else 2252 (n-arr) / (sumxp - sumyp) 2253 binit <- if (length( .ib )) rep_len( .ib , n) else 2254 (n-arr) / (sumx + sumyp) 2255 bpinit <- if (length( .ibp )) rep_len( .ibp , n) else 2256 arr / (sumy - sumx) 2257 2258 etastart <- 2259 cbind(theta2eta(rep_len(ainit, n), .la , earg = .ea ), 2260 theta2eta(rep_len(apinit, n), .lap , earg = .eap ), 2261 theta2eta(rep_len(binit, n), .lb , earg = .eb ), 2262 theta2eta(rep_len(bpinit, n), .lbp , earg = .ebp )) 2263 } 2264 }), list( .la = la, .lap = lap, .lb = lb, .lbp = lbp, 2265 .ea = ea, .eap = eap, .eb = eb, .ebp = ebp, 2266 .ia = ia, .iap = iap, .ib = ib, .ibp = ibp))), 2267 linkinv = eval(substitute(function(eta, extra = NULL) { 2268 NOS <- NCOL(eta) / c(M1 = 4) 2269 alpha <- eta2theta(eta[, 1], .la, earg = .ea ) 2270 alphap <- eta2theta(eta[, 2], .lap, earg = .eap ) 2271 beta <- eta2theta(eta[, 3], .lb, earg = .eb ) 2272 betap <- eta2theta(eta[, 4], .lbp, earg = .ebp ) 2273 fv.mat <- cbind((alphap + beta) / (alphap * (alpha + beta)), 2274 (alpha + betap) / (betap * (alpha + beta))) 2275 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 2276 }, list( .la = la, .lap = lap, .lb = lb, .lbp = lbp, 2277 .ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))), 2278 last = eval(substitute(expression({ 2279 misc$link <- c("a" = .la , "ap" = .lap , "b" = .lb , "bp" = .lbp ) 2280 misc$earg <- list("a" = .ea , "ap" = .eap , "b" = .eb , "bp" = .ebp ) 2281 2282 misc$multipleResponses <- FALSE 2283 }), list( .la = la, .lap = lap, .lb = lb, .lbp = lbp, 2284 .ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))), 2285 loglikelihood = eval(substitute( 2286 function(mu, y, w, residuals = FALSE, eta, 2287 extra = NULL, 2288 summation = TRUE) { 2289 alpha <- eta2theta(eta[, 1], .la, earg = .ea ) 2290 alphap <- eta2theta(eta[, 2], .lap, earg = .eap ) 2291 beta <- eta2theta(eta[, 3], .lb, earg = .eb ) 2292 betap <- eta2theta(eta[, 4], .lbp, earg = .ebp ) 2293 if (residuals) { 2294 stop("loglikelihood residuals not implemented yet") 2295 } else { 2296 tmp88 <- extra$y1.lt.y2 2297 ell1 <- log(alpha[tmp88]) + log(betap[tmp88]) - 2298 betap[tmp88] * y[tmp88, 2] - 2299 (alpha+beta-betap)[tmp88] * y[tmp88, 1] 2300 ell2 <- log(beta[!tmp88]) + log(alphap[!tmp88]) - 2301 alphap[!tmp88] * y[!tmp88, 1] - 2302 (alpha+beta-alphap)[!tmp88] * y[!tmp88, 2] 2303 all.vec <- alpha 2304 all.vec[ tmp88] <- ell1 2305 all.vec[!tmp88] <- ell2 2306 ll.elts <- c(w) * all.vec 2307 if (summation) { 2308 sum(ll.elts) 2309 } else { 2310 ll.elts 2311 } 2312 } 2313 }, list( .la = la, .lap = lap, .lb = lb, .lbp = lbp, 2314 .ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))), 2315 vfamily = c("freund61"), 2316 validparams = eval(substitute(function(eta, y, extra = NULL) { 2317 alpha <- eta2theta(eta[, 1], .la, earg = .ea ) 2318 alphap <- eta2theta(eta[, 2], .lap, earg = .eap ) 2319 beta <- eta2theta(eta[, 3], .lb, earg = .eb ) 2320 betap <- eta2theta(eta[, 4], .lbp, earg = .ebp ) 2321 okay1 <- all(is.finite(alpha )) && all(0 < alpha ) && 2322 all(is.finite(alphap)) && all(0 < alphap) && 2323 all(is.finite(beta )) && all(0 < beta ) && 2324 all(is.finite(betap )) && all(0 < betap ) 2325 okay1 2326 }, list( .la = la, .lap = lap, .lb = lb, .lbp = lbp, 2327 .ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))), 2328 deriv = eval(substitute(expression({ 2329 tmp88 <- extra$y1.lt.y2 2330 alpha <- eta2theta(eta[, 1], .la, earg = .ea ) 2331 alphap <- eta2theta(eta[, 2], .lap, earg = .eap ) 2332 beta <- eta2theta(eta[, 3], .lb, earg = .eb ) 2333 betap <- eta2theta(eta[, 4], .lbp, earg = .ebp ) 2334 2335 dalpha.deta <- dtheta.deta(alpha, .la, earg = .ea ) 2336 dalphap.deta <- dtheta.deta(alphap, .lap, earg = .eap ) 2337 dbeta.deta <- dtheta.deta(beta, .lb, earg = .eb ) 2338 dbetap.deta <- dtheta.deta(betap, .lbp, earg = .ebp ) 2339 2340 d1 <- 1/alpha - y[, 1] 2341 d1[!tmp88] <- -y[!tmp88, 2] 2342 d2 <- 0 * alphap 2343 d2[!tmp88] <- 1/alphap[!tmp88] - y[!tmp88, 1] + y[!tmp88, 2] 2344 d3 <- -y[, 1] 2345 d3[!tmp88] <- 1/beta[!tmp88] - y[!tmp88, 2] 2346 d4 <- 1/betap - y[, 2] + y[, 1] 2347 d4[!tmp88] <- 0 2348 2349 c(w) * cbind(d1 * dalpha.deta, 2350 d2 * dalphap.deta, 2351 d3 * dbeta.deta, 2352 d4 * dbetap.deta) 2353 }), list( .la = la, .lap = lap, .lb = lb, .lbp = lbp, 2354 .ea = ea, .eap = eap, .eb = eb, .ebp = ebp ))), 2355 weight = eval(substitute(expression({ 2356 py1.lt.y2 <- alpha / (alpha+beta) 2357 d11 <- py1.lt.y2 / alpha^2 2358 d22 <- (1-py1.lt.y2) / alphap^2 2359 d33 <- (1-py1.lt.y2) / beta^2 2360 d44 <- py1.lt.y2 / betap^2 2361 2362 wz <- matrix(0, n, M) # diagonal 2363 wz[, iam(1, 1, M)] <- dalpha.deta^2 * d11 2364 wz[, iam(2, 2, M)] <- dalphap.deta^2 * d22 2365 wz[, iam(3, 3, M)] <- dbeta.deta^2 * d33 2366 wz[, iam(4, 4, M)] <- dbetap.deta^2 * d44 2367 2368 c(w) * wz 2369 }), list( .la = la, .lap = lap, .lb = lb, .lbp = lbp, 2370 .ea = ea, .eap = eap, .eb = eb, .ebp = ebp )))) 2371} 2372 2373 2374 2375 2376 2377 2378 2379 2380 bigamma.mckay <- function(lscale = "loglink", 2381 lshape1 = "loglink", 2382 lshape2 = "loglink", 2383 iscale = NULL, 2384 ishape1 = NULL, 2385 ishape2 = NULL, 2386 imethod = 1, 2387 zero = "shape") { 2388 lscale <- as.list(substitute(lscale)) 2389 escale <- link2list(lscale) 2390 lscale <- attr(escale, "function.name") 2391 2392 lshape1 <- as.list(substitute(lshape1)) 2393 eshape1 <- link2list(lshape1) 2394 lshape1 <- attr(eshape1, "function.name") 2395 2396 lshape2 <- as.list(substitute(lshape2)) 2397 eshape2 <- link2list(lshape2) 2398 lshape2 <- attr(eshape2, "function.name") 2399 2400 2401 if (!is.null(iscale)) 2402 if (!is.Numeric(iscale, positive = TRUE)) 2403 stop("argument 'iscale' must be positive or NULL") 2404 if (!is.null(ishape1)) 2405 if (!is.Numeric(ishape1, positive = TRUE)) 2406 stop("argument 'ishape1' must be positive or NULL") 2407 if (!is.null(ishape2)) 2408 if (!is.Numeric(ishape2, positive = TRUE)) 2409 stop("argument 'ishape2' must be positive or NULL") 2410 2411 if (!is.Numeric(imethod, length.arg = 1, 2412 integer.valued = TRUE, positive = TRUE) || 2413 imethod > 2.5) 2414 stop("argument 'imethod' must be 1 or 2") 2415 2416 2417 2418 new("vglmff", 2419 blurb = c("Bivariate gamma: McKay's distribution\n", 2420 "Links: ", 2421 namesof("scale", lscale, earg = escale ), ", ", 2422 namesof("shape1", lshape1, earg = eshape1), ", ", 2423 namesof("shape2", lshape2, earg = eshape2)), 2424 constraints = eval(substitute(expression({ 2425 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 2426 predictors.names = predictors.names, 2427 M1 = 3) 2428 }), list( .zero = zero ))), 2429 2430 2431 2432 infos = eval(substitute(function(...) { 2433 list(M1 = 3, 2434 Q1 = 2, 2435 expected = TRUE, 2436 multipleResponses = FALSE, 2437 parameters.names = c("scale", "shape1", "shape2"), 2438 lscale = .lscale , 2439 lshape1 = .lshape1 , 2440 lshape2 = .lshape2 , 2441 zero = .zero ) 2442 }, list( .zero = zero, 2443 .lscale = lscale , 2444 .lshape1 = lshape1, 2445 .lshape2 = lshape2 ))), 2446 2447 2448 initialize = eval(substitute(expression({ 2449 2450 temp5 <- 2451 w.y.check(w = w, y = y, 2452 ncol.w.max = 1, 2453 ncol.y.max = 2, 2454 ncol.y.min = 2, 2455 out.wy = TRUE, 2456 colsyperw = 2, 2457 maximize = TRUE) 2458 w <- temp5$w 2459 y <- temp5$y 2460 2461 2462 extra$colnames.y <- colnames(y) 2463 2464 if (any(y[, 1] >= y[, 2])) 2465 stop("the second column minus the first column must be a vector ", 2466 "of positive values") 2467 2468 2469 predictors.names <- 2470 c(namesof("scale", .lscale, .escale, short = TRUE), 2471 namesof("shape1", .lshape1, .eshape1, short = TRUE), 2472 namesof("shape2", .lshape2, .eshape2, short = TRUE)) 2473 2474 if (!length(etastart)) { 2475 momentsY <- if ( .imethod == 1) { 2476 cbind(median(y[, 1]), # This may not be monotonic 2477 median(y[, 2])) + 0.01 2478 } else { 2479 cbind(weighted.mean(y[, 1], w), 2480 weighted.mean(y[, 2], w)) 2481 } 2482 2483 mcg2.loglik <- function(thetaval, y, x, w, extraargs) { 2484 ainit <- a <- thetaval 2485 momentsY <- extraargs$momentsY 2486 p <- (1/a) * abs(momentsY[1]) + 0.01 2487 q <- (1/a) * abs(momentsY[2] - momentsY[1]) + 0.01 2488 sum(c(w) * (-(p+q)*log(a) - lgamma(p) - lgamma(q) + 2489 (p - 1)*log(y[, 1]) + 2490 (q - 1)*log(y[, 2]-y[, 1]) - y[, 2] / a )) 2491 } 2492 2493 a.grid <- if (length( .iscale )) c( .iscale ) else 2494 c(0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100) 2495 extraargs <- list(momentsY = momentsY) 2496 ainit <- grid.search(a.grid, objfun = mcg2.loglik, 2497 y = y, x = x, w = w, maximize = TRUE, 2498 extraargs = extraargs) 2499 ainit <- rep_len(if (is.Numeric( .iscale )) .iscale else ainit, n) 2500 pinit <- (1/ainit) * abs(momentsY[1]) + 0.01 2501 qinit <- (1/ainit) * abs(momentsY[2] - momentsY[1]) + 0.01 2502 2503 pinit <- rep_len(if (is.Numeric( .ishape1 )) .ishape1 else pinit, n) 2504 qinit <- rep_len(if (is.Numeric( .ishape2 )) .ishape2 else qinit, n) 2505 2506 etastart <- 2507 cbind(theta2eta(ainit, .lscale), 2508 theta2eta(pinit, .lshape1), 2509 theta2eta(qinit, .lshape2)) 2510 } 2511 }), list( .lscale = lscale, .lshape1 = lshape1, .lshape2 = lshape2, 2512 .escale = escale, .eshape1 = eshape1, .eshape2 = eshape2, 2513 .iscale = iscale, .ishape1 = ishape1, .ishape2 = ishape2, 2514 .imethod = imethod ))), 2515 linkinv = eval(substitute(function(eta, extra = NULL) { 2516 NOS <- NCOL(eta) / c(M1 = 3) 2517 a <- eta2theta(eta[, 1], .lscale , .escale ) 2518 p <- eta2theta(eta[, 2], .lshape1 , .eshape1 ) 2519 q <- eta2theta(eta[, 3], .lshape2 , .eshape2 ) 2520 fv.mat <- cbind("y1" = p*a, 2521 "y2" = (p+q)*a) # Overwrite the colnames: 2522 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 2523 }, list( .lscale = lscale, .lshape1 = lshape1, .lshape2 = lshape2, 2524 .escale = escale, .eshape1 = eshape1, .eshape2 = eshape2 ))), 2525 last = eval(substitute(expression({ 2526 misc$link <- c("scale" = .lscale , 2527 "shape1" = .lshape1 , 2528 "shape2" = .lshape2 ) 2529 2530 misc$earg <- list("scale" = .escale , 2531 "shape1" = .eshape1 , 2532 "shape2" = .eshape2 ) 2533 2534 misc$ishape1 <- .ishape1 2535 misc$ishape2 <- .ishape2 2536 misc$iscale <- .iscale 2537 misc$expected <- TRUE 2538 misc$multipleResponses <- FALSE 2539 }), list( .lscale = lscale, .lshape1 = lshape1, .lshape2 = lshape2, 2540 .escale = escale, .eshape1 = eshape1, .eshape2 = eshape2, 2541 .iscale = iscale, .ishape1 = ishape1, .ishape2 = ishape2, 2542 .imethod = imethod ))), 2543 loglikelihood = eval(substitute( 2544 function(mu, y, w, residuals = FALSE, eta, 2545 extra = NULL, 2546 summation = TRUE) { 2547 a <- eta2theta(eta[, 1], .lscale , .escale ) 2548 p <- eta2theta(eta[, 2], .lshape1 , .eshape1 ) 2549 q <- eta2theta(eta[, 3], .lshape2 , .eshape2 ) 2550 2551 if (residuals) { 2552 stop("loglikelihood residuals not implemented yet") 2553 } else { 2554 ll.elts <- 2555 c(w) * (-(p+q)*log(a) - lgamma(p) - lgamma(q) + 2556 (p - 1)*log(y[, 1]) + (q - 1)*log(y[, 2]-y[, 1]) - 2557 y[, 2] / a) 2558 if (summation) { 2559 sum(ll.elts) 2560 } else { 2561 ll.elts 2562 } 2563 } 2564 }, list( .lscale = lscale, .lshape1 = lshape1, .lshape2 = lshape2, 2565 .escale = escale, .eshape1 = eshape1, .eshape2 = eshape2 ))), 2566 vfamily = c("bigamma.mckay"), 2567 validparams = eval(substitute(function(eta, y, extra = NULL) { 2568 aparam <- eta2theta(eta[, 1], .lscale , .escale ) 2569 shape1 <- eta2theta(eta[, 2], .lshape1 , .eshape1 ) 2570 shape2 <- eta2theta(eta[, 3], .lshape2 , .eshape2 ) 2571 okay1 <- all(is.finite(aparam)) && all(0 < aparam) && 2572 all(is.finite(shape1)) && all(0 < shape1) && 2573 all(is.finite(shape2)) && all(0 < shape2) 2574 okay1 2575 }, list( .lscale = lscale, .lshape1 = lshape1, .lshape2 = lshape2, 2576 .escale = escale, .eshape1 = eshape1, .eshape2 = eshape2 ))), 2577 deriv = eval(substitute(expression({ 2578 aparam <- eta2theta(eta[, 1], .lscale , .escale ) 2579 shape1 <- eta2theta(eta[, 2], .lshape1 , .eshape1 ) 2580 shape2 <- eta2theta(eta[, 3], .lshape2 , .eshape2 ) 2581 2582 dl.da <- (-(shape1+shape2) + y[, 2] / aparam) / aparam 2583 dl.dshape1 <- -log(aparam) - digamma(shape1) + log(y[, 1]) 2584 dl.dshape2 <- -log(aparam) - digamma(shape2) + log(y[, 2]-y[, 1]) 2585 2586 c(w) * cbind(dl.da * dtheta.deta(aparam, .lscale), 2587 dl.dshape1 * dtheta.deta(shape1, .lshape1), 2588 dl.dshape2 * dtheta.deta(shape2, .lshape2)) 2589 }), list( .lscale = lscale, .lshape1 = lshape1, .lshape2 = lshape2, 2590 .escale = escale, .eshape1 = eshape1, .eshape2 = eshape2 ))), 2591 weight = eval(substitute(expression({ 2592 d11 <- (shape1+shape2) / aparam^2 2593 d22 <- trigamma(shape1) 2594 d33 <- trigamma(shape2) 2595 d12 <- 1 / aparam 2596 d13 <- 1 / aparam 2597 d23 <- 0 2598 2599 wz <- matrix(0, n, dimm(M)) 2600 wz[, iam(1, 1, M)] <- dtheta.deta(aparam, .lscale )^2 * d11 2601 wz[, iam(2, 2, M)] <- dtheta.deta(shape1, .lshape1 )^2 * d22 2602 wz[, iam(3, 3, M)] <- dtheta.deta(shape2, .lshape2 )^2 * d33 2603 wz[, iam(1, 2, M)] <- dtheta.deta(aparam, .lscale ) * 2604 dtheta.deta(shape1, .lshape1 ) * d12 2605 wz[, iam(1, 3, M)] <- dtheta.deta(aparam, .lscale ) * 2606 dtheta.deta(shape2, .lshape2 ) * d13 2607 wz[, iam(2, 3, M)] <- dtheta.deta(shape1, .lshape1 ) * 2608 dtheta.deta(shape2, .lshape2 ) * d23 2609 2610 c(w) * wz 2611 }), list( .lscale = lscale, .lshape1 = lshape1, 2612 .lshape2 = lshape2 )))) 2613} 2614 2615 2616 2617 2618 2619 2620 2621 2622 2623 2624 2625rbifrankcop <- function(n, apar) { 2626 use.n <- if ((length.n <- length(n)) > 1) length.n else 2627 if (!is.Numeric(n, integer.valued = TRUE, 2628 length.arg = 1, positive = TRUE)) 2629 stop("bad input for argument 'n'") else n 2630 if (!is.Numeric(apar, positive = TRUE)) 2631 stop("bad input for argument 'apar'") 2632 if (length(apar) != use.n) 2633 apar <- rep_len(apar, use.n) 2634 U <- runif(use.n) 2635 V <- runif(use.n) 2636 2637 T <- apar^U + (apar - apar^U) * V 2638 X <- U 2639 index <- (abs(apar - 1) < .Machine$double.eps) 2640 Y <- U 2641 if (any(!index)) 2642 Y[!index] <- logb(T[!index] / (T[!index] + 2643 (1 - apar[!index]) * V[!index]), 2644 base = apar[!index]) 2645 ans <- matrix(c(X, Y), nrow = use.n, ncol = 2) 2646 if (any(index)) { 2647 ans[index, 1] <- runif(sum(index)) # Uniform density for apar == 1 2648 ans[index, 2] <- runif(sum(index)) 2649 } 2650 ans 2651} 2652 2653 2654pbifrankcop <- function(q1, q2, apar) { 2655 if (!is.Numeric(q1)) stop("bad input for 'q1'") 2656 if (!is.Numeric(q2)) stop("bad input for 'q2'") 2657 if (!is.Numeric(apar, positive = TRUE)) stop("bad input for 'apar'") 2658 2659 L <- max(length(q1), length(q2), length(apar)) 2660 if (length(apar ) != L) apar <- rep_len(apar, L) 2661 if (length(q1 ) != L) q1 <- rep_len(q1, L) 2662 if (length(q2 ) != L) q2 <- rep_len(q2, L) 2663 2664 x <- q1; y <- q2 2665 index <- (x >= 1 & y < 1) | (y >= 1 & x < 1) | 2666 (x <= 0 | y <= 0) | (x >= 1 & y >= 1) | 2667 (abs(apar - 1) < .Machine$double.eps) 2668 ans <- as.numeric(index) 2669 if (any(!index)) 2670 ans[!index] <- logb(1 + ((apar[!index])^(x[!index]) - 1)* 2671 ((apar[!index])^(y[!index]) - 1)/(apar[!index] - 1), 2672 base = apar[!index]) 2673 ind2 <- (abs(apar - 1) < .Machine$double.eps) 2674 ans[ind2] <- x[ind2] * y[ind2] 2675 ans[x >= 1 & y < 1] <- y[x >= 1 & y < 1] # P(Y2 < q2) = q2 2676 ans[y >= 1 & x < 1] <- x[y >= 1 & x < 1] # P(Y1 < q1) = q1 2677 ans[x <= 0 | y <= 0] <- 0 2678 ans[x >= 1 & y >= 1] <- 1 2679 ans 2680} 2681 2682 2683 2684 2685 2686if (FALSE) 2687dbifrank <- function(x1, x2, apar, log = FALSE) { 2688 if (!is.logical(log.arg <- log) || length(log) != 1) 2689 stop("bad input for argument 'log'") 2690 rm(log) 2691 logdens <- (x1+x2)*log(apar) + log(apar-1) + log(log(apar)) - 2692 2 * log(apar - 1 + (apar^x1 - 1) * (apar^x2 - 1)) 2693 2694 if (log.arg) logdens else exp(logdens) 2695} 2696 2697 2698 2699 2700dbifrankcop <- function(x1, x2, apar, log = FALSE) { 2701 if (!is.logical(log.arg <- log) || length(log) != 1) 2702 stop("bad input for argument 'log'") 2703 rm(log) 2704 2705 2706 if (!is.Numeric(x1)) stop("bad input for 'x1'") 2707 if (!is.Numeric(x2)) stop("bad input for 'x2'") 2708 if (!is.Numeric(apar, positive = TRUE)) stop("bad input for 'apar'") 2709 2710 L <- max(length(x1), length(x2), length(apar)) 2711 if (length(apar ) != L) apar <- rep_len(apar, L) 2712 if (length(x1 ) != L) x1 <- rep_len(x1, L) 2713 if (length(x2 ) != L) x2 <- rep_len(x2, L) 2714 2715 if (log.arg) { 2716 denom <- apar-1 + (apar^x1 - 1) * (apar^x2 - 1) 2717 denom <- abs(denom) 2718 log((apar - 1) * log(apar)) + (x1+x2)*log(apar) - 2 * log(denom) 2719 } else { 2720 temp <- (apar - 1) + (apar^x1 - 1) * (apar^x2 - 1) 2721 index <- (abs(apar - 1) < .Machine$double.eps) 2722 ans <- x1 2723 if (any(!index)) 2724 ans[!index] <- (apar[!index] - 1) * log(apar[!index]) * 2725 (apar[!index])^(x1[!index] + 2726 x2[!index]) / (temp[!index])^2 2727 ans[x1 <= 0 | x2 <= 0 | x1 >= 1 | x2 >= 1] <- 0 2728 ans[index] <- 1 2729 ans 2730 } 2731} 2732 2733 2734 2735 2736bifrankcop.control <- function(save.weights = TRUE, ...) { 2737 list(save.weights = save.weights) 2738} 2739 2740 2741 2742 2743 2744 2745 bifrankcop <- function(lapar = "loglink", iapar = 2, nsimEIM = 250) { 2746 2747 lapar <- as.list(substitute(lapar)) 2748 eapar <- link2list(lapar) 2749 lapar <- attr(eapar, "function.name") 2750 2751 2752 if (!is.Numeric(iapar, positive = TRUE)) 2753 stop("argument 'iapar' must be positive") 2754 2755 2756 if (length(nsimEIM) && 2757 (!is.Numeric(nsimEIM, length.arg = 1, 2758 integer.valued = TRUE) || 2759 nsimEIM <= 50)) 2760 stop("argument 'nsimEIM' should be an integer greater than 50") 2761 2762 2763 new("vglmff", 2764 blurb = c("Frank's bivariate copula\n", 2765 "Links: ", 2766 namesof("apar", lapar, earg = eapar )), 2767 initialize = eval(substitute(expression({ 2768 2769 if (any(y <= 0) || any(y >= 1)) 2770 stop("the response must have values between 0 and 1") 2771 2772 temp5 <- 2773 w.y.check(w = w, y = y, 2774 Is.positive.y = TRUE, 2775 ncol.w.max = 1, 2776 ncol.y.max = 2, 2777 ncol.y.min = 2, 2778 out.wy = TRUE, 2779 colsyperw = 2, 2780 maximize = TRUE) 2781 w <- temp5$w 2782 y <- temp5$y 2783 2784 2785 predictors.names <- 2786 c(namesof("apar", .lapar , earg = .eapar, short = TRUE)) 2787 2788 extra$colnames.y <- colnames(y) 2789 2790 if (!length(etastart)) { 2791 apar.init <- rep_len(.iapar , n) 2792 etastart <- cbind(theta2eta(apar.init, .lapar , earg = .eapar )) 2793 } 2794 }), list( .lapar = lapar, .eapar = eapar, .iapar = iapar))), 2795 linkinv = eval(substitute(function(eta, extra = NULL) { 2796 NOS <- NCOL(eta) / c(M1 = 1) 2797 Q1 <- 2 2798 fv.mat <- matrix(0.5, NROW(eta), NOS * Q1) 2799 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 2800 }, list( .lapar = lapar, .eapar = eapar ))), 2801 last = eval(substitute(expression({ 2802 misc$link <- c("apar" = .lapar ) 2803 2804 misc$earg <- list("apar" = .eapar ) 2805 2806 misc$expected <- TRUE 2807 misc$nsimEIM <- .nsimEIM 2808 misc$pooled.weight <- pooled.weight 2809 misc$multipleResponses <- FALSE 2810 }), list( .lapar = lapar, .eapar = eapar, .nsimEIM = nsimEIM ))), 2811 loglikelihood = eval(substitute( 2812 function(mu, y, w, residuals = FALSE, eta, 2813 extra = NULL, 2814 summation = TRUE) { 2815 apar <- eta2theta(eta, .lapar , earg = .eapar ) 2816 if (residuals) { 2817 stop("loglikelihood residuals not implemented yet") 2818 } else { 2819 ll.elts <- c(w) * dbifrankcop(x1 = y[, 1], x2 = y[, 2], 2820 apar = apar, log = TRUE) 2821 if (summation) { 2822 sum(ll.elts) 2823 } else { 2824 ll.elts 2825 } 2826 } 2827 }, list( .lapar = lapar, .eapar = eapar ))), 2828 vfamily = c("bifrankcop"), 2829 validparams = eval(substitute(function(eta, y, extra = NULL) { 2830 apar <- eta2theta(eta, .lapar , earg = .eapar ) 2831 okay1 <- all(is.finite(apar)) && all(0 < apar) 2832 okay1 2833 }, list( .lapar = lapar, .eapar = eapar ))), 2834 2835 2836 2837 simslot = eval(substitute( 2838 function(object, nsim) { 2839 2840 pwts <- if (length(pwts <- object@prior.weights) > 0) 2841 pwts else weights(object, type = "prior") 2842 if (any(pwts != 1)) 2843 warning("ignoring prior weights") 2844 eta <- predict(object) 2845 apar <- eta2theta(eta, .lapar , earg = .eapar ) 2846 rbifrankcop(nsim * length(apar), apar = c(apar)) 2847 }, list( .lapar = lapar, .eapar = eapar ))), 2848 2849 2850 2851 2852 deriv = eval(substitute(expression({ 2853 apar <- eta2theta(eta, .lapar , earg = .eapar ) 2854 dapar.deta <- dtheta.deta(apar, .lapar , earg = .eapar ) 2855 2856 de3 <- deriv3(~ (log((apar - 1) * log(apar)) + (y1+y2)*log(apar) - 2857 2 * log(apar-1 + (apar^y1 - 1) * (apar^y2 - 1))), 2858 name = "apar", hessian = TRUE) 2859 2860 denom <- apar-1 + (apar^y[, 1] - 1) * (apar^y[, 2] - 1) 2861 tmp700 <- 2*apar^(y[, 1]+y[, 2]) - apar^y[, 1] - apar^y[, 2] 2862 numerator <- 1 + y[, 1] * apar^(y[, 1] - 1) * (apar^y[, 2] - 1) + 2863 y[, 2] * apar^(y[, 2] - 1) * (apar^y[, 1] - 1) 2864 Dl.dapar <- 1/(apar - 1) + 1/(apar*log(apar)) + 2865 (y[, 1]+y[, 2])/apar - 2 * numerator / denom 2866 c(w) * Dl.dapar * dapar.deta 2867 }), list( .lapar = lapar, 2868 .eapar = eapar, .nsimEIM = nsimEIM ))), 2869 weight = eval(substitute(expression({ 2870 if ( is.Numeric( .nsimEIM)) { 2871 2872 pooled.weight <- FALSE # For @last 2873 2874 2875 run.mean <- 0 2876 for (ii in 1:( .nsimEIM )) { 2877 ysim <- rbifrankcop(n, apar = apar) 2878 y1 <- ysim[, 1]; y2 <- ysim[, 2]; 2879 eval.de3 <- eval(de3) 2880 d2l.dthetas2 <- attr(eval.de3, "hessian") 2881 rm(ysim) 2882 temp3 <- -d2l.dthetas2[, 1, 1] # M = 1 2883 run.mean <- ((ii - 1) * run.mean + temp3) / ii 2884 } 2885 wz <- if (intercept.only) 2886 matrix(mean(run.mean), n, dimm(M)) else run.mean 2887 2888 wz <- wz * dapar.deta^2 2889 c(w) * wz 2890 } else { 2891 nump <- apar^(y[, 1]+y[, 2]-2) * (2 * y[, 1] * y[, 2] + 2892 y[, 1]*(y[, 1] - 1) + y[, 2]*(y[, 2] - 1)) - 2893 y[, 1]*(y[, 1] - 1) * apar^(y[, 1]-2) - 2894 y[, 2]*(y[, 2] - 1) * apar^(y[, 2]-2) 2895 D2l.dapar2 <- 1/(apar - 1)^2 + (1+log(apar))/(apar*log(apar))^2 + 2896 (y[, 1]+y[, 2])/apar^2 + 2 * 2897 (nump / denom - (numerator/denom)^2) 2898 d2apar.deta2 <- d2theta.deta2(apar, .lapar , earg = .eapar ) 2899 wz <- c(w) * (dapar.deta^2 * D2l.dapar2 - Dl.dapar * d2apar.deta2) 2900 if (TRUE && intercept.only) { 2901 wz <- cbind(wz) 2902 sumw <- sum(w) 2903 for (iii in 1:ncol(wz)) 2904 wz[,iii] <- sum(wz[, iii]) / sumw 2905 pooled.weight <- TRUE 2906 wz <- c(w) * wz # Put back the weights 2907 } else { 2908 pooled.weight <- FALSE 2909 } 2910 wz 2911 } 2912 }), list( .lapar = lapar, 2913 .eapar = eapar, .nsimEIM = nsimEIM )))) 2914} 2915 2916 2917 2918 2919 2920 gammahyperbola <- 2921 function(ltheta = "loglink", itheta = NULL, expected = FALSE) { 2922 2923 ltheta <- as.list(substitute(ltheta)) 2924 etheta <- link2list(ltheta) 2925 ltheta <- attr(etheta, "function.name") 2926 2927 if (!is.logical(expected) || length(expected) != 1) 2928 stop("argument 'expected' must be a single logical") 2929 2930 2931 new("vglmff", 2932 blurb = c("Gamma hyperbola bivariate distribution\n", 2933 "Links: ", 2934 namesof("theta", ltheta, etheta)), 2935 initialize = eval(substitute(expression({ 2936 if (any(y[, 1] <= 0) || any(y[, 2] <= 1)) 2937 stop("the response has values that are out of range") 2938 2939 temp5 <- 2940 w.y.check(w = w, y = y, 2941 Is.positive.y = TRUE, 2942 ncol.w.max = 1, 2943 ncol.y.max = 2, 2944 ncol.y.min = 2, 2945 out.wy = TRUE, 2946 colsyperw = 2, 2947 maximize = TRUE) 2948 w <- temp5$w 2949 y <- temp5$y 2950 2951 extra$colnames.y <- colnames(y) 2952 2953 2954 predictors.names <- 2955 c(namesof("theta", .ltheta , .etheta , short = TRUE)) 2956 2957 if (!length(etastart)) { 2958 theta.init <- if (length( .itheta)) { 2959 rep_len( .itheta , n) 2960 } else { 2961 1 / (y[, 2] - 1 + 0.01) 2962 } 2963 etastart <- 2964 cbind(theta2eta(theta.init, .ltheta , .etheta )) 2965 } 2966 }), list( .ltheta = ltheta, .etheta = etheta, .itheta = itheta))), 2967 linkinv = eval(substitute(function(eta, extra = NULL) { 2968 NOS <- NCOL(eta) / c(M1 = 1) 2969 theta <- eta2theta(eta, .ltheta , .etheta ) 2970 fv.mat <- cbind(theta * exp(theta), 1 + 1 / theta) 2971 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 2972 }, list( .ltheta = ltheta, .etheta = etheta ))), 2973 last = eval(substitute(expression({ 2974 misc$link <- c("theta" = .ltheta ) 2975 2976 misc$earg <- list("theta" = .etheta ) 2977 2978 misc$expected <- .expected 2979 misc$multipleResponses <- FALSE 2980 }), list( .ltheta = ltheta, 2981 .etheta = etheta, .expected = expected ))), 2982 2983 loglikelihood = eval(substitute( 2984 function(mu, y, w, residuals = FALSE, eta, 2985 extra = NULL, 2986 summation = TRUE) { 2987 theta <- eta2theta(eta, .ltheta , .etheta ) 2988 if (residuals) { 2989 stop("loglikelihood residuals not implemented yet") 2990 } else { 2991 ll.elts <- c(w) * (-exp(-theta) * y[, 1] / theta - theta * y[, 2]) 2992 if (summation) { 2993 sum(ll.elts) 2994 } else { 2995 ll.elts 2996 } 2997 } 2998 }, list( .ltheta = ltheta, .etheta = etheta ))), 2999 vfamily = c("gammahyperbola"), 3000 validparams = eval(substitute(function(eta, y, extra = NULL) { 3001 theta <- eta2theta(eta, .ltheta , .etheta ) 3002 okay1 <- all(is.finite(theta)) && all(0 < theta) 3003 okay1 3004 }, list( .ltheta = ltheta, .etheta = etheta ))), 3005 deriv = eval(substitute(expression({ 3006 theta <- eta2theta(eta, .ltheta , .etheta ) 3007 Dl.dtheta <- exp(-theta) * y[, 1] * (1+theta) / theta^2 - y[, 2] 3008 DTHETA.deta <- dtheta.deta(theta, .ltheta , .etheta ) 3009 c(w) * Dl.dtheta * DTHETA.deta 3010 }), list( .ltheta = ltheta, .etheta = etheta ))), 3011 weight = eval(substitute(expression({ 3012 temp300 <- 2 + theta * (2 + theta) 3013 if ( .expected ) { 3014 D2l.dtheta2 <- temp300 / theta^2 3015 wz <- c(w) * DTHETA.deta^2 * D2l.dtheta2 3016 } else { 3017 D2l.dtheta2 <- temp300 * y[, 1] * exp(-theta) / theta^3 3018 D2theta.deta2 <- d2theta.deta2(theta, .ltheta ) 3019 wz <- c(w) * (DTHETA.deta^2 * D2l.dtheta2 - 3020 Dl.dtheta * D2theta.deta2) 3021 } 3022 wz 3023 }), list( .ltheta = ltheta, 3024 .etheta = etheta, .expected = expected )))) 3025} 3026 3027 3028 3029 bifgmexp <- 3030 function(lapar = "rhobitlink", 3031 iapar = NULL, tola0 = 0.01, 3032 imethod = 1) { 3033 lapar <- as.list(substitute(lapar)) 3034 earg <- link2list(lapar) 3035 lapar <- attr(earg, "function.name") 3036 3037 if (length(iapar) && 3038 (!is.Numeric(iapar, length.arg = 1) || 3039 abs(iapar) >= 1)) 3040 stop("argument 'iapar' must be a single number between -1 and 1") 3041 3042 if (!is.Numeric(tola0, length.arg = 1, positive = TRUE)) 3043 stop("argument 'tola0' must be a single positive number") 3044 3045 if (length(iapar) && abs(iapar) <= tola0) 3046 stop("argument 'iapar' must not be between -tola0 and tola0") 3047 if (!is.Numeric(imethod, length.arg = 1, 3048 integer.valued = TRUE, positive = TRUE) || 3049 imethod > 2.5) 3050 stop("argument 'imethod' must be 1 or 2") 3051 3052 3053 new("vglmff", 3054 blurb = c("Bivariate Farlie-Gumbel-Morgenstern ", 3055 "exponential distribution\n", # Morgenstern's 3056 "Links: ", 3057 namesof("apar", lapar, earg = earg )), 3058 initialize = eval(substitute(expression({ 3059 temp5 <- 3060 w.y.check(w = w, y = y, 3061 Is.nonnegative.y = TRUE, 3062 ncol.w.max = 1, 3063 ncol.y.max = 2, 3064 ncol.y.min = 2, 3065 out.wy = TRUE, 3066 colsyperw = 2, 3067 maximize = TRUE) 3068 w <- temp5$w 3069 y <- temp5$y 3070 3071 3072 3073 predictors.names <- 3074 c(namesof("apar", .lapar , earg = .earg , short = TRUE)) 3075 3076 extra$colnames.y <- colnames(y) 3077 3078 if (!length(etastart)) { 3079 ainit <- if (length(.iapar)) rep_len( .iapar , n) else { 3080 mean1 <- if ( .imethod == 1) median(y[, 1]) else mean(y[, 1]) 3081 mean2 <- if ( .imethod == 1) median(y[, 2]) else mean(y[, 2]) 3082 Finit <- 0.01 + mean(y[, 1] <= mean1 & y[, 2] <= mean2) 3083 ((Finit+expm1(-mean1)+exp(-mean2)) / exp(-mean1-mean2) - 1) / ( 3084 expm1(-mean1) * expm1(-mean2)) 3085 } 3086 etastart <- 3087 theta2eta(rep_len(ainit, n), .lapar , earg = .earg ) 3088 } 3089 }), list( .iapar = iapar, .lapar = lapar, .earg = earg, 3090 .imethod = imethod ))), 3091 linkinv = eval(substitute(function(eta, extra = NULL) { 3092 NOS <- NCOL(eta) / c(M1 = 1) 3093 Q1 <- 2 3094 fv.mat <- matrix(1, NROW(eta), NOS * Q1) 3095 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 3096 }, list( .lapar = lapar, .earg = earg ))), 3097 last = eval(substitute(expression({ 3098 misc$link <- c("apar" = .lapar ) 3099 3100 misc$earg <- list("apar" = .earg ) 3101 3102 misc$expected <- FALSE 3103 misc$pooled.weight <- pooled.weight 3104 misc$multipleResponses <- FALSE 3105 }), list( .lapar = lapar, .earg = earg ))), 3106 loglikelihood = eval(substitute( 3107 function(mu, y, w, residuals = FALSE, eta, 3108 extra = NULL, 3109 summation = TRUE) { 3110 alpha <- eta2theta(eta, .lapar , earg = .earg ) 3111 alpha[abs(alpha) < .tola0 ] <- .tola0 3112 if (residuals) { 3113 stop("loglikelihood residuals not implemented yet") 3114 } else { 3115 denom <- (1 + alpha - 2*alpha*(exp(-y[, 1]) + exp(-y[, 2])) + 3116 4*alpha*exp(-y[, 1] - y[, 2])) 3117 ll.elts <- c(w) * (-y[, 1] - y[, 2] + log(denom)) 3118 if (summation) { 3119 sum(ll.elts) 3120 } else { 3121 ll.elts 3122 } 3123 } 3124 }, list( .lapar = lapar, .earg = earg, .tola0 = tola0 ))), 3125 vfamily = c("bifgmexp"), # morgenstern 3126 validparams = eval(substitute(function(eta, y, extra = NULL) { 3127 alpha <- eta2theta(eta, .lapar , earg = .earg ) 3128 okay1 <- all(is.finite(alpha)) && all(abs(alpha) < 1) 3129 okay1 3130 }, list( .lapar = lapar, .earg = earg, .tola0 = tola0 ))), 3131 deriv = eval(substitute(expression({ 3132 alpha <- eta2theta(eta, .lapar , earg = .earg ) 3133 alpha[abs(alpha) < .tola0 ] <- .tola0 3134 numerator <- 1 - 2*(exp(-y[, 1]) + exp(-y[, 2])) + 3135 4*exp(-y[, 1] - y[, 2]) 3136 denom <- (1 + alpha - 2*alpha*(exp(-y[, 1]) + exp(-y[, 2])) + 3137 4 *alpha*exp(-y[, 1] - y[, 2])) 3138 dl.dalpha <- numerator / denom 3139 3140 dalpha.deta <- dtheta.deta(alpha, .lapar , earg = .earg ) 3141 3142 c(w) * cbind(dl.dalpha * dalpha.deta) 3143 }), list( .lapar = lapar, .earg = earg, .tola0 = tola0 ))), 3144 weight = eval(substitute(expression({ 3145 d2l.dalpha2 <- dl.dalpha^2 3146 d2alpha.deta2 <- d2theta.deta2(alpha, .lapar , earg = .earg ) 3147 wz <- c(w) * (dalpha.deta^2 * d2l.dalpha2 - 3148 d2alpha.deta2 * dl.dalpha) 3149 if (TRUE && 3150 intercept.only) { 3151 wz <- cbind(wz) 3152 sumw <- sum(w) 3153 for (iii in 1:ncol(wz)) 3154 wz[, iii] <- sum(wz[, iii]) / sumw 3155 pooled.weight <- TRUE 3156 wz <- c(w) * wz # Put back the weights 3157 } else { 3158 pooled.weight <- FALSE 3159 } 3160 wz 3161 }), list( .lapar = lapar, .earg = earg )))) 3162} 3163 3164 3165 3166 3167rbifgmcop <- function(n, apar) { 3168 use.n <- if ((length.n <- length(n)) > 1) length.n else 3169 if (!is.Numeric(n, integer.valued = TRUE, 3170 length.arg = 1, positive = TRUE)) 3171 stop("bad input for argument 'n'") else n 3172 3173 if (!is.Numeric(apar)) 3174 stop("bad input for argument 'apar'") 3175 if (any(abs(apar) > 1)) 3176 stop("argument 'apar' has values out of range") 3177 3178 y1 <- V1 <- runif(use.n) 3179 V2 <- runif(use.n) 3180 temp <- 2*y1 - 1 3181 A <- apar * temp - 1 3182 B <- sqrt(1 - 2 * apar * temp + (apar*temp)^2 + 4 * apar * V2 * temp) 3183 y2 <- 2 * V2 / (B - A) 3184 matrix(c(y1, y2), nrow = use.n, ncol = 2) 3185} 3186 3187 3188 3189dbifgmcop <- function(x1, x2, apar, log = FALSE) { 3190 if (!is.logical(log.arg <- log) || 3191 length(log) != 1) 3192 stop("bad input for argument 'log'") 3193 rm(log) 3194 3195 if (!is.Numeric(apar)) 3196 stop("bad input for 'apar'") 3197 if (any(abs(apar) > 1)) 3198 stop("'apar' values out of range") 3199 if ( !is.logical( log.arg ) || 3200 length( log.arg ) != 1 ) 3201 stop("bad input for argument 'log'") 3202 3203 L <- max(length(x1), length(x2), length(apar)) 3204 if (length(x1) != L) x1 <- rep_len(x1, L) 3205 if (length(x2) != L) x2 <- rep_len(x2, L) 3206 if (length(apar) != L) apar <- rep_len(apar, L) 3207 ans <- 0 * x1 3208 xnok <- (x1 <= 0) | (x1 >= 1) | (x2 <= 0) | (x2 >= 1) 3209 if ( log.arg ) { 3210 ans[!xnok] <- log1p(apar[!xnok] * (1-2*x1[!xnok]) * (1-2*x2[!xnok])) 3211 ans[xnok] <- log(0) 3212 } else { 3213 ans[!xnok] <- 1 + apar[!xnok] * (1-2*x1[!xnok]) * (1-2*x2[!xnok]) 3214 ans[xnok] <- 0 3215 if (any(ans < 0)) 3216 stop("negative values in the density (apar out of range)") 3217 } 3218 ans 3219} 3220 3221 3222pbifgmcop <- function(q1, q2, apar) { 3223 if (!is.Numeric(q1)) stop("bad input for 'q1'") 3224 if (!is.Numeric(q2)) stop("bad input for 'q2'") 3225 if (!is.Numeric(apar)) stop("bad input for 'apar'") 3226 if (any(abs(apar) > 1)) stop("'apar' values out of range") 3227 3228 L <- max(length(q1), length(q2), length(apar)) 3229 if (length(q1) != L) q1 <- rep_len(q1, L) 3230 if (length(q2) != L) q2 <- rep_len(q2, L) 3231 if (length(apar) != L) apar <- rep_len(apar, L) 3232 3233 x <- q1 3234 y <- q2 3235 index <- (x >= 1 & y < 1) | 3236 (y >= 1 & x < 1) | 3237 (x <= 0 | y <= 0) | 3238 (x >= 1 & y >= 1) 3239 ans <- as.numeric(index) 3240 if (any(!index)) { 3241 ans[!index] <- q1[!index] * q2[!index] * (1 + apar[!index] * 3242 (1-q1[!index])*(1-q2[!index])) 3243 } 3244 ans[x >= 1 & y<1] <- y[x >= 1 & y<1] # P(Y2 < q2) = q2 3245 ans[y >= 1 & x<1] <- x[y >= 1 & x<1] # P(Y1 < q1) = q1 3246 ans[x <= 0 | y <= 0] <- 0 3247 ans[x >= 1 & y >= 1] <- 1 3248 ans 3249} 3250 3251 3252 3253 3254 3255 3256 bifgmcop <- function(lapar = "rhobitlink", iapar = NULL, 3257 imethod = 1) { 3258 3259 lapar <- as.list(substitute(lapar)) 3260 earg <- link2list(lapar) 3261 lapar <- attr(earg, "function.name") 3262 3263 3264 if (!is.Numeric(imethod, length.arg = 1, 3265 integer.valued = TRUE, positive = TRUE) || 3266 imethod > 3.5) 3267 stop("argument 'imethod' must be 1 or 2 or 3") 3268 3269 if (length(iapar) && 3270 (abs(iapar) >= 1)) 3271 stop("'iapar' should be less than 1 in absolute value") 3272 3273 3274 new("vglmff", 3275 blurb = c("Farlie-Gumbel-Morgenstern copula \n", # distribution 3276 "Links: ", 3277 namesof("apar", lapar, earg = earg )), 3278 initialize = eval(substitute(expression({ 3279 if (any(y < 0) || any(y > 1)) 3280 stop("the response must have values in the unit square") 3281 3282 temp5 <- 3283 w.y.check(w = w, y = y, 3284 Is.nonnegative.y = TRUE, 3285 ncol.w.max = 1, 3286 ncol.y.max = 2, 3287 ncol.y.min = 2, 3288 out.wy = TRUE, 3289 colsyperw = 2, 3290 maximize = TRUE) 3291 w <- temp5$w 3292 y <- temp5$y 3293 3294 3295 predictors.names <- 3296 namesof("apar", .lapar , earg = .earg , short = TRUE) 3297 3298 extra$colnames.y <- colnames(y) 3299 3300 if (!length(etastart)) { 3301 ainit <- if (length( .iapar )) .iapar else { 3302 3303 3304 if ( .imethod == 1) { 3305 3 * cor(y[, 1], y[, 2], method = "spearman") 3306 } else if ( .imethod == 2) { 3307 9 * kendall.tau(y[, 1], y[, 2]) / 2 3308 } else { 3309 mean1 <- if ( .imethod == 1) weighted.mean(y[, 1], w) else 3310 median(y[, 1]) 3311 mean2 <- if ( .imethod == 1) weighted.mean(y[, 2], w) else 3312 median(y[, 2]) 3313 Finit <- weighted.mean(y[, 1] <= mean1 & y[, 2] <= mean2, w) 3314 (Finit / (mean1 * mean2) - 1) / ((1 - mean1) * (1 - mean2)) 3315 } 3316 } 3317 3318 ainit <- min(0.95, max(ainit, -0.95)) 3319 3320 etastart <- theta2eta(rep_len(ainit, n), .lapar , earg = .earg ) 3321 } 3322 }), list( .iapar = iapar, .lapar = lapar, .earg = earg, 3323 .imethod = imethod ))), 3324 linkinv = eval(substitute(function(eta, extra = NULL) { 3325 NOS <- NCOL(eta) / c(M1 = 1) 3326 Q1 <- 2 3327 fv.mat <- matrix(0.5, NROW(eta), NOS * Q1) 3328 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 3329 }, list( .lapar = lapar, .earg = earg ))), 3330 last = eval(substitute(expression({ 3331 misc$link <- c("apar" = .lapar ) 3332 misc$earg <- list("apar" = .earg ) 3333 3334 misc$expected <- FALSE 3335 misc$multipleResponses <- FALSE 3336 }), list( .lapar = lapar, .earg = earg))), 3337 loglikelihood = eval(substitute( 3338 function(mu, y, w, residuals = FALSE, eta, 3339 extra = NULL, 3340 summation = TRUE) { 3341 alpha <- eta2theta(eta, .lapar , earg = .earg ) 3342 if (residuals) { 3343 stop("loglikelihood residuals not implemented yet") 3344 } else { 3345 ll.elts <- c(w) * dbifgmcop(x1 = y[, 1], 3346 x2 = y[, 2], apar = alpha, log = TRUE) 3347 if (summation) { 3348 sum(ll.elts) 3349 } else { 3350 ll.elts 3351 } 3352 } 3353 }, list( .lapar = lapar, .earg = earg ))), 3354 vfamily = c("bifgmcop"), 3355 validparams = eval(substitute(function(eta, y, extra = NULL) { 3356 alpha <- eta2theta(eta, .lapar , earg = .earg ) 3357 okay1 <- all(is.finite(alpha)) && all(abs(alpha) < 1) 3358 okay1 3359 }, list( .lapar = lapar, .earg = earg ))), 3360 3361 3362 simslot = eval(substitute( 3363 function(object, nsim) { 3364 3365 pwts <- if (length(pwts <- object@prior.weights) > 0) 3366 pwts else weights(object, type = "prior") 3367 if (any(pwts != 1)) 3368 warning("ignoring prior weights") 3369 eta <- predict(object) 3370 alpha <- eta2theta(eta, .lapar , earg = .earg ) 3371 rbifgmcop(nsim * length(alpha), apar = c(alpha)) 3372 }, list( .lapar = lapar, .earg = earg ))), 3373 3374 3375 3376 deriv = eval(substitute(expression({ 3377 alpha <- eta2theta(eta, .lapar , earg = .earg ) 3378 3379 dalpha.deta <- dtheta.deta(alpha, .lapar , earg = .earg ) 3380 3381 numerator <- (1 - 2 * y[, 1]) * (1 - 2 * y[, 2]) 3382 denom <- 1 + alpha * numerator 3383 3384 mytolerance <- .Machine$double.eps 3385 bad <- (denom <= mytolerance) # Range violation 3386 if (any(bad)) { 3387 cat("There are some range violations in @deriv\n") 3388 flush.console() 3389 denom[bad] <- 2 * mytolerance 3390 } 3391 dl.dalpha <- numerator / denom 3392 c(w) * cbind(dl.dalpha * dalpha.deta) 3393 }), list( .lapar = lapar, .earg = earg))), 3394 3395 weight = eval(substitute(expression({ 3396 wz <- lerch(alpha^2, 2, 1.5) / 4 # Checked and correct 3397 wz <- wz * dalpha.deta^2 3398 c(w) * wz 3399 }), list( .lapar = lapar, .earg = earg)))) 3400} 3401 3402 3403 3404 3405 3406 bigumbelIexp <- 3407 function(lapar = "identitylink", iapar = NULL, imethod = 1) { 3408 3409 lapar <- as.list(substitute(lapar)) 3410 earg <- link2list(lapar) 3411 lapar <- attr(earg, "function.name") 3412 3413 3414 if (length(iapar) && 3415 !is.Numeric(iapar, length.arg = 1)) 3416 stop("'iapar' must be a single number") 3417 if (!is.Numeric(imethod, length.arg = 1, 3418 integer.valued = TRUE, positive = TRUE) || 3419 imethod > 2.5) 3420 stop("argument 'imethod' must be 1 or 2") 3421 3422 3423 new("vglmff", 3424 blurb = c("Gumbel's Type I bivariate exponential distribution\n", 3425 "Links: ", 3426 namesof("apar", lapar, earg = earg )), 3427 initialize = eval(substitute(expression({ 3428 3429 temp5 <- 3430 w.y.check(w = w, y = y, 3431 Is.nonnegative.y = TRUE, 3432 ncol.w.max = 1, 3433 ncol.y.max = 2, 3434 ncol.y.min = 2, 3435 out.wy = TRUE, 3436 colsyperw = 2, 3437 maximize = TRUE) 3438 w <- temp5$w 3439 y <- temp5$y 3440 3441 extra$colnames.y <- colnames(y) 3442 3443 3444 3445 predictors.names <- 3446 c(namesof("apar", .lapar , earg = .earg , short = TRUE)) 3447 3448 if (!length(etastart)) { 3449 ainit <- if (length( .iapar )) rep_len( .iapar, n) else { 3450 mean1 <- if ( .imethod == 1) median(y[, 1]) else mean(y[, 1]) 3451 mean2 <- if ( .imethod == 1) median(y[, 2]) else mean(y[, 2]) 3452 Finit <- 0.01 + mean(y[, 1] <= mean1 & y[, 2] <= mean2) 3453 (log(Finit+expm1(-mean1)+exp(-mean2))+mean1+mean2)/(mean1*mean2) 3454 } 3455 etastart <- 3456 theta2eta(rep_len(ainit, n), .lapar , earg = .earg ) 3457 } 3458 }), list( .iapar = iapar, .lapar = lapar, .earg = earg, 3459 .imethod = imethod ))), 3460 linkinv = eval(substitute(function(eta, extra = NULL) { 3461 NOS <- NCOL(eta) / c(M1 = 1) 3462 Q1 <- 2 3463 alpha <- eta2theta(eta, .lapar , earg = .earg ) 3464 fv.mat <- matrix(1, NROW(eta), NOS * Q1) 3465 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 3466 }, list( .lapar = lapar, .earg = earg ))), 3467 last = eval(substitute(expression({ 3468 misc$link <- c("apar" = .lapar ) 3469 3470 misc$earg <- list("apar" = .earg ) 3471 3472 misc$expected <- FALSE 3473 misc$pooled.weight <- pooled.weight 3474 misc$multipleResponses <- FALSE 3475 }), list( .lapar = lapar, .earg = earg ))), 3476 loglikelihood = eval(substitute( 3477 function(mu, y, w, residuals = FALSE, eta, 3478 extra = NULL, 3479 summation = TRUE) { 3480 alpha <- eta2theta(eta, .lapar , earg = .earg ) 3481 if (residuals) { 3482 stop("loglikelihood residuals not implemented yet") 3483 } else { 3484 denom <- (alpha*y[, 1] - 1) * (alpha*y[, 2] - 1) + alpha 3485 mytolerance <- .Machine$double.xmin 3486 bad <- (denom <= mytolerance) # Range violation 3487 if (any(bad)) { 3488 cat("There are some range violations in @deriv\n") 3489 flush.console() 3490 } 3491 3492 3493 3494 3495 if (summation) { 3496 sum(bad) * (-1.0e10) + 3497 sum(w[!bad] * (-y[!bad, 1] - y[!bad, 2] + 3498 alpha[!bad] * y[!bad, 1] * y[!bad, 2] + log(denom[!bad]))) 3499 } else { 3500 stop("argument 'summation = FALSE' does not work yet") 3501 } 3502 } 3503 }, list( .lapar = lapar, .earg = earg ))), 3504 vfamily = c("bigumbelIexp"), 3505 validparams = eval(substitute(function(eta, y, extra = NULL) { 3506 alpha <- eta2theta(eta, .lapar , earg = .earg ) 3507 okay1 <- all(is.finite(alpha)) 3508 okay1 3509 }, list( .lapar = lapar, .earg = earg ))), 3510 deriv = eval(substitute(expression({ 3511 alpha <- eta2theta(eta, .lapar , earg = .earg ) 3512 numerator <- (alpha * y[, 1] - 1) * y[, 2] + 3513 (alpha * y[, 2] - 1) * y[, 1] + 1 3514 denom <- (alpha * y[, 1] - 1) * (alpha * y[, 2] - 1) + alpha 3515 denom <- abs(denom) 3516 3517 dl.dalpha <- numerator / denom + y[, 1] * y[, 2] 3518 3519 dalpha.deta <- dtheta.deta(alpha, .lapar , earg = .earg ) 3520 3521 c(w) * cbind(dl.dalpha * dalpha.deta) 3522 }), list( .lapar = lapar, .earg = earg ))), 3523 weight = eval(substitute(expression({ 3524 d2l.dalpha2 <- (numerator/denom)^2 - 2*y[, 1]*y[, 2] / denom 3525 d2alpha.deta2 <- d2theta.deta2(alpha, .lapar , earg = .earg ) 3526 wz <- c(w) * (dalpha.deta^2 * d2l.dalpha2 - 3527 d2alpha.deta2 * dl.dalpha) 3528 if (TRUE && 3529 intercept.only) { 3530 wz <- cbind(wz) 3531 sumw <- sum(w) 3532 for (iii in 1:ncol(wz)) 3533 wz[, iii] <- sum(wz[, iii]) / sumw 3534 pooled.weight <- TRUE 3535 wz <- c(w) * wz # Put back the weights 3536 } else { 3537 pooled.weight <- FALSE 3538 } 3539 wz 3540 }), list( .lapar = lapar, .earg = earg )))) 3541} 3542 3543 3544 3545 3546 3547 3548 3549pbiplackcop <- function(q1, q2, oratio) { 3550 if (!is.Numeric(q1)) stop("bad input for 'q1'") 3551 if (!is.Numeric(q2)) stop("bad input for 'q2'") 3552 if (!is.Numeric(oratio, positive = TRUE)) 3553 stop("bad input for 'oratio'") 3554 3555 L <- max(length(q1), length(q2), length(oratio)) 3556 if (length(q1) != L) q1 <- rep_len(q1, L) 3557 if (length(q2) != L) q2 <- rep_len(q2, L) 3558 if (length(oratio) != L) oratio <- rep_len(oratio, L) 3559 3560 x <- q1; y <- q2 3561 index <- (x >= 1 & y < 1) | (y >= 1 & x < 1) | 3562 (x <= 0 | y <= 0) | (x >= 1 & y >= 1) | 3563 (abs(oratio - 1) < 1.0e-6) # .Machine$double.eps 3564 ans <- as.numeric(index) 3565 if (any(!index)) { 3566 temp1 <- 1 + (oratio[!index] - 1) * (q1[!index] + q2[!index]) 3567 temp2 <- temp1 - sqrt(temp1^2 - 4 * oratio[!index] * 3568 (oratio[!index] - 1) * q1[!index] * q2[!index]) 3569 ans[!index] <- 0.5 * temp2 / (oratio[!index] - 1) 3570 } 3571 3572 ind2 <- (abs(oratio - 1) < 1.0e-6) # .Machine$double.eps 3573 ans[ind2] <- x[ind2] * y[ind2] 3574 ans[x >= 1 & y<1] <- y[x >= 1 & y<1] # P(Y2 < q2) = q2 3575 ans[y >= 1 & x<1] <- x[y >= 1 & x<1] # P(Y1 < q1) = q1 3576 ans[x <= 0 | y <= 0] <- 0 3577 ans[x >= 1 & y >= 1] <- 1 3578 ans 3579} 3580 3581 3582 3583rbiplackcop <- function(n, oratio) { 3584 use.n <- if ((length.n <- length(n)) > 1) length.n else 3585 if (!is.Numeric(n, integer.valued = TRUE, 3586 length.arg = 1, positive = TRUE)) 3587 stop("bad input for argument 'n'") else n 3588 3589 3590 y1 <- U <- runif(use.n) 3591 V <- runif(use.n) 3592 Z <- V * (1-V) 3593 y2 <- (2*Z*(y1*oratio^2 + 1 - y1) + oratio * (1 - 2 * Z) - 3594 (1 - 2 * V) * 3595 sqrt(oratio * (oratio + 4*Z*y1*(1-y1)*(1-oratio)^2))) / (oratio + 3596 Z*(1-oratio)^2) 3597 matrix(c(y1, 0.5 * y2), nrow = use.n, ncol = 2) 3598} 3599 3600 3601 3602dbiplackcop <- function(x1, x2, oratio, log = FALSE) { 3603 if (!is.logical(log.arg <- log) || length(log) != 1) 3604 stop("bad input for argument 'log'") 3605 rm(log) 3606 3607 3608 ans <- log(oratio) + log1p((oratio - 1) * 3609 (x1+x2 - 2*x1*x2)) - 1.5 * 3610 log((1 + (x1+x2)*(oratio - 1))^2 - 3611 4 * oratio * (oratio - 1)*x1*x2) 3612 ans[ # !is.na(x1) & !is.na(x2) & !is.na(oratio) & 3613 ((x1 < 0) | (x1 > 1) | (x2 < 0) | (x2 > 1))] <- log(0) 3614 3615 3616 if (log.arg) ans else exp(ans) 3617} 3618 3619 3620 3621biplackettcop.control <- function(save.weights = TRUE, ...) { 3622 list(save.weights = save.weights) 3623} 3624 3625 3626 3627 biplackettcop <- function(link = "loglink", ioratio = NULL, 3628 imethod = 1, nsimEIM = 200) { 3629 3630 link <- as.list(substitute(link)) 3631 earg <- link2list(link) 3632 link <- attr(earg, "function.name") 3633 3634 3635 if (length(ioratio) && (!is.Numeric(ioratio, positive = TRUE))) 3636 stop("'ioratio' must be positive") 3637 3638 if (!is.Numeric(imethod, length.arg = 1, 3639 integer.valued = TRUE, positive = TRUE) || 3640 imethod > 2) 3641 stop("argument 'imethod' must be 1 or 2") 3642 3643 3644 new("vglmff", 3645 blurb = c("Plackett distribution (bivariate copula)\n", 3646 "Links: ", 3647 namesof("oratio", link, earg = earg )), 3648 initialize = eval(substitute(expression({ 3649 if (any(y < 0) || any(y > 1)) 3650 stop("the response must have values in the unit square") 3651 3652 temp5 <- 3653 w.y.check(w = w, y = y, 3654 Is.nonnegative.y = TRUE, 3655 ncol.w.max = 1, 3656 ncol.y.max = 2, 3657 ncol.y.min = 2, 3658 out.wy = TRUE, 3659 colsyperw = 2, 3660 maximize = TRUE) 3661 w <- temp5$w 3662 y <- temp5$y 3663 3664 3665 predictors.names <- 3666 namesof("oratio", .link , earg = .earg, short = TRUE) 3667 3668 extra$colnames.y <- colnames(y) 3669 3670 if (!length(etastart)) { 3671 orinit <- if (length( .ioratio )) .ioratio else { 3672 if ( .imethod == 2) { 3673 scorp <- cor(y)[1, 2] 3674 if (abs(scorp) <= 0.1) 1 else 3675 if (abs(scorp) <= 0.3) 3^sign(scorp) else 3676 if (abs(scorp) <= 0.6) 5^sign(scorp) else 3677 if (abs(scorp) <= 0.8) 20^sign(scorp) else 40^sign(scorp) 3678 } else { 3679 y10 <- weighted.mean(y[, 1], w) 3680 y20 <- weighted.mean(y[, 2], w) 3681 (0.5 + sum(w[(y[, 1] < y10) & (y[, 2] < y20)])) * 3682 (0.5 + sum(w[(y[, 1] >= y10) & (y[, 2] >= y20)])) / ( 3683 ((0.5 + sum(w[(y[, 1] < y10) & (y[, 2] >= y20)])) * 3684 (0.5 + sum(w[(y[, 1] >= y10) & (y[, 2] < y20)])))) 3685 } 3686 } 3687 etastart <- theta2eta(rep_len(orinit, n), .link , earg = .earg ) 3688 } 3689 }), list( .ioratio = ioratio, .link = link, .earg = earg, 3690 .imethod = imethod ))), 3691 linkinv = eval(substitute(function(eta, extra = NULL) { 3692 NOS <- NCOL(eta) / c(M1 = 1) 3693 Q1 <- 2 3694 fv.mat <- matrix(0.5, NROW(eta), NOS * Q1) 3695 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 3696 }, list( .link = link, .earg = earg ))), 3697 last = eval(substitute(expression({ 3698 misc$link <- c(oratio = .link) 3699 3700 misc$earg <- list(oratio = .earg) 3701 3702 misc$expected <- FALSE 3703 misc$nsimEIM <- .nsimEIM 3704 misc$multipleResponses <- FALSE 3705 }), list( .link = link, .earg = earg, 3706 .nsimEIM = nsimEIM ))), 3707 loglikelihood = eval(substitute( 3708 function(mu, y, w, residuals = FALSE, eta, 3709 extra = NULL, 3710 summation = TRUE) { 3711 oratio <- eta2theta(eta, .link , earg = .earg ) 3712 if (residuals) { 3713 stop("loglikelihood residuals not implemented yet") 3714 } else { 3715 ll.elts <- c(w) * dbiplackcop(x1 = y[, 1], x2 = y[, 2], 3716 oratio = oratio, log = TRUE) 3717 if (summation) { 3718 sum(ll.elts) 3719 } else { 3720 ll.elts 3721 } 3722 } 3723 }, list( .link = link, .earg = earg ))), 3724 vfamily = c("biplackettcop"), 3725 validparams = eval(substitute(function(eta, y, extra = NULL) { 3726 oratio <- eta2theta(eta, .link , earg = .earg ) 3727 okay1 <- all(is.finite(oratio)) && all(0 < oratio) 3728 okay1 3729 }, list( .link = link, .earg = earg ))), 3730 3731 3732 simslot = eval(substitute( 3733 function(object, nsim) { 3734 3735 pwts <- if (length(pwts <- object@prior.weights) > 0) 3736 pwts else weights(object, type = "prior") 3737 if (any(pwts != 1)) 3738 warning("ignoring prior weights") 3739 eta <- predict(object) 3740 oratio <- eta2theta(eta, .link , earg = .earg ) 3741 rbiplackcop(nsim * length(oratio), oratio = c(oratio)) 3742 }, list( .link = link, .earg = earg ))), 3743 3744 3745 3746 deriv = eval(substitute(expression({ 3747 oratio <- eta2theta(eta, .link , earg = .earg ) 3748 doratio.deta <- dtheta.deta(oratio, .link , earg = .earg ) 3749 y1 <- y[, 1] 3750 y2 <- y[, 2] 3751 de3 <- deriv3(~ (log(oratio) + log(1+(oratio - 1) * 3752 (y1+y2-2*y1*y2)) - 1.5 * 3753 log((1 + (y1+y2)*(oratio - 1))^2 - 3754 4 * oratio * (oratio - 1)*y1*y2)), 3755 name = "oratio", hessian = FALSE) 3756 eval.de3 <- eval(de3) 3757 3758 dl.doratio <- attr(eval.de3, "gradient") 3759 3760 c(w) * dl.doratio * doratio.deta 3761 }), list( .link = link, .earg = earg ))), 3762 weight = eval(substitute(expression({ 3763 sd3 <- deriv3(~ (log(oratio) + log(1+(oratio - 1) * 3764 (y1sim+y2sim-2*y1sim*y2sim)) - 1.5 * 3765 log((1 + (y1sim+y2sim)*(oratio - 1))^2 - 3766 4 * oratio * (oratio - 1)*y1sim*y2sim)), 3767 name = "oratio", hessian = FALSE) 3768 run.var <- 0 3769 for (ii in 1:( .nsimEIM )) { 3770 ysim <- rbiplackcop(n, oratio = oratio) 3771 y1sim <- ysim[, 1] 3772 y2sim <- ysim[, 1] 3773 eval.sd3 <- eval(sd3) 3774 dl.doratio <- attr(eval.sd3, "gradient") 3775 rm(ysim, y1sim, y2sim) 3776 temp3 <- dl.doratio 3777 run.var <- ((ii - 1) * run.var + temp3^2) / ii 3778 } 3779 wz <- if (intercept.only) 3780 matrix(colMeans(cbind(run.var)), 3781 n, dimm(M), byrow = TRUE) else cbind(run.var) 3782 3783 wz <- wz * doratio.deta^2 3784 c(w) * wz 3785 }), list( .link = link, .earg = earg, .nsimEIM = nsimEIM )))) 3786} 3787 3788 3789 3790 3791 3792dbiamhcop <- function(x1, x2, apar, log = FALSE) { 3793 if (!is.logical(log.arg <- log) || length(log) != 1) 3794 stop("bad input for argument 'log'") 3795 rm(log) 3796 3797 3798 3799 L <- max(length(x1), length(x2), length(apar)) 3800 if (length(apar) != L) apar <- rep_len(apar, L) 3801 if (length(x1) != L) x1 <- rep_len(x1, L) 3802 if (length(x2) != L) x2 <- rep_len(x2, L) 3803 temp <- 1 - apar*(1-x1)*(1-x2) 3804 3805 if (log.arg) { 3806 ans <- log1p(-apar+2*apar*x1*x2/temp) - 2*log(temp) 3807 ans[(x1 <= 0) | (x1 >= 1) | (x2 <= 0) | (x2 >= 1)] <- log(0) 3808 } else { 3809 ans <- (1-apar+2*apar*x1*x2/temp) / (temp^2) 3810 ans[(x1 <= 0) | (x1 >= 1) | (x2 <= 0) | (x2 >= 1)] <- 0 3811 } 3812 ans[abs(apar) > 1] <- NA 3813 ans 3814} 3815 3816 3817pbiamhcop <- function(q1, q2, apar) { 3818 if (!is.Numeric(q1)) stop("bad input for 'q1'") 3819 if (!is.Numeric(q2)) stop("bad input for 'q2'") 3820 if (!is.Numeric(apar)) stop("bad input for 'apar'") 3821 3822 L <- max(length(q1), length(q2), length(apar)) 3823 if (length(q1) != L) q1 <- rep_len(q1, L) 3824 if (length(q2) != L) q2 <- rep_len(q2, L) 3825 if (length(apar) != L) apar <- rep_len(apar, L) 3826 3827 x <- q1 3828 y <- q2 3829 index <- (x >= 1 & y < 1) | (y >= 1 & x < 1) | 3830 (x <= 0 | y<= 0) | (x >= 1 & y >= 1) 3831 ans <- as.numeric(index) 3832 if (any(!index)) { 3833 ans[!index] <- (q1[!index] * q2[!index]) / (1 - 3834 apar[!index] * (1-q1[!index]) * (1-q2[!index])) 3835 } 3836 ans[x >= 1 & y < 1] <- y[x >= 1 & y < 1] # P(Y2 < q2) = q2 3837 ans[y >= 1 & x < 1] <- x[y >= 1 & x < 1] # P(Y1 < q1) = q1 3838 ans[x <= 0 | y <= 0] <- 0 3839 ans[x >= 1 & y >= 1] <- 1 3840 ans[abs(apar) > 1] <- NA 3841 ans 3842} 3843 3844 3845rbiamhcop <- function(n, apar) { 3846 use.n <- if ((length.n <- length(n)) > 1) length.n else 3847 if (!is.Numeric(n, integer.valued = TRUE, 3848 length.arg = 1, positive = TRUE)) 3849 stop("bad input for argument 'n'") else n 3850 3851 3852 3853 3854 3855 3856 if (any(abs(apar) > 1)) 3857 stop("'apar' values out of range") 3858 3859 U1 <- V1 <- runif(use.n) 3860 V2 <- runif(use.n) 3861 b <- 1-V1 3862 A <- -apar*(2*b*V2+1)+2*apar^2*b^2*V2+1 3863 B <- apar^2*(4*b^2*V2-4*b*V2+1)+apar*(4*V2-4*b*V2-2)+1 3864 U2 <- (2*V2*(apar*b - 1)^2)/(A+sqrt(B)) 3865 matrix(c(U1, U2), nrow = use.n, ncol = 2) 3866} 3867 3868 3869biamhcop.control <- function(save.weights = TRUE, ...) { 3870 list(save.weights = save.weights) 3871} 3872 3873 3874 biamhcop <- function(lapar = "rhobitlink", iapar = NULL, 3875 imethod = 1, nsimEIM = 250) { 3876 lapar <- as.list(substitute(lapar)) 3877 eapar <- link2list(lapar) 3878 lapar <- attr(eapar, "function.name") 3879 3880 3881 3882 if (length(iapar) && (abs(iapar) > 1)) 3883 stop("'iapar' should be less than or equal to 1 in absolute value") 3884 if (!is.Numeric(imethod, length.arg = 1, 3885 integer.valued = TRUE, positive = TRUE) || 3886 imethod > 2) 3887 stop("imethod must be 1 or 2") 3888 3889 if (length(nsimEIM) && 3890 (!is.Numeric(nsimEIM, length.arg = 1, 3891 integer.valued = TRUE) || 3892 nsimEIM <= 50)) 3893 stop("'nsimEIM' should be an integer greater than 50") 3894 3895 3896 new("vglmff", 3897 blurb = c("Ali-Mikhail-Haq distribution\n", 3898 "Links: ", 3899 namesof("apar", lapar, earg = eapar )), 3900 initialize = eval(substitute(expression({ 3901 if (any(y < 0) || any(y > 1)) 3902 stop("the response must have values in the unit square") 3903 3904 temp5 <- 3905 w.y.check(w = w, y = y, 3906 Is.nonnegative.y = TRUE, 3907 ncol.w.max = 1, 3908 ncol.y.max = 2, 3909 ncol.y.min = 2, 3910 out.wy = TRUE, 3911 colsyperw = 2, 3912 maximize = TRUE) 3913 w <- temp5$w 3914 y <- temp5$y 3915 3916 3917 predictors.names <- 3918 c(namesof("apar", .lapar, earg = .eapar, short = TRUE)) 3919 3920 extra$colnames.y <- colnames(y) 3921 3922 if (!length(etastart)) { 3923 ainit <- if (length( .iapar )) .iapar else { 3924 mean1 <- if ( .imethod == 1) weighted.mean(y[, 1], w) else 3925 median(y[, 1]) 3926 mean2 <- if ( .imethod == 1) weighted.mean(y[, 2], w) else 3927 median(y[, 2]) 3928 Finit <- weighted.mean(y[, 1] <= mean1 & y[, 2] <= mean2, w) 3929 (1 - (mean1 * mean2 / Finit)) / ((1-mean1) * (1-mean2)) 3930 } 3931 ainit <- min(0.95, max(ainit, -0.95)) 3932 etastart <- theta2eta(rep_len(ainit, n), .lapar , earg = .eapar ) 3933 } 3934 }), list( .lapar = lapar, .eapar = eapar, .iapar = iapar, 3935 .imethod = imethod))), 3936 linkinv = eval(substitute(function(eta, extra = NULL) { 3937 NOS <- NCOL(eta) / c(M1 = 1) 3938 Q1 <- 2 3939 fv.mat <- matrix(0.5, NROW(eta), NOS * Q1) 3940 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 3941 }, list( .lapar = lapar, .eapar = eapar ))), 3942 last = eval(substitute(expression({ 3943 misc$link <- c("apar" = .lapar ) 3944 3945 misc$earg <- list("apar" = .eapar ) 3946 3947 misc$expected <- TRUE 3948 misc$nsimEIM <- .nsimEIM 3949 misc$multipleResponses <- FALSE 3950 }), list( .lapar = lapar, 3951 .eapar = eapar, .nsimEIM = nsimEIM ))), 3952 loglikelihood = eval(substitute( 3953 function(mu, y, w, residuals = FALSE, eta, 3954 extra = NULL, 3955 summation = TRUE) { 3956 apar <- eta2theta(eta, .lapar, earg = .eapar ) 3957 if (residuals) { 3958 stop("loglikelihood residuals not implemented yet") 3959 } else { 3960 ll.elts <- c(w) * dbiamhcop(x1 = y[, 1], x2 = y[, 2], 3961 apar = apar, log = TRUE) 3962 if (summation) { 3963 sum(ll.elts) 3964 } else { 3965 ll.elts 3966 } 3967 } 3968 }, list( .lapar = lapar, .eapar = eapar ))), 3969 vfamily = c("biamhcop"), 3970 validparams = eval(substitute(function(eta, y, extra = NULL) { 3971 apar <- eta2theta(eta, .lapar, earg = .eapar ) 3972 okay1 <- all(is.finite(apar)) && all(abs(apar) < 1) 3973 okay1 3974 }, list( .lapar = lapar, .eapar = eapar ))), 3975 3976 3977 3978 simslot = eval(substitute( 3979 function(object, nsim) { 3980 3981 pwts <- if (length(pwts <- object@prior.weights) > 0) 3982 pwts else weights(object, type = "prior") 3983 if (any(pwts != 1)) 3984 warning("ignoring prior weights") 3985 eta <- predict(object) 3986 apar <- eta2theta(eta, .lapar , earg = .eapar ) 3987 rbiamhcop(nsim * length(apar), apar = c(apar)) 3988 }, list( .lapar = lapar, .eapar = eapar ))), 3989 3990 3991 3992 deriv = eval(substitute(expression({ 3993 apar <- eta2theta(eta, .lapar, earg = .eapar ) 3994 3995 dapar.deta <- dtheta.deta(apar, .lapar, earg = .eapar ) 3996 3997 y1 <- y[, 1] 3998 y2 <- y[, 2] 3999 de3 <- deriv3(~ (log(1 - apar+ 4000 (2 * apar*y1*y2/(1-apar*(1-y1)*(1-y2)))) - 4001 2 * log(1 - apar*(1-y1)*(1-y2))) , 4002 name = "apar", hessian = FALSE) 4003 eval.de3 <- eval(de3) 4004 4005 dl.dapar <- attr(eval.de3, "gradient") 4006 4007 c(w) * dl.dapar * dapar.deta 4008 }), list( .lapar = lapar, .eapar = eapar ))), 4009 weight = eval(substitute(expression({ 4010 sd3 <- deriv3(~ (log(1 - apar + 4011 (2 * apar * y1sim * y2sim / (1 - apar * 4012 (1 - y1sim) * (1-y2sim)))) - 4013 2 * log(1-apar*(1-y1sim)*(1-y2sim))), 4014 name = "apar", hessian = FALSE) 4015 run.var <- 0 4016 for (ii in 1:( .nsimEIM )) { 4017 ysim <- rbiamhcop(n, apar = apar) 4018 y1sim <- ysim[, 1] 4019 y2sim <- ysim[, 1] 4020 eval.sd3 <- eval(sd3) 4021 dl.apar <- attr(eval.sd3, "gradient") 4022 rm(ysim, y1sim, y2sim) 4023 temp3 <- dl.dapar 4024 run.var <- ((ii - 1) * run.var + temp3^2) / ii 4025 } 4026 4027 wz <- if (intercept.only) 4028 matrix(colMeans(cbind(run.var)), 4029 n, dimm(M), byrow = TRUE) else cbind(run.var) 4030 4031 wz <- wz * dapar.deta^2 4032 4033 c(w) * wz 4034 }), list( .lapar = lapar, 4035 .eapar = eapar, .nsimEIM = nsimEIM )))) 4036} 4037 4038 4039 4040 4041 4042 4043 4044 4045 4046 4047 4048 4049 4050 4051 4052dbinorm <- function(x1, x2, mean1 = 0, mean2 = 0, 4053 var1 = 1, var2 = 1, cov12 = 0, 4054 log = FALSE) { 4055 if (!is.logical(log.arg <- log) || length(log) != 1) 4056 stop("bad input for argument 'log'") 4057 rm(log) 4058 4059 4060 sd1 <- sqrt(var1) 4061 sd2 <- sqrt(var2) 4062 rho <- cov12 / (sd1 * sd2) 4063 4064 4065 temp5 <- 1 - rho^2 4066 zedd1 <- (x1 - mean1) / sd1 4067 zedd2 <- (x2 - mean2) / sd2 4068 logpdf <- -log(2 * pi) - log(sd1) - log(sd2) - 4069 0.5 * log1p(-rho^2) + 4070 -(0.5 / temp5) * (zedd1^2 + (-2 * rho * zedd1 + zedd2) * zedd2) 4071 4072 logpdf[is.infinite(x1) | is.infinite(x2)] <- log(0) # 20141216 KaiH 4073 4074 if (log.arg) logpdf else exp(logpdf) 4075} 4076 4077 4078 4079rbinorm <- function(n, mean1 = 0, mean2 = 0, 4080 var1 = 1, var2 = 1, cov12 = 0) { 4081 4082 Y1 <- rnorm(n) 4083 Y2 <- rnorm(n) 4084 X1 <- sqrt(var1) * Y1 + mean1 4085 delta <- sqrt(var2 - (cov12^2) / var1) 4086 X2 <- cov12 * Y1 / sqrt(var1) + delta * Y2 + mean2 4087 4088 ans <- cbind(X1, X2) 4089 ans[is.na(delta), ] <- NA 4090 4091 ans 4092} 4093 4094 4095 4096 4097 binormal <- function(lmean1 = "identitylink", 4098 lmean2 = "identitylink", 4099 lsd1 = "loglink", 4100 lsd2 = "loglink", 4101 lrho = "rhobitlink", 4102 imean1 = NULL, imean2 = NULL, 4103 isd1 = NULL, isd2 = NULL, 4104 irho = NULL, imethod = 1, 4105 eq.mean = FALSE, eq.sd = FALSE, 4106 zero = c("sd", "rho")) { 4107 4108 lmean1 <- as.list(substitute(lmean1)) 4109 emean1 <- link2list(lmean1) 4110 lmean1 <- attr(emean1, "function.name") 4111 4112 lmean2 <- as.list(substitute(lmean2)) 4113 emean2 <- link2list(lmean2) 4114 lmean2 <- attr(emean2, "function.name") 4115 4116 lsd1 <- as.list(substitute(lsd1)) 4117 esd1 <- link2list(lsd1) 4118 lsd1 <- attr(esd1, "function.name") 4119 4120 lsd2 <- as.list(substitute(lsd2)) 4121 esd2 <- link2list(lsd2) 4122 lsd2 <- attr(esd2, "function.name") 4123 4124 lrho <- as.list(substitute(lrho)) 4125 erho <- link2list(lrho) 4126 lrho <- attr(erho, "function.name") 4127 4128 4129 4130 4131 trivial1 <- is.logical(eq.mean) && length(eq.mean) == 1 && !eq.mean 4132 trivial2 <- is.logical(eq.sd ) && length(eq.sd ) == 1 && !eq.sd 4133 4134 if (!is.Numeric(imethod, length.arg = 1, 4135 integer.valued = TRUE, positive = TRUE) || 4136 imethod > 2) 4137 stop("argument 'imethod' must be 1 or 2") 4138 4139 new("vglmff", 4140 blurb = c("Bivariate normal distribution\n", 4141 "Links: ", 4142 namesof("mean1", lmean1, earg = emean1 ), ", ", 4143 namesof("mean2", lmean2, earg = emean2 ), ", ", 4144 namesof("sd1", lsd1, earg = esd1 ), ", ", 4145 namesof("sd2", lsd2, earg = esd2 ), ", ", 4146 namesof("rho", lrho, earg = erho )), 4147 constraints = eval(substitute(expression({ 4148 4149 4150 constraints.orig <- constraints 4151 M1 <- 5 4152 NOS <- M / M1 4153 4154 cm1.m <- 4155 cmk.m <- kronecker(diag(NOS), rbind(diag(2), matrix(0, 3, 2))) 4156 con.m <- cm.VGAM(kronecker(diag(NOS), rbind(1, 1, 0, 0, 0)), 4157 x = x, 4158 bool = .eq.mean , # 4159 constraints = constraints.orig, 4160 apply.int = TRUE, 4161 cm.default = cmk.m, 4162 cm.intercept.default = cm1.m) 4163 4164 4165 cm1.s <- 4166 cmk.s <- kronecker(diag(NOS), 4167 rbind(matrix(0, 2, 2), diag(2), matrix(0, 1, 2))) 4168 con.s <- cm.VGAM(kronecker(diag(NOS), rbind(0, 0, 1, 1, 0)), 4169 x = x, 4170 bool = .eq.sd , # 4171 constraints = constraints.orig, 4172 apply.int = TRUE, 4173 cm.default = cmk.s, 4174 cm.intercept.default = cm1.s) 4175 4176 4177 con.use <- con.m 4178 for (klocal in seq_along(con.m)) { 4179 con.use[[klocal]] <- 4180 cbind(con.m[[klocal]], 4181 con.s[[klocal]], 4182 kronecker(matrix(1, NOS, 1), diag(5)[, 5])) 4183 4184 } 4185 4186 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 4187 predictors.names = predictors.names, 4188 M1 = 5) 4189 }), list( .zero = zero, 4190 .eq.sd = eq.sd, 4191 .eq.mean = eq.mean ))), 4192 4193 infos = eval(substitute(function(...) { 4194 list(M1 = 5, 4195 Q1 = 2, 4196 expected = TRUE, 4197 multipleResponses = FALSE, 4198 parameters.names = c("mean1", "mean2", "sd1", "sd2", "rho"), 4199 eq.mean = .eq.mean , 4200 eq.sd = .eq.sd , 4201 zero = .zero ) 4202 }, list( .zero = zero, 4203 .eq.mean = eq.mean, 4204 .eq.sd = eq.sd ))), 4205 4206 initialize = eval(substitute(expression({ 4207 Q1 <- 2 4208 4209 temp5 <- 4210 w.y.check(w = w, y = y, 4211 ncol.w.max = 1, 4212 ncol.y.max = Q1, 4213 ncol.y.min = Q1, 4214 out.wy = TRUE, 4215 colsyperw = Q1, 4216 maximize = TRUE) 4217 w <- temp5$w 4218 y <- temp5$y 4219 4220 4221 4222 predictors.names <- c( 4223 namesof("mean1", .lmean1 , earg = .emean1 , short = TRUE), 4224 namesof("mean2", .lmean2 , earg = .emean2 , short = TRUE), 4225 namesof("sd1", .lsd1 , earg = .esd1 , short = TRUE), 4226 namesof("sd2", .lsd2 , earg = .esd2 , short = TRUE), 4227 namesof("rho", .lrho , earg = .erho , short = TRUE)) 4228 4229 extra$colnames.y <- colnames(y) 4230 4231 if (!length(etastart)) { 4232 imean1 <- rep_len(if (length( .imean1 )) .imean1 else 4233 weighted.mean(y[, 1], w = w), n) 4234 imean2 <- rep_len(if (length( .imean2 )) .imean2 else 4235 weighted.mean(y[, 2], w = w), n) 4236 isd1 <- rep_len(if (length( .isd1 )) .isd1 else sd(y[, 1]), n) 4237 isd2 <- rep_len(if (length( .isd2 )) .isd2 else sd(y[, 2]), n) 4238 irho <- rep_len(if (length( .irho )) .irho else 4239 cor(y[, 1], y[, 2]), n) 4240 4241 if ( .imethod == 2) { 4242 imean1 <- abs(imean1) + 0.01 4243 imean2 <- abs(imean2) + 0.01 4244 } 4245 etastart <- 4246 cbind(theta2eta(imean1, .lmean1 , earg = .emean1 ), 4247 theta2eta(imean2, .lmean2 , earg = .emean2 ), 4248 theta2eta(isd1, .lsd1 , earg = .esd1 ), 4249 theta2eta(isd2, .lsd2 , earg = .esd2 ), 4250 theta2eta(irho, .lrho , earg = .erho )) 4251 } 4252 }), list( .lmean1 = lmean1, .lmean2 = lmean2, 4253 .emean1 = emean1, .emean2 = emean2, 4254 .lsd1 = lsd1 , .lsd2 = lsd2 , .lrho = lrho, 4255 .esd1 = esd1 , .esd2 = esd2 , .erho = erho, 4256 .imethod = imethod, 4257 .imean1 = imean1, .imean2 = imean2, 4258 .isd1 = isd1, .isd2 = isd2, 4259 .irho = irho ))), 4260 linkinv = eval(substitute(function(eta, extra = NULL) { 4261 NOS <- ncol(eta) / c(M1 = 5) 4262 mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 ) 4263 mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 ) 4264 fv.mat <- cbind(mean1, mean2) 4265 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 4266 } , list( .lmean1 = lmean1, .lmean2 = lmean2, 4267 .emean1 = emean1, .emean2 = emean2, 4268 .lsd1 = lsd1 , .lsd2 = lsd2 , .lrho = lrho, 4269 .esd1 = esd1 , .esd2 = esd2 , .erho = erho ))), 4270 4271 last = eval(substitute(expression({ 4272 misc$link <- c("mean1" = .lmean1 , 4273 "mean2" = .lmean2 , 4274 "sd1" = .lsd1 , 4275 "sd2" = .lsd2 , 4276 "rho" = .lrho ) 4277 4278 misc$earg <- list("mean1" = .emean1 , 4279 "mean2" = .emean2 , 4280 "sd1" = .esd1 , 4281 "sd2" = .esd2 , 4282 "rho" = .erho ) 4283 4284 misc$expected <- TRUE 4285 misc$multipleResponses <- FALSE 4286 }) , list( .lmean1 = lmean1, .lmean2 = lmean2, 4287 .emean1 = emean1, .emean2 = emean2, 4288 .lsd1 = lsd1 , .lsd2 = lsd2 , .lrho = lrho, 4289 .esd1 = esd1 , .esd2 = esd2 , .erho = erho ))), 4290 loglikelihood = eval(substitute( 4291 function(mu, y, w, residuals = FALSE, eta, 4292 extra = NULL, 4293 summation = TRUE) { 4294 mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 ) 4295 mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 ) 4296 sd1 <- eta2theta(eta[, 3], .lsd1 , earg = .esd1 ) 4297 sd2 <- eta2theta(eta[, 4], .lsd2 , earg = .esd2 ) 4298 Rho <- eta2theta(eta[, 5], .lrho , earg = .erho ) 4299 4300 if (residuals) { 4301 stop("loglikelihood residuals not implemented yet") 4302 } else { 4303 ll.elts <- 4304 c(w) * dbinorm(x1 = y[, 1], x2 = y[, 2], 4305 mean1 = mean1, mean2 = mean2, 4306 var1 = sd1^2, var2 = sd2^2, 4307 cov12 = Rho * sd1 * sd2, 4308 log = TRUE) 4309 if (summation) { 4310 sum(ll.elts) 4311 } else { 4312 ll.elts 4313 } 4314 } 4315 } , list( .lmean1 = lmean1, .lmean2 = lmean2, 4316 .emean1 = emean1, .emean2 = emean2, 4317 .lsd1 = lsd1 , .lsd2 = lsd2 , .lrho = lrho, 4318 .esd1 = esd1 , .esd2 = esd2 , .erho = erho, 4319 .imethod = imethod ))), 4320 vfamily = c("binormal"), 4321 validparams = eval(substitute(function(eta, y, extra = NULL) { 4322 mean1 <- eta2theta(eta[, 1], .lmean1, earg = .emean1) 4323 mean2 <- eta2theta(eta[, 2], .lmean2, earg = .emean2) 4324 sd1 <- eta2theta(eta[, 3], .lsd1 , earg = .esd1 ) 4325 sd2 <- eta2theta(eta[, 4], .lsd2 , earg = .esd2 ) 4326 Rho <- eta2theta(eta[, 5], .lrho , earg = .erho ) 4327 okay1 <- all(is.finite(mean1)) && 4328 all(is.finite(mean2)) && 4329 all(is.finite(sd1 )) && all(0 < sd1) && 4330 all(is.finite(sd2 )) && all(0 < sd2) && 4331 all(is.finite(Rho )) && all(abs(Rho) < 1) 4332 okay1 4333 } , list( .lmean1 = lmean1, .lmean2 = lmean2, 4334 .emean1 = emean1, .emean2 = emean2, 4335 .lsd1 = lsd1 , .lsd2 = lsd2 , .lrho = lrho, 4336 .esd1 = esd1 , .esd2 = esd2 , .erho = erho, 4337 .imethod = imethod ))), 4338 4339 4340 4341 simslot = eval(substitute( 4342 function(object, nsim) { 4343 4344 pwts <- if (length(pwts <- object@prior.weights) > 0) 4345 pwts else weights(object, type = "prior") 4346 if (any(pwts != 1)) 4347 warning("ignoring prior weights") 4348 eta <- predict(object) 4349 mean1 <- eta2theta(eta[, 1], .lmean1 , earg = .emean1 ) 4350 mean2 <- eta2theta(eta[, 2], .lmean2 , earg = .emean2 ) 4351 sd1 <- eta2theta(eta[, 3], .lsd1 , earg = .esd1 ) 4352 sd2 <- eta2theta(eta[, 4], .lsd2 , earg = .esd2 ) 4353 Rho <- eta2theta(eta[, 5], .lrho , earg = .erho ) 4354 rbinorm(nsim * length(sd1), 4355 mean1 = mean1, mean2 = mean2, 4356 var1 = sd1^2, var2 = sd2^2, cov12 = Rho * sd1 * sd2) 4357 } , list( .lmean1 = lmean1, .lmean2 = lmean2, 4358 .emean1 = emean1, .emean2 = emean2, 4359 .lsd1 = lsd1 , .lsd2 = lsd2 , .lrho = lrho, 4360 .esd1 = esd1 , .esd2 = esd2 , .erho = erho ))), 4361 4362 4363 4364 4365 deriv = eval(substitute(expression({ 4366 mean1 <- eta2theta(eta[, 1], .lmean1, earg = .emean1) 4367 mean2 <- eta2theta(eta[, 2], .lmean2, earg = .emean2) 4368 sd1 <- eta2theta(eta[, 3], .lsd1 , earg = .esd1 ) 4369 sd2 <- eta2theta(eta[, 4], .lsd2 , earg = .esd2 ) 4370 Rho <- eta2theta(eta[, 5], .lrho , earg = .erho ) 4371 4372 zedd1 <- (y[, 1] - mean1) / sd1 4373 zedd2 <- (y[, 2] - mean2) / sd2 4374 temp5 <- 1 - Rho^2 4375 4376 SigmaInv <- matrix(0, n, dimm(2)) 4377 SigmaInv[, iam(1, 1, M = 2)] <- 1 / ((sd1^2) * temp5) 4378 SigmaInv[, iam(2, 2, M = 2)] <- 1 / ((sd2^2) * temp5) 4379 SigmaInv[, iam(1, 2, M = 2)] <- -Rho / (sd1 * sd2 * temp5) 4380 dl.dmeans <- mux22(t(SigmaInv), y - cbind(mean1, mean2), M = 2, 4381 as.matrix = TRUE) 4382 dl.dsd1 <- -1 / sd1 + zedd1 * (zedd1 - Rho * zedd2) / (sd1 * temp5) 4383 dl.dsd2 <- -1 / sd2 + zedd2 * (zedd2 - Rho * zedd1) / (sd2 * temp5) 4384 dl.drho <- -Rho * (zedd1^2 - 2 * Rho * zedd1 * zedd2 + 4385 zedd2^2) / temp5^2 + 4386 zedd1 * zedd2 / temp5 + 4387 Rho / temp5 4388 4389 dmean1.deta <- dtheta.deta(mean1, .lmean1) 4390 dmean2.deta <- dtheta.deta(mean2, .lmean2) 4391 dsd1.deta <- dtheta.deta(sd1 , .lsd1 ) 4392 dsd2.deta <- dtheta.deta(sd2 , .lsd2 ) 4393 drho.deta <- dtheta.deta(Rho , .lrho ) 4394 dthetas.detas <- cbind(dmean1.deta, 4395 dmean2.deta, 4396 dsd1.deta, 4397 dsd2.deta, 4398 drho.deta) 4399 4400 c(w) * cbind(dl.dmeans[, 1], 4401 dl.dmeans[, 2], 4402 dl.dsd1, 4403 dl.dsd2, 4404 dl.drho) * dthetas.detas 4405 }), list( .lmean1 = lmean1, .lmean2 = lmean2, 4406 .emean1 = emean1, .emean2 = emean2, 4407 .lsd1 = lsd1 , .lsd2 = lsd2 , .lrho = lrho, 4408 .esd1 = esd1 , .esd2 = esd2 , .erho = erho, 4409 .imethod = imethod ))), 4410 4411 weight = eval(substitute(expression({ 4412 wz <- matrix(0.0, n, dimm(M)) 4413 wz[, iam(1, 1, M)] <- SigmaInv[, iam(1, 1, M = 2)] 4414 wz[, iam(2, 2, M)] <- SigmaInv[, iam(2, 2, M = 2)] 4415 wz[, iam(1, 2, M)] <- SigmaInv[, iam(1, 2, M = 2)] 4416 wz[, iam(3, 3, M)] <- (1 + 1 / temp5) / sd1^2 4417 wz[, iam(4, 4, M)] <- (1 + 1 / temp5) / sd2^2 4418 wz[, iam(3, 4, M)] <- -(Rho^2) / (temp5 * sd1 * sd2) 4419 wz[, iam(5, 5, M)] <- (1 + Rho^2) / temp5^2 4420 wz[, iam(3, 5, M)] <- -Rho / (sd1 * temp5) 4421 wz[, iam(4, 5, M)] <- -Rho / (sd2 * temp5) 4422 for (ilocal in 1:M) 4423 for (jlocal in ilocal:M) 4424 wz[, iam(ilocal, jlocal, M)] <- wz[, iam(ilocal, jlocal, M)] * 4425 dthetas.detas[, ilocal] * 4426 dthetas.detas[, jlocal] 4427 c(w) * wz 4428 }), list( .lmean1 = lmean1, .lmean2 = lmean2, 4429 .emean1 = emean1, .emean2 = emean2, 4430 .lsd1 = lsd1 , .lsd2 = lsd2 , .lrho = lrho, 4431 .esd1 = esd1 , .esd2 = esd2 , .erho = erho, 4432 .imethod = imethod )))) 4433} 4434 4435 4436 4437 4438 4439 4440 4441 4442 4443gumbelI <- 4444 function(la = "identitylink", earg = list(), ia = NULL, imethod = 1) { 4445 4446 la <- as.list(substitute(la)) 4447 earg <- link2list(la) 4448 la <- attr(earg, "function.name") 4449 4450 4451 4452 if (length(ia) && !is.Numeric(ia, length.arg = 1)) 4453 stop("'ia' must be a single number") 4454 4455 if (!is.Numeric(imethod, length.arg = 1, 4456 integer.valued = TRUE, positive = TRUE) || 4457 imethod > 2.5) 4458 stop("argument 'imethod' must be 1 or 2") 4459 4460 4461 new("vglmff", 4462 blurb = c("Gumbel's Type I Bivariate Distribution\n", 4463 "Links: ", 4464 namesof("a", la, earg = earg )), 4465 initialize = eval(substitute(expression({ 4466 if (!is.matrix(y) || ncol(y) != 2) 4467 stop("the response must be a 2 column matrix") 4468 4469 if (any(y < 0)) 4470 stop("the response must have non-negative values only") 4471 4472 4473 extra$colnames.y <- colnames(y) 4474 4475 predictors.names <- 4476 c(namesof("a", .la, earg = .earg , short = TRUE)) 4477 if (!length(etastart)) { 4478 ainit <- if (length( .ia )) rep_len( .ia , n) else { 4479 mean1 <- if ( .imethod == 1) median(y[,1]) else mean(y[,1]) 4480 mean2 <- if ( .imethod == 1) median(y[,2]) else mean(y[,2]) 4481 Finit <- 0.01 + mean(y[,1] <= mean1 & y[,2] <= mean2) 4482 (log(Finit+expm1(-mean1)+exp(-mean2))+mean1+mean2)/(mean1*mean2) 4483 } 4484 etastart <- theta2eta(rep_len(ainit, n), .la , earg = .earg ) 4485 } 4486 }), list( .ia = ia, 4487 .la = la, .earg = earg, .imethod = imethod ))), 4488 linkinv = eval(substitute(function(eta, extra = NULL) { 4489 NOS <- NCOL(eta) / c(M1 = 1) 4490 Q1 <- 2 4491 fv.mat <- matrix(1, NROW(eta), NOS * Q1) 4492 label.cols.y(fv.mat, colnames.y = extra$colnames.y, NOS = NOS) 4493 }, list( .la = la ))), 4494 last = eval(substitute(expression({ 4495 misc$link <- c("a" = .la ) 4496 misc$earg <- list("a" = .earg ) 4497 4498 misc$expected <- FALSE 4499 misc$pooled.weight <- pooled.weight 4500 }), list( .la = la, .earg = earg ))), 4501 loglikelihood = eval(substitute( 4502 function(mu, y, w, residuals = FALSE, eta, 4503 extra = NULL, 4504 summation = TRUE) { 4505 alpha <- eta2theta(eta, .la, earg = .earg ) 4506 if (residuals) { 4507 stop("loglikelihood residuals not implemented yet") 4508 } else { 4509 denom <- (alpha*y[,1] - 1) * (alpha*y[,2] - 1) + alpha 4510 mytolerance <- .Machine$double.xmin 4511 bad <- (denom <= mytolerance) # Range violation 4512 if (any(bad)) { 4513 cat("There are some range violations in @deriv\n") 4514 flush.console() 4515 denom[bad] <- 2 * mytolerance 4516 } 4517 ll.elts <- c(w) * (-y[,1] - y[,2] + alpha*y[,1]*y[,2] + log(denom)) 4518 if (summation) { 4519 sum(ll.elts) 4520 } else { 4521 ll.elts 4522 } 4523 } 4524 }, list( .la = la, .earg = earg ))), 4525 vfamily = c("gumbelI"), 4526 validparams = eval(substitute(function(eta, y, extra = NULL) { 4527 alpha <- eta2theta(eta, .la , earg = .earg ) 4528 okay1 <- all(is.finite(alpha)) 4529 okay1 4530 } , list( .la = la, .earg = earg ))), 4531 deriv = eval(substitute(expression({ 4532 alpha <- eta2theta(eta, .la, earg = .earg ) 4533 numerator <- (alpha*y[,1] - 1)*y[,2] + (alpha*y[,2] - 1)*y[,1] + 1 4534 denom <- (alpha*y[,1] - 1) * (alpha*y[,2] - 1) + alpha 4535 denom <- abs(denom) 4536 dl.dalpha <- numerator / denom + y[,1]*y[,2] 4537 dalpha.deta <- dtheta.deta(alpha, .la, earg = .earg ) 4538 c(w) * cbind(dl.dalpha * dalpha.deta) 4539 }), list( .la = la, .earg = earg ))), 4540 weight = eval(substitute(expression({ 4541 d2l.dalpha2 <- (numerator/denom)^2 - 2*y[,1]*y[,2] / denom 4542 d2alpha.deta2 <- d2theta.deta2(alpha, .la, earg = .earg ) 4543 wz <- w * (dalpha.deta^2 * d2l.dalpha2 - d2alpha.deta2 * dl.dalpha) 4544 if (TRUE && 4545 intercept.only) { 4546 wz <- cbind(wz) 4547 sumw <- sum(w) 4548 for (iii in 1:ncol(wz)) 4549 wz[,iii] <- sum(wz[,iii]) / sumw 4550 pooled.weight <- TRUE 4551 wz <- c(w) * wz # Put back the weights 4552 } else 4553 pooled.weight <- FALSE 4554 wz 4555 }), list( .la = la, .earg = earg )))) 4556} 4557 4558 4559 4560 4561kendall.tau <- function(x, y, exact = FALSE, max.n = 3000) { 4562 4563 if ((N <- length(x)) != length(y)) 4564 stop("arguments 'x' and 'y' do not have equal lengths") 4565 4566 NN <- if (!exact && N > max.n) { 4567 cindex <- sample.int(n = N, size = max.n, replace = FALSE) 4568 x <- x[cindex] 4569 y <- y[cindex] 4570 max.n 4571 } else { 4572 N 4573 } 4574 4575 4576 ans3 <- 4577 c( .C("VGAM_C_kend_tau", 4578 as.double(x), as.double(y), 4579 as.integer(NN), ans = double(3), 4580 NAOK = TRUE)$ans) 4581 4582 con <- ans3[1] + ans3[2] / 2 # Ties put half and half 4583 dis <- ans3[3] + ans3[2] / 2 4584 (con - dis) / (con + dis) 4585} 4586 4587 4588 4589 4590if (FALSE) 4591kendall.tau <- function(x, y, exact = TRUE, max.n = 1000) { 4592 4593 if ((N <- length(x)) != length(y)) 4594 stop("arguments 'x' and 'y' do not have equal lengths") 4595 index <- iam(NA, NA, M = N, both = TRUE) 4596 4597 index$row.index <- index$row.index[-(1:N)] 4598 index$col.index <- index$col.index[-(1:N)] 4599 4600 NN <- if (!exact && N > max.n) { 4601 cindex <- sample.int(n = N, size = max.n, replace = FALSE) 4602 index$row.index <- index$row.index[cindex] 4603 index$col.index <- index$col.index[cindex] 4604 max.n 4605 } else{ 4606 choose(N, 2) 4607 } 4608 4609 con <- sum((x[index$row.index] - x[index$col.index]) * 4610 (y[index$row.index] - y[index$col.index]) > 0) 4611 dis <- NN - con 4612 (con - dis) / (con + dis) 4613} 4614 4615 4616 4617 4618dbistudenttcop <- function(x1, x2, df, rho = 0, log = FALSE) { 4619 4620 if (!is.logical(log.arg <- log) || length(log) != 1) 4621 stop("bad input for argument 'log'") 4622 rm(log) 4623 4624 u1 <- qt(x1, df = df) 4625 u2 <- qt(x2, df = df) 4626 4627 logdensity <- 4628 -(df/2 + 1) * log1p( 4629 (u1^2 + u2^2 - 2 * rho * u1 * u2) / (df * (1 - rho^2))) - 4630 log(2*pi) - 0.5 * log1p(-rho^2) - 4631 dt(u1, df = df, log = TRUE) - 4632 dt(u2, df = df, log = TRUE) 4633 4634 if (log.arg) logdensity else exp(logdensity) 4635} 4636 4637 4638 4639 4640 4641