> 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) } >