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)
	for(j in 2:(n-1)) {
		M1[j] <- max(k(yy[j],(j-1)/n),k(yy[j],j/n))
	}
	M1[1] <- k(yy[1],1/n)
	M1[n] <- k(yy[n],(n-1)/n)
	return(max(M1))
}

K <- 500
lambda.sim <- rep(NA,K)
S <- c(seq(600,1000,100))
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=""))
}