Thermodynamic functions of crystalline solids can be derived from heat capacity data. Some sets of model parameters used to fit low temperature heat capacity measurements are purely empirical, other parameters are motivated by theory, e.g. Debye and Einstein-temperatures. Determining the amount and the significance of model parameters is a challenging task and is attacked by frequentist and Bayesian statistics using concepts of recently developed machine learning tools. In a first calculation a comprehensive model is introduced (with a large set of parameters) and transformed to a model with a reduced parameter set.  Unknown systematic errors and fluctuations are estimated from the difference of the calculated values to the measured values. In this context, model parameters and the variances corresponding to the measuring points are simultaneously optimized. In case of purely empirical polynomial fits the optimized fitting function for a certain temperature range is determined avoiding a high correlation of the model parameters. The thermodynamically motivated Debye-Einstein approach is also validated. A measure for the amount of Einstein integrals to be used in order to specify the temperature dependent heat capacity function of the specific crystalline solid is provided. The calculated model parameters and uncertainties are used to derive enthalpies, entropies and Gibbs energies including their uncertainties.