##The command
##source("http://parker.ad.siu.edu/Olive/lspack.txt")
##is an easy way to get these functions into R.


# lspack is a collection of R functions
#for Math 582  Large Sample Theory.


belimsim<-function(n = 100, p = 5, k = 1, nruns = 100,  
eps = 0.1, shift = 9, type = 1, psi=0.0, BB=1000, alph = 0.05){
##Gets proportion of times backwards elimination for MLR picks the full model. 
#needs library(leaps), n > 5p, p > 2
#Uses five iid error distributions: 
# type = 1 for N(0,1) errors, 2 for t3 errors, 3 for exp(1) - 1 errors
# 4 for uniform(-1,1) errors, 5 for (1-eps) N(0,1) + eps N(0,(1+shift)^2) errors.
# constant = 1 so there are p = q+1 coefficients
#1 <= k <= p-1, k is the number of nonnoise variables
#Need k < p-1 so there are some zeroes in the model. 
#need p > 2, beta = (1, 1, ..., 1, 0, ..., 0) with k+1 ones p-k-1  zeroes
# Multiply x by A: for MVN data this results
# in a covariance matrix with eigenvector c(1, ..., 1)^T
# corresponding to the largest eigenvalue. As psi gets
# close to 1, the data clusters about the line in the
# direction of (1, ..., 1)^T. 
# cor(X_i,X_j) = [2 psi +(q-2)psi^2]/[1 + (q-1)psi^2], i not = j
# when the correlation exists.
        q <- p-1
       	b <- 0 * 1:q
        b[1:k] <- 1 #b[1:0] acts like b[1:1] = b[1]
        beta <- c(1,b)       
        rho <- (2*psi + (q-2)*psi^2)/(1 + (q-1)*psi^2)
        A <- matrix(psi,nrow=q,ncol=q)
        diag(A) <- 1
        one <- as.vector(0*1:(k+1) + 1)
        vmax <- min(p,as.integer(n/5))
        vars <- as.vector(1:(p-1))
        fullct=0
	for(i in 1:nruns) {
		x <- matrix(rnorm(n * q), nrow = n, ncol = q)
                x <- x %*% A
		if(type == 1) {
			y <- 1 + x %*% b + rnorm(n)
					}
		if(type == 2) {
			y <- 1 + x %*% b + rt(n, df = 3)
					}
		if(type == 3) {
			y <- 1 + x %*% b + rexp(n) - 1
					}
		if(type == 4) {
			y <- 1 + x %*% b + runif(n, min = -1, max = 1)
					}
		if(type == 5) {
			err <- rnorm(n, sd = 1 + rbinom(n, 1, eps) * shift)
			y <- 1 + x %*% b + err
					} #make an MLR data set
#get the minimum Cp submodel
tem<-regsubsets(x,y,nvmax=vmax,method="backward")
out<-summary(tem)
mincp <- out$which[out$cp==min(out$cp)]
#do not need the constant in vin
vin <- vars[mincp[-1]]
if(length(vin)==(p-1))
 fullct=fullct+1		
           }
fullct<-fullct/nruns
list(fullprop=fullct)}



