High Damping Rubber Bearings (HDRBs) are widely used in seismic isolation of structures, and their force-deformation constitutive behavior controls the lateral structural response during an earthquake. However, given the many components and steps involved in the preparation of the rubber compound and the manufacturing process itself, the HDRB force-deformation behavior presents significant variability with respect to nominal design values. Consequently, this paper proposes a framework for parametric probabilistic modeling of the uncertain behavior of these devices. For the validation of this framework, a recently proposed force-deformation phenomenologically based numerical model is considered. Calibration of the HDRB model is based on experimental results of 125 devices and uses the Covariance Matrix Adaptative Evolution Strategy (CMA-ES) algorithm to fit the model parameters. Since the structural response is controlled by the isolator behavior, this parametric probabilistic model provides an efficient tool to propagate uncertainty from the isolation devices to the final dynamic response of the structure. This framework can also be used to evaluate the validity of the approximate code methodology of the upper and lower bound isolator property design limits, which are applied for all isolators simultaneously.