The accurate description of the structural and thermodynamic properties of ferroelectrics has been one of the most remarkable achievements of density functional theory (DFT). However, running large simulation cells with DFT is computationally demanding, while simulations of small cells are often plagued with nonphysical effects that are a consequence of the system's finite size. To avoid these finite-size effects one is thus often forced to use empirical models that describe the physics of the material in terms of effective interaction terms, that are fitted using the results from DFT. In this study we use a machine-learning (ML) potential trained on DFT, in combination with accelerated sampling techniques, to converge the thermodynamic properties of barium titanate (BTO) with first-principles accuracy and a full atomistic description. Our results indicate that the predicted Curie temperature depends strongly on the choice of DFT functional and system size, because of emergent long-range directional correlations in the local dipole fluctuations. Our findings demonstrate how the combination of ML models and traditional bottom-up modeling allow one to investigate emergent phenomena with the accuracy of first-principles calculations over the large size and time scales afforded by empirical models.