mlrsim<-function(n = 100, p = 4, k = 1, nruns = 100,  
eps = 0.1, shift = 9, type = 1, psi=0.0, indices=c(1,2), alph = 0.05){
#Simulates CIs for OLS MLR for unscaled and scaled predictors.
#Uses five iid error distributions: 
# type = 1 for N(0,1) errors, 2 for t3 errors, 3 for exp(1) - 1 errors
# 4 for uniform(-1,1) errors, 5 for (1-eps) N(0,1) + eps N(0,(1+shift)^2) errors.
# constant = 1 so there are p = q+1 coefficients
#1 <= k <= p-1, k is the number of nonnoise variables
#Need k < p-1 so there are some zeroes in the model. 
#need p > 2, beta = (1, 1, ..., 1, 0, ..., 0) with k+1 ones p-k-1  zeroes
# Multiply x by A: for MVN data this results
# in a covariance matrix with eigenvector c(1, ..., 1)^T
# corresponding to the largest eigenvalue. As psi gets
# close to 1, the data clusters about the line in the
# direction of (1, ..., 1)^T. 
# cor(X_i,X_j) = [2 psi +(q-2)psi^2]/[1 + (q-1)psi^2], i not = j
# when the correlation exists.
        q <- p-1
       	b <- 0 * 1:q
        b[1:k] <- 1 #b[1:0] acts like b[1:1] = b[1]
        beta <- c(1,b)
        cicov <- 0*(1:p)
        avelen <- 0*(1:p)
        cicovs <- 0*(1:p)
        avelens <- 0*(1:p)
        rho <- (2*psi + (q-2)*psi^2)/(1 + (q-1)*psi^2)
        A <- matrix(psi,nrow=q,ncol=q)
        diag(A) <- 1
        covx <- A %*% t(A)
        sds <- sqrt(diag(covx)) #pop sds for the predictors
        sds<- c(1,sds)
        betas <- beta*sds #pop beta for OLS using scaled predictors
        val <- 1-alph/2
        valf <- 1 - alph
        dfd <- n-p
        dfn <- length(indices)
        cut <- qt(val,df=dfd)
        fcut <- dfn*qf(valf,df1=dfn,df2=dfd)
        Up <- 1:p
        Lp <- Up
        Ups <- Up
        Lps <- Up
        one <- 1 + 0*1:n
        testcov <- 0
        testcovs <- 0
	for(i in 1:nruns) {
		x <- matrix(rnorm(n * q), nrow = n, ncol = q)
                x <- x %*% A
		if(type == 1) {
			y <- 1 + x %*% b + rnorm(n)
					}
		if(type == 2) {
			y <- 1 + x %*% b + rt(n, df = 3)
					}
		if(type == 3) {
			y <- 1 + x %*% b + rexp(n) - 1
					}
		if(type == 4) {
			y <- 1 + x %*% b + runif(n, min = -1, max = 1)
					}
		if(type == 5) {
			err <- rnorm(n, sd = 1 + rbinom(n, 1, eps) * shift)
			y <- 1 + x %*% b + err
					} #make an MLR data set
		out <-lsfit(x,y) #
                temp <- ls.diag(out)
                ses <- temp$std.err
                mse <- temp$std.dev^2
                bhat <- out$coef
                tem <- cut*ses
                Up <- bhat+tem
                Lp <- bhat-tem
                w<-scale(x, center = FALSE, scale = apply(x, 2, sd, na.rm = TRUE))
                outs <- lsfit(w,y)
                temps <- ls.diag(outs)
                sess <- temps$std.err
                mses <- temps$std.dev^2
                bhats <- outs$coef
                tems <- cut*sess
                Ups <- bhats+tems
                Lps <- bhats-tems
                for (j in 1:p){                   
                   if(beta[j] >= Lp[j] && beta[j] <= Up[j])
                    cicov[j] <- cicov[j] + 1
                    avelen[j] <- avelen[j] + Up[j] - Lp[j]
                  #"CIs for OLS" using scaled predictors
                   if(betas[j] >= Lps[j] && betas[j] <= Ups[j])
                    cicovs[j] <- cicovs[j] + 1
                    avelens[j] <- avelens[j] + Ups[j] - Lps[j]
                    } 
                #test whether beta_I = 0 with I = indices
               x1 <- cbind(one,x)
               xtx <- t(x1)%*%x1
               cinv <- solve(xtx)                             
               diff <- bhat[indices] - beta[indices]
               dsq <- t(diff) %*% solve(cinv[indices,indices]) %*% diff
               dsq<-dsq[1,1]  #make dsq a real instead of a matrix
               if(dsq < mse*fcut) 
               testcov <- testcov +1  
                #test whether beta_I = 0 with I = indices and scaled predictors
               w1 <- cbind(one,w)
               wtw <- t(w1)%*%w1
               cinv <- solve(wtw)                             
               diff <- bhats[indices] - betas[indices]
               dsq <- t(diff) %*% solve(cinv[indices,indices]) %*% diff
               dsq<-dsq[1,1]  #make dsq a real instead of a matrix
               if(dsq < mses*fcut) 
               testcovs <- testcovs +1              
           }
cicov <- cicov/nruns
avelen <- sqrt(n)*avelen/nruns
cicovs <- cicovs/nruns
avelens <- sqrt(n)*avelens/nruns
testcov <- testcov/nruns
testcovs <- testcovs/nruns
list(cicov=cicov,avelen=avelen,testcov=testcov,beta=beta,
cicovs=cicovs,avelens=avelens,testcovs=testcovs,betas=betas,k=k)}