accisimf<- function(n = 50, N = 500, nruns = 500, p = 0.5, alpha = 0.05){ #Simulates the classical, Agresti Coull type and modified CIs #for p for finite samples. The AC CI seems to work well for #n < 0.4 N, N large. The classical CI seems to work well for #p = 0.01 if n > 0.3 N and N is large. # Changes CI [L,U] to [max(0,L),min(1,U)] to get shorter lengths ccov <- 0 accov <- 0 mcov <- 0 cut <- 1 - alpha/2 clow <- 1:nruns cup <- clow aclow <- clow acup <- clow mlow <- clow mup <- cup cp <- clow acp <- clow val <- qnorm(cut) valsq <- val^2 nt <- n + valsq sqrtfpc <- sqrt((N - n)/N) pop <- 0*1:N # make p = pop p approx nominal p pp <- as.integer(p*N) pop[1:pp] <- 1 p <- mean(pop) for(i in 1:nruns) { x <- sample(pop, n) # get classical CI phat <- mean(x) se <- sqrt((phat * (1 - phat))/(n - 1)) * sqrtfpc tem <- val * se clow[i] <- phat - tem cup[i] <- phat + tem clow[i] <- max(0, clow[i]) cup[i] <- min(1, cup[i]) if(clow[i] <= p && cup[i] >= p) ccov <- ccov + 1 cp[i] <- phat #get AC CI pt <- (n * phat + 0.5 * valsq)/nt se <- sqrt((pt * (1 - pt))/nt) * sqrtfpc tem <- val * se aclow[i] <- pt - tem acup[i] <- pt + tem aclow[i] <- max(0, aclow[i]) acup[i] <- min(1,acup[i]) if(aclow[i] <= p && acup[i] >= p) accov <- accov + 1 acp[i] <- pt # get modified CI mlow[i] <- max(0,min(clow[i],aclow[i])) mup[i] <- min(1,max(cup[i],acup[i])) if(mlow[i] <= p && mup[i] >= p) mcov <- mcov + 1 } ccov <- ccov/nruns clen <- sqrt(n) * mean(cup - clow) accov <- accov/nruns aclen <- sqrt(n) * mean(acup - aclow) mcov <- mcov/nruns mlen <- sqrt(n) * mean(mup - mlow) print("p,mean(phat),mean(ptilde)") print(p) print(mean(cp)) print(mean(acp)) list(ccov = ccov, clen = clen, accov = accov, aclen = aclen, mcov = mcov, mlen = mlen) } bcisim<- function(n = 100, nruns = 500, p = 0.5, alpha = 0.05) { # Simulates the classical, Agresti Coull and exact CIs for p. # Changes CI (L,U) to (max(0,L),min(1,U)) to get shorter lengths. ccov <- 0 accov <- 0 ecov <- 0 cut <- 1 - alpha/2 cut2 <- alpha/2 clow <- 1:nruns cup <- clow aclow <- clow acup <- clow elow <- clow eup <- clow val <- qnorm(cut) valsq <- val^2 nt <- n + valsq zup <- 1/(1 + n * qf(alpha, 2 * n, 2)) zlow <- n/(n + qf((1 - alpha), 2, 2 * n)) for(i in 1:nruns) { x <- rbinom(n, 1, p) # get classical CI phat <- mean(x) se <- sqrt((phat * (1 - phat))/n) tem <- val * se clow[i] <- max(0,phat - tem) cup[i] <- min(1,phat + tem) if(clow[i] < p && cup[i] > p) ccov <- ccov + 1 #get AC CI pt <- (n * phat + 0.5 * valsq)/nt se <- sqrt((pt * (1 - pt))/nt) tem <- val * se aclow[i] <- max(0,pt - tem) acup[i] <- min(1,pt + tem) if(aclow[i] < p && acup[i] > p) accov <- accov + 1 #get exact CI w <- sum(x) if(w == 0) { elow[i] <- 0 eup[i] <- zup } else if(w == n) { elow[i] <- zlow eup[i] <- 1 } else { val2 <- w + (n - w + 1) * qf(cut, (2 * (n - w + 1)), (2 * w)) elow[i] <- w/val2 val2 <- w + 1 + (n - w) * qf(cut2, (2 * (n - w)), (2 * ( w + 1))) eup[i] <- (w + 1)/val2 } if(elow[i] < p && eup[i] > p) ecov <- ecov + 1 } ccov <- ccov/nruns clen <- sqrt(n) * mean(cup - clow) accov <- accov/nruns aclen <- sqrt(n) * mean(acup - aclow) ecov <- ecov/nruns elen <- sqrt(n) * mean(eup - elow) list(ccov = ccov, clen = clen, accov = accov, aclen = aclen, ecov = ecov, elen = elen) } ducisim<- function(n = 100, nruns = 5000, alpha = 0.05, eta = 1000) { # Simulates 100(1-alpha)% CI for eta for discrete uniform eta # data. Coverage and n CI length are also given. The CI is # [max(Y), max(Y)/alpha^(1/n)). Notice that max(Y) <= eta wp1. cov <- 0 low <- 1:nruns up <- low cut <- alpha^(1/n) pop <- 1:eta for(i in 1:nruns) { y <- sample(x = pop, size = n, replace = T) maxy <- max(y) up[i] <- maxy/cut low[i] <- maxy if(low[i] <= eta && up[i] > eta) cov <- cov + 1 } cov <- cov/nruns nlen <- n * mean(up - low) ups <- sum(up < eta)/nruns list(ups = ups, cov = cov, nlen = nlen) } expsim<- function(n = 10, nruns = 5000, theta = 0, lambda = 1, alpha = 0.05, p = 1) { #Simulates exp 100(1-alpha)% CI for lambda, CI for theta #and CI for lambda with theta known. scov <- 0 ccov <- 0 mcov <- 0 slow <- 1:nruns sup <- slow mlow <- slow mup <- slow clow <- slow cup <- slow ucut <- alpha/2 lcut <- 1 - ucut d <- 2 * (n - p) d2 <- 2 * n lval <- log(alpha/2)/n uval <- log(1 - alpha/2)/n mlval <- alpha^(-1/(n - 1)) - 1 for(i in 1:nruns) { y <- theta + lambda * rexp(n) miny <- min(y) dn <- sum((y - miny)) #get CI for lambda lamhat <- dn/n num <- 2 * dn slow[i] <- num/qchisq(lcut, df = d) sup[i] <- num/qchisq(ucut, df = d) if(slow[i] < lambda && sup[i] > lambda) scov <- scov + 1 #get Mann et al CI for theta mlow[i] <- miny - lamhat * mlval mup[i] <- miny if(mlow[i] < theta) mcov <- mcov + 1 #get CI for lambda when theta is known tn <- sum((y - theta)) num <- 2 * tn clow[i] <- num/qchisq(lcut, df = d2) cup[i] <- num/qchisq(ucut, df = d2) if(clow[i] < lambda && cup[i] > lambda) ccov <- ccov + 1 } scov <- scov/nruns slen <- sqrt(n) * mean(sup - slow) mcov <- mcov/nruns mlen <- n * mean(mup - mlow) ccov <- ccov/nruns clen <- sqrt(n) * mean(cup - clow) list(d = d, scov = scov, slen = slen, mcov = mcov, mlen = mlen, ccov = ccov, clen = clen) } hnsim<- function(n = 10, nruns = 5000, mu = 0, sigma = 1, alpha = 0.05, p = 1) { #Simulates Pewsey HN 100(1-alpha)% CI for sigma^2, CI for mu, #and CI for sigma^2 with mu known. scov <- 0 lcov <- 0 ccov <- 0 slow <- 1:nruns sup <- slow llow <- slow lup <- slow clow <- slow cup <- sup ucut <- alpha/2 lcut <- 1 - ucut sigsq <- sigma^2 d <- n - p phiinv <- qnorm((0.5 + 1/(2 * n))) lval <- log(alpha) * phiinv * (1 + 13/n^2) for(i in 1:nruns) { y <- mu + sigma * abs(rnorm(n)) miny <- min(y) dn <- sum((y - miny)^2) #get CI for sigma^2 slow[i] <- dn/qchisq(lcut, df = d) sup[i] <- dn/qchisq(ucut, df = d) if(slow[i] < sigsq && sup[i] > sigsq) scov <- scov + 1 #get CI for mu shat <- sqrt(dn/n) llow[i] <- miny + shat * lval lup[i] <- miny if(llow[i] < mu) lcov <- lcov + 1 #get CI for sigma^2 if mu is known tn <- sum((y - mu)^2) clow[i] <- tn/qchisq(lcut, df = n) cup[i] <- tn/qchisq(ucut, df = n) if(clow[i] < sigsq && cup[i] > sigsq) ccov <- ccov + 1 } scov <- scov/nruns slen <- sqrt(n) * mean(sup - slow) lcov <- lcov/nruns llen <- n * mean(lup - llow) ccov <- ccov/nruns clen <- sqrt(n) * mean(cup - clow) list(d = d, scov = scov, slen = slen, lcov = lcov, llen = llen, ccov = ccov, clen = clen) } pcisim<- function(nruns = 1000, n1 = 10, n2 = 10, mu1 = 0, mu2 = 0, var1 = 1, var2 = 1, dist = 1, alpha = 0.05) { #gets the pooled, modified pooled, and Welch 100 (1-alpha)% CIs for mu1 - mu2 #defaults are alpha = .05, lengths make the most sense if n1 = n2 #dist = 1 for x ~ N(mu1, var1), y ~ N(mu2, var2) #dist = 2 for x ~ EXP with mean mu1 variance var1, # y ~ EXP with mean mu2 variance var2 #dist = 3 for x ~ N(mu1, var1), y EXP with mean mu2 variance var2 diff <- mu1 - mu2 up <- 1 - alpha/2 pdf <- n1 + n2 - 2 rhohat <- n1/(n1 + n2) pcov <- 0 mpcov <- 0 wcov <- 0 plow <- 1:nruns pup <- plow mplow <- plow mpup <- plow wlow <- plow wup <- plow tpcut <- qt(up, pdf) mtpcut <- qt(up, n1 + n2 - 4) for(i in 1:nruns) { if(dist == 1) { x <- mu1 + sqrt(var1) * rnorm(n1) y <- mu2 + sqrt(var2) * rnorm(n2) } if(dist == 2) { x <- mu1 - sqrt(var1) + sqrt(var1) * rexp(n1) y <- mu2 - sqrt(var2) + sqrt(var2) * rexp(n2) } if(dist == 3) { x <- mu1 + sqrt(var1) * rnorm(n1) y <- mu2 - sqrt(var2) + sqrt(var2) * rexp(n2) } xbar <- mean(x) s1sq <- var(x) ybar <- mean(y) s2sq <- var(y) xmybar <- xbar - ybar sp <- (n1 - 1) * s1sq + (n2 - 1) * s2sq sp <- sqrt(sp/pdf) val <- sp * sqrt(1/n1 + 1/n2) #get pooled CI tem <- tpcut * val plow[i] <- xmybar - tem pup[i] <- xmybar + tem if(plow[i] < diff && pup[i] > diff) pcov <- pcov + 1 #get modified pooled CI thetahat <- s2sq/s1sq tauhatsq <- (1 - rhohat + rhohat * thetahat)/(rhohat + (1 - rhohat) * thetahat) tauhat <- sqrt(tauhatsq) tem <- mtpcut * val * tauhat mplow[i] <- xmybar - tem mpup[i] <- xmybar + tem if(mplow[i] < diff && mpup[i] > diff) mpcov <- mpcov + 1 #get Welch CI t1 <- s1sq/n1 t2 <- s2sq/n2 do <- (t1 + t2)^2/(t1^2/(n1 - 1) + t2^2/(n2 - 1)) d <- max(1, floor(do)) wtcut <- qt(up, d) val <- sqrt(s1sq/n1 + s2sq/n2) tem <- wtcut * val wlow[i] <- xmybar - tem wup[i] <- xmybar + tem if(wlow[i] < diff && wup[i] > diff) wcov <- wcov + 1 } pcov <- pcov/nruns plen <- sqrt(n1) * mean(pup - plow) mpcov <- mpcov/nruns mplen <- sqrt(n1) * mean(mpup - mplow) wcov <- wcov/nruns wlen <- sqrt(n1) * mean(wup - wlow) #print(mean(x)) #print(var(x)) #print(mean(y)) #print(var(y)) list(pcov = pcov, plen = plen, mpcov = mpcov, mplen = mplen, wcov = wcov, wlen = wlen ) } poiscisim<- function(n = 100, nruns = 500, theta = 1, alpha = 0.05) { # Simulates the classical, modified and ``exact" CIs for theta. ccov <- 0 mcov <- 0 ecov <- 0 cut <- 1 - alpha/2 cut2 <- alpha/2 clow <- 1:nruns cup <- clow mlow <- clow mup <- clow elow <- clow eup <- clow zup <- qchisq(1 - alpha, 2)/(2 * n) val <- qnorm(cut) valsq <- val^2 for(i in 1:nruns) { x <- rpois(n, theta) # get classical CI xbar <- mean(x) se <- sqrt(xbar/n) tem <- val * se clow[i] <- xbar - tem cup[i] <- xbar + tem if(clow[i] < theta && cup[i] > theta) ccov <- ccov + 1 #get modified CI xsum <- sum(x) tem <- xsum - 0.5 + 0.5 * valsq - val * sqrt(xsum - 0.5 + valsq/4) mlow[i] <- tem/n tem <- xsum + 0.5 + 0.5 * valsq + val * sqrt(xsum + 0.5 + valsq/4) mup[i] <- tem/n if(mlow[i] < theta && mup[i] > theta) mcov <- mcov + 1 #get Grosh exact CI w <- as.integer(xsum) if(w == 0) { elow[i] <- 0 eup[i] <- zup } else { elow[i] <- qchisq(cut2, 2 * w)/(2 * n) eup[i] <- qchisq(cut, (2 * w + 2))/(2 * n) } if(elow[i] < theta && eup[i] > theta) ecov <- ecov + 1 } ccov <- ccov/nruns clen <- sqrt(n) * mean(cup - clow) mcov <- mcov/nruns mlen <- sqrt(n) * mean(mup - mlow) ecov <- ecov/nruns elen <- sqrt(n) * mean(eup - elow) list(ccov = ccov, clen = clen, mcov = mcov, mlen = mlen, ecov = ecov, elen = elen) } raysim<- function(n = 100, mu = 1, sigma = 1, iter = 100, runs = 100, alpha = 0.05) { #For simulation of Rayleigh MLEs. # Note that mum and sigm are the initial estimators. muhat <- 1:runs sighat <- muhat mum <- muhat sigm <- muhat muo <- muhat muo2 <- muo sigo <- muhat sigo2 <- sigo mlow <- 1:runs mup <- mlow slow <- mlow sup <- mlow mcov <- 0 scov <- 0 lcut <- qnorm(1 - alpha/2) for(i in 1:runs) { w <- 2 * sigma^2 * rexp(n) y <- sqrt(w) + mu shat <- sqrt(var(y)/0.429204) mhat <- mean(y) - 1.253314 * shat mum[i] <- min(mhat, (2 * min(y) - mhat)) sigm[i] <- sqrt(sum((y - mum[i])^2)/(2 * n)) muold <- mum[i] sigold <- sigm[i] #use Newton's method to find the MLE for(j in 1:iter) { g1 <- - sum((y - muold)^(-1)) + sum(y - muold)/sigold^2 g2 <- (-2 * n)/sigold + sum((y - muold)^2)/sigold^3 d11 <- - sum((y - muold)^(-2)) - n/sigold^2 d12 <- (-2 * sum((y - muold)))/sigold^3 d21 <- d12 d22 <- (2 * n)/sigold^2 - (3 * sum((y - muold)^2))/sigold^4 D <- rbind(c(d11, d12), c(d21, d22)) told <- cbind(muold, sigold) g <- cbind(g1, g2) tnew <- t(told) - solve(D) %*% t(g) muold <- tnew[1] sigold <- tnew[2] } #muhat[i] <- tnew[1] or try muhat[i] <- min(tnew[1], (2 * min(y) - tnew[1])) #sighat[i] <- tnew[2] or try sighat[i] <- sqrt(sum((y - muhat[i])^2)/(2 * n)) muo[i] <- told[1] muo2[i] <- tnew[1] - told[1] sigo[i] <- told[2] sigo2[i] <- tnew[2] - told[2] #get CIs for mu and sigma d11 <- - sum((y - muhat[i])^(-2)) - n/sighat[i]^2 d12 <- (-2 * sum((y - muhat[i])))/sighat[i]^3 d21 <- d12 d22 <- (2 * n)/sighat[i]^2 - (3 * sum((y - muhat[i])^2))/sighat[i]^4 D <- rbind(c(d11, d12), c(d21, d22)) C <- solve(D) tem <- lcut * sqrt( - C[1, 1]) mlow[i] <- muhat[i] - tem mup[i] <- muhat[i] + tem if(mlow[i] < mu && mup[i] > mu) mcov <- mcov + 1 tem <- lcut * sqrt( - C[2, 2]) slow[i] <- sighat[i] - tem sup[i] <- sighat[i] + tem if(slow[i] < sigma && sup[i] > sigma) scov <- scov + 1 } mummn <- mean(mum) sdm <- sqrt(n * var(mum)) sigmmn <- mean(sigm) sds <- sqrt(n * var(sigm)) muhatmn <- mean(muhat) sdmuhat <- sqrt(n * var(muhat)) sighatmn <- mean(sighat) sdsighat <- sqrt(n * var(sighat)) muconv <- max(abs(muhat - muo)) muconv2 <- max(abs(muo2)) sigconv <- max(abs(sighat - sigo)) sigconv2 <- max(abs(sigo2)) mcov <- mcov/runs smlen <- sqrt(n) * mean(mup - mlow) scov <- scov/runs sslen <- sqrt(n) * mean(sup - slow) list(muconv = muconv, muconv2 = muconv2, sigconv = sigconv, sigconv2 = sigconv2, mummn = mummn, sdm = sdm, sigmmn = sigmmn, sds = sds, muhatmn = muhatmn, sdmuhat = sdmuhat, sighatmn = sighatmn, sdsighat = sdsighat, mcov = mcov, smlen = smlen, scov = scov, sslen = sslen) } shorth<-function(y, alpha = 0.05){ # computes the shorth(c) interval [Ln,Un] for c = ceiling[n(1-alpha)]. n <- length(y) cc <- ceiling(n * (1 - alpha)) sy <- sort(y) rup <- sy[cc] rlow <- sy[1] olen <- rup - rlow if(cc < n){ for(j in (cc + 1):n){ zlen <- sy[j] - sy[j - cc + 1] if(zlen < olen) { olen <- zlen rup <- sy[j] rlow <- sy[j - cc + 1] } } } list(Ln = rlow, Un = rup) } varcisim<- function(n = 100, nruns = 1000, alpha = 0.05, type = 1) { #Simulates 100(1-alpha)% CI for the variance sigsq. vcov <- 0 vlow <- 1:nruns vup <- vlow cut <- qt(1 - alpha/2, df = n - 1) sigsqln <- exp(1) * (exp(1) - 1) tav <- vcov for(i in 1:nruns) { if(type == 1) { y <- rnorm(n) sigsq <- 1 } if(type == 2) { y <- rexp(n) sigsq <- 1 } if(type == 3) { y <- exp(rnorm(n)) sigsq <- sigsqln } if(type == 4) { y <- runif(n) sigsq <- 1/12 } if(type == 5) { y <- rbinom(n, 1, 0.5) sigsq <- 0.25 } if(type == 6) { y <- rpois(n, 1) sigsq <- 1 } ssq <- var(y) tem <- (y - mean(y))^4 tem <- mean(tem) tauhat <- tem/(ssq^2) - 1 if(tauhat < 0) tauhat <- 0 val <- sqrt(n * tauhat) * cut if(tauhat < 2) { numl <- (n - 2) * ssq numu <- (n + 1) * ssq } if(tauhat >= 2 && tauhat < 4) { numl <- (n - 5) * ssq numu <- (n + 5) * ssq } if(tauhat >= 4 && tauhat < 5) { numl <- (n - 20) * ssq numu <- (n + 25) * ssq } if(tauhat >= 5 && tauhat < 6) { numl <- (n - 20) * ssq numu <- (n + 50) * ssq } if(tauhat >= 6 && tauhat < 9) { numl <- (n - 20) * ssq numu <- (n + 55) * ssq } if(tauhat >= 9) { numl <- (n - 25) * ssq numu <- (n + 55) * ssq } if(n - val < 0) { vup[i] <- 2 * max(ssq, sqrt(ssq)) } else vup[i] <- numu/(n - val) vlow[i] <- numl/(n + val) if(vlow[i] < sigsq && vup[i] > sigsq) vcov <- vcov + 1 tav[i] <- tauhat } vcov <- vcov/nruns vlen <- sqrt(n) * mean(vup - vlow) lows <- sum(vlow > sigsq)/nruns ups <- sum(vup < sigsq)/nruns tauave <- mean(tav) #could list zz <- sum(tav > 6)/nruns list(lows = lows, ups = ups, tauave = tauave, vcov = vcov, vlen = vlen) } wcisim<- function(n = 100, phi = 1, lam = 1, iter = 100, runs = 100, alpha = 0.05) { #For simulation of Weibull CIs for phi and lambda. # Note that phir and lamr are the robust initial estimators. phihat <- 1:runs lamhat <- phihat lamr <- phihat phir <- phihat phio <- phihat lamo <- phihat plow <- 1:runs pup <- plow llow <- plow lup <- plow pcov <- 0 lcov <- 0 lcut <- qnorm(1 - alpha/2) pcut <- sqrt(0.608/n) * qnorm(1 - alpha/2) for(i in 1:runs) { weib <- (lam * rexp(n))^(1/phi) lw <- log(weib) tem <- mad(lw, constant = 1) phir[i] <- 0.767049/tem ahat <- median(lw) - log(log(2))/phir[i] lamr[i] <- exp(ahat * phir[i]) #get the MLEs y <- weib suml <- sum(log(y)) pold <- phir[i] lold <- lamr[i] for(j in 1:iter) { lnew <- mean(y^pold) lvo <- lold lold <- lnew den <- (1/lold) * sum(y^pold * log(y)) - suml pnew <- n/den pvo <- pold pold <- pnew } phihat[i] <- pnew lamhat[i] <- lnew phio[i] <- pvo lamo[i] <- lvo #get CI for phi plow[i] <- phihat[i] * (1 - pcut) pup[i] <- phihat[i] * (1 + pcut) if(plow[i] < phi && pup[i] > phi) pcov <- pcov + 1 #get CI for lambda logl <- log(lamhat[i]) tem <- 1 + 0.4635 * logl + 0.5824 * (logl)^2 val <- lcut * lamhat[i] * sqrt((1.109 * tem)/n) llow[i] <- lamhat[i] - val lup[i] <- lamhat[i] + val if(llow[i] < lam && lup[i] > lam) lcov <- lcov + 1 } mnphir <- mean(phir) sdphir <- sqrt(var(phir)) mnlamr <- mean(lamr) sdlamr <- sqrt(var(lamr)) pconv <- max(abs(phihat - phio)) lconv <- max(abs(lamhat - lamo)) mnphihat <- mean(phihat) sdphihat <- sqrt(n * var(phihat)) mnlamhat <- mean(lamhat) sdlamhat <- sqrt(n * var(lamhat)) pcov <- pcov/runs lcov <- lcov/runs phiasd <- phi * sqrt(0.608) lamasd <- sqrt(1.109 * lam^2 * (1 + 0.4635 * log(lam) + 0.5282 * (log(lam))^2)) list(pconv = pconv, lconv = lconv, mnphir = mnphir, sdphir = sdphir, mnlamr = mnlamr, sdlamr = sdlamr, mnphihat = mnphihat, sdphihat = sdphihat, phiasd = phiasd, mnlamhat = mnlamhat, sdlamhat = sdlamhat, lamasd = lamasd, pcov = pcov, lcov = lcov) }