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 12fill <- 13fill1 <- fill2 <- fill3 <- 14 function(x, values = 0, ncolx = ncol(x)) { 15 x <- as.matrix(x) 16 matrix(values, nrow = nrow(x), ncol = ncolx, byrow = TRUE) 17} 18 19 20 21 22 23extract.arg <- function(a) { 24 s <- substitute(a) 25 as.character(s) 26} 27 28 29 30 31 32 33remove.arg <- function(string) { 34 35 nc <- nchar(string) 36 bits <- substring(string, 1:nc, 1:nc) 37 b1 <- (1:nc)[bits == "("] 38 b1 <- if (length(b1)) b1[1]-1 else nc 39 if (b1 == 0) 40 return("") 41 string <- paste(bits[1:b1], collapse = "") 42 string 43} 44 45 46 47 48 49 50add.arg <- function(string, arg.string) { 51 52 if (arg.string == "") 53 return(string) 54 nc <- nchar(string) 55 lastc <- substring(string, nc, nc) 56 if (lastc == ")") { 57 if (substring(string, nc-1, nc-1) == "(") { 58 paste(substring(string, 1, nc-2), "(", arg.string, ")", 59 sep = "") 60 } else { 61 paste(substring(string, 1, nc-1), ", ", arg.string, ")", 62 sep = "") 63 } 64 } else { 65 paste(string, "(", arg.string, ")", sep = "") 66 } 67} 68 69 70 71 72 73 74get.arg <- function(string) { 75 76 nc <- nchar(string) 77 bits <- substring(string, 1:nc, 1:nc) 78 b1 <- (1:nc)[bits == "("] 79 b2 <- (1:nc)[bits == ")"] 80 b1 <- if (length(b1)) min(b1) else return("") 81 b2 <- if (length(b2)) max(b2) else return("") 82 if (b2-b1 == 1) "" else paste(bits[(1+b1):(b2-1)], collapse = "") 83} 84 85 86 87 88 89 90 91 eifun <- function(i, n) 92 cbind(as.numeric((1:n) == i)) 93 94 95 96 eifun <- 97 I.col <- function(i, n) 98 diag(n)[, i, drop = FALSE] 99 100 101 102 eijfun <- function(i, n) { 103 temp <- matrix(0, n, 1) 104 if (length(i)) 105 temp[i, ] <- 1 106 temp 107} 108 109 110 111 112 113 114tapplymat1 <- 115 function(mat, function.arg = c("cumsum", "diff", "cumprod")) { 116 117 118 if (!missing(function.arg)) 119 function.arg <- as.character(substitute(function.arg)) 120 function.arg <- match.arg(function.arg, 121 c("cumsum", "diff", "cumprod"))[1] 122 123 type <- switch(function.arg, cumsum = 1, diff = 2, cumprod = 3, 124 stop("argument 'function.arg' not matched")) 125 126 if (!is.matrix(mat)) 127 mat <- as.matrix(mat) 128 NR <- nrow(mat) 129 NC <- ncol(mat) 130 fred <- .C("tapply_mat1", mat = as.double(mat), as.integer(NR), 131 as.integer(NC), as.integer(type)) # , PACKAGE = "VGAM" 132 dim(fred$mat) <- c(NR, NC) 133 dimnames(fred$mat) <- dimnames(mat) 134 switch(function.arg, 135 cumsum = fred$mat, 136 diff = fred$mat[, -1, drop = FALSE], 137 cumprod = fred$mat) 138} 139 140 141 142 143 144 145matrix.power <- function(wz, M, power, fast = TRUE) { 146 147 148 149 150 n <- nrow(wz) 151 index <- iam(NA, NA, M, both = TRUE, diag = TRUE) 152 dimm.value <- if (is.matrix(wz)) ncol(wz) else 1 153 if (dimm.value > M*(M+1)/2) 154 stop("too many columns") 155 156 157 if (M == 1 || dimm.value == M) { 158 WW <- wz^power # May contain NAs 159 return(t(WW)) 160 } 161 162 if (fast) { 163 k <- veigen(t(wz), M = M) # matrix.arg) 164 evals <- k$values # M x n 165 evects <- k$vectors # M x M x n 166 } else { 167 stop("sorry, cannot handle matrix-band form yet") 168 k <- unlist(apply(wz, 3, eigen), use.names = FALSE) 169 dim(k) <- c(M, M+1, n) 170 evals <- k[, 1, , drop = TRUE] # M x n 171 evects <- k[, -1, , drop = TRUE] # M x M x n 172 } 173 174 temp <- evals^power # Some values may be NAs 175 176 177 index <- as.vector( matrix(1, 1, M) %*% is.na(temp) ) 178 179 180 index <- (index == 0) 181 if (!all(index)) { 182 warning("Some weight matrices have negative ", 183 "eigenvalues. They will be assigned NAs") 184 temp[,!index] <- 1 185 } 186 187 WW <- mux55(evects, temp, M = M) 188 WW[,!index] <- NA 189 WW 190} 191 192 193 194 195 196 197ResSS.vgam <- function(z, wz, M) { 198 199 200 if (M == 1) 201 return(sum(c(wz) * c(z^2))) 202 203 wz.z <- mux22(t(wz), z, M = M, as.matrix = TRUE) 204 sum(wz.z * z) 205} 206 207 208 209 210 211 212wweighted.mean <- function(y, w = NULL, matrix.arg = TRUE) { 213 if (!matrix.arg) 214 stop("currently, argument 'matrix.arg' must be TRUE") 215 y <- as.matrix(y) 216 M <- ncol(y) 217 n <- nrow(y) 218 if (M == 1) { 219 if (missing(w)) mean(y) else sum(w * y) / sum(w) 220 } else { 221 if (missing(w)) y %*% rep(1, n) else { 222 numer <- mux22(t(w), y, M, as.matrix = TRUE) 223 numer <- t(numer) %*% rep(1, n) 224 denom <- t(w) %*% rep(1, n) 225 denom <- matrix(denom, 1, length(denom)) 226 if (matrix.arg) 227 denom <- m2a(denom, M = M)[, , 1] 228 c(solve(denom, numer)) 229 } 230 } 231} 232 233 234 235 236 237 238veigen <- function(x, M) { 239 240 241 n <- ncol(x) 242 index <- iam(NA, NA, M = M, both = TRUE, diag = TRUE) 243 dimm.value <- nrow(x) # usually M or M(M+1)/2 244 245 z <- .Fortran("veigenf", 246 as.integer(M), 247 as.integer(n), 248 as.double(x), 249 values = double(M * n), 250 as.integer(1), 251 vectors = double(M*M*n), 252 double(M), 253 double(M), 254 wk = double(M*M), 255 as.integer(index$row), as.integer(index$col), 256 as.integer(dimm.value), 257 error.code = integer(1)) 258 259 if (z$error.code) 260 stop("eigen algorithm (rs) returned error code ", z$error.code) 261 ord <- M:1 262 dim(z$values) <- c(M, n) 263 z$values <- z$values[ord, , drop = FALSE] 264 dim(z$vectors) <- c(M, M, n) 265 z$vectors <- z$vectors[, ord, , drop = FALSE] 266 return(list(values = z$values, 267 vectors = z$vectors)) 268} 269 270 271 272 273 274 275ima <- function(j, k, M) { 276 if (length(M) > 1 || M <= 0 || j <= 0 || k <= 0 || 277 j > M || k > M) 278 stop("input wrong in ima()") 279 m <- diag(M) 280 m[col(m) <= row(m)] <- 1:(M*(M+1)/2) 281 if (j >= k) m[j, k] else m[k, j] 282} 283 284 285 286 287 288 289checkwz <- function(wz, M, trace = FALSE, 290 wzepsilon = .Machine$double.eps^0.75) { 291 if (wzepsilon > 0.5) 292 warning("argument 'wzepsilon' is probably too large") 293 if (!is.matrix(wz)) 294 wz <- as.matrix(wz) 295 wzsubset <- wz[, 1:M, drop = FALSE] 296 if (any(is.na(wzsubset))) 297 stop("NAs found in the working weights variable 'wz'") 298 if (any(!is.finite(wzsubset))) 299 stop("Some elements in the working weights variable 'wz' are ", 300 "not finite") 301 302 if ((temp <- sum(wzsubset < wzepsilon))) 303 warning(temp, " diagonal elements of the working weights variable ", 304 "'wz' have been replaced by ", signif(wzepsilon, 5)) 305 wz[, 1:M] <- pmax(wzepsilon, wzsubset) 306 wz 307} 308 309 310 311 312 313 314 315 316label.cols.y <- 317 function(answer, 318 colnames.y = NULL, 319 NOS = 1, 320 percentiles = c(25, 50, 75), 321 one.on.one = TRUE, 322 byy = TRUE) { 323 if (!is.matrix(answer)) 324 answer <- as.matrix(answer) 325 326 if (one.on.one) { 327 colnames(answer) <- 328 if (length(colnames.y) == ncol(answer)) 329 colnames.y else NULL 330 return(answer) 331 } 332 333 334 335 if (is.null(percentiles)) 336 percentiles <- c(25, 50, 75) # Restore to the default 337 338 if (!is.Numeric(percentiles) || 339 min(percentiles) <= 0 || 340 max(percentiles) >= 100) 341 stop("values of 'percentiles' should be in [0, 100]") 342 343 percentiles <- signif(percentiles, digits = 5) 344 345 ab1 <- rep(as.character(percentiles), length = ncol(answer)) 346 ab1 <- paste(ab1, "%", sep = "") 347 if (NOS > 1) { 348 suffix.char <- if (length(colnames.y) == NOS) 349 colnames.y else as.character(1:NOS) 350 ab1 <- paste(ab1, rep(suffix.char, each = length(percentiles)), 351 sep = "") 352 } 353 colnames(answer) <- ab1 354 355 356 if (byy) { 357 answer <- 358 answer[, interleave.VGAM(.M = NCOL(answer), 359 M1 = NOS), # length(percentiles)), 360 drop = FALSE] 361 } 362 answer 363} 364 365 366 367 368 369 370 371 372