question on fitting a function to an actual distribution

Hello:

I have been pushing forward, using a function in gp-pari script written 
by Max Alekseyev, which finds how many sets of square triads might sum 
to a given number. For example n=194 results in [[0, 5, 13], [1, 7, 12], 
[3, 4, 13], [3, 8, 11], [5, 5, 12], [7, 8, 9]] or six triads, the number 
194 is the smallest for 6 triads. His program works very well, and I 
verified it for the values of 1 <= n <= 1e7. (first 240 values of rel 
max highs)

I have a collection of 41,374 triad counts for the 1 <= count <= 41374 
with the n value given.

This distribution seems to follow:

(1)  n =  5.32459104391076 * count ^ 1.93934231682073

which is a linear graphic when plotted on a log-log scale.

But at the very first the distribution curiously drops below (1) for the 
first dozen or two dozen values.

Would anyone be curious to find a better fit for the value of n versus 
the count value of square triads? I am happy to send you the 41374 row 
file of currently determined values, for testing purposes.

Randall