The Baryonic Tully-Fisher relation (BTFR) is a clear manifestation of the underlying physics of galaxy formation. As such, it is used to constrain and test galaxy formation and evolution models. Of particular interest, apart from the slope of the relation, is its intrinsic scatter. In this paper, we use the EAGLE simulation to study the dependence of the BTFR on the size of the simulated galaxy sample. The huge number of datapoint available in the simulation is indeed not available with current observations. Observational studies that computed the BTFR used various (small) size samples with the only obligation to have galaxies spanning over a large range of masses and rotation rates. Accordingly, to compare observational and theoretical results, we build a large number of various size datasets using the same criterion and derive the BTFR for all of them. Unmistakably, their is an effect of the number of galaxies used to derive the relation. The smaller the number, the larger the standard deviation around the average slope and intrinsic scatter of a given size sample of galaxies. This observation allows us to alleviate the tensions between observational measurements and LCDM predictions. Namely, the size of the observational samples adds up to the complexity in comparing observed and simulated relations to discredit or confirm LCDM. Similarly, samples, even large, that do not reflect the galaxy distribution give on average biased results. Large size samples reproducing the underlying distribution of galaxies constitute a supplementary necessity to compare efficiently observations and simulations.