The description of the behavior of the geomagnetic field in the geological past is greatly hampered by the paucity of the data both in space and time. In paleomagnetic studies, this circumstance is partly overcome through the use of the Geomagnetic Axial Dipole (GAD) model, which however becomes unsatisfactory when multipolar terms with non-zero time averages are introduced. Under these conditions, the only way to describe the spatio-temporal evolution of the paleofield is to investigate the temporal evolution of its statistical characteristics. An applicable method for such a description is the use of the Giant Gaussian Process model which does not determine the individual parameters but the probabilistic compatibility of a given model with empirical data. However, the specificity of the data entails many technical difficulties for the implementation of this method when applied to paleomagnetic and paleointensity data. An example of such an analysis applied to the Brunhes epoch clearly revealed the necessity of introducing the significant non-dipole terms in the ancient geomagnetic field configuration for this epoch. In addition, the problem of the discrepancy between models constructed separately for paleointensity and paleoinclination data was discovered. Hypotheses are to explain this discrepancy.