Hello.
Interesting problem 
I don’t have any certainties about the cause, but to me, the second formula looks a bit more ‘unstable’ to compute from a purely numerical point of view. I suppose that if your numerator and denominator are two big numbers (and as you have probably a lot of pixels, they are), you will loose precision in the final result, compared to the fraction done using already averaged numbers. I may be completely wrong though, but I already had some comparable issues in the past with the computation of fractions.
Also, just in case you could be interested by this. I use to compute my linear regressions using the command -solve that computes the least square solution to a linear system.
The idea is to build the system [X,1].[a;b] = [Y], where X and Y are the column vectors composed of the values of images X and Y you want to match, and [a;b] the vector of the coefficients you are looking for.
This small piece of code compute a and b for two input images [0] = X and [1] = Y such that ||Y - (a.X + b)||^2 is minimum:
linear_reg :
--l -y 1,100%,1,1,1 -a[0,-1] x -r 100%,{min(h,10000)} -solve. [-2] -k. a={[0]} b={[1]} -rm -endl # Get vector [a,b]
--*[0] $a -+. $b # Compute [2] = aX+b (should match [1]).
I hope this helps.