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 double.cens.normal <- 16 function(r1 = 0, r2 = 0, 17 lmu = "identitylink", 18 lsd = "loglink", 19 imu = NULL, isd = NULL, 20 zero = "sd") { 21 if (!is.Numeric(r1, length.arg = 1, integer.valued = TRUE) || 22 r1 < 0) 23 stop("bad input for 'r1'") 24 if (!is.Numeric(r2, length.arg = 1, integer.valued = TRUE) || 25 r2 < 0) 26 stop("bad input for 'r2'") 27 28 lmu <- as.list(substitute(lmu)) 29 emu <- link2list(lmu) 30 lmu <- attr(emu, "function.name") 31 32 lsd <- as.list(substitute(lsd)) 33 esd <- link2list(lsd) 34 lsd <- attr(esd, "function.name") 35 36 37 new("vglmff", 38 blurb = c("Univariate normal distribution with double censoring\n\n", 39 "Links: ", 40 namesof("mu", lmu, earg = emu, tag = TRUE), ", ", 41 namesof("sd", lsd, earg = esd, tag = TRUE), 42 "\n", 43 "Variance: sd^2"), 44 constraints = eval(substitute(expression({ 45 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 46 predictors.names = predictors.names, 47 M1 = 2) 48 }) , list( .zero = zero))), 49 50 infos = eval(substitute(function(...) { 51 list(M1 = 2, 52 Q1 = 1, 53 expected = TRUE, 54 multipleResponses = FALSE, 55 parameters.names = c("mu", "sd"), 56 lmu = .lmu , 57 lsd = .lsd , 58 zero = .zero ) 59 }, list( .zero = zero, .lmu = lmu, .lsd = lsd ))), 60 61 initialize = eval(substitute(expression({ 62 predictors.names <- 63 c(namesof("mu", .lmu , earg = .emu , tag = FALSE), 64 namesof("sd", .lsd , earg = .esd , tag = FALSE)) 65 66 if (ncol(y <- cbind(y)) != 1) 67 stop("the response must be a vector or a one-column matrix") 68 69 if (length(w) != n || 70 !is.Numeric(w, integer.valued = TRUE, positive = TRUE)) 71 stop("the argument 'weights' must be a vector ", 72 "of positive integers") 73 74 sumw <- sum(w) 75 extra$bign <- sumw + .r1 + .r2 # Tot num of censored & uncensored obsns 76 77 if (!length(etastart)) { 78 yyyy.est <- if (length( .imu )) .imu else median(y) 79 sd.y.est <- if (length( .isd )) .isd else { 80 junk <- lm.wfit(x = x, y = c(y), w = c(w)) 81 1.25 * sqrt( sum(w * junk$resid^2) / junk$df.residual ) 82 } 83 yyyy.est <- rep_len(yyyy.est , n) 84 sd.y.est <- rep_len(sd.y.est , n) 85 etastart <- cbind(mu = theta2eta(yyyy.est, .lmu , earg = .emu ), 86 sd = theta2eta(sd.y.est, .lsd , earg = .esd )) 87 } 88 }) , list( .lmu = lmu, .lsd = lsd, 89 .emu = emu, .esd = esd, 90 .imu = imu, .isd = isd, 91 .r1 = r1, .r2 = r2 ))), 92 linkinv = function(eta, extra = NULL) eta[, 1], 93 last = eval(substitute(expression({ 94 misc$link <- c(mu = .lmu , sd = .lsd ) 95 96 misc$earg <- list(mu = .emu , sd = .esd ) 97 98 misc$multipleResponses <- FALSE 99 misc$expected <- TRUE 100 misc$r1 <- .r1 101 misc$r2 <- .r2 102 }) , list( .lmu = lmu, .lsd = lsd, 103 .emu = emu, .esd = esd, 104 .r1 = r1, .r2 = r2 ))), 105 loglikelihood = eval(substitute( 106 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 107 summation = TRUE) { 108 sd <- eta2theta(eta[, 2], .lsd , earg = .esd ) 109 110 if (!summation) 111 stop("cannot handle 'summation = FALSE' yet") 112 113 if (residuals) { 114 stop("loglikelihood residuals not implemented yet") 115 } else { 116 sum(w * dnorm(y, m = mu, sd = sd, log = TRUE)) + 117 (if ( .r1 == 0) 0 else { 118 z1 <- min((y - mu) / sd) 119 Fz1 <- pnorm(z1) 120 .r1 * log(Fz1)}) + 121 (if ( .r2 == 0) 0 else { 122 z2 <- max((y - mu) / sd) 123 Fz2 <- pnorm(z2) 124 .r2 * log1p(-Fz2)}) 125 } 126 } , list( .lmu = lmu, .lsd = lsd, 127 .emu = emu, .esd = esd, 128 .r1 = r1, .r2 = r2 ))), 129 vfamily = c("double.cens.normal"), 130 validparams = eval(substitute(function(eta, y, extra = NULL) { 131 mu <- eta[, 1] 132 sd <- eta2theta(eta[, 2], .lsd , earg = .esd ) 133 okay1 <- all(is.finite(mu)) && 134 all(is.finite(sd)) && all(0 < sd) 135 okay1 136 }, list( .lmu = lmu, .lsd = lsd, 137 .emu = emu, .esd = esd, 138 .r1 = r1, .r2 = r2 ))), 139 deriv = eval(substitute(expression({ 140 sd <- eta2theta(eta[, 2], .lsd, earg =.esd) 141 142 q1 <- .r1 / extra$bign 143 q2 <- .r2 / extra$bign 144 pee <- 1 - q1 - q2 # 1 if r1==r2==0 145 z1 <- if ( .r1 == 0) - 100 else min((y - mu) / sd) # 100==Inf 146 z2 <- if ( .r2 == 0) + 100 else max((y - mu) / sd) # 100==Inf 147 fz1 <- if ( .r1 == 0) 0 else dnorm(z1) 148 fz2 <- if ( .r2 == 0) 0 else dnorm(z2) 149 Fz1 <- if ( .r1 == 0) 0.02 else pnorm(z1) # 0/0 undefined 150 Fz2 <- if ( .r2 == 0) 0.99 else pnorm(z2) 151 152 dl.dmu <- (y - mu) / sd^2 + 153 ((- .r1 * fz1/Fz1 + .r2 * fz2/(1-Fz2)) / sd) / (n*w) 154 dl.dsd <- -1/sd + (y-mu)^2 / sd^3 + 155 ((- .r1 * z1*fz1/Fz1 + .r2 * z2*fz2/(1-Fz2)) / sd) / (n*w) 156 157 dmu.deta <- dtheta.deta(mu, .lmu , earg =.emu ) 158 dsd.deta <- dtheta.deta(sd, .lsd , earg =.esd ) 159 160 c(w) * cbind(dl.dmu * dmu.deta, dl.dsd * dsd.deta) 161 }) , list( .lmu = lmu, .lsd = lsd, 162 .emu = emu, .esd = esd, 163 .r1 = r1, .r2 = r2 ))), 164 weight = expression({ 165 wz <- matrix(NA_real_, n, dimm(M)) 166 167 Q.1 <- ifelse(q1 == 0, 1, q1) # Saves division by 0 below; not elegant 168 Q.2 <- ifelse(q2 == 0, 1, q2) # Saves division by 0 below; not elegant 169 170 ed2l.dmu2 <- 1 / (sd^2) + 171 ((fz1*(z1+fz1/Q.1) - fz2*(z2-fz2/Q.2)) / sd^2) / (pee*w) 172 ed2l.dmusd <- ((fz1-fz2 + z1*fz1*(z1+fz1/Q.1) - 173 z2*fz2*(z2-fz2/Q.2)) / sd^2) / (pee*w) 174 ed2l.dsd2 <- 2 / (sd^2) + 175 ((z1*fz1-z2*fz2 + z1^2 *fz1 *(z1+fz1/Q.1) - 176 z2^2 *fz2*(z2-fz2/Q.2)) / sd^2) / (pee*w) 177 178 wz[,iam(1,1,M)] <- w * ed2l.dmu2 * dmu.deta^2 179 wz[,iam(2,2,M)] <- w * ed2l.dsd2 * dsd.deta^2 180 wz[,iam(1,2,M)] <- w * ed2l.dmusd * dsd.deta * dmu.deta 181 wz 182 })) 183} 184 185 186 187 188 189 190dbisa <- function(x, scale = 1, shape, log = FALSE) { 191 if (!is.logical(log.arg <- log) || length(log) != 1) 192 stop("bad input for argument 'log'") 193 rm(log) 194 195 196 L <- max(length(x), length(shape), length(scale)) 197 if (length(x) != L) x <- rep_len(x, L) 198 if (length(shape) != L) shape <- rep_len(shape, L) 199 if (length(scale) != L) scale <- rep_len(scale, L) 200 201 logdensity <- rep_len(log(0), L) 202 203 xok <- (x > 0) 204 xifun <- function(x) { 205 temp <- sqrt(x) 206 temp - 1/temp 207 } 208 logdensity[xok] <- 209 dnorm(xifun(x[xok] / scale[xok]) / shape[xok], log = TRUE) + 210 log1p(scale[xok]/x[xok]) - log(2) - log(shape[xok]) - 211 0.5 * log(x[xok]) - 0.5 * log(scale[xok]) 212 logdensity[scale <= 0] <- NaN 213 logdensity[shape <= 0] <- NaN 214 if (log.arg) logdensity else exp(logdensity) 215} 216 217 218 219pbisa <- function(q, scale = 1, shape, 220 lower.tail = TRUE, log.p = FALSE) { 221 222 223 ans <- pnorm(((temp <- sqrt(q/scale)) - 1/temp) / shape, 224 lower.tail = lower.tail, log.p = log.p) 225 ans[scale < 0 | shape < 0] <- NaN 226 ans[q <= 0] <- if (lower.tail) ifelse(log.p, log(0), 0) else 227 ifelse(log.p, log(1), 1) 228 ans 229} 230 231 232 233qbisa <- function(p, scale = 1, shape, 234 lower.tail = TRUE, log.p = FALSE) { 235 236 if (!is.logical(lower.tail) || length(lower.tail ) != 1) 237 stop("bad input for argument 'lower.tail'") 238 239 if (!is.logical(log.p) || length(log.p) != 1) 240 stop("bad input for argument 'log.p'") 241 242 243 A <- qnorm(p, lower.tail = lower.tail, log.p = log.p) 244 temp1 <- A * shape * sqrt(4 + A^2 * shape^2) 245 ans1 <- (2 + A^2 * shape^2 + temp1) * scale / 2 246 ans2 <- (2 + A^2 * shape^2 - temp1) * scale / 2 247 248 249 250 if (lower.tail) { 251 if (log.p) { 252 ln.p <- p 253 ans <- ifelse(exp(p) < 0.5, pmin(ans1, ans2), pmax(ans1, ans2)) 254 ans[ln.p == -Inf] <- 0 255 ans[ln.p == 0] <- Inf 256 #ans[ln.p > 0] <- NaN 257 } else { 258 ans <- ifelse(p < 0.5, pmin(ans1, ans2), pmax(ans1, ans2)) 259 #ans[p < 0] <- NaN 260 ans[p == 0] <- 0 261 ans[p == 1] <- Inf 262 #ans[p > 1] <- NaN 263 } 264 } else { 265 if (log.p) { 266 ln.p <- p 267 ans <- ifelse(-expm1(p) < 0.5, pmin(ans1, ans2), pmax(ans1, ans2)) 268 ans[ln.p == -Inf] <- Inf 269 ans[ln.p == 0] <- 0 270 #ans[ln.p > 0] <- NaN 271 } else { 272 ans <- ifelse(p > 0.5, pmin(ans1, ans2), pmax(ans1, ans2)) 273 #ans[p < 0] <- NaN 274 ans[p == 0] <- Inf 275 ans[p == 1] <- 0 276 #ans[p > 1] <- NaN 277 } 278 } 279 280 ans[scale < 0 | shape < 0] <- NaN 281 ans 282} 283 284 285 286rbisa <- function(n, scale = 1, shape) { 287 288 A <- rnorm(n) 289 temp1 <- A * shape 290 temp1 <- temp1 * sqrt(4 + temp1^2) 291 ans1 <- (2 + A^2 * shape^2 + temp1) * scale / 2 292 ans2 <- (2 + A^2 * shape^2 - temp1) * scale / 2 293 294 295 296 ans <- ifelse(A < 0, pmin(ans1, ans2), pmax(ans1, ans2)) 297 ans[shape <= 0] <- NaN 298 ans[scale <= 0] <- NaN 299 ans 300} 301 302 303 304 305 306 307 308 309 bisa <- function(lscale = "loglink", lshape = "loglink", 310 iscale = 1, ishape = NULL, 311 imethod = 1, 312 zero = "shape", 313 nowarning = FALSE) { 314 315 316 317 lshape <- as.list(substitute(lshape)) 318 eshape <- link2list(lshape) 319 lshape <- attr(eshape, "function.name") 320 321 lscale <- as.list(substitute(lscale)) 322 escale <- link2list(lscale) 323 lscale <- attr(escale, "function.name") 324 325 326 if (length(ishape) && !is.Numeric(ishape, positive = TRUE)) 327 stop("bad input for argument 'ishape'") 328 if (!is.Numeric(iscale, positive = TRUE)) 329 stop("bad input for argument 'iscale'") 330 if (!is.Numeric(imethod, length.arg = 1, 331 integer.valued = TRUE, positive = TRUE) || 332 imethod > 3) 333 stop("argument 'imethod' must be 1 or 2 or 3") 334 335 336 new("vglmff", 337 blurb = c("Birnbaum-Saunders distribution\n\n", 338 "Links: ", 339 namesof("scale", lscale, earg = escale, tag = TRUE), "; ", 340 namesof("shape", lshape, earg = eshape, tag = TRUE)), 341 constraints = eval(substitute(expression({ 342 constraints <- cm.zero.VGAM(constraints, x = x, .zero , M = M, 343 predictors.names = predictors.names, 344 M1 = 2) 345 }) , list( .zero = zero))), 346 347 infos = eval(substitute(function(...) { 348 list(M1 = 2, 349 Q1 = 1, 350 expected = TRUE, 351 multipleResponses = FALSE, 352 parameters.names = c("scale", "shape"), 353 lscale = .lscale , 354 lshape = .lshape , 355 zero = .zero ) 356 }, list( .zero = zero, .lscale = lscale, .lshape = lshape 357 ))), 358 359 360 initialize = eval(substitute(expression({ 361 if (ncol(y <- cbind(y)) != 1) 362 stop("the response must be a vector or a one-column matrix") 363 364 predictors.names <- 365 c(namesof("scale", .lscale , earg = .escale , tag = FALSE), 366 namesof("shape", .lshape , earg = .eshape , tag = FALSE)) 367 368 if (!length(etastart)) { 369 scale.init <- rep_len( .iscale , n) 370 shape.init <- if (is.Numeric( .ishape)) rep_len( .ishape , n) else { 371 if ( .imethod == 1) { 372 ybar <- rep_len(weighted.mean(y, w), n) 373 ybarr <- rep_len(1 / weighted.mean(1/y, w), n) # Reqrs y > 0 374 sqrt(ybar / scale.init + scale.init / ybarr - 2) 375 } else if ( .imethod == 2) { 376 sqrt(2*( pmax(y, scale.init+0.1) / scale.init - 1)) 377 } else { 378 ybar <- rep_len(weighted.mean(y, w), n) 379 sqrt(2*(pmax(ybar, scale.init + 0.1) / scale.init - 1)) 380 } 381 } 382 etastart <- cbind(theta2eta(scale.init, .lscale , earg = .escale ), 383 theta2eta(shape.init, .lshape , earg = .eshape )) 384 } 385 }) , list( .lshape = lshape, .lscale = lscale, 386 .ishape = ishape, .iscale = iscale, 387 .eshape = eshape, .escale = escale, 388 .imethod = imethod ))), 389 linkinv = eval(substitute(function(eta, extra = NULL) { 390 sc <- eta2theta(eta[, 1], .lscale , earg = .escale ) 391 sh <- eta2theta(eta[, 2], .lshape , earg = .eshape ) 392 sc * (1 + sh^2 / 2) 393 }, list( .lshape = lshape, .lscale = lscale, 394 .eshape = eshape, .escale = escale ))), 395 last = eval(substitute(expression({ 396 misc$link <- c(scale = .lscale , shape = .lshape ) 397 398 misc$earg <- list(scale = .escale , shape = .eshape ) 399 400 misc$expected <- TRUE 401 misc$multipleResponses <- FALSE 402 }) , list( .lshape = lshape, .lscale = lscale, 403 .eshape = eshape, .escale = escale ))), 404 loglikelihood = eval(substitute( 405 function(mu, y, w, residuals = FALSE, eta, 406 extra = NULL, 407 summation = TRUE) { 408 sc <- eta2theta(eta[, 1], .lscale , earg = .escale ) 409 sh <- eta2theta(eta[, 2], .lshape , earg = .eshape ) 410 if (residuals) { 411 stop("loglikelihood residuals not implemented yet") 412 } else { 413 ll.elts <- c(w) * dbisa(x = y, scale = sc, shape = sh, log = TRUE) 414 if (summation) { 415 sum(ll.elts) 416 } else { 417 ll.elts 418 } 419 } 420 }, list( .lshape = lshape, .lscale = lscale, 421 .eshape = eshape, .escale = escale ))), 422 423 vfamily = c("bisa"), 424 validparams = eval(substitute(function(eta, y, extra = NULL) { 425 sc <- eta2theta(eta[, 1], .lscale , earg = .escale ) 426 sh <- eta2theta(eta[, 2], .lshape , earg = .eshape ) 427 okay1 <- all(is.finite(sc)) && all(0 < sc) && 428 all(is.finite(sh)) && all(0 < sh) 429 okay1 430 }, list( .lshape = lshape, .lscale = lscale, 431 .eshape = eshape, .escale = escale ))), 432 deriv = eval(substitute(expression({ 433 sc <- eta2theta(eta[, 1], .lscale , earg = .escale ) 434 sh <- eta2theta(eta[, 2], .lshape , earg = .eshape ) 435 436 dl.dsh <- ((y/sc - 2 + sc/y) / sh^2 - 1) / sh 437 dl.dsc <- -0.5 / sc + 1/(y+sc) + sqrt(y) * ((y+sc)/y) * 438 (sqrt(y/sc) - sqrt(sc/y)) / (2 * sh^2 * sc^1.5) 439 440 dsh.deta <- dtheta.deta(sh, .lshape , earg = .eshape ) 441 dsc.deta <- dtheta.deta(sc, .lscale , earg = .escale ) 442 443 c(w) * cbind(dl.dsc * dsc.deta, 444 dl.dsh * dsh.deta) 445 }) , list( .lshape = lshape, .lscale = lscale, 446 .eshape = eshape, .escale = escale ))), 447 weight = eval(substitute(expression({ 448 wz <- matrix(NA_real_, n, M) # Diagonal!! 449 wz[, iam(2, 2, M)] <- 2 * dsh.deta^2 / sh^2 450 hfunction <- function(alpha) 451 alpha * sqrt(pi/2) - pi * exp(2/alpha^2) * 452 pnorm(2/alpha, lower.tail = FALSE) 453 wz[, iam(1, 1, M)] <- dsc.deta^2 * 454 (sh * hfunction(sh) / sqrt(2*pi) + 1) / (sh*sc)^2 455 c(w) * wz 456 }), list( .zero = zero )))) 457} 458 459 460