function(x, s = rep(1, lx), mu = rep(0, lx), p = NULL, conf = 0.95, dist = 
	"norm", df = NULL)
{
# Last edited 20 March 2011.
#
# This function calculates the simultaneous confidence intervals
# that are less likely to contain 0, as described in 
# Benjamini and Stark JASA 1996 paper.
# x is a vector of estimates of the parameter for which the CI are needed.
# s is a vector of standard errors of the estimators. (1 by default)
# mu is a vector of parameters which is tested using the CI. (0 by default)
# p is an optional vector of probabilities, satisfying
#   (1-p1)*(1-p2)*...*(1-pn)=1-p=simultaneous confidence level desired.
#   If the vector is not supplied, 
#   geometrically decreasing probabilities are used.
# conf is the simultaneous confidence level desired.
# dist is the distribution of the estimates : either "norm" or "t". 
# df is a vector of degrees of freedom for each estimator that should
#   be specified if "t" distribution is used. 
#   (We make now use only of min(df))
#
#Checking for consistency of arguments
#
	lx <- length(x)
	if(lx == 0)
		stop("No Data!")
	if(length(s) != lx)
		stop("The number of standard errors specified is not the same as the number of estimates"
			)
	if(length(mu) != lx)
		stop("The number of parameters specified is not the same as the number of estimates"
			)
	if(dist == "t") if(length(df) != lx)
			stop("The number of degrres of freedom  specified is not the same as the number of estimates"
				)	#
#Defining functions
#
	crit <- function(i, conf, dist, df)
	{
		a <- (conf)^(2/(i * (i + 1)))
		b <- 1:i
		e <- (1 - a^b)/2
		if(dist == "norm")
			crit <- qnorm(1 - e)
		if(dist == "t")
			crit <- qt(1 - e, min(df))
		crit
	}
#Equation (19) in paper
	Rlw <- function(l, w, oax, d)
	{
		max((oax - d[w])[w >= l])
	}
#Equation (20) in paper
	Llw <- function(l, w, oax, d)
	{
		min((oax + d[w])[w <= l])
	}
#Equation (21) in paper
	fi <- function(l, j, lx)
	{
		fi <- 1:lx
		fi[j] <- l
		if(min(l, j) > 1)
			for(k in 1:(min(l, j) - 1))
				fi[k] <- l - k
		if((l > j) & (j + 1 <= l))
			for(k in (j + 1):l)
				fi[k] <- l - (k - 1)
		if((l < j) & (l <= (j - 1)))
			for(k in l:(j - 1))
				fi[k] <- k + 1
		fi
	}
#
#Defining vector of probabilities
#
	if(length(p) != lx) d <- crit(lx, conf, dist)	#Stndardizing
	oldx <- x
	x <- (x - mu)/s	#Begin main loop
#
	absx <- abs(x)
	orx <- order(absx)[lx:1]
	oax <- absx[orx]
	t <- oax	#lower ends
	u <- oax + d	#upper ends
	for(j in lx:1) {
		for(l in lx:1) {
			L <- Llw(l, fi(l, j, lx), oax, d)
			R <- Rlw(l, fi(l, j, lx), oax, d)
			if(L >= R)
				u[j] <- max(u[j], L)
			if(((oax[j] - d[l]) < 0) & ((oax[j] - d[l]) < (.Uminus(
				R))))
				t[j] <- min(t[j], oax[j] - d[j])
			else t[j] <- min(t[j], R)
		}
	}
#End of main loop
#
#Arranging the results accoding to original order and signs
#
	scu <- cbind(oldx, oldx, oldx)
	for(i in 1:lx)
		if(x[orx[i]] >= 0) {
			scu[orx[i], 1] <- t[i] * s[orx[i]] + mu[orx[i]]
			scu[orx[i], 3] <- u[i] * s[orx[i]] + mu[orx[i]]
		}
		else {
			scu[orx[i], 3] <- .Uminus(t[i]) * s[orx[i]] + mu[orx[i]
				]
			scu[orx[i], 1] <- .Uminus(u[i]) * s[orx[i]] + mu[orx[i]
				]
		}
	dimnames(scu)[[2]] <- c("Lower End", "Estimate of Center", "Upper End")
	scu
}