slope.ci <- function(x, y, alpha=.05, epsilon=.01) {

compute.k <- function(n, p, alpha = 0.05)
{
k <- qbinom(alpha/2, n, p)
if(k > 0) {
k <- k - 1
}
k0 <- k
beta <- 0
while(beta < alpha) {
gama <- beta
beta <- pbinom(k, n, p) + 1 - pbinom(n - k - 1, n, p) 
k <- k + 1
}
k <-  k - 1
d1 <- beta - alpha
d2 <- alpha - gama
if(d1 > d2) {
k <- k - 1
beta <- gama
}
return(k)
}

compute.c <- function(eps) {
b2 <- qnorm(.5/(1-eps))
d1 <- 0
if(eps==0.010) d1 <- .019
if(eps==0.025) d1 <- .046
if(eps==0.050) d1 <- .096
if(eps==0.100) d1 <- .198
if(eps==0.150) d1 <- .339
if(eps==0.200) d1 <- .505
if(eps==0.250) d1 <- .73
if(d1 == 0) stop("can't use that epsilon")
b1 <- qnorm(.5/(1-eps))*(d1 + sqrt(1 + d1^2))
# b1 = max bias centercept of RMS
# b2 = max bias median
p1 <- pnorm(b1)
p2 <- pnorm(b2)
return( p1 + p2 - 2 * p1 * p2 )
}

p <- (1-epsilon)*compute.c(epsilon)
n <- length(x)
k <- compute.k(n=n, p=p, alpha=alpha)
slope <- repeated.r(x=x, y=y)
inte <- median( y - slope * (x - median(x)) )
w <- sort((y - inte)/(x - median(x)))
return( c(w[k+1], w[n-k]) )
}



repeated.r <- function(x, y) {

n <- length(y)
a <- rep(0,n)
for(i in 1:n) {
	tmp <- (x[-i] != x[i] )
	tmp.x <- (x[-i])[tmp]
	tmp.y <- (y[-i])[tmp]
	a[i] <- median( (y[i] - tmp.y) / (x[i] - tmp.x))
}
return(median(a))

}