I guess by "closer fit" you mean that the values for lower avalanche sizes (where you have higher frequencies) should be closer to the fitted function. You may want to add weights to your fit, as the values for higher frequencies should have a higher fidelity. Consider the code
q = 1.16
s = 100
m = 100
Q(x)=(1+(1-q)*((s+x)/m))**(1/(1-q))
fit Q(x) 'tableavalanchesizeGSA' using 1:2:(1/sqrt($2)) via s, m, q
Here, the third column of the using statement is considered as the standard deviation of the value in the second column (see gnuplot documentation). This should give you a better fit, but you may need to use the correct standard deviations.
See also the gnuplot fit demos for examples with weighting.