It seems to me (check it first) that instead of

expSum = Exp[sum - Max[sum]];
q = expSum/Total[expSum];
lnQ = Log[q];
logL = Total[Take[lnQ, n]]/n;

you can write

logL = Total[Take[sum, n]]/n - Log[Total[Exp[sum]]];

which is symbolically simpler. I don't know if this helps with the optimization.

It's also a slightly better fit again. Still, I would like to know how it compares to a purely numerical fit to a compiled version of the code as given, since the main calculation is a matrix multiplication, which should be very fast. (Surely the symbolic expression generated internally by FindMaximum loses the matrix form?)

But thanks so much!

Ha! You are right it doesn't change the output, and removing it makes it 10x faster again. So the lesson here folks is improvements in the underlying function can make the biggest difference. In my defense I was following a template made by others, and now I have to figure out why, for example, the 'divide by the max value' step was in there. Maybe it is used later. But even if it is, it can be left out of the fitting process!

Having said this, I still want to know how one would code this as a direct numerical optimization, ideally specifying the 12 parameters as a single vector (list). It is possible?