Y.data <- c(4,3,4,4,8,6,5,9,3,4,2,4,6,6,9); num.categories <- 10

binom.ll <- function(Y,k,p)  {
	n <- length(Y)
	return(n*mean(Y)*log(p) + (n*k - n*mean(Y))*log(1-p))	
}

ruler <- seq(0,1,length=100)

for (i in 1:length(ruler))  {
	if (i == 1)  {
		ll.pos <- ruler[i]
		ll.max <- binom.ll(Y.data,num.categories,ruler[i])
	}
	else  {
		test.val <- binom.ll(Y.data,num.categories,ruler[i]) 
		if (test.val > ll.max)  {
			ll.pos <- ruler[i]
			ll.max <- test.val
		}
	}
	
}

> ll.max
[1] -103.9197
> ll.pos    
[1] 0.5151515


round(ll.pos,3)
round(mean(Y.data)/num.categories,3)

ll.vec <- binom.ll(Y.data,num.categories,ruler)
par(mfrow=c(1,1))
plot(ruler,ll.vec,type="l",col="purple",lwd=2)
abline(h=max(ll.vec)); abline(v=ll.pos)
text(0.58,-140,paste("[",round(ll.pos,3),round(max(ll.vec),3),"]",sep=""))