In this paper, we introduce a statistical estimation framework for magnetic resonance fingerprinting (MRF), a recently proposed quantitative imaging paradigm. Within this framework, we present a maximum likelihood formulation to simultaneously estimate multiple parameter maps from highly undersampled, noisy k-space data. A novel iterative algorithm, based on variable splitting, the alternating direction method of multipliers, and the variable projection method, is proposed to solve the resulting optimization problem. Representative results demonstrate that compared to the conventional MRF reconstruction, the proposed method yields improved accuracy and/or reduced acquisition time. Moreover, the proposed formulation enables theoretical analysis of MRF. For example, we show that with the gridding reconstruction as an initialization, the first iteration of the proposed method exactly produces the conventional MRF reconstruction.