A key component of Earthquake Early Warning Systems is the development of an accurate and robust predictive model relating the recorded waveform to seismicity and/or ground shaking characteristics. This has been traditionally established by adopting deterministic calibration methods and aggregating or even binning the available data within a region. Unfortunately, these practices ignore the fact that the use of limited information to establish these models within a real-time estimation setting, typically corresponding to the first few seconds after the identification of the seismic wave arrival, introduces significant uncertainty (variability) in the predictions and a strong dependence of the results on regional wave propagation effects. For developing robust predictive models, the resultant variability needs to be explicitly and formally considered in the model calibration stage.  This dissertation presents a probabilistic Bayesian inference methodology to address this challenge. To better focus on the research effort, the calibration of earthquake magnitude regression models is investigated, while the application to the Sichuan region of Southwestern China is examined. Specifically, the records from the 2008 Wenchuan earthquake and the 2013 Lushan earthquake sequences are used to develop a relationship between the seismic maximum predominant period and the earthquake magnitude. Both classical and hierarchical Bayesian calibration formulations are examined. The classical Bayesian learning accommodates the following improvements: a model class selection is established, comparing across different candidate models to promote the most appropriate from accuracy and robustness perspectives; the full posterior distribution of the model parameters is identified, quantifying relevant uncertainties in their values given the available data; a heteroscedastic model is considered for the estimation error variance; the dependence of the predictive model robustness on binning practices for the available data is investigated through a comprehensive inclusion of the regression error on the analysis. The hierarchical Bayesian formulation extends these efforts, exploring how regional wave propagation effects impact the predictive models, while simultaneously allowing for the use of a sufficient amount of data for the calibration scheme. Within this hierarchical updating framework, the predictive models for each station are separately calibrated (first level), while sharing information across the stations by selecting common calibration hyper-parameters (second level). A computational framework relying on Markov Chain Monte Carlo techniques is established for accommodating the developments for both the classical and hierarchical calibration. This pertains to both sampling from the posterior distribution of the updated parameters, as well as the estimation of model evidence for the modal class selection and averaging. Results demonstrate the value of the formal implementation of Bayesian model updating for the calibration of predictive models for Earthquake Early Warning applications, as it accommodates a more comprehensive incorporation of the different sources of uncertainty/error impacting the development and implementation of these models. This ultimately improves the robustness of predictions. Additionally, results demonstrate the value of hierarchical calibration formulations, contributing to a significant reduction of the prediction variability that originates from regional wave propagation effects.