MLE for beta binomial estimation

Why not use the method of moments estimates? They are more stable in small samples anyways.

(cut and pasted from my R code, should be easy enough to port if you choose to go this way:

xb <- mean(na.omit(x))
s2 <- var(x)
a <- xb * (((xb*(1-xb))/s2)-1)
b <- (1-xb) * (((xb*(1-xb))/s2)-1)

If you must use the MLE, look into Simas (2010) for bias correction approaches to beta MLEs.

What's your exact setup? Do you have n trials and for each you draw theta_i from Beta(alpha,beta) and have observed X_i? And you want to find point alpha, beta to maximize marginal likelihood of the X_is? i.e.

L(alpha,beta) = P(X_1:n | alpha, beta) = int_{theta_{1:n}} P(X_i | theta_i) P(theta_i | alpha,beta) d theta_{1:n}?

I believe what your looking for is described in a Kevin Murphy paper: Empirical Bayes for Beta-Binomial model (2007). If you have some variant of this problem (like a shared theta amongst trials), its not too difficult to get there from this paper?

Update: Actually he only does a method-of-moments estimate, I'm pretty sure its possible to numerically optimize the empirical bayes likelihood, I'll think about it.

You only have two parameters to estimate and you can compute and store the gradient and full Hessian. So you can simply use Newton's method. You'll need a numerical library for evaluating the Gamma function and its derivatives. The GNU Scientific Library has everything you need. Its relatively easy to call it from Java via JNI.