The equation for N(k) is the Poisson Distribution. Stirling's Formula was used to simplify the computations. The deviations of the data points from the solid blue line appear to statistical in nature. The largest earthquakes were omitted because of their low probabilities. The fit indicates the the peak magnitude is at approximately M 1.175.
(edit: In the plot what was referred to as the variance (Var) was not the average but the weighted sum of the square of the deviations of the natural logarithms of the number of earthquakes in each interval. It is the sum of the square of the errors for all the events which was the function that was optimized in the fit.)