Edit: Ok, I was going to take a stab at getting Newton to work, but I'm a little confused-- I don't see it computing terms of the Weiner prior or the Hessian matrix.

Edit 2: Optimization wise, you binary search the entire rating space for every player every main loop-- it seems like you could decrease the space you binary search with each loop. For example, take the maximum change from the last loop, and search a space twice that big instead of going from MAX_RATING all the way to MIN_RATING. This ought to speed up later iterations a fair amount without changing the result. I bet you could also vary CLOSE_ENOUGH, start at (say) 5 or 10 on the first iteration and decrease it logarithmically every iteration until it gets down to .1-- that ought to decrease time spent in the inner loop in early iterations.

I'm not really much of a mathematician, so I was blindly following Remi's example. But from my limited understanding, the Hessian matrix is just used to execute Newton's method, so when I switched to bisection I ditched it. Similarly the Weiner prior is implemented by making virtual draws between a players previous/next vpds.

Where is the sigma parameter buried in the program? I don't see it, and you need a value even if you're only doing a bisection method.

I got as far as implementing a test version in C++ some time ago. If you want to cross-check test data, feel free to drop me a line.

sigma? To me that implies an algorithm that takes a player's previous rating, has an uncertainty value (this is sigma), and then adds in new game results. This algorithm recalculates using all results so you don't really have the concept of the old rating.