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