Nothing
# Automatically generated from the noweb directory
survexpmsetup <- function(rmat) {
# check the validity of the transition matrix, and determine if it
# is acyclic, i.e., can be reordered into an upper triangular matrix.
if (!is.matrix(rmat) || nrow(rmat) != ncol(rmat) || any(diag(rmat) > 0) ||
any(rmat[row(rmat) != col(rmat)] < 0))
stop ("input is not a transition matrix")
if (!is.logical(all.equal(rowSums(rmat), rep(0, ncol(rmat)))))
stop("input is not a transition matrix")
nc <- ncol(rmat)
lower <- row(rmat) > col(rmat)
if (all(rmat[lower] ==0)) return(0) # already in order
# score each state by (number of states it follows) - (number it precedes)
temp <- 1*(rmat >0) # 0/1 matrix
indx <- order(colSums(temp) - rowSums(temp))
temp <- rmat[indx, indx] # try that ordering
if (all(temp[lower]== 0)) indx # it worked!
else -1 # there is a loop in the states
}
survexpm <- function(rmat, time=1.0, setup, eps=1e-6) {
# rmat is a transition matrix, so the diagonal elements are 0 or negative
if (length(rmat)==1) exp(rmat[1]*time) #failsafe -- should never be called
else {
nonzero <- (diag(rmat) != 0)
if (sum(nonzero ==0)) diag(nrow(rmat)) # expm(0 matrix) = identity
if (sum(nonzero) ==1) { # only one state had departures
j <- which(nonzero)
emat <- diag(nrow(rmat))
temp <- exp(rmat[j,j] * time)
emat[j,j] <- temp
emat[j, -j] <- (1-temp)* rmat[j, -j]/sum(rmat[j,-j])
emat
}
else if (missing(setup) || setup[1] < 0 ||
any(diff(sort(diag(rmat)))< eps)) pade(rmat*time)
else {
if (setup[1]==0) .Call(Ccdecomp, rmat, time)$P
else {
temp <- rmat
temp[setup, setup] <- .Call(Ccdecomp, rmat[setup, setup], time)
temp$P
}
}
}
}
derivative <- function(rmat, time, dR, setup, eps=1e-8) {
if (missing(setup) || setup[1] <0 || any(diff(sort(diag(rmat)))< eps))
return (pade(rmat*time, dR*time))
if (setup==0) dlist <- .Call(Ccdecomp, rmat, time)
else dlist <- .Call(Ccdecomp, rmat[setup, setup], time)
ncoef <- dim(dR)[3]
nstate <- nrow(rmat)
dmat <- array(0.0, dim=c(nstate, nstate, ncoef))
vtemp <- outer(dlist$d, dlist$d,
function(a, b) {
ifelse(abs(a-b)< eps, time* exp(time* (a+b)/2),
(exp(a*time) - exp(b*time))/(a-b))})
# two transitions can share a coef, but only for the same X variable
for (i in 1:ncoef) {
G <- dlist$Ainv %*% dR[,,i] %*% dlist$A
V <- G*vtemp
dmat[,,i] <- dlist$A %*% V %*% dlist$Ainv
}
dlist$dmat <- dmat
# undo the reordering, if needed
if (setup[1] >0) {
indx <- order(setup)
dlist <- list(P = dlist$P[indx, indx],
dmat = apply(dmat,1:2, function(x) x[indx, indx]))
}
dlist
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.