In many scientific applications, it is important to investigate how certain properties scale with the parameters of the systems. The experimental studies of scalings have traditionally been addressed with log regression, which limits the results to power laws and to theoretical and not data-driven dimensionless quantities. This has also been the case in nuclear fusion, in which the scaling of the energy confinement time is a crucial aspect in understanding the physics of transport and in the design of future devices. Traditionally two main assumptions are at the basis of the most widely accepted empirical scaling laws for the confinement time: (a) the dimensionless variables used are the ones derived from the symmetries of the Vlasov equation; (b) the final scalings have the mathematical form of power laws. In this paper, it is shown how symbolic regression (SR), implemented with genetic programming (GP) techniques, can be used to test these hypotheses. Neither assumption is confirmed by the available data of the multi-machine International Tokamak Physics Activity (ITPA) of validated tokamak discharges. The statistically soundest expressions are not power laws and cannot be formulated in terms of the traditional dimensionless quantities. The consequences of the data-driven scaling laws obtained are both practical and theoretical: the confinement time for the ITER can be significantly shorter than foreseen by power laws and different dimensionless variables should be considered for theoretical investigations. On the other hand, higher quality databases should be built to reduce the uncertainties in the extrapolations. It is also worth emphasising that the proposed methodology is fully general and therefore can be applied to any field of science.

The most widely used forms of model selection criteria, the Bayesian Information Criterion (BIC) and the Akaike Information Criterion (AIC), are expressed in terms of synthetic indicators of the residual distribution: the variance and the mean-squared error of the residuals respectively. In many applications in science, the noise affecting the data can be expected to have a Gaussian distribution. Therefore, at the same level of variance and mean-squared error, models, whose residuals are more uniformly distributed, should be favoured. The degree of uniformity of the residuals can be quantified by the Shannon entropy. Including the Shannon entropy in the BIC and AIC expressions improves significantly these criteria. The better performances have been demonstrated empirically with a series of simulations for various classes of functions and for different levels and statistics of the noise. In presence of outliers, a better treatment of the errors, using the Geodesic Distance, has proved essential.