k <- function(x,y) {
 	return(x*log(x/y)+(1-x)*log((1-x)/(1-y)))
}

LRStat <- function(n) {
 	x <- runif(n)		#generates random numbers
	yy <- sort(x)
 	M1 <- rep(NA,n+2)
 	for(j in 2:(n-1)) {
 		M1[j] <- max(k((j-1)/n,yy[j]),k(j/n, yy[j]))
 	}
 	M1[1] <- k(1/n, yy[1])
 	M1[n] <- k((n-1)/n,yy[n])
         M1[n+1]<--log(1-yy[1])
          M1[n+2]<- -log(yy[n])
 	return(max(M1))
 }
 
K <- 500
lambda.sim <- rep(NA,K)
S <- c(seq(10,100,10),seq(125,500,25))
for (s in S) {
m <- 100000
results <- rep(NA,m)
for (i in 1:m){
	results[i] <- LRStat(s)
}
zz <- sort(results)
lambda.sim[s] <- zz[0.95*m]
print(paste("lambda[",s,"]=",lambda.sim[s],sep=""))
}