Estimating the distribution of hail sizes is crucial for assessing related weather hazards and potential damage to buildings, vehicles and agriculture. In this study, we present a novel technique for estimating the hail size number distribution (HSND) using polarimetric C-band radar data. A generalized additive model (GAM) is employed to estimate two empirical moments of the HSND, which is then reconstructed using double-moment normalization. This approach capitalizes on the relative invariance of the double-moment normalized HSND. The model is trained on data from the Swiss network of automatic hail sensors, spanning from September 2018 to August 2024 and covering three regions of Switzerland particularly prone to hail. Several polarimetric features are extracted from a 3D radar composite that combines observations from all operational Swiss radars. Among the various extracted features, the model selects the echo-top height of 50 dBZ reflectivity value at vertical polarization and the volume of the region with a cross-correlation coefficient rhoHV below 0.97, as these provided the best predictive performance. Radar-derived HSND estimates show good agreement with independent hail sensor observations. Additionally, the model is evaluated through comparisons with photogrammetric drone surveys and crowd-sourced reports of hail. This technique enables high spatio-temporal resolution (1 km and 5 minutes) retrievals of HSND and related metrics, such as kinetic energy. Further ground observations, particularly drone-based, are essential for more comprehensive evaluation of the retrieved HSND.