LISTINGS OF SPLUS FUNCTIONS (written and tested in ver. 3.4) REFERENCED IN LECTURE NOTES AND HANDOUTS, Stat 798C ============================================================ > TrialSim function(Delta, nsiz = 60, Nrep = 2500, lam0 = 0.7, cut1 = 2.4, cut2 = 1.7) { Utim <- runif(Nrep * nsiz) Ttim <- rexp(Nrep * nsiz, rate = lam0) Zgp <- rbinom(Nrep * nsiz, 1, 0.5) Ttim <- ifelse(Zgp == 1, Ttim * exp(Delta), Ttim) Timtmp1 <- matrix(ifelse(Ttim < 1 - Utim, Ttim, 1 - Utim), ncol = nsiz) Dthtmp <- matrix(ifelse(Ttim < 1 - Utim, 1, 0), ncol = nsiz) Onevec <- rep(1, nsiz) Tim1A <- c((Timtmp1 * Zgp) %*% Onevec) Tim1 <- c(Timtmp1 %*% Onevec) Dth1A <- c((Dthtmp * Zgp) %*% Onevec) Dth1 <- c(Dthtmp %*% Onevec) Timtmp2 <- matrix(ifelse(Ttim < 2 - Utim, Ttim, 2 - Utim), ncol = nsiz) Dthtmp <- matrix(ifelse(Ttim < 2 - Utim, 1, 0), ncol = nsiz) Tim2A <- c((Timtmp2 * Zgp) %*% Onevec) Tim2 <- c(Timtmp2 %*% Onevec) Dth2A <- c((Dthtmp * Zgp) %*% Onevec) Dth2 <- c(Dthtmp %*% Onevec) S1 <- log(((Dth1A + 0.0001) * Tim1)/((Dth1 + 0.0001) * Tim1A)) * sqrt( Dth1) S2 <- log((Dth2A * Tim2)/(Dth2 * Tim2A)) * sqrt(Dth2) Rej <- ifelse(S1 > cut1 | S2 > cut2, 1, 0) TT1 <- c(Timtmp1 %*% Onevec) TT2 <- c(Timtmp2 %*% Onevec) TimeTest <- ifelse(S1 <= cut1, TT2, TT1) ET1T2 <- function(lam) { Int1 <- (1 - (1 - exp( - lam))/lam)/lam c(Int1, exp( - lam) * Int1 + (1 - exp( - lam))/lam) } e12 <- 0.5 * (ET1T2(lam0) + ET1T2(lam0/exp(Delta))) TT1 <- TT1/nsiz - e12[1] TT2 <- TT2/nsiz - e12[2] Bmat <- matrix(c(mean(TT1^2), rep(mean(TT1 * TT2), 2), mean(TT2^2)), ncol = 2) ET <- mean(TimeTest)/nsiz mu1 <- mean(TT1) mu2 <- mean(TT2) cvec <- c(mean(TT1 * TimeTest), mean(TT2 * TimeTest))/nsiz - ET * c(mu1, mu2) VarT <- var(TimeTest)/(nsiz^2) ContrVar <- ifelse(S1 > cut1, 1, 0) TimEst <- ET - sum(solve(Bmat, cvec) * c(mu1, mu2)) VarFac <- 1 - (sum(solve(Bmat, cvec) * cvec)/VarT) # Note that in final output, Time-on-Test for each replicated trial is # calculated on a per-person basis list(Power = mean(Rej), ET = ET, VarT = VarT, S1sq = c(mean(S1^2), var( S1^2)), S2sq = c(mean(S2^2), var(S2^2)), EY = mean(ContrVar), EYT = mean(ContrVar * TimeTest)/nsiz, TT1mean = mu1, TT2mean = mu2, TTvar = Bmat, Et12Tim = cvec, TimEst = TimEst, VarFac = VarFac) } > SmExpTst <- function(lamsamp, lamval, nsiz = 40, Nrep = 10000, Upper = T) { # Function to simulate the rejection probability (and confidence # interval for this probability) for a one-sided hypothesis test of # lam0=1 vs lam1<1, based upon nominal sample size nsiz and number # Nrep of replications , where lamval <=1 is the actual exponential # hazard-rate parameter for which the probability is desired, and # lamsamp is the exponential rate parameter under which the # simulation is conducted before importance-sampling adjustment cutoff <- nsiz + qnorm(0.95) * sqrt(nsiz) largarr <- matrix(rexp(Nrep * nsiz, rate = lamsamp), ncol = nsiz) sstat <- c(largarr %*% rep(1, nsiz)) modstat <- if(Upper) ifelse(sstat > cutoff, exp(sstat * (lamsamp - lamval) - nsiz * log(lamsamp/lamval)), 0) else ifelse( sstat < cutoff, exp(sstat * (lamsamp - lamval) - nsiz * log(lamsamp/lamval)), 0) aa <- if(Upper) c(mean(modstat), sum(rep(1, Nrep)[sstat > cutoff])) else c(1 - mean(modstat), sum(rep(1, Nrep)[sstat < cutoff])) ## Outputs are respectively: estimated probability, variance of ## estimator, and number of nonzero values averaged to obtain ### the estimated probability. c(aa[1], var(modstat), aa[2]) } =================================================================== > Gradmat function(parvec, infcn, eps = 1e-06) { # Function to calculate the difference-quotient approximate gradient # (matrix) of an arbitrary input (vector) function infcn # Now recoded to use central differences ! dd <- length(parvec) aa <- length(infcn(parvec)) epsmat <- (diag(dd) * eps)/2 gmat <- array(0, dim = c(aa, dd)) for(i in 1:dd) gmat[, i] <- (infcn(parvec + epsmat[, i]) - infcn(parvec - epsmat[, i]))/eps if(aa > 1) gmat else c(gmat) } NRroot function(inipar, infcn, nmax = 25, stoptol = 1e-05, eps = 1e-06, gradfunc = NULL) { ### This function implements a simple Newton-Raphson iteration ## to find the root x of the (vector) function infcn(x), ## but for the solution to be well-defined and for NRroot to ## work properly, the dimension q of the function values infcn(x) ## and of the parameter x should both be the same. In this case, ## if the gradient gradfunc is supplied, it should be the qxq ## Jacobian matrix of the function infcn with respect to its ## arguments. Function checked/updated as of 3/31/03 assign("Infcn", infcn, frame = 0) assign("Eps", eps, frame = 0) if(is.null(gradfunc)) gradfunc <- function(x) Gradmat(x, Infcn, Eps) ctr <- 0 newpar <- inipar oldpar <- inipar - 1 while(ctr < nmax & sqrt(sum((newpar - oldpar)^2)) > stoptol) { oldpar <- newpar newpar <- oldpar - solve(gradfunc(oldpar), infcn(oldpar)) ctr <- ctr + 1 } list(nstep = ctr, initial = inipar, final = newpar, funcval = infcn(newpar)) } > GradSrch function(inipar, infcn, step, nmax = 25, stoptol = 1e-05, unitfac = F, eps = 1e-06, gradfunc = NULL) { # Function to implement Steepest-descent. The unitfac condition # indicates whether or not the supplied step-length factor(s) multiply # the negative gradient itself, or the unit vector in the same direction. assign("Infcn", infcn, frame = 0) assign("Eps", eps, frame = 0) if(is.null(gradfunc)) gradfunc <- function(x) Gradmat(x, Infcn, Eps) steps <- if(length(step) > 1) step else rep(step, nmax) newpar <- inipar oldpar <- newpar - 1 ctr <- 0 while(ctr < nmax & sqrt(sum((newpar - oldpar)^2)) > stoptol) { ctr <- ctr + 1 oldpar <- newpar newstep <- gradfunc(oldpar) newstep <- if(unitfac) newstep/sum(newstep^2) else newstep newpar <- oldpar - steps[ctr] * newstep } list(nstep = ctr, initial = inipar, final = newpar, funcval = infcn(newpar)) }