In this paper, a novel approach to estimate fractional snow cover (FSC) from MODIS data in a complex and heterogeneous Alpine terrain is represented by using a state-of-the-art nonparametric spline regression method, namely, multivariate adaptive regression splines (MARS). For this purpose, twenty MODIS - Landsat 8 image pairs acquired between April 2013 and December 2016 over European Alps are used. Fifteen of the image pairs are employed during model training and five images are reserved as an independent test dataset. MARS models are trained by using MODIS top-of-atmosphere reflectance values of bands 1-7, normalized difference snow index, normalized difference vegetation index and land cover class as predictor variables. Reference FSC maps are generated from higher spatial resolution Landsat 8 binary snow cover maps. Multilayer feedforward artificial neural network (ANN) models are also trained by using the same input data. During the training and the testing, the effects of the training data size and the sampling type on the predictive performance of ANN and MARS models are investigated. An additional search is also conducted to reveal whether the choice of the transfer function used in the output layer of ANN has a significant contribution to the network's FSC mapping performance. The final ANN and MARS FSC products are at 500 m spatial resolution. The results on the independent test scenes indicate that the developed ANN models with linear and hyperbolic tangent transfer functions in the output layer and the MARS models are in good agreement with reference FSC data with the same average values of R = 0.93. In contrast, the standard MODIS snow fraction product, namely, MOD10 FSC, exhibits slightly poorer performance with average R = 0.88. The proposed MARS approach is statistically proven to have the same performance with ANN, yet it is computationally more efficient in model building.