> FitNtupl
function(ncoord, nquant, indata)
{
## assumes block of pseudo-random uniform[0,1)
## numbers in indata; to be tested for fit based on
## empirically generated contingency-table of nquant
## equal-length intervals in each of ncoord
## consecutive coordinates
ntup = length(indata) %/% ncoord
idata = c(matrix(trunc(nquant * indata - 1e-11), ncol =
ncoord) %*% nquant^(0:(ncoord - 1))) + 1
cellexp = ntup/(nquant^ncoord)
cells = table(idata)
diagind = 1 + (0:(nquant - 1)) * sum(nquant^(0:(ncoord - 1)))
chistat = sum((cells - cellexp)^2)/cellexp
diagstat = (sum(cells[diagind]) - nquant * cellexp)^2/(nquant *
cellexp * (1 - ((cellexp * nquant)/ntup)))
list(chisq = chistat, pval = 1 - pchisq(chistat,
nquant^ncoord - 1), diagstat = diagstat, diagPval =
1-pchisq(diagstat, 1), CountTbl = cells)
}

> FitNtupl(5, 5, runif(1e7))
$chisq
[1] 3115.641

$pval
[1] 0.5388012

$diagstat
[1] 0.6059696

$diagPval
[1] 0.4363094

##--------------------------------------------------------------------

> 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 
   #   exponential hazard-rate (`alternative') parameter 
   #   for which the probability is desired, and
   # lamnew  is the exponential rate parameter under 
   #   which the simulation is conducted before 
   #   importance-sampling adjustment.
   # Upper=T means that the probabilities are simulated as
   #   empirical averages of modified indicators that the
   #   sum-statistic is larger than its cutoff; Upper=F 
   #   means the indicators used are that the sum statistic
   #   falls below its cutoff.
    largarr = matrix(rexp(Nrep * nsiz, rate = lamsamp), ncol = nsiz)
    cutoff = nsiz + qnorm(.95)*sqrt(nsiz)
    sstat = c(largarr %*% rep(1, nsiz))
    del = lamsamp-lamval
    modstat = if (Upper) 
        ifelse(sstat > cutoff, 
            exp(sstat * del - nsiz * log(lamsamp/lamval)), 0)
      else ifelse(sstat < cutoff, 
            exp(sstat * del - 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]))
    c(aa[1], var(modstat), aa[2])
}

##-----------------------------------------------------------
(Older, superseded)


> TrialSim
function (Delta = 0, 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)
    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)
    Dth2 <- c(Dthtmp %*% Onevec)
    S1 <- log((`(`)/(`(`)) * sqrt(Dth1)
    S2 <- log((`(`)/(`(`)) * 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 <- (`(`)/lam
        c(Int1, exp(-lam) * Int1 + (`(`)/lam)
    }
    e12 <- 0.5 * (`(`)
    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)/(`(`)
    ContrVar <- ifelse(S1 > cut1, 1, 0)
    TimEst <- ET - sum(solve(Bmat, cvec) * c(mu1, mu2))
    VarFac <- 1 - (`(`)
    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)
}
>